SUBROUTINE CALERF(ARG,RESULT,JINT)
C------------------------------------------------------------------
C
C This packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x)
C for a real argument x. It contains three FUNCTION type
C subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX),
C and one SUBROUTINE type subprogram, CALERF. The calling
C statements for the primary entries are:
C
C Y=ERF(X) (or Y=DERF(X)),
C
C Y=ERFC(X) (or Y=DERFC(X)),
C and
C Y=ERFCX(X) (or Y=DERFCX(X)).
C
C The routine CALERF is intended for internal packet use only,
C all computations within the packet being concentrated in this
C routine. The function subprograms invoke CALERF with the
C statement
C
C CALL CALERF(ARG,RESULT,JINT)
C
C where the parameter usage is as follows
C
C Function Parameters for CALERF
C call ARG Result JINT
C
C ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0
C ERFC(ARG) ABS(ARG) .LT. XBIG ERFC(ARG) 1
C ERFCX(ARG) XNEG .LT. ARG .LT. XMAX ERFCX(ARG) 2
C
C The main computation evaluates near-minimax approximations
C from "Rational Chebyshev approximations for the error function"
C by W. J. Cody, Math. Comp., 1969, PP. 631-638. This
C transportable program uses rational functions that theoretically
C approximate erf(x) and erfc(x) to at least 18 significant
C decimal digits. The accuracy achieved depends on the arithmetic
C system, the compiler, the intrinsic functions, and proper
C selection of the machine-dependent constants.
C
C*******************************************************************
C*******************************************************************
C
C Explanation of machine-dependent constants
C
C XMIN = the smallest positive floating-point number.
C XINF = the largest positive finite floating-point number.
C XNEG = the largest negative argument acceptable to ERFCX;
C the negative of the solution to the equation
C 2*exp(x*x) = XINF.
C XSMALL = argument below which erf(x) may be represented by
C 2*x/sqrt(pi) and above which x*x will not underflow.
C A conservative value is the largest machine number X
C such that 1.0 + X = 1.0 to machine precision.
C XBIG = largest argument acceptable to ERFC; solution to
C the equation: W(x) * (1-0.5/x**2) = XMIN, where
C W(x) = exp(-x*x)/[x*sqrt(pi)].
C XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to
C machine precision. A conservative value is
C 1/[2*sqrt(XSMALL)]
C XMAX = largest acceptable argument to ERFCX; the minimum
C of XINF and 1/[sqrt(pi)*XMIN].
C
C Approximate values for some important machines are:
C
C XMIN XINF XNEG XSMALL
C
C CDC 7600 (S.P.) 3.13E-294 1.26E+322 -27.220 7.11E-15
C CRAY-1 (S.P.) 4.58E-2467 5.45E+2465 -75.345 7.11E-15
C IEEE (IBM/XT,
C SUN, etc.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8
C IEEE (IBM/XT,
C SUN, etc.) (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16
C IBM 195 (D.P.) 5.40D-79 7.23E+75 -13.190 1.39D-17
C UNIVAC 1108 (D.P.) 2.78D-309 8.98D+307 -26.615 1.73D-18
C VAX D-Format (D.P.) 2.94D-39 1.70D+38 -9.345 1.39D-17
C VAX G-Format (D.P.) 5.56D-309 8.98D+307 -26.615 1.11D-16
C
C
C XBIG XHUGE XMAX
C
C CDC 7600 (S.P.) 25.922 8.39E+6 1.80X+293
C CRAY-1 (S.P.) 75.326 8.39E+6 5.45E+2465
C IEEE (IBM/XT,
C SUN, etc.) (S.P.) 9.194 2.90E+3 4.79E+37
C IEEE (IBM/XT,
C SUN, etc.) (D.P.) 26.543 6.71D+7 2.53D+307
C IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75
C UNIVAC 1108 (D.P.) 26.582 5.37D+8 8.98D+307
C VAX D-Format (D.P.) 9.269 1.90D+8 1.70D+38
C VAX G-Format (D.P.) 26.569 6.71D+7 8.98D+307
C
C*******************************************************************
C*******************************************************************
C
C Error returns
C
C The program returns ERFC = 0 for ARG .GE. XBIG;
C
C ERFCX = XINF for ARG .LT. XNEG;
C and
C ERFCX = 0 for ARG .GE. XMAX.
C
C
C Intrinsic functions required are:
C
C ABS, AINT, EXP
C
C
C Author: W. J. Cody
C Mathematics and Computer Science Division
C Argonne National Laboratory
C Argonne, IL 60439
C
C Latest modification: March 19, 1990
C
C------------------------------------------------------------------
INTEGER I,JINT
CS REAL
CD DOUBLE PRECISION
1 A,ARG,B,C,D,DEL,FOUR,HALF,P,ONE,Q,RESULT,SIXTEN,SQRPI,
2 TWO,THRESH,X,XBIG,XDEN,XHUGE,XINF,XMAX,XNEG,XNUM,XSMALL,
3 Y,YSQ,ZERO
DIMENSION A(5),B(4),C(9),D(8),P(6),Q(5)
C------------------------------------------------------------------
C Mathematical constants
C------------------------------------------------------------------
CS DATA FOUR,ONE,HALF,TWO,ZERO/4.0E0,1.0E0,0.5E0,2.0E0,0.0E0/,
CS 1 SQRPI/5.6418958354775628695E-1/,THRESH/0.46875E0/,
CS 2 SIXTEN/16.0E0/
CD DATA FOUR,ONE,HALF,TWO,ZERO/4.0D0,1.0D0,0.5D0,2.0D0,0.0D0/,
CD 1 SQRPI/5.6418958354775628695D-1/,THRESH/0.46875D0/,
CD 2 SIXTEN/16.0D0/
C------------------------------------------------------------------
C Machine-dependent constants
C------------------------------------------------------------------
CS DATA XINF,XNEG,XSMALL/3.40E+38,-9.382E0,5.96E-8/,
CS 1 XBIG,XHUGE,XMAX/9.194E0,2.90E3,4.79E37/
CD DATA XINF,XNEG,XSMALL/1.79D308,-26.628D0,1.11D-16/,
CD 1 XBIG,XHUGE,XMAX/26.543D0,6.71D7,2.53D307/
C------------------------------------------------------------------
C Coefficients for approximation to erf in first interval
C------------------------------------------------------------------
CS DATA A/3.16112374387056560E00,1.13864154151050156E02,
CS 1 3.77485237685302021E02,3.20937758913846947E03,
CS 2 1.85777706184603153E-1/
CS DATA B/2.36012909523441209E01,2.44024637934444173E02,
CS 1 1.28261652607737228E03,2.84423683343917062E03/
CD DATA A/3.16112374387056560D00,1.13864154151050156D02,
CD 1 3.77485237685302021D02,3.20937758913846947D03,
CD 2 1.85777706184603153D-1/
CD DATA B/2.36012909523441209D01,2.44024637934444173D02,
CD 1 1.28261652607737228D03,2.84423683343917062D03/
C------------------------------------------------------------------
C Coefficients for approximation to erfc in second interval
C------------------------------------------------------------------
CS DATA C/5.64188496988670089E-1,8.88314979438837594E0,
CS 1 6.61191906371416295E01,2.98635138197400131E02,
CS 2 8.81952221241769090E02,1.71204761263407058E03,
CS 3 2.05107837782607147E03,1.23033935479799725E03,
CS 4 2.15311535474403846E-8/
CS DATA D/1.57449261107098347E01,1.17693950891312499E02,
CS 1 5.37181101862009858E02,1.62138957456669019E03,
CS 2 3.29079923573345963E03,4.36261909014324716E03,
CS 3 3.43936767414372164E03,1.23033935480374942E03/
CD DATA C/5.64188496988670089D-1,8.88314979438837594D0,
CD 1 6.61191906371416295D01,2.98635138197400131D02,
CD 2 8.81952221241769090D02,1.71204761263407058D03,
CD 3 2.05107837782607147D03,1.23033935479799725D03,
CD 4 2.15311535474403846D-8/
CD DATA D/1.57449261107098347D01,1.17693950891312499D02,
CD 1 5.37181101862009858D02,1.62138957456669019D03,
CD 2 3.29079923573345963D03,4.36261909014324716D03,
CD 3 3.43936767414372164D03,1.23033935480374942D03/
C------------------------------------------------------------------
C Coefficients for approximation to erfc in third interval
C------------------------------------------------------------------
CS DATA P/3.05326634961232344E-1,3.60344899949804439E-1,
CS 1 1.25781726111229246E-1,1.60837851487422766E-2,
CS 2 6.58749161529837803E-4,1.63153871373020978E-2/
CS DATA Q/2.56852019228982242E00,1.87295284992346047E00,
CS 1 5.27905102951428412E-1,6.05183413124413191E-2,
CS 2 2.33520497626869185E-3/
CD DATA P/3.05326634961232344D-1,3.60344899949804439D-1,
CD 1 1.25781726111229246D-1,1.60837851487422766D-2,
CD 2 6.58749161529837803D-4,1.63153871373020978D-2/
CD DATA Q/2.56852019228982242D00,1.87295284992346047D00,
CD 1 5.27905102951428412D-1,6.05183413124413191D-2,
CD 2 2.33520497626869185D-3/
C------------------------------------------------------------------
X = ARG
Y = ABS(X)
IF (Y .LE. THRESH) THEN
C------------------------------------------------------------------
C Evaluate erf for |X| <= 0.46875
C------------------------------------------------------------------
YSQ = ZERO
IF (Y .GT. XSMALL) YSQ = Y * Y
XNUM = A(5)*YSQ
XDEN = YSQ
DO 20 I = 1, 3
XNUM = (XNUM + A(I)) * YSQ
XDEN = (XDEN + B(I)) * YSQ
20 CONTINUE
RESULT = X * (XNUM + A(4)) / (XDEN + B(4))
IF (JINT .NE. 0) RESULT = ONE - RESULT
IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT
GO TO 800
C------------------------------------------------------------------
C Evaluate erfc for 0.46875 <= |X| <= 4.0
C------------------------------------------------------------------
ELSE IF (Y .LE. FOUR) THEN
XNUM = C(9)*Y
XDEN = Y
DO 120 I = 1, 7
XNUM = (XNUM + C(I)) * Y
XDEN = (XDEN + D(I)) * Y
120 CONTINUE
RESULT = (XNUM + C(8)) / (XDEN + D(8))
IF (JINT .NE. 2) THEN
YSQ = AINT(Y*SIXTEN)/SIXTEN
DEL = (Y-YSQ)*(Y+YSQ)
RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
END IF
C------------------------------------------------------------------
C Evaluate erfc for |X| > 4.0
C------------------------------------------------------------------
ELSE
RESULT = ZERO
IF (Y .GE. XBIG) THEN
IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 300
IF (Y .GE. XHUGE) THEN
RESULT = SQRPI / Y
GO TO 300
END IF
END IF
YSQ = ONE / (Y * Y)
XNUM = P(6)*YSQ
XDEN = YSQ
DO 240 I = 1, 4
XNUM = (XNUM + P(I)) * YSQ
XDEN = (XDEN + Q(I)) * YSQ
240 CONTINUE
RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5))
RESULT = (SQRPI - RESULT) / Y
IF (JINT .NE. 2) THEN
YSQ = AINT(Y*SIXTEN)/SIXTEN
DEL = (Y-YSQ)*(Y+YSQ)
RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT
END IF
END IF
C------------------------------------------------------------------
C Fix up for negative argument, erf, etc.
C------------------------------------------------------------------
300 IF (JINT .EQ. 0) THEN
RESULT = (HALF - RESULT) + HALF
IF (X .LT. ZERO) RESULT = -RESULT
ELSE IF (JINT .EQ. 1) THEN
IF (X .LT. ZERO) RESULT = TWO - RESULT
ELSE
IF (X .LT. ZERO) THEN
IF (X .LT. XNEG) THEN
RESULT = XINF
ELSE
YSQ = AINT(X*SIXTEN)/SIXTEN
DEL = (X-YSQ)*(X+YSQ)
Y = EXP(YSQ*YSQ) * EXP(DEL)
RESULT = (Y+Y) - RESULT
END IF
END IF
END IF
800 RETURN
C---------- Last card of CALERF ----------
END
CS REAL FUNCTION ERF(X)
CD DOUBLE PRECISION FUNCTION DERF(X)
C--------------------------------------------------------------------
C
C This subprogram computes approximate values for erf(x).
C (see comments heading CALERF).
C
C Author/date: W. J. Cody, January 8, 1985
C
C--------------------------------------------------------------------
INTEGER JINT
CS REAL X, RESULT
CD DOUBLE PRECISION X, RESULT
C------------------------------------------------------------------
JINT = 0
CALL CALERF(X,RESULT,JINT)
CS ERF = RESULT
CD DERF = RESULT
RETURN
C---------- Last card of DERF ----------
END
CS REAL FUNCTION ERFC(X)
CD DOUBLE PRECISION FUNCTION DERFC(X)
C--------------------------------------------------------------------
C
C This subprogram computes approximate values for erfc(x).
C (see comments heading CALERF).
C
C Author/date: W. J. Cody, January 8, 1985
C
C--------------------------------------------------------------------
INTEGER JINT
CS REAL X, RESULT
CD DOUBLE PRECISION X, RESULT
C------------------------------------------------------------------
JINT = 1
CALL CALERF(X,RESULT,JINT)
CS ERFC = RESULT
CD DERFC = RESULT
RETURN
C---------- Last card of DERFC ----------
END
CS REAL FUNCTION ERFCX(X)
CD DOUBLE PRECISION FUNCTION DERFCX(X)
C------------------------------------------------------------------
C
C This subprogram computes approximate values for exp(x*x) * erfc(x).
C (see comments heading CALERF).
C
C Author/date: W. J. Cody, March 30, 1987
C
C------------------------------------------------------------------
INTEGER JINT
CS REAL X, RESULT
CD DOUBLE PRECISION X, RESULT
C------------------------------------------------------------------
JINT = 2
CALL CALERF(X,RESULT,JINT)
CS ERFCX = RESULT
CD DERFCX = RESULT
RETURN
C---------- Last card of DERFCX ----------
END