2
SUBROUTINE EXINT (X, N, KODE, M, TOL, EN, NZ, IERR)
3
C***BEGIN PROLOGUE EXINT
4
C***PURPOSE Compute an M member sequence of exponential integrals
5
C E(N+K,X), K=0,1,...,M-1 for N .GE. 1 and X .GE. 0.
8
C***TYPE SINGLE PRECISION (EXINT-S, DEXINT-D)
9
C***KEYWORDS EXPONENTIAL INTEGRAL, SPECIAL FUNCTIONS
10
C***AUTHOR Amos, D. E., (SNLA)
13
C EXINT computes M member sequences of exponential integrals
14
C E(N+K,X), K=0,1,...,M-1 for N .GE. 1 and X .GE. 0. The
15
C exponential integral is defined by
17
C E(N,X)=integral on (1,infinity) of EXP(-XT)/T**N
19
C where X=0.0 and N=1 cannot occur simultaneously. Formulas
20
C and notation are found in the NBS Handbook of Mathematical
23
C The power series is implemented for X .LE. XCUT and the
24
C confluent hypergeometric representation
26
C E(A,X) = EXP(-X)*(X**(A-1))*U(A,A,X)
28
C is computed for X .GT. XCUT. Since sequences are computed in
29
C a stable fashion by recurring away from X, A is selected as
30
C the integer closest to X within the constraint N .LE. A .LE.
31
C N+M-1. For the U computation, A is further modified to be the
32
C nearest even integer. Indices are carried forward or
33
C backward by the two term recursion relation
35
C K*E(K+1,X) + X*E(K,X) = EXP(-X)
37
C once E(A,X) is computed. The U function is computed by means
38
C of the backward recursive Miller algorithm applied to the
39
C three term contiguous relation for U(A+K,A,X), K=0,1,...
40
C This produces accurate ratios and determines U(A+K,A,X), and
41
C hence E(A,X), to within a multiplicative constant C.
42
C Another contiguous relation applied to C*U(A,A,X) and
43
C C*U(A+1,A,X) gets C*U(A+1,A+1,X), a quantity proportional to
44
C E(A+1,X). The normalizing constant C is obtained from the
45
C two term recursion relation above with K=A.
47
C Description of Arguments
50
C X X .GT. 0.0 for N=1 and X .GE. 0.0 for N .GE. 2
51
C N order of the first member of the sequence, N .GE. 1
52
C (X=0.0 and N=1 is an error)
53
C KODE a selection parameter for scaled values
54
C KODE=1 returns E(N+K,X), K=0,1,...,M-1.
55
C =2 returns EXP(X)*E(N+K,X), K=0,1,...,M-1.
56
C M number of exponential integrals in the sequence,
58
C TOL relative accuracy wanted, ETOL .LE. TOL .LE. 0.1
59
C ETOL = single precision unit roundoff = R1MACH(4)
62
C EN a vector of dimension at least M containing values
63
C EN(K) = E(N+K-1,X) or EXP(X)*E(N+K-1,X), K=1,M
65
C NZ underflow indicator
66
C NZ=0 a normal return
67
C NZ=M X exceeds XLIM and an underflow occurs.
68
C EN(K)=0.0E0 , K=1,M returned on KODE=1
70
C IERR=0, normal return, computation completed
71
C IERR=1, input error, no computation
72
C IERR=2, error, no computation
73
C algorithm termination condition not met
75
C***REFERENCES M. Abramowitz and I. A. Stegun, Handbook of
76
C Mathematical Functions, NBS AMS Series 55, U.S. Dept.
78
C D. E. Amos, Computation of exponential integrals, ACM
79
C Transactions on Mathematical Software 6, (1980),
80
C pp. 365-377 and pp. 420-428.
81
C***ROUTINES CALLED I1MACH, PSIXN, R1MACH
82
C***REVISION HISTORY (YYMMDD)
84
C 890531 Changed all specific intrinsics to generic. (WRB)
85
C 890531 REVISION DATE from Version 3.2
86
C 891214 Prologue converted to Version 4.0 format. (BAB)
87
C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
88
C 900326 Removed duplicate information from DESCRIPTION section.
90
C 910408 Updated the REFERENCES section. (WRB)
91
C 920207 Updated with code with a revision date of 880811 from
92
C D. Amos. Included correction of argument list. (WRB)
93
C 920501 Reformatted the REFERENCES section. (WRB)
94
C***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,
99
INTEGER I,IC,ICASE,ICT,IERR,IK,IND,IX,I1M,JSET,K,KK,KN,KODE,KS,M,
102
DIMENSION EN(*), A(99), B(99), Y(2)
103
C***FIRST EXECUTABLE STATEMENT EXINT
106
ETOL = MAX(R1MACH(4),0.5E-18)
107
IF (X.LT.0.0E0) IERR = 1
109
IF (KODE.LT.1 .OR. KODE.GT.2) 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
115
PT = 2.3026E0*R1MACH(5)*I1M
116
XLIM = PT - 6.907755E0
118
IF (BT.GT.1000.0E0) XLIM = PT - LOG(BT)
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
124
C-----------------------------------------------------------------------
125
C SERIES FOR E(N,X) FOR X.LE.XCUT
126
C-----------------------------------------------------------------------
129
C-----------------------------------------------------------------------
130
C ICASE=1 MEANS INTEGER CLOSEST TO X IS 2 AND N=1
131
C ICASE=2 MEANS INTEGER CLOSEST TO X IS 0,1, OR 2 AND N.GE.2
132
C-----------------------------------------------------------------------
134
IF (IX.GT.N) ICASE = 1
144
IF (ND.EQ.1) GO TO 10
151
IF (X.LT.ETOL) IC = 1
154
IF (I.EQ.NM) GO TO 30
156
IF (ABS(AA).LE.XTOL*ABS(S)) GO TO 20
161
IF (ND-2.GT.I .OR. I.GT.ND-1) GO TO 60
164
30 S = S + AA*(-LOG(X)+PSIXN(ND))
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)
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)
184
EN(I) = 1.0E0/(N+I-2)
187
C-----------------------------------------------------------------------
188
C BACKWARD RECURSIVE MILLER ALGORITHM FOR
189
C E(N,X)=EXP(-X)*(X**(N-1))*U(N,N,X)
190
C WITH RECURSION AWAY FROM N=INTEGER CLOSEST TO X.
191
C U(A,B,X) IS THE SECOND CONFLUENT HYPERGEOMETRIC FUNCTION
192
C-----------------------------------------------------------------------
195
IF (KODE.EQ.2) GO TO 130
196
IF (X.LE.XLIM) GO TO 120
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
215
IF (KN.GT.1) GO TO 180
223
IF (N.GT.1) GO TO 180
224
IF (KN.EQ.1) GO TO 150
234
JSET = 1 + KS - (IK+IK)
235
C-----------------------------------------------------------------------
236
C START COMPUTATION FOR
237
C EN(IND) = C*U( A , A ,X) JSET=1
238
C EN(IND) = C*U(A+1,A+1,X) JSET=2
239
C FOR AN EVEN INTEGER A.
240
C-----------------------------------------------------------------------
249
IF (TOL.LE.1.0E-3) XTOL = 20.0E0*TOL
251
EM = (AH+1.0E0)/((X+AA)*XTOL*SQRT(CT))
254
C-----------------------------------------------------------------------
255
C FORWARD RECURSION FOR P(IC),P(IC+1) AND INDEX IC FOR BACKWARD
257
C-----------------------------------------------------------------------
261
IF (IC.EQ.99) GO TO 340
264
AT = BK/(BK+AK+CC+IC)
267
BT = (AK+AK+X)/(AK+1.0E0)
273
EM = EM*AT*(1.0E0-TX/CT)
274
IF (EM*(AK+1.0E0).GT.P1*P1) GO TO 190
278
Y2 = (BK/(BK+CC+KK))*(P1/P2)*(1.0E0-BT+0.375E0*BT*BT)
280
C-----------------------------------------------------------------------
281
C BACKWARD RECURRENCE FOR
283
C Y2= C*(A/(1+A/2))*U(A+1,A,X)
284
C-----------------------------------------------------------------------
288
Y1 = (B(KK)*Y1-Y2)/A(KK)
291
C-----------------------------------------------------------------------
292
C THE CONTIGUOUS RELATION
293
C X*U(B,C+1,X)=(C-B)*U(B,C,X)+U(B-1,C,X)
294
C WITH B=A+1 , C=A IS USED FOR
295
C Y(2) = C * U(A+1,A+1,X)
296
C X IS INCORPORATED INTO THE NORMALIZING RELATION
297
C-----------------------------------------------------------------------
299
CNORM = 1.0E0 - PT*(AH+1.0E0)/AA
300
Y(1) = 1.0E0/(CNORM*AA+X)
302
IF (ICASE.EQ.3) GO TO 210
303
EN(IND) = EMX*Y(JSET)
306
GO TO (220, 240), ICASE
307
C-----------------------------------------------------------------------
308
C RECURSION SECTION N*E(N+1,X) + X*E(N,X)=EMX
309
C-----------------------------------------------------------------------
310
210 EN(1) = EMX*(1.0E0-Y(1))/X
315
EN(K) = (EMX-AA*EN(K+1))/X
322
EN(K+1) = (EMX-X*EN(K))/AA
added
removed
deps/openlibm/slatec/exint.f
