double precision function  dtanh (x)
*     (Hyperbolic Tangent)
*     Implement the hyperbolic tangent using the Cody-Waite algoritm.
*     The computation is divided into 4 regions:
*
*     0 <= x < XSMALL                 tanh(x) = x
*     XSMALL <= x < XMEDIUM           tanh(x) = rational polynomial
*     XMEDIUM <= x < XLARGE           tanh(x) = 1 - 2/(1 + exp(2x))
*     XLARGE <= x <= infinity         tanh(x) = 1
*
*     (07-Feb-1991)

*     In this implementation, t = 53, B = 2.

*     XLARGE = (ln(2) + (t + 1)ln(B))/2
*          = smallest value for which tanh(x) = 1

      double precision XLARGE
      parameter (XLARGE = 19.06154 74653 98496 00897 D+00)

*     XMED = ln(3)/2
*          = cutoff for second approximation

      double precision XMED
      parameter (XMED = 0.54930 61443 34054 84570 D+00)

*     XSMALL = sqrt(3) * B**((-t-1)/2)
*            = cutoff for third approximation

      double precision XSMALL
      parameter (XSMALL = 1.29047 84139 75892 43466D-08)

      double precision f, g, r, result, x
      double precision p(0:2), q(0:3)

      data p / -0.16134 11902 39962 28053 D+04,
     x         -0.99225 92967 22360 83313 D+02,
     x         -0.96437 49277 72254 69787 D+00 /
      data q /  0.48402 35707 19886 88686 D+04,
     x          0.22337 72071 89623 12926 D+04,
     x          0.11274 47438 05349 49335 D+03,
     x          1.00000 00000 00000 00000 D+00 /

      f = dabs(x)
      if (f .ne. f) then
*        x is a NaN, so tanh is too--generate a run-time trap
         result = f + f
      else if (f .ge. XLARGE) then
         result = 1.0D+00
      else if (f .ge. XMED) then
         result = 0.5D+00 - 1.0D+00 / (dexp(f + f) + 1.0D+00)
         result = result + result
      else if (f .ge. XSMALL) then
         g = f * f
         r = ((p(2) * g + p(1)) * g + p(0)) * g /
     X     (((g + q(2))*g + q(1)) * g + q(0))
         result = f + f * r
      else
         result = f
      end if
      if (x .lt. 0.0D+00) result = -result
      dtanh = result
      end