Fcalva/OpenLibm
1
0
Fork 0
mirror of https://git.planet-casio.com/Lephenixnoir/OpenLibm.git synced 2025-01-01 06:23:39 +01:00
Code Issues Projects Releases Packages Wiki Activity Actions
OpenLibm/slatec/dexint.f
Viral B. Shah c977aa998f Add Makefile.extras to build libopenlibm-extras. 
Replace amos with slatec
2012-12-31 16:37:05 -05:00

336 lines
11 KiB
Fortran
Raw Blame History

*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
Reference in a new issue View git blame Copy permalink