double precision function dtanh (x)
c june 1977 edition.   w. fullerton, c3, los alamos scientific lab.
      double precision x, tanhcs(31), sqeps, xmax, y, yrec,
     1  dcsevl, d1mach, dexp, dlog, dsqrt
      external d1mach, dcsevl, dexp, dlog, dsqrt, initds
c
c series for tanh       on the interval  0.          to  1.00000e+00
c                                        with weighted error   9.92e-33
c                                         log weighted error  32.00
c                               significant figures required  31.25
c                                    decimal places required  32.75
c
      data tanhcs(  1) / -.2582875664 3634710438 3381514506 05 d+0     /
      data tanhcs(  2) / -.1183610633 0053496535 3836719402 04 d+0     /
      data tanhcs(  3) / +.9869442648 0063988762 8273079996 81 d-2     /
      data tanhcs(  4) / -.8357986623 4458257836 1636903986 38 d-3     /
      data tanhcs(  5) / +.7090432119 8943582626 7780343634 13 d-4     /
      data tanhcs(  6) / -.6016424318 1207040390 7434790010 10 d-5     /
      data tanhcs(  7) / +.5105241908 0064402965 1362977234 11 d-6     /
      data tanhcs(  8) / -.4332072907 7584087216 5454673871 92 d-7     /
      data tanhcs(  9) / +.3675999055 3445306144 9300762337 14 d-8     /
      data tanhcs( 10) / -.3119284961 2492011117 2156514809 53 d-9     /
      data tanhcs( 11) / +.2646882819 9718962579 3777584453 81 d-10    /
      data tanhcs( 12) / -.2246023930 7504140621 8709970061 96 d-11    /
      data tanhcs( 13) / +.1905873376 8288196054 3194683961 39 d-12    /
      data tanhcs( 14) / -.1617237144 6432292391 3307692797 01 d-13    /
      data tanhcs( 15) / +.1372313614 2294289632 8977612893 86 d-14    /
      data tanhcs( 16) / -.1164482687 0554194634 4396472937 81 d-15    /
      data tanhcs( 17) / +.9881268497 1669738285 5405143381 33 d-17    /
      data tanhcs( 18) / -.8384793367 7744865122 2692290559 99 d-18    /
      data tanhcs( 19) / +.7114952886 9124351310 7235061760 00 d-19    /
      data tanhcs( 20) / -.6037424222 9442045413 2888371199 99 d-20    /
      data tanhcs( 21) / +.5123082587 7768084883 4046634666 66 d-21    /
      data tanhcs( 22) / -.4347214015 7782110106 0478293333 33 d-22    /
      data tanhcs( 23) / +.3688847363 9031328479 4231466666 66 d-23    /
      data tanhcs( 24) / -.3130187477 4939399883 3254399999 99 d-24    /
      data tanhcs( 25) / +.2656134200 6551994468 4885333333 33 d-25    /
      data tanhcs( 26) / -.2253874230 4145029883 4943999999 99 d-26    /
      data tanhcs( 27) / +.1912534782 7973995102 2080000000 00 d-27    /
      data tanhcs( 28) / -.1622889709 6543663117 6533333333 33 d-28    /
      data tanhcs( 29) / +.1377110122 9854738786 9866666666 66 d-29    /
      data tanhcs( 30) / -.1168552784 0188950118 3999999999 99 d-30    /
      data tanhcs( 31) / +.9915805538 4640389120 0000000000 00 d-32    /
c
      data nterms, sqeps, xmax / 0, 2*0.d0 /
c
      if (nterms.ne.0) go to 10
      nterms = initds (tanhcs, 31, 0.1*sngl(d1mach(3)) )
      sqeps = dsqrt (3.0d0*d1mach(3))
      xmax = -0.5d0*dlog (d1mach(3))
c
 10   y = dabs(x)
      if (y.gt.1.d0) go to 20
c
      dtanh = x
      if (y.gt.sqeps) dtanh = x*(1.d0 + dcsevl(2.d0*x*x-1.d0, tanhcs,
     1  nterms) )
      return
c
 20   if (y.gt.xmax) go to 30
      y = dexp(y)
      yrec = 1.d0/y
      dtanh = dsign ((y-yrec)/(y+yrec), x)
      return
c
 30   dtanh = dsign (1.d0, x)
      return
      end