1*DECK EXINT
2      SUBROUTINE EXINT (X, N, KODE, M, TOL, EN, NZ, IERR)
3C***BEGIN PROLOGUE  EXINT
4C***PURPOSE  Compute an M member sequence of exponential integrals
5C            E(N+K,X), K=0,1,...,M-1 for N .GE. 1 and X .GE. 0.
6C***LIBRARY   SLATEC
7C***CATEGORY  C5
8C***TYPE      SINGLE PRECISION (EXINT-S, DEXINT-D)
9C***KEYWORDS  EXPONENTIAL INTEGRAL, SPECIAL FUNCTIONS
10C***AUTHOR  Amos, D. E., (SNLA)
11C***DESCRIPTION
12C
13C         EXINT computes M member sequences of exponential integrals
14C         E(N+K,X), K=0,1,...,M-1 for N .GE. 1 and X .GE. 0.  The
15C         exponential integral is defined by
16C
17C         E(N,X)=integral on (1,infinity) of EXP(-XT)/T**N
18C
19C         where X=0.0 and N=1 cannot occur simultaneously.  Formulas
20C         and notation are found in the NBS Handbook of Mathematical
21C         Functions (ref. 1).
22C
23C         The power series is implemented for X .LE. XCUT and the
24C         confluent hypergeometric representation
25C
26C                     E(A,X) = EXP(-X)*(X**(A-1))*U(A,A,X)
27C
28C         is computed for X .GT. XCUT.  Since sequences are computed in
29C         a stable fashion by recurring away from X, A is selected as
30C         the integer closest to X within the constraint N .LE. A .LE.
31C         N+M-1.  For the U computation, A is further modified to be the
32C         nearest even integer.  Indices are carried forward or
33C         backward by the two term recursion relation
34C
35C                     K*E(K+1,X) + X*E(K,X) = EXP(-X)
36C
37C         once E(A,X) is computed.  The U function is computed by means
38C         of the backward recursive Miller algorithm applied to the
39C         three term contiguous relation for U(A+K,A,X), K=0,1,...
40C         This produces accurate ratios and determines U(A+K,A,X), and
41C         hence E(A,X), to within a multiplicative constant C.
42C         Another contiguous relation applied to C*U(A,A,X) and
43C         C*U(A+1,A,X) gets C*U(A+1,A+1,X), a quantity proportional to
44C         E(A+1,X).  The normalizing constant C is obtained from the
45C         two term recursion relation above with K=A.
46C
47C     Description of Arguments
48C
49C         Input
50C           X       X .GT. 0.0 for N=1 and  X .GE. 0.0 for N .GE. 2
51C           N       order of the first member of the sequence, N .GE. 1
52C                   (X=0.0 and N=1 is an error)
53C           KODE    a selection parameter for scaled values
54C                   KODE=1   returns        E(N+K,X), K=0,1,...,M-1.
55C                       =2   returns EXP(X)*E(N+K,X), K=0,1,...,M-1.
56C           M       number of exponential integrals in the sequence,
57C                   M .GE. 1
58C           TOL     relative accuracy wanted, ETOL .LE. TOL .LE. 0.1
59C                   ETOL = single precision unit roundoff = R1MACH(4)
60C
61C         Output
62C           EN      a vector of dimension at least M containing values
63C                   EN(K) = E(N+K-1,X) or EXP(X)*E(N+K-1,X), K=1,M
64C                   depending on KODE
65C           NZ      underflow indicator
66C                   NZ=0   a normal return
67C                   NZ=M   X exceeds XLIM and an underflow occurs.
68C                          EN(K)=0.0E0 , K=1,M returned on KODE=1
69C           IERR    error flag
70C                   IERR=0, normal return, computation completed
71C                   IERR=1, input error,   no computation
72C                   IERR=2, error,         no computation
73C                           algorithm termination condition not met
74C
75C***REFERENCES  M. Abramowitz and I. A. Stegun, Handbook of
76C                 Mathematical Functions, NBS AMS Series 55, U.S. Dept.
77C                 of Commerce, 1955.
78C               D. E. Amos, Computation of exponential integrals, ACM
79C                 Transactions on Mathematical Software 6, (1980),
80C                 pp. 365-377 and pp. 420-428.
81C***ROUTINES CALLED  I1MACH, PSIXN, R1MACH
82C***REVISION HISTORY  (YYMMDD)
83C   800501  DATE WRITTEN
84C   890531  Changed all specific intrinsics to generic.  (WRB)
85C   890531  REVISION DATE from Version 3.2
86C   891214  Prologue converted to Version 4.0 format.  (BAB)
87C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
88C   900326  Removed duplicate information from DESCRIPTION section.
89C           (WRB)
90C   910408  Updated the REFERENCES section.  (WRB)
91C   920207  Updated with code with a revision date of 880811 from
92C           D. Amos.  Included correction of argument list.  (WRB)
93C   920501  Reformatted the REFERENCES section.  (WRB)
94C***END PROLOGUE  EXINT
95      REAL             A,AA,AAMS,AH,AK,AT,B,BK,BT,CC,CNORM,CT,EM,EMX,EN,
96     1                 ETOL,FNM,FX,PT,P1,P2,S,TOL,TX,X,XCUT,XLIM,XTOL,Y,
97     2                 YT,Y1,Y2
98      REAL             R1MACH,PSIXN
99      INTEGER I,IC,ICASE,ICT,IERR,IK,IND,IX,I1M,JSET,K,KK,KN,KODE,KS,M,
100     1        ML,MU,N,ND,NM,NZ
101      INTEGER I1MACH
102      DIMENSION EN(*), A(99), B(99), Y(2)
103C***FIRST EXECUTABLE STATEMENT  EXINT
104      IERR = 0
105      NZ = 0
106      ETOL = MAX(R1MACH(4),0.5E-18)
107      IF (X.LT.0.0E0) IERR = 1
108      IF (N.LT.1) IERR = 1
109      IF (KODE.LT.1 .OR. KODE.GT.2) IERR = 1
110      IF (M.LT.1) IERR = 1
111      IF (TOL.LT.ETOL .OR. TOL.GT.0.1E0) IERR = 1
112      IF (X.EQ.0.0E0 .AND. N.EQ.1) IERR = 1
113      IF (IERR.NE.0) RETURN
114      I1M = -I1MACH(12)
115      PT = 2.3026E0*R1MACH(5)*I1M
116      XLIM = PT - 6.907755E0
117      BT = PT + (N+M-1)
118      IF (BT.GT.1000.0E0) XLIM = PT - LOG(BT)
119C
120      XCUT = 2.0E0
121      IF (ETOL.GT.2.0E-7) XCUT = 1.0E0
122      IF (X.GT.XCUT) GO TO 100
123      IF (X.EQ.0.0E0 .AND. N.GT.1) GO TO 80
124C-----------------------------------------------------------------------
125C     SERIES FOR E(N,X) FOR X.LE.XCUT
126C-----------------------------------------------------------------------
127      TX = X + 0.5E0
128      IX = TX
129C-----------------------------------------------------------------------
130C     ICASE=1 MEANS INTEGER CLOSEST TO X IS 2 AND N=1
131C     ICASE=2 MEANS INTEGER CLOSEST TO X IS 0,1, OR 2 AND N.GE.2
132C-----------------------------------------------------------------------
133      ICASE = 2
134      IF (IX.GT.N) ICASE = 1
135      NM = N - ICASE + 1
136      ND = NM + 1
137      IND = 3 - ICASE
138      MU = M - IND
139      ML = 1
140      KS = ND
141      FNM = NM
142      S = 0.0E0
143      XTOL = 3.0E0*TOL
144      IF (ND.EQ.1) GO TO 10
145      XTOL = 0.3333E0*TOL
146      S = 1.0E0/FNM
147   10 CONTINUE
148      AA = 1.0E0
149      AK = 1.0E0
150      IC = 35
151      IF (X.LT.ETOL) IC = 1
152      DO 50 I=1,IC
153        AA = -AA*X/AK
154        IF (I.EQ.NM) GO TO 30
155        S = S - AA/(AK-FNM)
156        IF (ABS(AA).LE.XTOL*ABS(S)) GO TO 20
157        AK = AK + 1.0E0
158        GO TO 50
159   20   CONTINUE
160        IF (I.LT.2) GO TO 40
161        IF (ND-2.GT.I .OR. I.GT.ND-1) GO TO 60
162        AK = AK + 1.0E0
163        GO TO 50
164   30   S = S + AA*(-LOG(X)+PSIXN(ND))
165        XTOL = 3.0E0*TOL
166   40   AK = AK + 1.0E0
167   50 CONTINUE
168      IF (IC.NE.1) GO TO 340
169   60 IF (ND.EQ.1) S = S + (-LOG(X)+PSIXN(1))
170      IF (KODE.EQ.2) S = S*EXP(X)
171      EN(1) = S
172      EMX = 1.0E0
173      IF (M.EQ.1) GO TO 70
174      EN(IND) = S
175      AA = KS
176      IF (KODE.EQ.1) EMX = EXP(-X)
177      GO TO (220, 240), ICASE
178   70 IF (ICASE.EQ.2) RETURN
179      IF (KODE.EQ.1) EMX = EXP(-X)
180      EN(1) = (EMX-S)/X
181      RETURN
182   80 CONTINUE
183      DO 90 I=1,M
184        EN(I) = 1.0E0/(N+I-2)
185   90 CONTINUE
186      RETURN
187C-----------------------------------------------------------------------
188C     BACKWARD RECURSIVE MILLER ALGORITHM FOR
189C              E(N,X)=EXP(-X)*(X**(N-1))*U(N,N,X)
190C     WITH RECURSION AWAY FROM N=INTEGER CLOSEST TO X.
191C     U(A,B,X) IS THE SECOND CONFLUENT HYPERGEOMETRIC FUNCTION
192C-----------------------------------------------------------------------
193  100 CONTINUE
194      EMX = 1.0E0
195      IF (KODE.EQ.2) GO TO 130
196      IF (X.LE.XLIM) GO TO 120
197      NZ = M
198      DO 110 I=1,M
199        EN(I) = 0.0E0
200  110 CONTINUE
201      RETURN
202  120 EMX = EXP(-X)
203  130 CONTINUE
204      IX = X+0.5E0
205      KN = N + M - 1
206      IF (KN.LE.IX) GO TO 140
207      IF (N.LT.IX .AND. IX.LT.KN) GO TO 170
208      IF (N.GE.IX) GO TO 160
209      GO TO 340
210  140 ICASE = 1
211      KS = KN
212      ML = M - 1
213      MU = -1
214      IND = M
215      IF (KN.GT.1) GO TO 180
216  150 KS = 2
217      ICASE = 3
218      GO TO 180
219  160 ICASE = 2
220      IND = 1
221      KS = N
222      MU = M - 1
223      IF (N.GT.1) GO TO 180
224      IF (KN.EQ.1) GO TO 150
225      IX = 2
226  170 ICASE = 1
227      KS = IX
228      ML = IX - N
229      IND = ML + 1
230      MU = KN - IX
231  180 CONTINUE
232      IK = KS/2
233      AH = IK
234      JSET = 1 + KS - (IK+IK)
235C-----------------------------------------------------------------------
236C     START COMPUTATION FOR
237C              EN(IND) = C*U( A , A ,X)    JSET=1
238C              EN(IND) = C*U(A+1,A+1,X)    JSET=2
239C     FOR AN EVEN INTEGER A.
240C-----------------------------------------------------------------------
241      IC = 0
242      AA = AH + AH
243      AAMS = AA - 1.0E0
244      AAMS = AAMS*AAMS
245      TX = X + X
246      FX = TX + TX
247      AK = AH
248      XTOL = TOL
249      IF (TOL.LE.1.0E-3) XTOL = 20.0E0*TOL
250      CT = AAMS + FX*AH
251      EM = (AH+1.0E0)/((X+AA)*XTOL*SQRT(CT))
252      BK = AA
253      CC = AH*AH
254C-----------------------------------------------------------------------
255C     FORWARD RECURSION FOR P(IC),P(IC+1) AND INDEX IC FOR BACKWARD
256C     RECURSION
257C-----------------------------------------------------------------------
258      P1 = 0.0E0
259      P2 = 1.0E0
260  190 CONTINUE
261      IF (IC.EQ.99) GO TO 340
262      IC = IC + 1
263      AK = AK + 1.0E0
264      AT = BK/(BK+AK+CC+IC)
265      BK = BK + AK + AK
266      A(IC) = AT
267      BT = (AK+AK+X)/(AK+1.0E0)
268      B(IC) = BT
269      PT = P2
270      P2 = BT*P2 - AT*P1
271      P1 = PT
272      CT = CT + FX
273      EM = EM*AT*(1.0E0-TX/CT)
274      IF (EM*(AK+1.0E0).GT.P1*P1) GO TO 190
275      ICT = IC
276      KK = IC + 1
277      BT = TX/(CT+FX)
278      Y2 = (BK/(BK+CC+KK))*(P1/P2)*(1.0E0-BT+0.375E0*BT*BT)
279      Y1 = 1.0E0
280C-----------------------------------------------------------------------
281C     BACKWARD RECURRENCE FOR
282C              Y1=             C*U( A ,A,X)
283C              Y2= C*(A/(1+A/2))*U(A+1,A,X)
284C-----------------------------------------------------------------------
285      DO 200 K=1,ICT
286        KK = KK - 1
287        YT = Y1
288        Y1 = (B(KK)*Y1-Y2)/A(KK)
289        Y2 = YT
290  200 CONTINUE
291C-----------------------------------------------------------------------
292C     THE CONTIGUOUS RELATION
293C              X*U(B,C+1,X)=(C-B)*U(B,C,X)+U(B-1,C,X)
294C     WITH  B=A+1 , C=A IS USED FOR
295C              Y(2) = C * U(A+1,A+1,X)
296C     X IS INCORPORATED INTO THE NORMALIZING RELATION
297C-----------------------------------------------------------------------
298      PT = Y2/Y1
299      CNORM = 1.0E0 - PT*(AH+1.0E0)/AA
300      Y(1) = 1.0E0/(CNORM*AA+X)
301      Y(2) = CNORM*Y(1)
302      IF (ICASE.EQ.3) GO TO 210
303      EN(IND) = EMX*Y(JSET)
304      IF (M.EQ.1) RETURN
305      AA = KS
306      GO TO (220, 240), ICASE
307C-----------------------------------------------------------------------
308C     RECURSION SECTION  N*E(N+1,X) + X*E(N,X)=EMX
309C-----------------------------------------------------------------------
310  210 EN(1) = EMX*(1.0E0-Y(1))/X
311      RETURN
312  220 K = IND - 1
313      DO 230 I=1,ML
314        AA = AA - 1.0E0
315        EN(K) = (EMX-AA*EN(K+1))/X
316        K = K - 1
317  230 CONTINUE
318      IF (MU.LE.0) RETURN
319      AA = KS
320  240 K = IND
321      DO 250 I=1,MU
322        EN(K+1) = (EMX-X*EN(K))/AA
323        AA = AA + 1.0E0
324        K = K + 1
325  250 CONTINUE
326      RETURN
327  340 CONTINUE
328      IERR = 2
329      RETURN
330      END
331