Error function and complementary error function

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>