dcot.f
3.6 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
double precision function dcot (x)
c may 1980 edition. w. fullerton, c3, los alamos scientific lab.
double precision x, cotcs(15), ainty, ainty2, pi2rec, sqeps,
1 xmax, xmin, xsml, y, yrem, prodbg, dint, dcsevl, d1mach,
2 dexp, dlog, dsqrt
external d1mach, dcsevl, initds
c
c series for cot on the interval 0. to 6.25000e-02
c with weighted error 5.52e-34
c log weighted error 33.26
c significant figures required 32.34
c decimal places required 33.85
c
data cot cs( 1) / +.2402591609 8295630250 9553617744 970 d+0 /
data cot cs( 2) / -.1653303160 1500227845 4746025255 758 d-1 /
data cot cs( 3) / -.4299839193 1724018935 6476228239 895 d-4 /
data cot cs( 4) / -.1592832233 2754104602 3490851122 445 d-6 /
data cot cs( 5) / -.6191093135 1293487258 8620579343 187 d-9 /
data cot cs( 6) / -.2430197415 0726460433 1702590579 575 d-11 /
data cot cs( 7) / -.9560936758 8000809842 7062083100 000 d-14 /
data cot cs( 8) / -.3763537981 9458058041 6291539706 666 d-16 /
data cot cs( 9) / -.1481665746 4674657885 2176794666 666 d-18 /
data cot cs( 10) / -.5833356589 0366657947 7984000000 000 d-21 /
data cot cs( 11) / -.2296626469 6464577392 8533333333 333 d-23 /
data cot cs( 12) / -.9041970573 0748332671 9999999999 999 d-26 /
data cot cs( 13) / -.3559885519 2060006400 0000000000 000 d-28 /
data cot cs( 14) / -.1401551398 2429866666 6666666666 666 d-30 /
data cot cs( 15) / -.5518004368 7253333333 3333333333 333 d-33 /
c
c pi2rec = 2/pi - 0.625
data pi2rec / .01161977236 7581343075 5350534900 57 d0 /
data nterms, xmax, xsml, xmin, sqeps /0, 4*0.d0 /
c
if (nterms.ne.0) go to 10
nterms = initds (cotcs, 15, 0.1*sngl(d1mach(3)) )
xmax = 1.0d0/d1mach(4)
xsml = dsqrt (3.0d0*d1mach(3))
xmin = dexp (dmax1(dlog(d1mach(1)), -dlog(d1mach(2))) + 0.01d0)
sqeps = dsqrt (d1mach(4))
c
10 y = dabs(x)
if (y.lt.xmin) call seteru (
1 50hdcot dabs(x) is zero or so small dcot overflows, 50, 2, 2)
if (y.gt.xmax) call seteru (
1 43hdcot no precision because dabs(x) is big, 43, 3, 2)
c
c carefully compute y * (2/pi) = (aint(y) + rem(y)) * (.625 + pi2rec)
c = aint(.625*y) + rem(.625*y) + y*pi2rec = aint(.625*y) + z
c = aint(.625*y) + aint(z) + rem(z)
c
ainty = dint (y)
yrem = y - ainty
prodbg = 0.625d0*ainty
ainty = dint (prodbg)
y = (prodbg-ainty) + 0.625d0*yrem + pi2rec*y
ainty2 = dint (y)
ainty = ainty + ainty2
y = y - ainty2
c
ifn = dmod (ainty, 2.0d0)
if (ifn.eq.1) y = 1.0d0 - y
c
if (dabs(x).gt.0.5d0 .and. y.lt.dabs(x)*sqeps) call seteru (
1 72hdcot answer lt half precision, abs(x) too big or x near n*
2pi (n.ne.0), 72, 1, 1)
c
if (y.gt.0.25d0) go to 20
if (y.eq.0.0d0) call seteru (29hdcot x is a multiple of pi,
1 29, 4, 2)
dcot = 1.0d0/y
if (y.gt.xsml) dcot = (0.5d0 + dcsevl (32.0d0*y*y-1.d0, cotcs,
1 nterms)) / y
go to 40
c
20 if (y.gt.0.5d0) go to 30
dcot = (0.5d0 + dcsevl (8.d0*y*y-1.d0, cotcs, nterms))/(0.5d0*y)
dcot = (dcot*dcot-1.d0)*0.5d0/dcot
go to 40
c
30 dcot = (0.5d0 + dcsevl (2.d0*y*y-1.d0, cotcs, nterms))/(.25d0*y)
dcot = (dcot*dcot-1.d0)*0.5d0/dcot
dcot = (dcot*dcot-1.d0)*0.5d0/dcot
c
40 if (x.ne.0.d0) dcot = dsign (dcot, x)
if (ifn.eq.1) dcot = -dcot
c
return
end