r9atn1.f
2.39 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
function r9atn1 (x)
c april 1978 edition. w. fullerton c3, los alamos scientific lab.
c
c evaluate atan(x) from first order, that is, evaluate
c (atan(x)-x)/x**3 with relative error accuracy so that
c atan(x) = x + x**3*r9atn1(x).
c
dimension atn1cs(21)
external csevl, inits, r1mach
c
c series for atn1 on the interval 0. to 1.00000d+00
c with weighted error 2.21e-17
c log weighted error 16.66
c significant figures required 15.44
c decimal places required 17.32
c
data atn1cs( 1) / -.0328399753 5355202e0 /
data atn1cs( 2) / .0583343234 3172412e0 /
data atn1cs( 3) / -.0074003696 9671964e0 /
data atn1cs( 4) / .0010097841 9933728e0 /
data atn1cs( 5) / -.0001439787 1635652e0 /
data atn1cs( 6) / .0000211451 2648992e0 /
data atn1cs( 7) / -.0000031723 2107425e0 /
data atn1cs( 8) / .0000004836 6203654e0 /
data atn1cs( 9) / -.0000000746 7746546e0 /
data atn1cs(10) / .0000000116 4800896e0 /
data atn1cs(11) / -.0000000018 3208837e0 /
data atn1cs(12) / .0000000002 9019082e0 /
data atn1cs(13) / -.0000000000 4623885e0 /
data atn1cs(14) / .0000000000 0740552e0 /
data atn1cs(15) / -.0000000000 0119135e0 /
data atn1cs(16) / .0000000000 0019240e0 /
data atn1cs(17) / -.0000000000 0003118e0 /
data atn1cs(18) / .0000000000 0000506e0 /
data atn1cs(19) / -.0000000000 0000082e0 /
data atn1cs(20) / .0000000000 0000013e0 /
data atn1cs(21) / -.0000000000 0000002e0 /
c
data ntatn1, xsml, xbig, xmax / 0, 3*0.0 /
c
if (ntatn1.ne.0) go to 10
eps = r1mach(3)
ntatn1 = inits (atn1cs, 21, 0.1*eps)
c
xsml = sqrt (0.1*eps)
xbig = 1.571/sqrt(eps)
xmax = 1.571/eps
c
10 y = abs(x)
if (y.gt.1.0) go to 20
c
if (y.le.xsml) r9atn1 = -1.0/3.0
if (y.le.xsml) return
c
r9atn1 = -0.25 + csevl (2.0*y*y-1., atn1cs, ntatn1)
return
c
20 if (y.gt.xmax) call seteru (
1 51hr9atn1 no precision in answer because x is too big, 51,2,2)
if (y.gt.xbig) call seteru (
1 53hr9atn1 answer lt half precision because x is too big, 53,
2 1, 1)
c
r9atn1 = (atan(x) - x) / x**3
return
c
end