Blame view

fvn_fnlib/datanh.f 3.13 KB
38581db0c   daniau   git-svn-id: https...
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