LetsBeRational-1.0.0.0: external/src/erf_cody.cpp
//
// Original Fortran code taken from http://www.netlib.org/specfun/erf, compiled with f2c, and adapted by hand.
//
// Created with command line f2c -C++ -c -a -krd -r8 cody_erf.f
//
// Translated by f2c (version 20100827).
//
//
// This source code resides at www.jaeckel.org/LetsBeRational.7z .
//
// ======================================================================================
// WARRANTY DISCLAIMER
// The Software is provided "as is" without warranty of any kind, either express or implied,
// including without limitation any implied warranties of condition, uninterrupted use,
// merchantability, fitness for a particular purpose, or non-infringement.
// ======================================================================================
//
#if defined( _DEBUG ) || defined( BOUNDS_CHECK_STL_ARRAYS )
#define _SECURE_SCL 1
#define _SECURE_SCL_THROWS 1
#define _SCL_SECURE_NO_WARNINGS
#define _HAS_ITERATOR_DEBUGGING 0
#else
#define _SECURE_SCL 0
#endif
#if defined(_MSC_VER)
# define NOMINMAX // to suppress MSVC's definitions of min() and max()
// These four pragmas are the equivalent to /fp:fast.
# pragma float_control( except, off )
# pragma float_control( precise, off )
# pragma fp_contract( on )
# pragma fenv_access( off )
#endif
#include "normaldistribution.h"
#include <math.h>
#include <float.h>
namespace {
inline double d_int(const double x){ return( (x>0) ? floor(x) : -floor(-x) ); }
}
/*< SUBROUTINE CALERF(ARG,RESULT,JINT) >*/
double calerf(double x, const int jint) {
static const double a[5] = { 3.1611237438705656,113.864154151050156,377.485237685302021,3209.37758913846947,.185777706184603153 };
static const double b[4] = { 23.6012909523441209,244.024637934444173,1282.61652607737228,2844.23683343917062 };
static const double c__[9] = { .564188496988670089,8.88314979438837594,66.1191906371416295,298.635138197400131,881.95222124176909,1712.04761263407058,2051.07837782607147,1230.33935479799725,2.15311535474403846e-8 };
static const double d__[8] = { 15.7449261107098347,117.693950891312499,537.181101862009858,1621.38957456669019,3290.79923573345963,4362.61909014324716,3439.36767414372164,1230.33935480374942 };
static const double p[6] = { .305326634961232344,.360344899949804439,.125781726111229246,.0160837851487422766,6.58749161529837803e-4,.0163153871373020978 };
static const double q[5] = { 2.56852019228982242,1.87295284992346047,.527905102951428412,.0605183413124413191,.00233520497626869185 };
static const double zero = 0.;
static const double half = .5;
static const double one = 1.;
static const double two = 2.;
static const double four = 4.;
static const double sqrpi = 0.56418958354775628695;
static const double thresh = .46875;
static const double sixten = 16.;
double y, del, ysq, xden, xnum, result;
/* ------------------------------------------------------------------ */
/* 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) ABS(ARG) .LT. XBIG ERFC(ARG) 1 */
/* ERFCX(ARG) XNEG .LT. ARG .LT. 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]. */
// The numbers below were preselected for IEEE .
static const double xinf = 1.79e308;
static const double xneg = -26.628;
static const double xsmall = 1.11e-16;
static const double xbig = 26.543;
static const double xhuge = 6.71e7;
static const double xmax = 2.53e307;
/* 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 .GE. XBIG; */
/* ERFCX = XINF for ARG .LT. XNEG; */
/* and */
/* ERFCX = 0 for ARG .GE. XMAX. */
/* Intrinsic functions required are: */
/* ABS, AINT, EXP */
/* Author: W. J. Cody */
/* Mathematics and Computer Science Division */
/* Argonne National Laboratory */
/* Argonne, IL 60439 */
/* Latest modification: March 19, 1990 */
/* ------------------------------------------------------------------ */
/*< INTEGER I,JINT >*/
/* S REAL */
/*< >*/
/*< DIMENSION A(5),B(4),C(9),D(8),P(6),Q(5) >*/
/* ------------------------------------------------------------------ */
/* Mathematical constants */
/* ------------------------------------------------------------------ */
/* S DATA FOUR,ONE,HALF,TWO,ZERO/4.0E0,1.0E0,0.5E0,2.0E0,0.0E0/, */
/* S 1 SQRPI/5.6418958354775628695E-1/,THRESH/0.46875E0/, */
/* S 2 SIXTEN/16.0E0/ */
/*< >*/
/* ------------------------------------------------------------------ */
/* Machine-dependent constants */
/* ------------------------------------------------------------------ */
/* S DATA XINF,XNEG,XSMALL/3.40E+38,-9.382E0,5.96E-8/, */
/* S 1 XBIG,XHUGE,XMAX/9.194E0,2.90E3,4.79E37/ */
/*< >*/
/* ------------------------------------------------------------------ */
/* Coefficients for approximation to erf in first interval */
/* ------------------------------------------------------------------ */
/* S DATA A/3.16112374387056560E00,1.13864154151050156E02, */
/* S 1 3.77485237685302021E02,3.20937758913846947E03, */
/* S 2 1.85777706184603153E-1/ */
/* S DATA B/2.36012909523441209E01,2.44024637934444173E02, */
/* S 1 1.28261652607737228E03,2.84423683343917062E03/ */
/*< >*/
/*< >*/
/* ------------------------------------------------------------------ */
/* Coefficients for approximation to erfc in second interval */
/* ------------------------------------------------------------------ */
/* S DATA C/5.64188496988670089E-1,8.88314979438837594E0, */
/* S 1 6.61191906371416295E01,2.98635138197400131E02, */
/* S 2 8.81952221241769090E02,1.71204761263407058E03, */
/* S 3 2.05107837782607147E03,1.23033935479799725E03, */
/* S 4 2.15311535474403846E-8/ */
/* S DATA D/1.57449261107098347E01,1.17693950891312499E02, */
/* S 1 5.37181101862009858E02,1.62138957456669019E03, */
/* S 2 3.29079923573345963E03,4.36261909014324716E03, */
/* S 3 3.43936767414372164E03,1.23033935480374942E03/ */
/*< >*/
/*< >*/
/* ------------------------------------------------------------------ */
/* Coefficients for approximation to erfc in third interval */
/* ------------------------------------------------------------------ */
/* S DATA P/3.05326634961232344E-1,3.60344899949804439E-1, */
/* S 1 1.25781726111229246E-1,1.60837851487422766E-2, */
/* S 2 6.58749161529837803E-4,1.63153871373020978E-2/ */
/* S DATA Q/2.56852019228982242E00,1.87295284992346047E00, */
/* S 1 5.27905102951428412E-1,6.05183413124413191E-2, */
/* S 2 2.33520497626869185E-3/ */
/*< >*/
/*< >*/
/* ------------------------------------------------------------------ */
/*< X = ARG >*/
// x = *arg;
/*< Y = ABS(X) >*/
y = fabs(x);
/*< IF (Y .LE. THRESH) THEN >*/
if (y <= thresh) {
/* ------------------------------------------------------------------ */
/* Evaluate erf for |X| <= 0.46875 */
/* ------------------------------------------------------------------ */
/*< YSQ = ZERO >*/
ysq = zero;
/*< IF (Y .GT. XSMALL) YSQ = Y * Y >*/
if (y > xsmall) {
ysq = y * y;
}
/*< XNUM = A(5)*YSQ >*/
xnum = a[4] * ysq;
/*< XDEN = YSQ >*/
xden = ysq;
/*< DO 20 I = 1, 3 >*/
for (int i__ = 1; i__ <= 3; ++i__) {
/*< XNUM = (XNUM + A(I)) * YSQ >*/
xnum = (xnum + a[i__ - 1]) * ysq;
/*< XDEN = (XDEN + B(I)) * YSQ >*/
xden = (xden + b[i__ - 1]) * ysq;
/*< 20 CONTINUE >*/
/* L20: */
}
/*< RESULT = X * (XNUM + A(4)) / (XDEN + B(4)) >*/
result = x * (xnum + a[3]) / (xden + b[3]);
/*< IF (JINT .NE. 0) RESULT = ONE - RESULT >*/
if (jint != 0) {
result = one - result;
}
/*< IF (JINT .EQ. 2) RESULT = EXP(YSQ) * RESULT >*/
if (jint == 2) {
result = exp(ysq) * result;
}
/*< GO TO 800 >*/
goto L800;
/* ------------------------------------------------------------------ */
/* Evaluate erfc for 0.46875 <= |X| <= 4.0 */
/* ------------------------------------------------------------------ */
/*< ELSE IF (Y .LE. FOUR) THEN >*/
} else if (y <= four) {
/*< XNUM = C(9)*Y >*/
xnum = c__[8] * y;
/*< XDEN = Y >*/
xden = y;
/*< DO 120 I = 1, 7 >*/
for (int i__ = 1; i__ <= 7; ++i__) {
/*< XNUM = (XNUM + C(I)) * Y >*/
xnum = (xnum + c__[i__ - 1]) * y;
/*< XDEN = (XDEN + D(I)) * Y >*/
xden = (xden + d__[i__ - 1]) * y;
/*< 120 CONTINUE >*/
/* L120: */
}
/*< RESULT = (XNUM + C(8)) / (XDEN + D(8)) >*/
result = (xnum + c__[7]) / (xden + d__[7]);
/*< IF (JINT .NE. 2) THEN >*/
if (jint != 2) {
/*< YSQ = AINT(Y*SIXTEN)/SIXTEN >*/
double d__1 = y * sixten;
ysq = d_int(d__1) / sixten;
/*< DEL = (Y-YSQ)*(Y+YSQ) >*/
del = (y - ysq) * (y + ysq);
/*< RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT >*/
d__1 = exp(-ysq * ysq) * exp(-del);
result = d__1 * result;
/*< END IF >*/
}
/* ------------------------------------------------------------------ */
/* Evaluate erfc for |X| > 4.0 */
/* ------------------------------------------------------------------ */
/*< ELSE >*/
} else {
/*< RESULT = ZERO >*/
result = zero;
/*< IF (Y .GE. XBIG) THEN >*/
if (y >= xbig) {
/*< IF ((JINT .NE. 2) .OR. (Y .GE. XMAX)) GO TO 300 >*/
if (jint != 2 || y >= xmax) {
goto L300;
}
/*< IF (Y .GE. XHUGE) THEN >*/
if (y >= xhuge) {
/*< RESULT = SQRPI / Y >*/
result = sqrpi / y;
/*< GO TO 300 >*/
goto L300;
/*< END IF >*/
}
/*< END IF >*/
}
/*< YSQ = ONE / (Y * Y) >*/
ysq = one / (y * y);
/*< XNUM = P(6)*YSQ >*/
xnum = p[5] * ysq;
/*< XDEN = YSQ >*/
xden = ysq;
/*< DO 240 I = 1, 4 >*/
for (int i__ = 1; i__ <= 4; ++i__) {
/*< XNUM = (XNUM + P(I)) * YSQ >*/
xnum = (xnum + p[i__ - 1]) * ysq;
/*< XDEN = (XDEN + Q(I)) * YSQ >*/
xden = (xden + q[i__ - 1]) * ysq;
/*< 240 CONTINUE >*/
/* L240: */
}
/*< RESULT = YSQ *(XNUM + P(5)) / (XDEN + Q(5)) >*/
result = ysq * (xnum + p[4]) / (xden + q[4]);
/*< RESULT = (SQRPI - RESULT) / Y >*/
result = (sqrpi - result) / y;
/*< IF (JINT .NE. 2) THEN >*/
if (jint != 2) {
/*< YSQ = AINT(Y*SIXTEN)/SIXTEN >*/
double d__1 = y * sixten;
ysq = d_int(d__1) / sixten;
/*< DEL = (Y-YSQ)*(Y+YSQ) >*/
del = (y - ysq) * (y + ysq);
/*< RESULT = EXP(-YSQ*YSQ) * EXP(-DEL) * RESULT >*/
d__1 = exp(-ysq * ysq) * exp(-del);
result = d__1 * result;
/*< END IF >*/
}
/*< END IF >*/
}
/* ------------------------------------------------------------------ */
/* Fix up for negative argument, erf, etc. */
/* ------------------------------------------------------------------ */
/*< 300 IF (JINT .EQ. 0) THEN >*/
L300:
if (jint == 0) {
/*< RESULT = (HALF - RESULT) + HALF >*/
result = (half - result) + half;
/*< IF (X .LT. ZERO) RESULT = -RESULT >*/
if (x < zero) {
result = -(result);
}
/*< ELSE IF (JINT .EQ. 1) THEN >*/
} else if (jint == 1) {
/*< IF (X .LT. ZERO) RESULT = TWO - RESULT >*/
if (x < zero) {
result = two - result;
}
/*< ELSE >*/
} else {
/*< IF (X .LT. ZERO) THEN >*/
if (x < zero) {
/*< IF (X .LT. XNEG) THEN >*/
if (x < xneg) {
/*< RESULT = XINF >*/
result = xinf;
/*< ELSE >*/
} else {
/*< YSQ = AINT(X*SIXTEN)/SIXTEN >*/
double d__1 = x * sixten;
ysq = d_int(d__1) / sixten;
/*< DEL = (X-YSQ)*(X+YSQ) >*/
del = (x - ysq) * (x + ysq);
/*< Y = EXP(YSQ*YSQ) * EXP(DEL) >*/
y = exp(ysq * ysq) * exp(del);
/*< RESULT = (Y+Y) - RESULT >*/
result = y + y - result;
/*< END IF >*/
}
/*< END IF >*/
}
/*< END IF >*/
}
/*< 800 RETURN >*/
L800:
return result;
/* ---------- Last card of CALERF ---------- */
/*< END >*/
} /* calerf_ */
/* S REAL FUNCTION ERF(X) */
/*< DOUBLE PRECISION FUNCTION DERF(X) >*/
double erf_cody(double x){
/* -------------------------------------------------------------------- */
/* This subprogram computes approximate values for erf(x). */
/* (see comments heading CALERF). */
/* Author/date: W. J. Cody, January 8, 1985 */
/* -------------------------------------------------------------------- */
/*< INTEGER JINT >*/
/* S REAL X, RESULT */
/*< DOUBLE PRECISION X, RESULT >*/
/* ------------------------------------------------------------------ */
/*< JINT = 0 >*/
/*< CALL CALERF(X,RESULT,JINT) >*/
return calerf(x, 0);
/* S ERF = RESULT */
/*< DERF = RESULT >*/
/*< RETURN >*/
/* ---------- Last card of DERF ---------- */
/*< END >*/
} /* derf_ */
/* S REAL FUNCTION ERFC(X) */
/*< DOUBLE PRECISION FUNCTION DERFC(X) >*/
double erfc_cody(double x) {
/* -------------------------------------------------------------------- */
/* This subprogram computes approximate values for erfc(x). */
/* (see comments heading CALERF). */
/* Author/date: W. J. Cody, January 8, 1985 */
/* -------------------------------------------------------------------- */
/*< INTEGER JINT >*/
/* S REAL X, RESULT */
/*< DOUBLE PRECISION X, RESULT >*/
/* ------------------------------------------------------------------ */
/*< JINT = 1 >*/
/*< CALL CALERF(X,RESULT,JINT) >*/
return calerf(x, 1);
/* S ERFC = RESULT */
/*< DERFC = RESULT >*/
/*< RETURN >*/
/* ---------- Last card of DERFC ---------- */
/*< END >*/
} /* derfc_ */
/* S REAL FUNCTION ERFCX(X) */
/*< DOUBLE PRECISION FUNCTION DERFCX(X) >*/
double erfcx_cody(double x) {
/* ------------------------------------------------------------------ */
/* This subprogram computes approximate values for exp(x*x) * erfc(x). */
/* (see comments heading CALERF). */
/* Author/date: W. J. Cody, March 30, 1987 */
/* ------------------------------------------------------------------ */
/*< INTEGER JINT >*/
/* S REAL X, RESULT */
/*< DOUBLE PRECISION X, RESULT >*/
/* ------------------------------------------------------------------ */
/*< JINT = 2 >*/
/*< CALL CALERF(X,RESULT,JINT) >*/
return calerf(x, 2);
/* S ERFCX = RESULT */
/*< DERFCX = RESULT >*/
/*< RETURN >*/
/* ---------- Last card of DERFCX ---------- */
/*< END >*/
} /* derfcx_ */