atanh.f
1.84 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
function atanh (x)
c april 1977 edition. w. fullerton, c3, los alamos scientific lab.
dimension atnhcs(15)
external csevl, inits, r1mach
c
c series for atnh on the interval 0. to 2.50000d-01
c with weighted error 6.70e-18
c log weighted error 17.17
c significant figures required 16.01
c decimal places required 17.76
c
data atnhcs( 1) / .0943951023 93195492e0 /
data atnhcs( 2) / .0491984370 55786159e0 /
data atnhcs( 3) / .0021025935 22455432e0 /
data atnhcs( 4) / .0001073554 44977611e0 /
data atnhcs( 5) / .0000059782 67249293e0 /
data atnhcs( 6) / .0000003505 06203088e0 /
data atnhcs( 7) / .0000000212 63743437e0 /
data atnhcs( 8) / .0000000013 21694535e0 /
data atnhcs( 9) / .0000000000 83658755e0 /
data atnhcs(10) / .0000000000 05370503e0 /
data atnhcs(11) / .0000000000 00348665e0 /
data atnhcs(12) / .0000000000 00022845e0 /
data atnhcs(13) / .0000000000 00001508e0 /
data atnhcs(14) / .0000000000 00000100e0 /
data atnhcs(15) / .0000000000 00000006e0 /
c
data nterms, dxrel, sqeps /0, 0., 0./
c
if (nterms.ne.0) go to 10
nterms = inits (atnhcs, 15, 0.1*r1mach(3))
dxrel = sqrt (r1mach(4))
sqeps = sqrt (3.0*r1mach(3))
c
10 y = abs(x)
if (y.ge.1.0) call seteru (19hatanh abs(x) ge 1, 19, 2, 2)
c
if (1.0-y.lt.dxrel) call seteru ( 58hatanh answer lt half precis
1ion because abs(x) too near 1, 58, 1, 1)
c
atanh = x
if (y.gt.sqeps .and. y.le.0.5) atanh = x*(1.0 + csevl (8.*x*x-1.,
1 atnhcs, nterms))
if (y.gt.0.5) atanh = 0.5*alog((1.0+x)/(1.0-x))
c
return
end