- Language
- C++
- Purpose
- Calculate the error function, <math>\mathbf{erf}(x)</math>, and complementary error function, <math>\mathbf{erfc}(x)</math>
- Author
- James R. Craig, as adapted from the VB code of R.J. Charbeneau
#include <math.h>
/***********************************************************************
erf & erfc
************************************************************************
Error Function (erf), Complementary Error Function (erfc)
Author: James R. Craig
Adapted from R. Charbeneau, Groundwater Hydraulics and Pollutant Transport, 2000
-----------------------------------------------------------------------*/
double erfc(const double &x)
{
double tmp(fabs(x)); //take abs so that we are always in positive quadrant.
static double fun;
static double f1;
static double tmp2;
static double tmp3;
if(tmp > 3.0)
{
f1 = (1.0 - 1.0/(2.0 * tmp * tmp)
+ 3.0/(4.0 * pow(tmp,4))
- 5.0/(6.0 * pow(tmp,6)));
fun = f1 * exp(-tmp * tmp) / (tmp * sqrt(PI));
}
else
{
tmp2 = 1.0 / (1.0 + (0.3275911 * tmp));
tmp3 = 0.254829592 * tmp2 //5th order polynomial interpolation
- (0.284496736 * tmp2 * tmp2)
+ (1.421413741 * pow(tmp2,3))
- (1.453152027 * pow(tmp2,4))
+ (1.061405429 * pow(tmp2,5));
fun = tmp3 * exp(-tmp * tmp);
}
if (tmp == x) {return fun;}
else {return (2.0-fun);}
}
//------------------------------------------------------------------------------
double erf(const double &x)
{
return 1.0-erfc(x);
}
Notes
- This algorithm seems to implement eq. 7.1.26 from Abramowitz and Stegun's handbook of mathematical functions
- The polynomial approximation used yields an error bound of <math>1.5\times 10^{-7}</math>