Algorithm 556: Exponential Integrals [S13]
Donald E. Amos · ACM Transactions on Mathematical Software · 1980
The Fortran subroutine EXPINT given here is an implementation of [1].EX-PINT has four machine-dependent parameters XCUT, XLIM, ETOL, and EU-LER which are set into DATA statements.XCUT is a breakpoint such that for x _ XCUT, the series is evaluated and for x > XCUT, the Miller algorithm is applied.XLIM is the approximate underflow limit for e -=, x _ 0, and ETOL is nominally set to 1.E-D where D is the number of base 10 digits in a word.EULER is the negative Euler constant, -7-Thus DATA XCUT, XLIM, ETOL/2.0,667.0, 1.E-12/ would be appropriate for CDC single-precision arithmetic, while DATA XCUT, XLIM, ETOL/1.0,172.0, 1.E-6/ would be appropriate for IBM single-precision arithmetic.The two choices for XCUT reflect the fact that there is a loss of up to two digits on 1 < x _< 2 with the series evaluation.This loss can be tolerated on longer word length machines, but not on shorter word length machines.Maximum accuracy can always be achieved with XCUT = 1.However, for longer word lengths where D = 14, the reduction in computation by moving XCUT from 1 to 2, achieved at the expense of two digits of accuracy, seems to be a worthwhile trade-off, and D ffi 12 reflects this modification for CDC machines.D = 12 is also more consistent with the accuracy attainable from E X P ( -X ) near the underflow limit X = 667.While the subroutine EXPINT is almost portable, the function DIGAM, which computes the psi function at integer arguments, is supplied as a CDC 6600-7600 Fortran function.Modifications for other machines can be made easily from eq. ( 7) of [1].