1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
|
*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
DOUBLE PRECISION A,AA,AAMS,AH,AK,AT,B,BK,BT,CC,CNORM,CT,EM,EMX,EN,
1 ETOL,FNM,FX,PT,P1,P2,S,TOL,TX,X,XCUT,XLIM,XTOL,Y,
2 YT,Y1,Y2
DOUBLE PRECISION D1MACH,DPSIXN
INTEGER I,IC,ICASE,ICT,IERR,IK,IND,IX,I1M,JSET,K,KK,KN,KODE,KS,M,
1 ML,MU,N,ND,NM,NZ
INTEGER I1MACH
DIMENSION EN(*), A(99), B(99), Y(2)
SAVE XCUT
DATA XCUT / 2.0D0 /
C***FIRST EXECUTABLE STATEMENT DEXINT
IERR = 0
NZ = 0
ETOL = MAX(D1MACH(4),0.5D-18)
IF (X.LT.0.0D0) IERR = 1
IF (N.LT.1) IERR = 1
IF (KODE.LT.1 .OR. KODE.GT.2) IERR = 1
IF (M.LT.1) IERR = 1
IF (TOL.LT.ETOL .OR. TOL.GT.0.1D0) IERR = 1
IF (X.EQ.0.0D0 .AND. N.EQ.1) IERR = 1
IF(IERR.NE.0) RETURN
I1M = -I1MACH(15)
PT = 2.3026D0*I1M*D1MACH(5)
XLIM = PT - 6.907755D0
BT = PT + (N+M-1)
IF (BT.GT.1000.0D0) XLIM = PT - LOG(BT)
C
IF (X.GT.XCUT) GO TO 100
IF (X.EQ.0.0D0 .AND. N.GT.1) GO TO 80
C-----------------------------------------------------------------------
C SERIES FOR E(N,X) FOR X.LE.XCUT
C-----------------------------------------------------------------------
TX = X + 0.5D0
IX = TX
C-----------------------------------------------------------------------
C ICASE=1 MEANS INTEGER CLOSEST TO X IS 2 AND N=1
C ICASE=2 MEANS INTEGER CLOSEST TO X IS 0,1, OR 2 AND N.GE.2
C-----------------------------------------------------------------------
ICASE = 2
IF (IX.GT.N) ICASE = 1
NM = N - ICASE + 1
ND = NM + 1
IND = 3 - ICASE
MU = M - IND
ML = 1
KS = ND
FNM = NM
S = 0.0D0
XTOL = 3.0D0*TOL
IF (ND.EQ.1) GO TO 10
XTOL = 0.3333D0*TOL
S = 1.0D0/FNM
10 CONTINUE
AA = 1.0D0
AK = 1.0D0
IC = 35
IF (X.LT.ETOL) IC = 1
DO 50 I=1,IC
AA = -AA*X/AK
IF (I.EQ.NM) GO TO 30
S = S - AA/(AK-FNM)
IF (ABS(AA).LE.XTOL*ABS(S)) GO TO 20
AK = AK + 1.0D0
GO TO 50
20 CONTINUE
IF (I.LT.2) GO TO 40
IF (ND-2.GT.I .OR. I.GT.ND-1) GO TO 60
AK = AK + 1.0D0
GO TO 50
30 S = S + AA*(-LOG(X)+DPSIXN(ND))
XTOL = 3.0D0*TOL
40 AK = AK + 1.0D0
50 CONTINUE
IF (IC.NE.1) GO TO 340
60 IF (ND.EQ.1) S = S + (-LOG(X)+DPSIXN(1))
IF (KODE.EQ.2) S = S*EXP(X)
EN(1) = S
EMX = 1.0D0
IF (M.EQ.1) GO TO 70
EN(IND) = S
AA = KS
IF (KODE.EQ.1) EMX = EXP(-X)
GO TO (220, 240), ICASE
70 IF (ICASE.EQ.2) RETURN
IF (KODE.EQ.1) EMX = EXP(-X)
EN(1) = (EMX-S)/X
RETURN
80 CONTINUE
DO 90 I=1,M
EN(I) = 1.0D0/(N+I-2)
90 CONTINUE
RETURN
C-----------------------------------------------------------------------
C BACKWARD RECURSIVE MILLER ALGORITHM FOR
C E(N,X)=EXP(-X)*(X**(N-1))*U(N,N,X)
C WITH RECURSION AWAY FROM N=INTEGER CLOSEST TO X.
C U(A,B,X) IS THE SECOND CONFLUENT HYPERGEOMETRIC FUNCTION
C-----------------------------------------------------------------------
100 CONTINUE
EMX = 1.0D0
IF (KODE.EQ.2) GO TO 130
IF (X.LE.XLIM) GO TO 120
NZ = M
DO 110 I=1,M
EN(I) = 0.0D0
110 CONTINUE
RETURN
120 EMX = EXP(-X)
130 CONTINUE
TX = X + 0.5D0
IX = TX
KN = N + M - 1
IF (KN.LE.IX) GO TO 140
IF (N.LT.IX .AND. IX.LT.KN) GO TO 170
IF (N.GE.IX) GO TO 160
GO TO 340
140 ICASE = 1
KS = KN
ML = M - 1
MU = -1
IND = M
IF (KN.GT.1) GO TO 180
150 KS = 2
ICASE = 3
GO TO 180
160 ICASE = 2
IND = 1
KS = N
MU = M - 1
IF (N.GT.1) GO TO 180
IF (KN.EQ.1) GO TO 150
IX = 2
170 ICASE = 1
KS = IX
ML = IX - N
IND = ML + 1
MU = KN - IX
180 CONTINUE
IK = KS/2
AH = IK
JSET = 1 + KS - (IK+IK)
C-----------------------------------------------------------------------
C START COMPUTATION FOR
C EN(IND) = C*U( A , A ,X) JSET=1
C EN(IND) = C*U(A+1,A+1,X) JSET=2
C FOR AN EVEN INTEGER A.
C-----------------------------------------------------------------------
IC = 0
AA = AH + AH
AAMS = AA - 1.0D0
AAMS = AAMS*AAMS
TX = X + X
FX = TX + TX
AK = AH
XTOL = TOL
IF (TOL.LE.1.0D-3) XTOL = 20.0D0*TOL
CT = AAMS + FX*AH
EM = (AH+1.0D0)/((X+AA)*XTOL*SQRT(CT))
BK = AA
CC = AH*AH
C-----------------------------------------------------------------------
C FORWARD RECURSION FOR P(IC),P(IC+1) AND INDEX IC FOR BACKWARD
C RECURSION
C-----------------------------------------------------------------------
P1 = 0.0D0
P2 = 1.0D0
190 CONTINUE
IF (IC.EQ.99) GO TO 340
IC = IC + 1
AK = AK + 1.0D0
AT = BK/(BK+AK+CC+IC)
BK = BK + AK + AK
A(IC) = AT
BT = (AK+AK+X)/(AK+1.0D0)
B(IC) = BT
PT = P2
P2 = BT*P2 - AT*P1
P1 = PT
CT = CT + FX
EM = EM*AT*(1.0D0-TX/CT)
IF (EM*(AK+1.0D0).GT.P1*P1) GO TO 190
ICT = IC
KK = IC + 1
BT = TX/(CT+FX)
Y2 = (BK/(BK+CC+KK))*(P1/P2)*(1.0D0-BT+0.375D0*BT*BT)
Y1 = 1.0D0
C-----------------------------------------------------------------------
C BACKWARD RECURRENCE FOR
C Y1= C*U( A ,A,X)
C Y2= C*(A/(1+A/2))*U(A+1,A,X)
C-----------------------------------------------------------------------
DO 200 K=1,ICT
KK = KK - 1
YT = Y1
Y1 = (B(KK)*Y1-Y2)/A(KK)
Y2 = YT
200 CONTINUE
C-----------------------------------------------------------------------
C THE CONTIGUOUS RELATION
C X*U(B,C+1,X)=(C-B)*U(B,C,X)+U(B-1,C,X)
C WITH B=A+1 , C=A IS USED FOR
C Y(2) = C * U(A+1,A+1,X)
C X IS INCORPORATED INTO THE NORMALIZING RELATION
C-----------------------------------------------------------------------
PT = Y2/Y1
CNORM = 1.0E0 - PT*(AH+1.0E0)/AA
Y(1) = 1.0E0/(CNORM*AA+X)
Y(2) = CNORM*Y(1)
IF (ICASE.EQ.3) GO TO 210
EN(IND) = EMX*Y(JSET)
IF (M.EQ.1) RETURN
AA = KS
GO TO (220, 240), ICASE
C-----------------------------------------------------------------------
C RECURSION SECTION N*E(N+1,X) + X*E(N,X)=EMX
C-----------------------------------------------------------------------
210 EN(1) = EMX*(1.0E0-Y(1))/X
RETURN
220 K = IND - 1
DO 230 I=1,ML
AA = AA - 1.0D0
EN(K) = (EMX-AA*EN(K+1))/X
K = K - 1
230 CONTINUE
IF (MU.LE.0) RETURN
AA = KS
240 K = IND
DO 250 I=1,MU
EN(K+1) = (EMX-X*EN(K))/AA
AA = AA + 1.0D0
K = K + 1
250 CONTINUE
RETURN
340 CONTINUE
IERR = 2
RETURN
END
|
