dexint.f

      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