DEXINT(3)

SYNOPSIS
SUBROUTINE DEXINT (X, N, KODE, M, TOL, EN, NZ, IERR)

NAME
dexint - Compute an M member sequence of exponential integrals E(N+K,X), K=0,1,...,M-1 for N .GE. 1 and X .GE. 0.

LIBRARY
SLATEC

CATEGORY
C5

TYPE
DOUBLE PRECISION (EXINT-S, DEXINT-D)

KEYWORDS
EXPONENTIAL INTEGRAL, SPECIAL FUNCTIONS

AUTHOR
Amos, D. E., (SNLA)

DESCRIPTION
DEXINT computes M member sequences of exponential integrals E(N+K,X), K=0,1,...,M-1 for N .GE. 1 and X .GE. 0. The exponential integral is defined by E(N,X)=integral on (1,infinity) of EXP(-XT)/T**N where X=0.0 and N=1 cannot occur simultaneously. Formulas and notation are found in the NBS Handbook of Mathematical Functions (ref. 1). The power series is implemented for X .LE. XCUT and the confluent hypergeometric representation E(A,X) = EXP(-X)*(X**(A-1))*U(A,A,X) is computed for X .GT. XCUT. Since sequences are computed in a stable fashion by recurring away from X, A is selected as the integer closest to X within the constraint N .LE. A .LE. N+M-1. For the U computation, A is further modified to be the nearest even integer. Indices are carried forward or backward by the two term recursion relation K*E(K+1,X) + X*E(K,X) = EXP(-X) once E(A,X) is computed. The U function is computed by means of the backward recursive Miller algorithm applied to the three term contiguous relation for U(A+K,A,X), K=0,1,... This produces accurate ratios and determines U(A+K,A,X), and hence E(A,X), to within a multiplicative constant C. Another contiguous relation applied to C*U(A,A,X) and C*U(A+1,A,X) gets C*U(A+1,A+1,X), a quantity proportional to E(A+1,X). The normalizing constant C is obtained from the two term recursion relation above with K=A. The maximum number of significant digits obtainable is the smaller of 14 and the number of digits carried in double precision arithmetic. Description of Arguments Input * X and TOL are double precision * X X .GT. 0.0 for N=1 and X .GE. 0.0 for N .GE. 2 N order of the first member of the sequence, N .GE. 1 (X=0.0 and N=1 is an error) KODE a selection parameter for scaled values KODE=1 returns E(N+K,X), K=0,1,...,M-1. =2 returns EXP(X)*E(N+K,X), K=0,1,...,M-1. M number of exponential integrals in the sequence, M .GE. 1 TOL relative accuracy wanted, ETOL .LE. TOL .LE. 0.1 ETOL is the larger of double precision unit roundoff = D1MACH(4) and 1.0D-18 Output * EN is a double precision vector * EN a vector of dimension at least M containing values EN(K) = E(N+K-1,X) or EXP(X)*E(N+K-1,X), K=1,M depending on KODE NZ underflow indicator NZ=0 a normal return NZ=M X exceeds XLIM and an underflow occurs. EN(K)=0.0D0 , K=1,M returned on KODE=1 IERR error flag IERR=0, normal return, computation completed IERR=1, input error, no computation IERR=2, error, no computation algorithm termination condition not met

REFERENCES
M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, NBS AMS Series 55, U.S. Dept. of Commerce, 1955. D. E. Amos, Computation of exponential integrals, ACM Transactions on Mathematical Software 6, (1980), pp. 365-377 and pp. 420-428.

ROUTINES CALLED
D1MACH, DPSIXN, I1MACH

REVISION HISTORY (YYMMDD)
800501 DATE WRITTEN
890531 Changed all specific intrinsics to generic. (WRB)
890531 REVISION DATE from Version 3.2
891214 Prologue converted to Version 4.0 format. (BAB)
900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
900326 Removed duplicate information from DESCRIPTION section. (WRB)
910408 Updated the REFERENCES section. (WRB)
920207 Updated with code with a revision date of 880811 from D. Amos. Included correction of argument list. (WRB)
920501 Reformatted the REFERENCES section. (WRB)