/* Double precision version of routines in ERF from the netlib SPECFUN */
/* library. Translated by f2c and modified. */
#include "xlisp.h"
#include "xlstat.h"
#include "linalg.h"
static VOID calerf P3C(double, arg, double *, result, int, jint)
{
/* ------------------------------------------------------------------ */
/* 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]. */
/* 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 */
/* ------------------------------------------------------------------ */
double xden, xnum;
int i;
double x, y, del, ysq;
/* ------------------------------------------------------------------ */
/* Mathematical constants */
/* ------------------------------------------------------------------ */
static double four = 4.;
static double one = 1.;
static double half = .5;
static double two = 2.;
static double zero = 0.;
static double sqrpi = .56418958354775628695;
static double thresh = .46875;
static double sixten = 16.;
/* ------------------------------------------------------------------ */
/* Machine-dependent constants */
/* ------------------------------------------------------------------ */
#ifdef IEEEFP
static double xinf = 1.79e308;
static double xneg = -26.628;
static double xsmall = 1.11e-16;
static double xbig = 26.543;
static double xhuge = 6.71e7;
static double xmax = 2.53e307;
#else
#ifdef CRAYCC
static double xinf = 5.45e2465;
static double xneg = -75.345;
static double xsmall = 7.11e-15;
static double xbig = 75.326;
static double xhuge = 8.39e6;
static double xmax = 5.45e2465;
#else /* use IBM 196 values */
static double xinf = 7.23e75;
static double xneg = -13.190;
static double xsmall = 1.39e-17;
static double xbig = 13.306;
static double xhuge = 1.90e8;
static double xmax = 7.23e75;
#endif /* CRAYCC */
#endif /* IEEEFP */
/* ------------------------------------------------------------------ */
/* Coefficients for approximation to erf in first interval */
/* ------------------------------------------------------------------ */
static double a[5] = { 3.1611237438705656,113.864154151050156,
377.485237685302021,3209.37758913846947,
.185777706184603153 };
static double b[4] = { 23.6012909523441209,244.024637934444173,
1282.61652607737228,2844.23683343917062 };
/* ------------------------------------------------------------------ */
/* Coefficients for approximation to erfc in second interval */
/* ------------------------------------------------------------------ */
static double c[9] = { .564188496988670089,8.88314979438837594,
66.1191906371416295,298.635138197400131,
881.95222124176909,1712.04761263407058,
2051.07837782607147,1230.33935479799725,
2.15311535474403846e-8 };
static double d[8] = { 15.7449261107098347,117.693950891312499,
537.181101862009858,1621.38957456669019,
3290.79923573345963,4362.61909014324716,
3439.36767414372164,1230.33935480374942 };
/* ------------------------------------------------------------------ */
/* Coefficients for approximation to erfc in third interval */
/* ------------------------------------------------------------------ */
static double p[6] = { .305326634961232344,.360344899949804439,
.125781726111229246,.0160837851487422766,
6.58749161529837803e-4,.0163153871373020978 };
static double q[5] = { 2.56852019228982242,1.87295284992346047,
.527905102951428412,.0605183413124413191,
.00233520497626869185 };
/* ------------------------------------------------------------------ */
x = arg;
y = abs(x);
if (y <= thresh) {
/* ------------------------------------------------------------------ */
/* Evaluate erf for |X| <= 0.46875 */
/* ------------------------------------------------------------------ */
ysq = zero;
if (y > xsmall) {
ysq = y * y;
}
xnum = a[4] * 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[3]) / (xden + b[3]);
if (jint != 0) {
*result = one - *result;
}
if (jint == 2) {
*result = exp(ysq) * *result;
}
return;
/* ------------------------------------------------------------------ */
/* Evaluate erfc for 0.46875 <= |X| <= 4.0 */
/* ------------------------------------------------------------------ */
}
else if (y <= four) {
xnum = c[8] * y;
xden = y;
for (i = 1; i <= 7; ++i) {
xnum = (xnum + c[i - 1]) * y;
xden = (xden + d[i - 1]) * y;
}
*result = (xnum + c[7]) / (xden + d[7]);
if (jint != 2) {
ysq = floor(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[5] * 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[4]) / (xden + q[4]);
*result = (sqrpi - *result) / y;
if (jint != 2) {
ysq = floor(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 = -floor(-(x * sixten)) / sixten;
del = (x - ysq) * (x + ysq);
y = exp(ysq * ysq) * exp(del);
*result = y + y - *result;
}
}
}
return;
}
#ifdef DEBUG
double derf_ P1C(double *, x)
{
int jint;
double result;
/* This subprogram computes approximate values for erf(x). */
/* (see comments heading CALERF). */
/* Author/date: W. J. Cody, January 8, 1985 */
jint = 0;
calerf(*x, &result, jint);
return result;
}
double derfc_ P1C(double *, x)
{
int jint;
double result;
/* This subprogram computes approximate values for erfc(x). */
/* (see comments heading CALERF). */
/* Author/date: W. J. Cody, January 8, 1985 */
jint = 1;
calerf(*x, &result, jint);
return result;
}
double derfcx_ (double *, x)
{
int jint;
double result;
/* This subprogram computes approximate values for exp(x*x) * erfc(x). */
/* (see comments heading CALERF). */
/* Author/date: W. J. Cody, March 30, 1987 */
jint = 2;
calerf(*x, &result, jint);
return result;
}
#endif /* DEBUG */
#define SQRT2 1.414213562373095049
VOID normbase P2C(double *, x, double *, phi)
{
double y;
y = *x / SQRT2;
if (y < 0) {
calerf(-y, phi, 1);
*phi = 0.5 * *phi;
}
else {
calerf(y, phi, 1);
*phi = 1.0 - 0.5 * *phi;
}
}
syntax highlighted by Code2HTML, v. 0.9.1