SUBROUTINE INERFC(X, NMAX, ACC, FZERO, F, IFLAG) INE 10
C THIS PROGRAM GENERATES FN(X) FOR N=1,2,...,NMAX, WHERE
C FN(X)=EXP(X**2)*INERFC(X), IF X.GE.0., AND FN(X)=INERFC(X), IF
C X.LT.0., INERFC(X) BEING THE N-TH REPEATED INTEGRAL OF THE COERROR
C FUNCTION ERFC(X), EXTENDED FROM X TO INFINITY. THE PROGRAM ALSO
C SUPPLIES EXP(X**2)*ERFC(X),IF X.GE.0., AND ERFC(X), IF X.LT.0..
C THE USER HAS TO SPECIFY THE DESIRED ACCURACY, AS WELL AS THE PRECISION
C OF THE HIGHEST PRECISION MODE HE IS PREPARED TO USE. THE FORMER IS
C INPUT VIA THE PARAMETER ACC IN THE CALL LIST, WHICH DESIGNATES THE
C NUMBER OF CORRECT SIGNIFICANT DECIMAL DIGITS DESIRED IN THE RESULTS.
C THE LATTER IS INPUT VIA THE PARAMETER PREC IN THE DATA STATEMENT. THE
C VALUE OF PREC SHOULD BE SET APPROXIMATELY EQUAL TO BETA*ALOG(2.)/
C ALOG(10.), WHERE BETA IS THE NUMBER OF BINARY DIGITS AVAILABLE IN THE
C MANTISSA OF THE HIGHEST PRECISION FLOATING-POINT WORD. ORDINARILY, THE
C HIGHEST PRECISION IS DOUBLE PRECISION, EITHER HARDWARE GENERATED OR
C SOFTWARE GENERATED, BUT IT MAY ALSO BE SINGLE PRECISION. IN ANY CASE,
C ACC HAS TO BE LESS THAN PREC, THE DIFFERENCE BEING AT LEAST EQUAL TO
C ONE, PREFERABLY SEVERAL UNITS LARGER. EVIDENTLY, ACC SHOULD NOT BE
C SPECIFIED TO EXCEED THE PRECISION OF THE LOWEST PRECISION MODE USED.
C THE DATA STATEMENT CONTAINS ANOTHER MACHINE-DEPENDABLE PARAMETER,
C BOTEXP, WHICH IS, APPROXIMATELY, THE SMALLEST NEGATIVE NUMBER SUCH
C THAT 10.**BOTEXP IS STILL REPRESENTABLE ON THE COMPUTER IN SINGLE PRE-
C CISION FLOATING-POINT ARITHMETIC. THE TYPE OF BOTEXP IS REAL. IN THE
C DATA STATEMENT BELOW THE PARAMETERS PREC AND BOTEXP ARE SET TO CORRE-
C SPOND TO THE MACHINE CHARACTERISTICS OF THE CDC 6500 COMPUTER.
C PARAMETER LIST
C X - - THE ARGUMENT OF INERFC. TYPE REAL.
C NMAX - - THE UPPER LIMIT ON N. TYPE INTEGER.
C ACC - - THE NUMBER OF CORRECT SIGNIFICANT DECIMAL DIGITS DESIRED
C IN THE RESULTS. TYPE REAL.
C FZERO - - AN OUTPUT VARIABLE RETURNING EXP(X**2)*ERFC(X), IF
C X.GE.0., AND ERFC(X), IF X.LT.0. TYPE REAL.
C F - - THE NAME OF A ONE-DIMENSIONAL ARRAY HOLDING THE RESULTS
C FN(X), N=1,2,...,NMAX. THE ARRAY HAS DIMENSION NMAX.
C IFLAG - - AN ERROR FLAG INDICATING A NUMBER OF FATAL ERROR CON-
C DITIONS UPON EXECUTION. TYPE INTEGER. THE VALUES OF THIS
C VARIABLE HAVE THE FOLLOWING MEANINGS.
C IFLAG=0 - NO ERROR CONDITION.
C IFLAG=1 - ILLEGAL NEGATIVE OR ZERO VALUE OF NMAX.
C IFLAG=2 - EXCESSIVELY LARGE VALUE OF NMAX. REDIMENSION THE
C ARRAYS R0 AND R1 TO HAVE DIMENSION .GE. NMAX, IF
C INDEED NMAX MUST BE THAT LARGE.
C IFLAG=3 - INAPPROPRIATE ACCURACY SPECIFICATION.
C IFLAG=4 - UNUSUAL UNDERFLOW OF EXP(-X**2). TRY REDUCING
C THE ERROR TOLERANCE PREC-ACC, EITHER BY DE-
C CREASING PREC OR BY INCREASING ACC.
C IFLAG=5 - TAYLOR SERIES FOR ERFC(X) DOES NOT CONVERGE
C WITHIN 200 TERMS FOR SOME UNKNOWN REASON.
C IFLAG=6 - NU IN THE BACKWARD RECURRENCE ALGORITHM EXCEEDS
C 5000, PROBABLY BECAUSE X IS A SMALL POSITIVE
C NUMBER, NMAX RELATIVELY LARGE, AND THE ERROR
C TOLERANCE SMALL. TRY INCREASING THE ERROR TOLER-
C ANCE PREC-ACC, EITHER BY INCREASING PREC OR BY
C DECREASING ACC.
C THE SUBROUTINE BELOW IS WRITTEN WITH THE ASSUMPTION THAT PART OF
C THE COMPUTATION IS DONE IN DOUBLE PRECISION AND THE REST IN SINGLE
C PRECISION. TO OBTAIN A VERSION OF THIS SUBROUTINE, THAT OPERATES
C ENTIRELY IN SINGLE PRECISION, THE FOLLOWING CHANGES MUST BE MADE.
C (1) DELETE THE DOUBLE PRECISION TYPE DECLARATION.
C (2) REPLACE THE SECTION OF PROGRAM FROM STATEMENT LABELED 30
C (INCLUSIVE) TO STATEMENT LABELED 60 (INCLUSIVE) BY THE FOLLOWING
C PIECE OF PROGRAM.
C 30 DEPS=.5*10.**(-PREC)
C TE=C*X
C SUM=1.-TE
C N=0
C N1=1
C 40 N=N+1
C N1=N1+2
C IF(N.GT.200) GOTO 80
C TE=-TE*XSQ/FLOAT(N)
C TERM=TE/FLOAT(N1)
C SUM=SUM-TERM
C IF(ABS(TERM).GT.DEPS*ABS(SUM)) GOTO 40
C F1=T*SUM
C THE ARRAYS R0 AND R1 IN THE DIMENSION STATEMENT ARE LOCAL TO THE
C SUBROUTINE. THEY BOTH NEED TO HAVE DIMENSION LARGER OR EQUAL TO NMAX.
C WE DECLARED THEIR DIMENSION TO BE 500, ASSUMING THAT NMAX WILL NOT
C EXCEED 500. IN THE EVENT IT DOES, THE SUBROUTINE EXITS WITH THE ERROR
C FLAG AT 2. IF THE USER EXPECTS HIS VALUES OF NMAX TO BE CONSISTENTLY
C MUCH LESS THAN 500, HE CAN SAVE STORAGE BY LOWERING THE DIMENSION OF
C R0 AND R1 AND AT THE SAME TIME ADJUSTING THE ERROR EXIT STATEMENT (THE
C THIRD EXECUTABLE STATEMENT IN THE PROGRAM).
C WITH THE EXCEPTION NOTED UNDER IFLAG=4, ANY VARIABLE THAT UNDER-
C FLOWS MAY BE SET EQUAL TO ZERO. NO PRECAUTION IS TAKEN AGAINST OVER-
C FLOW, WHICH MAY OCCUR IN THE DO-LOOP AFTER STATEMENT 15 IF X.LT.0.
C AND ABS(X) AND/OR NMAX IS LARGE.
C REFERENCE - W. GAUTSCHI, EVALUATION OF THE REPEATED INTEGRALS OF
C THE COERROR FUNCTION, TRANS. ON MATH. SOFTWARE.
DIMENSION F(NMAX), R0(500), R1(500)
DOUBLE PRECISION DEPS, DC, DX, DXSQ, DFM1, DF0, DF1, DN, DN1,
* DTE, DTERM, DSUM
LOGICAL FRSTTM
DATA FRSTTM, PREC, BOTEXP /.TRUE.,28.8989,-293./
IFLAG = 0
IF (NMAX.LT.1) GO TO 230
IF (NMAX.GT.500) GO TO 240
TOL = PREC - ACC
IF (TOL.LT.1.) GO TO 250
I = 1
EPS = .5*10.**(-ACC)
XSQ = X*X
IF (.NOT.FRSTTM) GO TO 10
P = ATAN(1.)
AL10 = ALOG(10.)
C = SQRT(1./P)
DC = DSQRT(1.D0/DATAN(1.D0))
FRSTTM = .FALSE.
10 S = 0.
IF (-XSQ.GT.AL10*BOTEXP) S = EXP(-XSQ)
B = .5*AL10*TOL + .25*ALOG(2.*P)
B1 = .5*AL10*(TOL-1.)
IF (XSQ.GT.B) GO TO 90
IF (X.LT.0.) GO TO 20
F0 = C
IF (S.EQ.0.) GO TO 260
T = 1./S
SQ = SQRT(2.*FLOAT(NMAX))
IF ((X.GE.1.12 .AND. XSQ+SQ*X.GT.B) .OR. (X.LT.1.12 .AND.
* SQ*X.GT.B1)) GO TO 100
GO TO 30
20 F0 = C*S
T = 1.
C
C EVALUATION OF ERFC(X) BY TAYLOR SERIES IN DOUBLE PRECISION
C
30 DEPS = .5D0*10.D0**(-PREC)
DX = DBLE(X)
DXSQ = DX*DX
DTE = DC*DX
DSUM = 1.D0 - DTE
DN = 0.D0
DN1 = 1.D0
40 DN = DN + 1.D0
DN1 = DN1 + 2.D0
IF (DN.GT.2.D2) GO TO 270
DTE = -DTE*DXSQ/DN
DTERM = DTE/DN1
DSUM = DSUM - DTERM
IF (DABS(DTERM).GT.DEPS*DABS(DSUM)) GO TO 40
IF (X.LT.0.) GO TO 60
C
C EVALUATION OF EXP(X**2)*INERFC(X),X.GE.0.,BY FORWARD RECURSION IN
C DOUBLE PRECISION
C
DF0 = DC
DF1 = DEXP(DXSQ)*DSUM
FZERO = SNGL(DF1)
DO 50 N=1,NMAX
DFM1 = DF0
DF0 = DF1
DF1 = (-DX*DF0+.5D0*DFM1)/DBLE(FLOAT(N))
F(N) = SNGL(DF1)
50 CONTINUE
RETURN
C
C EVALUATION OF INERFC(X),X.LT.0.,BY FORWARD RECURSION
C
60 F1 = SNGL(DSUM)
70 FZERO = F1
DO 80 N=1,NMAX
FM1 = F0
F0 = F1
F1 = (-X*F0+.5*FM1)/FLOAT(N)
F(N) = F1
80 CONTINUE
RETURN
90 IF (X.LT.0.) GO TO 110
100 N0 = NMAX
X0 = X
GO TO 130
110 IF (XSQ.LT.(ACC+1.)*AL10-.572) GO TO 120
F0 = C*S
F1 = 2.
GO TO 70
120 I = 2
N0 = 1
X0 = -X
C
C BACKWARD RECURRENCE ALGORITHM FOR EVALUATING FN(ABS(X)) FOR N=1,2,
C ...,N0, WHERE N0=NMAX IF X.GT.0. AND N0=1 IF X.LT.0.
C
130 DO 140 N=1,N0
R1(N) = 0.
140 CONTINUE
FN0 = FLOAT(N0)
NU = IFIX((SQRT(FN0)+(2.3026*ACC+1.3863)/(2.8284*X0))**2-5.)
150 NU = NU + 10
IF (NU.GT.5000) GO TO 280
N0P1 = N0 + 1
NUM = NU - N0P1
DO 160 N=1,N0
R0(N) = R1(N)
160 CONTINUE
R = 0.
DO 170 K=1,NUM
N = NU - K
R = .5/(X0+FLOAT(N+1)*R)
170 CONTINUE
DO 180 K=1,N0
N = N0P1 - K
R = .5/(X0+FLOAT(N+1)*R)
R1(N) = R
180 CONTINUE
DO 190 N=1,N0
IF (ABS(R1(N)-R0(N)).GT.EPS*ABS(R1(N))) GO TO 150
190 CONTINUE
F0 = .5*C/(X0+R1(1))
FZERO = F0
FF = F0
DO 200 N=1,N0
FF = R1(N)*FF
F(N) = FF
200 CONTINUE
GO TO (210, 220), I
210 RETURN
C
C STARTING VALUE FOR INERFC(X),X.LT.0.,PRIOR TO FORWARD RECURSION
C
220 F1 = 2. - S*F0
F0 = C*S
GO TO 70
230 IFLAG = 1
RETURN
240 IFLAG = 2
RETURN
250 IFLAG = 3
RETURN
260 IFLAG = 4
RETURN
270 IFLAG = 5
RETURN
280 IFLAG = 6
RETURN
END
.