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 .