datanh.f

Go to the documentation of this file.

*DECK DATANH
      DOUBLE PRECISION FUNCTION datanh (X)
C***BEGIN PROLOGUE  DATANH
C***PURPOSE  Compute the arc hyperbolic tangent.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C4C
C***TYPE      DOUBLE PRECISION (ATANH-S, DATANH-D, CATANH-C)
C***KEYWORDS  ARC HYPERBOLIC TANGENT, ATANH, ELEMENTARY FUNCTIONS,
C             FNLIB, INVERSE HYPERBOLIC TANGENT
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C DATANH(X) calculates the double precision arc hyperbolic
C tangent for double precision argument X.
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
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D1MACH, DCSEVL, INITDS, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C***END PROLOGUE  DATANH
      DOUBLE PRECISION x, atnhcs(27), dxrel, sqeps, y, dcsevl, d1mach
      LOGICAL first
      SAVE atnhcs, nterms, dxrel, sqeps, first
      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     /
      DATA first /.true./
C***FIRST EXECUTABLE STATEMENT  DATANH
      IF (first) THEN
         nterms = initds(atnhcs, 27, 0.1*REAL(D1MACH(3)) )
         dxrel = sqrt(d1mach(4))
         sqeps = sqrt(3.0d0*d1mach(3))
      ENDIF
      first = .false.
C
      y = abs(x)
      IF (y .GE. 1.d0) THEN
         IF (y .GT. 1.d0) THEN 
            datanh = (x - x) / (x - x)
         ELSE
            datanh = x / 0.d0
         ENDIF
         RETURN
      ENDIF
C
      IF (1.d0-y .LT. dxrel) CALL xermsg('SLATEC', 'DATANH',
     +   'ANSWER LT HALF PRECISION BECAUSE ABS(X) TOO NEAR 1', 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*log((1.0d0+x)/(1.0d0-x))
C
      RETURN
      END