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