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