00001 SUBROUTINE CALERF(ARG,RESULT,JINT) 00002 C------------------------------------------------------------------ 00003 C 00004 C This packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x) 00005 C for a real argument x. It contains three FUNCTION type 00006 C subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX), 00007 C and one SUBROUTINE type subprogram, CALERF. The calling 00008 C statements for the primary entries are: 00009 C 00010 C Y=ERF(X) (or Y=DERF(X)), 00011 C 00012 C Y=ERFC(X) (or Y=DERFC(X)), 00013 C and 00014 C Y=ERFCX(X) (or Y=DERFCX(X)). 00015 C 00016 C The routine CALERF is intended for internal packet use only, 00017 C all computations within the packet being concentrated in this 00018 C routine. The function subprograms invoke CALERF with the 00019 C statement 00020 C 00021 C CALL CALERF(ARG,RESULT,JINT) 00022 C 00023 C where the parameter usage is as follows 00024 C 00025 C Function Parameters for CALERF 00026 C call ARG Result JINT 00027 C 00028 C ERF(ARG) ANY REAL ARGUMENT ERF(ARG) 0 00029 C ERFC(ARG) ABS(ARG) .LT. XBIG ERFC(ARG) 1 00030 C ERFCX(ARG) XNEG .LT. ARG .LT. XMAX ERFCX(ARG) 2 00031 C 00032 C The main computation evaluates near-minimax approximations 00033 C from "Rational Chebyshev approximations for the error function" 00034 C by W. J. Cody, Math. Comp., 1969, PP. 631-638. This 00035 C transportable program uses rational functions that theoretically 00036 C approximate erf(x) and erfc(x) to at least 18 significant 00037 C decimal digits. The accuracy achieved depends on the arithmetic 00038 C system, the compiler, the intrinsic functions, and proper 00039 C selection of the machine-dependent constants. 00040 C 00041 C******************************************************************* 00042 C******************************************************************* 00043 C 00044 C Explanation of machine-dependent constants 00045 C 00046 C XMIN = the smallest positive floating-point number. 00047 C XINF = the largest positive finite floating-point number. 00048 C XNEG = the largest negative argument acceptable to ERFCX; 00049 C the negative of the solution to the equation 00050 C 2*exp(x*x) = XINF. 00051 C XSMALL = argument below which erf(x) may be represented by 00052 C 2*x/sqrt(pi) and above which x*x will not underflow. 00053 C A conservative value is the largest machine number X 00054 C such that 1.0 + X = 1.0 to machine precision. 00055 C XBIG = largest argument acceptable to ERFC; solution to 00056 C the equation: W(x) * (1-0.5/x**2) = XMIN, where 00057 C W(x) = exp(-x*x)/[x*sqrt(pi)]. 00058 C XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to 00059 C machine precision. A conservative value is 00060 C 1/[2*sqrt(XSMALL)] 00061 C XMAX = largest acceptable argument to ERFCX; the minimum 00062 C of XINF and 1/[sqrt(pi)*XMIN]. 00063 C 00064 C Approximate values for some important machines are: 00065 C 00066 C XMIN XINF XNEG XSMALL 00067 C 00068 C CDC 7600 (S.P.) 3.13E-294 1.26E+322 -27.220 7.11E-15 00069 C CRAY-1 (S.P.) 4.58E-2467 5.45E+2465 -75.345 7.11E-15 00070 C IEEE (IBM/XT, 00071 C SUN, etc.) (S.P.) 1.18E-38 3.40E+38 -9.382 5.96E-8 00072 C IEEE (IBM/XT, 00073 C SUN, etc.) (D.P.) 2.23D-308 1.79D+308 -26.628 1.11D-16 00074 C IBM 195 (D.P.) 5.40D-79 7.23E+75 -13.190 1.39D-17 00075 C UNIVAC 1108 (D.P.) 2.78D-309 8.98D+307 -26.615 1.73D-18 00076 C VAX D-Format (D.P.) 2.94D-39 1.70D+38 -9.345 1.39D-17 00077 C VAX G-Format (D.P.) 5.56D-309 8.98D+307 -26.615 1.11D-16 00078 C 00079 C 00080 C XBIG XHUGE XMAX 00081 C 00082 C CDC 7600 (S.P.) 25.922 8.39E+6 1.80X+293 00083 C CRAY-1 (S.P.) 75.326 8.39E+6 5.45E+2465 00084 C IEEE (IBM/XT, 00085 C SUN, etc.) (S.P.) 9.194 2.90E+3 4.79E+37 00086 C IEEE (IBM/XT, 00087 C SUN, etc.) (D.P.) 26.543 6.71D+7 2.53D+307 00088 C IBM 195 (D.P.) 13.306 1.90D+8 7.23E+75 00089 C UNIVAC 1108 (D.P.) 26.582 5.37D+8 8.98D+307 00090 C VAX D-Format (D.P.) 9.269 1.90D+8 1.70D+38 00091 C VAX G-Format (D.P.) 26.569 6.71D+7 8.98D+307 00092 C 00093 C******************************************************************* 00094 C******************************************************************* 00095 C 00096 C Error returns 00097 C 00098 C The program returns ERFC = 0 for ARG .GE. XBIG; 00099 C 00100 C ERFCX = XINF for ARG .LT. XNEG; 00101 C and 00102 C ERFCX = 0 for ARG .GE. XMAX. 00103 C 00104 C 00105 C Intrinsic functions required are: 00106 C 00107 C ABS, AINT, EXP 00108 C 00109 C 00110 C Author: W. J. Cody 00111 C Mathematics and Computer Science Division 00112 C Argonne National Laboratory 00113 C Argonne, IL 60439 00114 C 00115 C Latest modification: March 19, 1990 00116 C 00117 C------------------------------------------------------------------ 00118 INTEGER I,JINT 00119 CS REAL 00120 CD DOUBLE PRECISION 00121 1 A,ARG,B,C,D,DEL,FOUR,HALF,P,ONE,Q,RESULT,SIXTEN,SQRPI, 00122 2 TWO,THRESH,X,XBIG,XDEN,XHUGE,XINF,XMAX,XNEG,XNUM,XSMALL, 00123 3 Y,YSQ,ZERO 00124 DIMENSION A(5),B(4),C(9),D(8),P(6),Q(5) 00125 C------------------------------------------------------------------ 00126 C Mathematical constants 00127 C------------------------------------------------------------------ 00128 CS DATA FOUR,ONE,HALF,TWO,ZERO/4.0E0,1.0E0,0.5E0,2.0E0,0.0E0/, 00129 CS 1 SQRPI/5.6418958354775628695E-1/,THRESH/0.46875E0/, 00130 CS 2 SIXTEN/16.0E0/ 00131 CD DATA FOUR,ONE,HALF,TWO,ZERO/4.0D0,1.0D0,0.5D0,2.0D0,0.0D0/, 00132 CD 1 SQRPI/5.6418958354775628695D-1/,THRESH/0.46875D0/, 00133 CD 2 SIXTEN/16.0D0/ 00134 C------------------------------------------------------------------ 00135 C Machine-dependent constants 00136 C------------------------------------------------------------------ 00137 CS DATA XINF,XNEG,XSMALL/3.40E+38,-9.382E0,5.96E-8/, 00138 CS 1 XBIG,XHUGE,XMAX/9.194E0,2.90E3,4.79E37/ 00139 CD DATA XINF,XNEG,XSMALL/1.79D308,-26.628D0,1.11D-16/, 00140 CD 1 XBIG,XHUGE,XMAX/26.543D0,6.71D7,2.53D307/ 00141 C------------------------------------------------------------------ 00142 C Coefficients for approximation to erf in first interval 00143 C------------------------------------------------------------------ 00144 CS DATA A/3.16112374387056560E00,1.13864154151050156E02, 00145 CS 1 3.77485237685302021E02,3.20937758913846947E03, 00146 CS 2 1.85777706184603153E-1/ 00147 CS DATA B/2.36012909523441209E01,2.44024637934444173E02, 00148 CS 1 1.28261652607737228E03,2.84423683343917062E03/ 00149 CD DATA A/3.16112374387056560D00,1.13864154151050156D02, 00150 CD 1 3.77485237685302021D02,3.20937758913846947D03, 00151 CD 2 1.85777706184603153D-1/ 00152 CD DATA B/2.36012909523441209D01,2.44024637934444173D02, 00153 CD 1 1.28261652607737228D03,2.84423683343917062D03/ 00154 C------------------------------------------------------------------ 00155 C Coefficients for approximation to erfc in second interval 00156 C------------------------------------------------------------------ 00157 CS DATA C/5.64188496988670089E-1,8.88314979438837594E0, 00158 CS 1 6.61191906371416295E01,2.98635138197400131E02, 00159 CS 2 8.81952221241769090E02,1.71204761263407058E03, 00160 CS 3 2.05107837782607147E03,1.23033935479799725E03, 00161 CS 4 2.15311535474403846E-8/ 00162 CS DATA D/1.57449261107098347E01,1.17693950891312499E02, 00163 CS 1 5.37181101862009858E02,1.62138957456669019E03, 00164 CS 2 3.29079923573345963E03,4.36261909014324716E03, 00165 CS 3 3.43936767414372164E03,1.23033935480374942E03/ 00166 CD DATA C/5.64188496988670089D-1,8.88314979438837594D0, 00167 CD 1 6.61191906371416295D01,2.98635138197400131D02, 00168 CD 2 8.81952221241769090D02,1.71204761263407058D03, 00169 CD 3 2.05107837782607147D03,1.23033935479799725D03, 00170 CD 4 2.15311535474403846D-8/ 00171 CD DATA D/1.57449261107098347D01,1.17693950891312499D02, 00172 CD 1 5.37181101862009858D02,1.62138957456669019D03, 00173 CD 2 3.29079923573345963D03,4.36261909014324716D03, 00174 CD 3 3.43936767414372164D03,1.23033935480374942D03/ 00175 C------------------------------------------------------------------ 00176 C Coefficients for approximation to erfc in third interval 00177 C------------------------------------------------------------------ 00178 CS DATA P/3.05326634961232344E-1,3.60344899949804439E-1, 00179 CS 1 1.25781726111229246E-1,1.60837851487422766E-2, 00180 CS 2 6.58749161529837803E-4,1.63153871373020978E-2/ 00181 CS DATA Q/2.56852019228982242E00,1.87295284992346047E00, 00182 CS 1 5.27905102951428412E-1,6.05183413124413191E-2, 00183 CS 2 2.33520497626869185E-3/ 00184 CD DATA P/3.05326634961232344D-1,3.60344899949804439D-1, 00185 CD 1 1.25781726111229246D-1,1.60837851487422766D-2, 00186 CD 2 6.58749161529837803D-4,1.63153871373020978D-2/ 00187 CD DATA Q/2.56852019228982242D00,1.87295284992346047D00, 00188 CD 1 5.27905102951428412D-1,6.05183413124413191D-2, 00189 CD 2 2.33520497626869185D-3/ 00190 C------------------------------------------------------------------ 00191 X = ARG 00192 Y = ABS(X) 00193 IF (Y .LE. THRESH) THEN 00194 C------------------------------------------------------------------ 00195 C Evaluate erf for |X| <= 0.46875 00196 C------------------------------------------------------------------ 00197 YSQ = ZERO 00198 IF (Y .GT. XSMALL) YSQ = Y * Y 00199 XNUM = A(5)*YSQ 00200 XDEN = YSQ 00201 DO 20 I = 1, 3 00202 XNUM = (XNUM + A(I)) * YSQ 00203 XDEN = (XDEN + B(I)) * YSQ 00204 20 CONTINUE 00205 RESULT = X * (XNUM + A(4)) / (XDEN + B(4)) 00206 IF (JINT .NE. 0) RESULT = ONE - RESULT 00207 IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT 00208 GO TO 800 00209 C------------------------------------------------------------------ 00210 C Evaluate erfc for 0.46875 <= |X| <= 4.0 00211 C------------------------------------------------------------------ 00212 ELSE IF (Y .LE. FOUR) THEN 00213 XNUM = C(9)*Y 00214 XDEN = Y 00215 DO 120 I = 1, 7 00216 XNUM = (XNUM + C(I)) * Y 00217 XDEN = (XDEN + D(I)) * Y 00218 120 CONTINUE 00219 RESULT = (XNUM + C(8)) / (XDEN + D(8)) 00220 IF (JINT .NE. 2) THEN 00221 YSQ = AINT(Y*SIXTEN)/SIXTEN 00222 DEL = (Y-YSQ)*(Y+YSQ) 00223 RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT 00224 END IF 00225 C------------------------------------------------------------------ 00226 C Evaluate erfc for |X| > 4.0 00227 C------------------------------------------------------------------ 00228 ELSE 00229 RESULT = ZERO 00230 IF (Y .GE. XBIG) THEN 00231 IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 300 00232 IF (Y .GE. XHUGE) THEN 00233 RESULT = SQRPI / Y 00234 GO TO 300 00235 END IF 00236 END IF 00237 YSQ = ONE / (Y * Y) 00238 XNUM = P(6)*YSQ 00239 XDEN = YSQ 00240 DO 240 I = 1, 4 00241 XNUM = (XNUM + P(I)) * YSQ 00242 XDEN = (XDEN + Q(I)) * YSQ 00243 240 CONTINUE 00244 RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5)) 00245 RESULT = (SQRPI - RESULT) / Y 00246 IF (JINT .NE. 2) THEN 00247 YSQ = AINT(Y*SIXTEN)/SIXTEN 00248 DEL = (Y-YSQ)*(Y+YSQ) 00249 RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT 00250 END IF 00251 END IF 00252 C------------------------------------------------------------------ 00253 C Fix up for negative argument, erf, etc. 00254 C------------------------------------------------------------------ 00255 300 IF (JINT .EQ. 0) THEN 00256 RESULT = (HALF - RESULT) + HALF 00257 IF (X .LT. ZERO) RESULT = -RESULT 00258 ELSE IF (JINT .EQ. 1) THEN 00259 IF (X .LT. ZERO) RESULT = TWO - RESULT 00260 ELSE 00261 IF (X .LT. ZERO) THEN 00262 IF (X .LT. XNEG) THEN 00263 RESULT = XINF 00264 ELSE 00265 YSQ = AINT(X*SIXTEN)/SIXTEN 00266 DEL = (X-YSQ)*(X+YSQ) 00267 Y = EXP(YSQ*YSQ) * EXP(DEL) 00268 RESULT = (Y+Y) - RESULT 00269 END IF 00270 END IF 00271 END IF 00272 800 RETURN 00273 C---------- Last card of CALERF ---------- 00274 END 00275 CS REAL FUNCTION ERF(X) 00276 CD DOUBLE PRECISION FUNCTION DERF(X) 00277 C-------------------------------------------------------------------- 00278 C 00279 C This subprogram computes approximate values for erf(x). 00280 C (see comments heading CALERF). 00281 C 00282 C Author/date: W. J. Cody, January 8, 1985 00283 C 00284 C-------------------------------------------------------------------- 00285 INTEGER JINT 00286 CS REAL X, RESULT 00287 CD DOUBLE PRECISION X, RESULT 00288 C------------------------------------------------------------------ 00289 JINT = 0 00290 CALL CALERF(X,RESULT,JINT) 00291 CS ERF = RESULT 00292 CD DERF = RESULT 00293 RETURN 00294 C---------- Last card of DERF ---------- 00295 END 00296 CS REAL FUNCTION ERFC(X) 00297 CD DOUBLE PRECISION FUNCTION DERFC(X) 00298 C-------------------------------------------------------------------- 00299 C 00300 C This subprogram computes approximate values for erfc(x). 00301 C (see comments heading CALERF). 00302 C 00303 C Author/date: W. J. Cody, January 8, 1985 00304 C 00305 C-------------------------------------------------------------------- 00306 INTEGER JINT 00307 CS REAL X, RESULT 00308 CD DOUBLE PRECISION X, RESULT 00309 C------------------------------------------------------------------ 00310 JINT = 1 00311 CALL CALERF(X,RESULT,JINT) 00312 CS ERFC = RESULT 00313 CD DERFC = RESULT 00314 RETURN 00315 C---------- Last card of DERFC ---------- 00316 END 00317 CS REAL FUNCTION ERFCX(X) 00318 CD DOUBLE PRECISION FUNCTION DERFCX(X) 00319 C------------------------------------------------------------------ 00320 C 00321 C This subprogram computes approximate values for exp(x*x) * erfc(x). 00322 C (see comments heading CALERF). 00323 C 00324 C Author/date: W. J. Cody, March 30, 1987 00325 C 00326 C------------------------------------------------------------------ 00327 INTEGER JINT 00328 CS REAL X, RESULT 00329 CD DOUBLE PRECISION X, RESULT 00330 C------------------------------------------------------------------ 00331 JINT = 2 00332 CALL CALERF(X,RESULT,JINT) 00333 CS ERFCX = RESULT 00334 CD DERFCX = RESULT 00335 RETURN 00336 C---------- Last card of DERFCX ---------- 00337 END