# exint.f

          SUBROUTINE EXINT (X, N, KODE, M, TOL, EN, NZ, IERR)
    C***BEGIN PROLOGUE  EXINT
    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      SINGLE PRECISION (EXINT-S, DEXINT-D)
    C***KEYWORDS  EXPONENTIAL INTEGRAL, SPECIAL FUNCTIONS
    C***AUTHOR  Amos, D. E., (SNLA)
    C***DESCRIPTION
    C
    C         EXINT 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     Description of Arguments
    C
    C         Input
    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 = single precision unit roundoff = R1MACH(4)
    C
    C         Output
    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.0E0 , 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  I1MACH, PSIXN, R1MACH
    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  EXINT