cot.f
2.69 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
function cot (x)
c may 1980 edition. w. fullerton, c3, los alamos scientific lab.
dimension cotcs(8)
external csevl, inits, r1mach
c
c series for cot on the interval 0. to 6.25000d-02
c with weighted error 3.76e-17
c log weighted error 16.42
c significant figures required 15.51
c decimal places required 16.88
c
data cot cs( 1) / .2402591609 8295630e0 /
data cot cs( 2) / -.0165330316 01500228e0 /
data cot cs( 3) / -.0000429983 91931724e0 /
data cot cs( 4) / -.0000001592 83223327e0 /
data cot cs( 5) / -.0000000006 19109313e0 /
data cot cs( 6) / -.0000000000 02430197e0 /
data cot cs( 7) / -.0000000000 00009560e0 /
data cot cs( 8) / -.0000000000 00000037e0 /
c
c pi2rec = 2/pi - 0.625
data pi2rec / .01161977236 75813430 e0 /
data nterms, xmax, xsml, xmin, sqeps / 0, 4*0.0 /
c
if (nterms.ne.0) go to 10
nterms = inits (cotcs, 8, 0.1*r1mach(3))
xmax = 1.0/r1mach(4)
xsml = sqrt (3.0*r1mach(3))
xmin = exp ( amax1(alog(r1mach(1)), -alog(r1mach(2))) + 0.01)
sqeps = sqrt (r1mach(4))
c
10 y = abs(x)
if (abs(x).lt.xmin) call seteru (
1 48hcot abs(x) is zero or so small cot overflows, 48, 2, 2)
if (y.gt.xmax) call seteru (
1 42hcot no precision because abs(x) is big, 42, 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 = aint (y)
yrem = y - ainty
prodbg = 0.625*ainty
ainty = aint (prodbg)
y = (prodbg-ainty) + 0.625*yrem + y*pi2rec
ainty2 = aint (y)
ainty = ainty + ainty2
y = y - ainty2
c
ifn = amod (ainty, 2.)
if (ifn.eq.1) y = 1.0 - y
c
if (abs(x).gt.0.5 .and. y.lt.abs(x)*sqeps) call seteru (
1 72hcot answer lt half precision, abs(x) too big or x near n*
2pi (n.ne.0), 72, 1, 1)
c
if (y.gt.0.25) go to 20
if (y.eq.0.0) call seteru (29hcot x is a multiple of pi,
1 29, 4, 2)
cot = 1.0/y
if (y.gt.xsml) cot = (0.5 + csevl (32.0*y*y-1., cotcs, nterms)) /y
go to 40
c
20 if (y.gt.0.5) go to 30
cot = (0.5 + csevl (8.0*y*y-1., cotcs, nterms)) / (0.5*y)
cot = (cot**2 - 1.0) * 0.5 / cot
go to 40
c
30 cot = (0.5 + csevl (2.0*y*y-1., cotcs, nterms)) / (0.25*y)
cot = (cot**2 - 1.0) * 0.5 / cot
cot = (cot**2 - 1.0) * 0.5 / cot
c
40 if (x.ne.0.) cot = sign (cot, x)
if (ifn.eq.1) cot = -cot
c
return
end