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