relibc_openlibm/slatec/dexint.f

*DECK DEXINT
      SUBROUTINE DEXINT (X, N, KODE, M, TOL, EN, NZ, IERR)
C***BEGIN PROLOGUE  DEXINT
C***PURPOSE  Compute an M member sequence of exponential integrals
C            E(N+K,X), K=0,1,...,M-1 for N .GE. 1 and X .GE. 0.
C***LIBRARY   SLATEC
C***CATEGORY  C5
C***TYPE      DOUBLE PRECISION (EXINT-S, DEXINT-D)
C***KEYWORDS  EXPONENTIAL INTEGRAL, SPECIAL FUNCTIONS
C***AUTHOR  Amos, D. E., (SNLA)
C***DESCRIPTION
C
C         DEXINT computes M member sequences of exponential integrals
C         E(N+K,X), K=0,1,...,M-1 for N .GE. 1 and X .GE. 0.  The
C         exponential integral is defined by
C
C         E(N,X)=integral on (1,infinity) of EXP(-XT)/T**N
C
C         where X=0.0 and N=1 cannot occur simultaneously.  Formulas
C         and notation are found in the NBS Handbook of Mathematical
C         Functions (ref. 1).
C
C         The power series is implemented for X .LE. XCUT and the
C         confluent hypergeometric representation
C
C                     E(A,X) = EXP(-X)*(X**(A-1))*U(A,A,X)
C
C         is computed for X .GT. XCUT.  Since sequences are computed in
C         a stable fashion by recurring away from X, A is selected as
C         the integer closest to X within the constraint N .LE. A .LE.
C         N+M-1.  For the U computation, A is further modified to be the
C         nearest even integer.  Indices are carried forward or
C         backward by the two term recursion relation
C
C                     K*E(K+1,X) + X*E(K,X) = EXP(-X)
C
C         once E(A,X) is computed.  The U function is computed by means
C         of the backward recursive Miller algorithm applied to the
C         three term contiguous relation for U(A+K,A,X), K=0,1,...
C         This produces accurate ratios and determines U(A+K,A,X), and
C         hence E(A,X), to within a multiplicative constant C.
C         Another contiguous relation applied to C*U(A,A,X) and
C         C*U(A+1,A,X) gets C*U(A+1,A+1,X), a quantity proportional to
C         E(A+1,X).  The normalizing constant C is obtained from the
C         two term recursion relation above with K=A.
C
C         The maximum number of significant digits obtainable
C         is the smaller of 14 and the number of digits carried in
C         double precision arithmetic.
C
C     Description of Arguments
C
C         Input     * X and TOL are double precision *
C           X       X .GT. 0.0 for N=1 and  X .GE. 0.0 for N .GE. 2
C           N       order of the first member of the sequence, N .GE. 1
C                   (X=0.0 and N=1 is an error)
C           KODE    a selection parameter for scaled values
C                   KODE=1   returns        E(N+K,X), K=0,1,...,M-1.
C                       =2   returns EXP(X)*E(N+K,X), K=0,1,...,M-1.
C           M       number of exponential integrals in the sequence,
C                   M .GE. 1
C           TOL     relative accuracy wanted, ETOL .LE. TOL .LE. 0.1
C                   ETOL is the larger of double precision unit
C                   roundoff = D1MACH(4) and 1.0D-18
C
C         Output    * EN is a double precision vector *
C           EN      a vector of dimension at least M containing values
C                   EN(K) = E(N+K-1,X) or EXP(X)*E(N+K-1,X), K=1,M
C                   depending on KODE
C           NZ      underflow indicator
C                   NZ=0   a normal return
C                   NZ=M   X exceeds XLIM and an underflow occurs.
C                          EN(K)=0.0D0 , K=1,M returned on KODE=1
C           IERR    error flag
C                   IERR=0, normal return, computation completed
C                   IERR=1, input error,   no computation
C                   IERR=2, error,         no computation
C                           algorithm termination condition not met
C
C***REFERENCES  M. Abramowitz and I. A. Stegun, Handbook of
C                 Mathematical Functions, NBS AMS Series 55, U.S. Dept.
C                 of Commerce, 1955.
C               D. E. Amos, Computation of exponential integrals, ACM
C                 Transactions on Mathematical Software 6, (1980),
C                 pp. 365-377 and pp. 420-428.
C***ROUTINES CALLED  D1MACH, DPSIXN, I1MACH
C***REVISION HISTORY  (YYMMDD)
C   800501  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   900326  Removed duplicate information from DESCRIPTION section.
C           (WRB)
C   910408  Updated the REFERENCES section.  (WRB)
C   920207  Updated with code with a revision date of 880811 from
C           D. Amos.  Included correction of argument list.  (WRB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  DEXINT
      DOUBLE PRECISION A,AA,AAMS,AH,AK,AT,B,BK,BT,CC,CNORM,CT,EM,EMX,EN,
     1                 ETOL,FNM,FX,PT,P1,P2,S,TOL,TX,X,XCUT,XLIM,XTOL,Y,
     2                 YT,Y1,Y2
      DOUBLE PRECISION D1MACH,DPSIXN
      INTEGER I,IC,ICASE,ICT,IERR,IK,IND,IX,I1M,JSET,K,KK,KN,KODE,KS,M,
     1        ML,MU,N,ND,NM,NZ
      INTEGER I1MACH
      DIMENSION EN(*), A(99), B(99), Y(2)
      SAVE XCUT
      DATA XCUT / 2.0D0 /
C***FIRST EXECUTABLE STATEMENT  DEXINT
      IERR = 0
      NZ = 0
      ETOL = MAX(D1MACH(4),0.5D-18)
      IF (X.LT.0.0D0) IERR = 1
      IF (N.LT.1) IERR = 1
      IF (KODE.LT.1 .OR. KODE.GT.2) IERR = 1
      IF (M.LT.1) IERR = 1
      IF (TOL.LT.ETOL .OR. TOL.GT.0.1D0) IERR = 1
      IF (X.EQ.0.0D0 .AND. N.EQ.1) IERR = 1
      IF(IERR.NE.0) RETURN
      I1M = -I1MACH(15)
      PT = 2.3026D0*I1M*D1MACH(5)
      XLIM = PT - 6.907755D0
      BT = PT + (N+M-1)
      IF (BT.GT.1000.0D0) XLIM = PT - LOG(BT)
C
      IF (X.GT.XCUT) GO TO 100
      IF (X.EQ.0.0D0 .AND. N.GT.1) GO TO 80
C-----------------------------------------------------------------------
C     SERIES FOR E(N,X) FOR X.LE.XCUT
C-----------------------------------------------------------------------
      TX = X + 0.5D0
      IX = TX
C-----------------------------------------------------------------------
C     ICASE=1 MEANS INTEGER CLOSEST TO X IS 2 AND N=1
C     ICASE=2 MEANS INTEGER CLOSEST TO X IS 0,1, OR 2 AND N.GE.2
C-----------------------------------------------------------------------
      ICASE = 2
      IF (IX.GT.N) ICASE = 1
      NM = N - ICASE + 1
      ND = NM + 1
      IND = 3 - ICASE
      MU = M - IND
      ML = 1
      KS = ND
      FNM = NM
      S = 0.0D0
      XTOL = 3.0D0*TOL
      IF (ND.EQ.1) GO TO 10
      XTOL = 0.3333D0*TOL
      S = 1.0D0/FNM
   10 CONTINUE
      AA = 1.0D0
      AK = 1.0D0
      IC = 35
      IF (X.LT.ETOL) IC = 1
      DO 50 I=1,IC
        AA = -AA*X/AK
        IF (I.EQ.NM) GO TO 30
        S = S - AA/(AK-FNM)
        IF (ABS(AA).LE.XTOL*ABS(S)) GO TO 20
        AK = AK + 1.0D0
        GO TO 50
   20   CONTINUE
        IF (I.LT.2) GO TO 40
        IF (ND-2.GT.I .OR. I.GT.ND-1) GO TO 60
        AK = AK + 1.0D0
        GO TO 50
   30   S = S + AA*(-LOG(X)+DPSIXN(ND))
        XTOL = 3.0D0*TOL
   40   AK = AK + 1.0D0
   50 CONTINUE
      IF (IC.NE.1) GO TO 340
   60 IF (ND.EQ.1) S = S + (-LOG(X)+DPSIXN(1))
      IF (KODE.EQ.2) S = S*EXP(X)
      EN(1) = S
      EMX = 1.0D0
      IF (M.EQ.1) GO TO 70
      EN(IND) = S
      AA = KS
      IF (KODE.EQ.1) EMX = EXP(-X)
      GO TO (220, 240), ICASE
   70 IF (ICASE.EQ.2) RETURN
      IF (KODE.EQ.1) EMX = EXP(-X)
      EN(1) = (EMX-S)/X
      RETURN
   80 CONTINUE
      DO 90 I=1,M
        EN(I) = 1.0D0/(N+I-2)
   90 CONTINUE
      RETURN
C-----------------------------------------------------------------------
C     BACKWARD RECURSIVE MILLER ALGORITHM FOR
C              E(N,X)=EXP(-X)*(X**(N-1))*U(N,N,X)
C     WITH RECURSION AWAY FROM N=INTEGER CLOSEST TO X.
C     U(A,B,X) IS THE SECOND CONFLUENT HYPERGEOMETRIC FUNCTION
C-----------------------------------------------------------------------
  100 CONTINUE
      EMX = 1.0D0
      IF (KODE.EQ.2) GO TO 130
      IF (X.LE.XLIM) GO TO 120
      NZ = M
      DO 110 I=1,M
        EN(I) = 0.0D0
  110 CONTINUE
      RETURN
  120 EMX = EXP(-X)
  130 CONTINUE
      TX = X + 0.5D0
      IX = TX
      KN = N + M - 1
      IF (KN.LE.IX) GO TO 140
      IF (N.LT.IX .AND. IX.LT.KN) GO TO 170
      IF (N.GE.IX) GO TO 160
      GO TO 340
  140 ICASE = 1
      KS = KN
      ML = M - 1
      MU = -1
      IND = M
      IF (KN.GT.1) GO TO 180
  150 KS = 2
      ICASE = 3
      GO TO 180
  160 ICASE = 2
      IND = 1
      KS = N
      MU = M - 1
      IF (N.GT.1) GO TO 180
      IF (KN.EQ.1) GO TO 150
      IX = 2
  170 ICASE = 1
      KS = IX
      ML = IX - N
      IND = ML + 1
      MU = KN - IX
  180 CONTINUE
      IK = KS/2
      AH = IK
      JSET = 1 + KS - (IK+IK)
C-----------------------------------------------------------------------
C     START COMPUTATION FOR
C              EN(IND) = C*U( A , A ,X)    JSET=1
C              EN(IND) = C*U(A+1,A+1,X)    JSET=2
C     FOR AN EVEN INTEGER A.
C-----------------------------------------------------------------------
      IC = 0
      AA = AH + AH
      AAMS = AA - 1.0D0
      AAMS = AAMS*AAMS
      TX = X + X
      FX = TX + TX
      AK = AH
      XTOL = TOL
      IF (TOL.LE.1.0D-3) XTOL = 20.0D0*TOL
      CT = AAMS + FX*AH
      EM = (AH+1.0D0)/((X+AA)*XTOL*SQRT(CT))
      BK = AA
      CC = AH*AH
C-----------------------------------------------------------------------
C     FORWARD RECURSION FOR P(IC),P(IC+1) AND INDEX IC FOR BACKWARD
C     RECURSION
C-----------------------------------------------------------------------
      P1 = 0.0D0
      P2 = 1.0D0
  190 CONTINUE
      IF (IC.EQ.99) GO TO 340
      IC = IC + 1
      AK = AK + 1.0D0
      AT = BK/(BK+AK+CC+IC)
      BK = BK + AK + AK
      A(IC) = AT
      BT = (AK+AK+X)/(AK+1.0D0)
      B(IC) = BT
      PT = P2
      P2 = BT*P2 - AT*P1
      P1 = PT
      CT = CT + FX
      EM = EM*AT*(1.0D0-TX/CT)
      IF (EM*(AK+1.0D0).GT.P1*P1) GO TO 190
      ICT = IC
      KK = IC + 1
      BT = TX/(CT+FX)
      Y2 = (BK/(BK+CC+KK))*(P1/P2)*(1.0D0-BT+0.375D0*BT*BT)
      Y1 = 1.0D0
C-----------------------------------------------------------------------
C     BACKWARD RECURRENCE FOR
C              Y1=             C*U( A ,A,X)
C              Y2= C*(A/(1+A/2))*U(A+1,A,X)
C-----------------------------------------------------------------------
      DO 200 K=1,ICT
        KK = KK - 1
        YT = Y1
        Y1 = (B(KK)*Y1-Y2)/A(KK)
        Y2 = YT
  200 CONTINUE
C-----------------------------------------------------------------------
C     THE CONTIGUOUS RELATION
C              X*U(B,C+1,X)=(C-B)*U(B,C,X)+U(B-1,C,X)
C     WITH  B=A+1 , C=A IS USED FOR
C              Y(2) = C * U(A+1,A+1,X)
C     X IS INCORPORATED INTO THE NORMALIZING RELATION
C-----------------------------------------------------------------------
      PT = Y2/Y1
      CNORM = 1.0E0 - PT*(AH+1.0E0)/AA
      Y(1) = 1.0E0/(CNORM*AA+X)
      Y(2) = CNORM*Y(1)
      IF (ICASE.EQ.3) GO TO 210
      EN(IND) = EMX*Y(JSET)
      IF (M.EQ.1) RETURN
      AA = KS
      GO TO (220, 240), ICASE
C-----------------------------------------------------------------------
C     RECURSION SECTION  N*E(N+1,X) + X*E(N,X)=EMX
C-----------------------------------------------------------------------
  210 EN(1) = EMX*(1.0E0-Y(1))/X
      RETURN
  220 K = IND - 1
      DO 230 I=1,ML
        AA = AA - 1.0D0
        EN(K) = (EMX-AA*EN(K+1))/X
        K = K - 1
  230 CONTINUE
      IF (MU.LE.0) RETURN
      AA = KS
  240 K = IND
      DO 250 I=1,MU
        EN(K+1) = (EMX-X*EN(K))/AA
        AA = AA + 1.0D0
        K = K + 1
  250 CONTINUE
      RETURN
  340 CONTINUE
      IERR = 2
      RETURN
      END