datanh.f
3.13 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
double precision function datanh (x)
c june 1977 edition. w. fullerton, c3, los alamos scientific lab.
double precision x, atnhcs(27), dxrel, sqeps, y, dcsevl,
1 d1mach, dlog, dsqrt
external d1mach, dcsevl, initds
c
c series for atnh on the interval 0. to 2.50000e-01
c with weighted error 6.86e-32
c log weighted error 31.16
c significant figures required 30.00
c decimal places required 31.88
c
data atnhcs( 1) / +.9439510239 3195492308 4289221863 3 d-1 /
data atnhcs( 2) / +.4919843705 5786159472 0003457666 8 d-1 /
data atnhcs( 3) / +.2102593522 4554327634 7932733175 2 d-2 /
data atnhcs( 4) / +.1073554449 7761165846 4073104527 6 d-3 /
data atnhcs( 5) / +.5978267249 2930314786 4278751787 2 d-5 /
data atnhcs( 6) / +.3505062030 8891348459 6683488620 0 d-6 /
data atnhcs( 7) / +.2126374343 7653403508 9621931443 1 d-7 /
data atnhcs( 8) / +.1321694535 7155271921 2980172305 5 d-8 /
data atnhcs( 9) / +.8365875501 1780703646 2360405295 9 d-10 /
data atnhcs( 10) / +.5370503749 3110021638 8143458777 2 d-11 /
data atnhcs( 11) / +.3486659470 1571079229 7124578429 0 d-12 /
data atnhcs( 12) / +.2284549509 6034330155 2402411972 2 d-13 /
data atnhcs( 13) / +.1508407105 9447930448 7422906755 8 d-14 /
data atnhcs( 14) / +.1002418816 8041091261 3699572283 7 d-15 /
data atnhcs( 15) / +.6698674738 1650695397 1552688298 6 d-17 /
data atnhcs( 16) / +.4497954546 4949310830 8332762453 3 d-18 /
data atnhcs( 17) / +.3032954474 2794535416 8236714666 6 d-19 /
data atnhcs( 18) / +.2052702064 1909368264 6386141866 6 d-20 /
data atnhcs( 19) / +.1393848977 0538377131 9301461333 3 d-21 /
data atnhcs( 20) / +.9492580637 2245769719 5895466666 6 d-23 /
data atnhcs( 21) / +.6481915448 2423076049 8244266666 6 d-24 /
data atnhcs( 22) / +.4436730205 7236152726 3232000000 0 d-25 /
data atnhcs( 23) / +.3043465618 5431616389 1200000000 0 d-26 /
data atnhcs( 24) / +.2091881298 7923934740 4799999999 9 d-27 /
data atnhcs( 25) / +.1440445411 2340505613 6533333333 3 d-28 /
data atnhcs( 26) / +.9935374683 1416404650 6666666666 6 d-30 /
data atnhcs( 27) / +.6863462444 3582600533 3333333333 3 d-31 /
c
data nterms, dxrel, sqeps / 0, 2*0.d0 /
c
if (nterms.ne.0) go to 10
nterms = initds (atnhcs, 27, 0.1*sngl(d1mach(3)) )
dxrel = dsqrt (d1mach(4))
sqeps = dsqrt (3.0d0*d1mach(3))
c
10 y = dabs(x)
if (y.ge.1.d0) call seteru (20hdatanh dabs(x) ge 1, 20, 2, 2)
c
if (1.d0-y.lt.dxrel) call seteru (
1 59hdatanh answer lt half precision because dabs(x) too near 1,
2 59, 1, 1)
c
datanh = x
if (y.gt.sqeps .and. y.le.0.5d0) datanh = x*(1.0d0 +
1 dcsevl (8.d0*x*x-1.d0, atnhcs, nterms) )
if (y.gt.0.5d0) datanh = 0.5d0*dlog ((1.0d0+x)/(1.0d0-x) )
c
return
end