SUBROUTINE EXPINT (RN,X,ENX,EXPENX,IERR) C C LANGUAGE. AMERICAN NATIONAL STANDARD FORTRAN C C DEFINITIONS. C EN(X)= INTEGRAL (EXP(-X*T)DT/(T**N)), FROM 1 TO INFINITY C RN(=N), POSITIVE INTEGER C X, REAL AND POSITIVE C C SPECIAL CASES C EN(0)=INFINITY(=RINF, MAXIMUM MACHINE VALUE) N .LE. 1 C EN(0)= 1/N-1 N .GT. 1 C C E0(X)= EXP(-X)/X X .GT. 1/RINF C E0(X)= INFINITY X .LE. 1/RINF C C USAGE. CALL EXPINT (RN,X,ENX,EXPENX,IERR) C C FORMAL PARAMETERS C RN REAL OR DOUBLE PRECISION TYPE INPUT C FOR A POSITIVE INTEGER N C X (SAME TYPE AS RN) INPUT C ENX= EN(X) (SAME TYPE AS X) OUTPUT C EXPENX= EXP(X)*EN(X) (SAME TYPE AS X) OUTPUT C IERR INTEGER VARIABLE OUTPUT C NORMAL RETURN IERR= 0 C ERROR RETURN IERR= 1, X AND/OR N NEGATIVE C ENX=EXPENX=-INFINITY (IMPOSSIBLE VALUE) C IERR= 2, N NON-INTEGER C ENX=EXPENX=INFINITY C C MODIFICATIONS C DOUBLE PRECISION UNIVAC 1108 RESULTS ARE OBTAINED IF AS C SET UP BELOW WHERE C C NBM=ACCURACY DESIRED OR MAXIMUM NUMBER OF BINARY C DIGITS IN THE MANTISSA OF A FLOATING POINT NUMBER C C WITH C (1) THE DOUBLE PRECISION TYPE STATEMENT C (2) THE MAXIMUM MACHINE VALUE AND THE MAXIMUM INTEGER C (=RMAXI) CONVERTIBLE TO A FLOATING POINT NUMBER C GIVEN AS DOUBLE PRECISION CONSTANTS C (3) DOUBLE PRECISION DECIMAL CONSTANTS C (4) DATA STATEMENT NBM=60 FOR THE CONTROL VARIABLE C (5) FUNCTION TYPE STATEMENTS - DEXP, DLOG. C C SINGLE PRECISION UNIVAC 1108 RESULTS ARE OBTAINED BY C (1) DELETING THE DOUBLE PRECISION TYPE STATEMENT C (2) ADJUSTING MAXIMUM MACHINE VALUE C (3) CHANGING THE D'S TO E'S ON THE DATA CARDS FOR ALL C DECIMAL CONSTANTS C (4) SETTING NBM=27 IN THE CONTROL VARIABLE DATA C (5) CHANGING FUNCTION TYPE - EXP, ALOG. C C FOR OTHER COMPUTERS THE APPROPRIATE VALUE OF NBM MUST C BE INSERTED AND ALL VALUES ADJUSTED ACCORDINGLY. C C IF A PRECOMPUTED VALUE OF TOLER(=2**(-NBM)) IS C INCLUDED IN A DATA STATEMENT, COMPUTATION OF THE CONTROL C VARIABLE MAY BE OMITTED AND THE DATA STATEMENT FOR NBM C DELETED. C C CAUTION - THE SUBROUTINE CANNOT READILY BE ADAPTED TO C COMPUTE THE EXPONENTIAL INTEGRAL FOR A COMPLEX C ARGUMENT AS THE CONTINUED FRACTION IS INVALID C ALONG THE NEGATIVE REAL AXIS. IN ADDITION MANY C OF THE COMPARISONS BECOME MEANINGLESS. C C METHOD. C POWER SERIES, X .LE. 1(=ULPS, UPPER LIMIT FOR POWER C SERIES) C ENX = SUM(TM), M=0,1,2,...,K C TM =-((-X)**M/1*2*3...M)/(M-N+1=D) M .NE. N-1 C = PTERM/D C TM = PTERM*(LOG(X)-PSI(N)) M .EQ. N-1 C PTERM(0)=-1 C PTERM(M+1)= PTERM(M)*(-X)/(M+1) C PSI(N)= -EULER+1+1/2+ ... +1/(N-1) C IF R.E.(=ABS(TM/SUM)) .LE. TOLER (=2**(-NBM)), M=K C C CONTINUED FRACTION, X .GT. 1 C EN(X)=EXP(-X)*(1I/I(X+N)- 1*N I/I(X+N+2)- C 2*(N+1)I/I(X+N+4)-... C EN(X)=EXP(-X)*II(AM I/I BM) M=1,2,...,K C AM(1)=1 AM(M)=-(M-1)*(N+M-2) C BM(M)=X+N+2*(M-1) C C EN(X)=EXP(-X)*FM(K)/GM(K)=EXP(-X)*F(K) C IF R.E. (=ABS(1-PREV/F)) .LE. TOLER(=2**(-NBM)),M=K C F=FM/GM PREV=FMM1/GMM1 C C RANGE. C FOR EXP(X)*EN(X) (N+X) .LE. MAXIMUM MACHINE VALUE C FOR EN(X) X=APPROXIMATELY 85.0, SINGLE PRECISION C 704.0, DOUBLE PRECISION C BEYOND THIS RANGE EN(X)=0 C C ACCURACY. THE NUMBER OF ACCURATE BINARY DIGITS IS ESSEN- C TIALLY THE LESSER OF C NBM - SQUARE ROOT OF NBM C NBM - I(NUMBER OF BINARY DIGITS REPRESENT- C ING THE INTEGER PART OF X) C (ACCURACY OF THE EXPONENTIAL ROUTINE) C C PRECISION. VARIABLE - BY SETTING THE DESIRED NBM OR TOLER. C C MAXIMUM TIME. UNIVAC 1108-EXEC II C (SECONDS) NBM=27 NBM=60 C .0025 .015 C C STORAGE. COMPILED BY THE UNIVAC 1108, EXEC 8/FORTRAN V, C THIS SUBROUTINE REQUIRES 478 WORDS OF STORAGE. C C C C C C C C C MACHINE DEPENDENT STATEMENTS C DOUBLE PRECISION AM,ASUM,BM,D,ENX,EULER, 1 EXPENX,EXPNX,F,FM,FMM1,FMM2,GM,GMM1, 2 GMM2,HALF,ONE,ONPTFV,PREV,PSI,PTERM, 3 RINF,RM,RMAXI,RN,RNM1,RRN,SUM,TEMP, 4 TM,TOLER,TWO,ULPS,X,XLOG,ZERO C C CONSTANTS C MAXIMUM MACHINE VALUE C MAXIMUM CONVERTIBLE INTEGER C DATA RINF / .898846567431157954D308 / DATA RMAXI / 134217727.D0 / C C ADDITIONAL CONSTANTS C DATA EULER / .577215664901532861D0 / DATA HALF,ONE,ONPTFV,TWO,ULPS,ZERO / 1 .5D0,1.D0,1.5D0,2.D0,1.D0,0.0D0 / C C CONTROL VARIABLE C C NBM=ACCURACY DESIRED OR MAXIMUM NUMBER OF BINARY C DIGITS IN THE MANTISSA OF A FLOATING POINT NUMBER C DATA NBM / 60 / TOLER=TWO**(-NBM) C C VALIDITY TEST FOR INPUT PARAMETERS C ERROR CONDITION, IMPOSSIBLE VALUES RETURNED C C NEGATIVE ZERO CHECKS C IERR=0 IF (RN .LT. ZERO) GO TO 3 IF (RN .GT. RMAXI) GO TO 10 C C VALIDITY TEST FOR N INTEGER C N=RN RRN=N IF (RRN) 2,1,2 1 IF ((RN-RRN)-TOLER) 10,10,2 2 IF ((ONE-(RRN/RN)) .LE. TOLER) GO TO 10 RRN=RRN+ONE IF ((ONE-(RN/RRN)) .LE. TOLER) GO TO 10 C C ERROR RETURN C N NON-INTEGER C IERR=2 ENX=RINF EXPENX=ENX RETURN C 3 IF (-RN) 12,10,12 C 10 IF (X) 11,30,20 11 IF (-X) 12,30,12 C C ERROR RETURN C X AND/OR N NEGATIVE C 12 IERR=1 ENX=-RINF EXPENX=ENX RETURN C C FUNCTION TYPE STATEMENTS C 20 EXPNX=DEXP(-X) XLOG=DLOG(X) C IF (RN .GE. HALF) GO TO 60 C C SPECIAL CASES C IF (X .LE. ONE/RINF) GO TO 40 EXPENX=ONE/X ENX=EXPNX*EXPENX RETURN C 30 IF (RN .LT. ONPTFV) GO TO 40 ENX=ONE/(RN-ONE) GO TO 50 C 40 ENX=RINF 50 EXPENX=ENX RETURN C 60 IF (RN .GT. RMAXI) GO TO 70 IF (X .LE. ULPS) GO TO 100 70 IF (X .LE. (RINF-RN)) GO TO 200 EXPENX=ZERO ENX=EXPENX RETURN C C METHOD --- POWER SERIES C 100 RM=0 PTERM=-ONE SUM=0 PSI=-EULER D=-(RN-ONE) C 110 IF (D .GE. HALF) GO TO 130 IF (-D .GE. HALF) GO TO 120 C C COMPUTE TM FOR M .EQ. N-1 C SUM=PTERM*(XLOG-PSI)+SUM GO TO 170 C C COMPUTE PSI(N) C 120 PSI=PSI+(ONE/(RM+ONE)) C C COMPUTE TM FOR M .NE. N-1 C 130 TM=PTERM/D SUM=TM+SUM C C TOLERANCE CHECK C C ZERO CHECKS C IF (SUM .LT. ZERO) GO TO 140 ASUM=SUM GO TO 150 C 140 ASUM=-SUM 150 IF (ASUM) 160,170,160 160 IF (TM .LT. ZERO) TM=-TM IF ( TM/ASUM .GT. TOLER) GO TO 170 C ENX= SUM EXPENX=ENX/EXPNX RETURN C C ADDITIONAL TERMS C 170 RM=RM+ONE D=D+ONE PTERM=-(X*PTERM)/RM GO TO 110 C C METHOD --- CONTINUED FRACTION C 200 RM=ONE FMM2=ONE GMM2=0 FMM1=0 GMM1=ONE PREV=FMM1/GMM1 AM=ONE BM=X+RN RNM1=RN-ONE C 210 FM=BM*FMM1 + AM*FMM2 GM=BM*GMM1 + AM*GMM2 F=FM/GM C C TOLERANCE CHECK C TEMP=ONE-(PREV/F) IF (TEMP .GT. ZERO) GO TO 220 C EXPENX=PREV GO TO 230 C 220 IF (TEMP .GT. TOLER) GO TO 240 C EXPENX=F 230 ENX=EXPNX*EXPENX RETURN C 240 IF( GM .LT. RINF/BM) GO TO 250 C C SCALING C C BOTH FM AND GM MUST BE TESTED IF N=1 AND X .LT. .44 C SCALING SHOULD NOT BE DELETED AS THE VALUES OF FM AND GM C MAY OVERFLOW FOR SOME CHOICES OF THE PARAMETERS. C FMM1=FMM1/BM GMM1=GMM1/BM FM=FM/BM GM=GM/BM C C ADDITIONAL CONVERGENTS C 250 AM=-RM*(RNM1+RM) RM=RM+ONE BM=BM+TWO FMM2=FMM1 GMM2=GMM1 FMM1=FM GMM1=GM PREV=F GO TO 210 C END