SLATEC Routines --- DEXINT ---
*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