erf.c

操作系统SunOS 4.1.3版本的源码

#ifndef lint
static	char sccsid[] = "@(#)erf.c 1.1 92/07/30 SMI";
#endif

/*
 * Copyright (c) 1989 by Sun Microsystems, Inc.
 */

/* double erf,erfc(double x)
 * K.C. Ng, March, 1989.
 *			     x
 *		      2      |\
 *     erf(x)  =  ---------  | exp(-t*t)dt
 *	 	   sqrt(pi) \| 
 *			     0
 *
 *     erfc(x) =  1-erf(x)
 *
 * Part of the algorithm is based on cody's erf,erfc program.
 */

extern double exp();
static double
zero	    = 0.0,
tiny	    = 1e-300,
one	    = 1.0,
invsqrtpi_2 = 1.1283791670955125738961589031,
invsqrtpi   = 5.6418958354775628695e-1,
/*
 * Coefficients for approximation to  erf  in first interval
 */
a0 =	3.20937758913846947e03,
a1 =	3.77485237685302021e02,
a2 =	1.13864154151050156e02,
a3 =	3.16112374387056560e00,
a4 =	1.85777706184603153e-1,
b0 =	2.84423683343917062e03,
b1 =	1.28261652607737228e03,
b2 =	2.44024637934444173e02,
b3 = 	2.36012909523441209e01,
/*
 * Coefficients for approximation to  erfc  in second interval
 */
c7 =	5.64188496988670089e-1,
c6 =	8.88314979438837594e0,
c5 =	6.61191906371416295e01,
c4 =	2.98635138197400131e02,
c3 =	8.81952221241769090e02,
c2 =	1.71204761263407058e03,
c1 =	2.05107837782607147e03,
c0 =	1.23033935479799725e03,
c8 =	2.15311535474403846e-8,
d7 =	1.57449261107098347e01,
d6 =	1.17693950891312499e02,
d5 =	5.37181101862009858e02,
d4 =	1.62138957456669019e03,
d3 =	3.29079923573345963e03,
d2 =	4.36261909014324716e03,
d1 =	3.43936767414372164e03,
d0 =	1.23033935480374942e03,
/*
 * Coefficients for approximation to  erfc  in third interval
 */
p4 =	3.05326634961232344e-1,
p3 =	3.60344899949804439e-1,
p2 =	1.25781726111229246e-1,
p1 =	1.60837851487422766e-2,
p0 =	6.58749161529837803e-4,
p5 =	1.63153871373020978e-2,
q4 =	2.56852019228982242e00,
q3 =	1.87295284992346047e00,
q2 =	5.27905102951428412e-1,
q1 =	6.05183413124413191e-2,
q0 =	2.33520497626869185e-3;


#include <math.h>
double erf(x) 
double x;
{
	double erfc(),s,y,t;
	long *px = (long*)&x;
	int ix,iy,i0;

	if(!finite(x)) {
	    if(x!=x) return x+x; 	/* NaN */
	    return copysign(1.0,x);	/* Inf */
	}

	y = fabs(x);
	if(y <= 0.46875){
	    t = one+y;
	    if(one==t) return invsqrtpi_2*x;
	    s = y*y;
	    y = (a0+s*(a1+s*(a2+s*(a3+s*a4))))/(b0+s*(b1+s*(b2+s*(b3+s))));
	    return x*y;
	}

	if(y>10.0) t=tiny; else t = erfc(y);
	return (x>zero)? one-t:t-one;
}

double erfc(x) 
double x;
{
	double erf(),y,s,t,z,r;
	int i;

	if(!finite(x)) {
	    if(x!=x) return x+x; 	/* NaN */
	    return (x>zero)? zero:2.0;	/* +-Inf */
	}
	y = fabs(x);
	if(y <= 0.46875) return one-erf(x);
	if(y <= 4.0) {
	    z  = d0+y*(d1+y*(d2+y*(d3+y*(d4+y*(d5+y*(d6+y*(d7+y)))))));
	    t  = c0+y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*(c6+y*(c7+y*c8)))))));
	    s  = y*y;
	    r  = t/z;
	} else {	/* erfc for |x| > 4.0 */
	    if(x<zero&&y>10.0) return 2.0-tiny;
	    if(y>28.0) return tiny*tiny;	/* underflow */
	    s  = one/(y*y);
	    z  = q0+s*(q1+s*(q2+s*(q3+s*(q4+s))));
	    t  = p0+s*(p1+s*(p2+s*(p3+s*(p4+s*p5))));
	    r  = t/z;
	    r  = (invsqrtpi-s*r)/y;
	}
	i  = (int)(y*16.0); t=(double)i*0.0625;
	r *= exp(-t*t)*exp(-(y-t)*(y+t));
	return (x>zero)? r: 2.0-r;
}