The docs for the Boost Math Toolkit have a long list of links , including Abramowitz and Stegun. The erf-function interface contains a policy parameter that can be used to control numerical accuracy (and therefore its complexity at run time).
#include <boost/math/special_functions/erf.hpp> namespace boost{ namespace math{ template <class T> calculated-result-type erf(T z); template <class T, class Policy> calculated-result-type erf(T z, const Policy&); template <class T> calculated-result-type erfc(T z); template <class T, class Policy> calculated-result-type erfc(T z, const Policy&); }} // namespaces
UPDATE
Below is a shorthand copy of the "Implementation" section of the previously provided link to the erf function:
Implementation
All versions of these functions first use the usual reflection formulas to make their arguments positive:
erf(-z) = 1 - erf(z); erfc(-z) = 2 - erfc(z); // preferred when -z < -0.5 erfc(-z) = 1 + erf(z); // preferred when -0.5 <= -z < 0
General versions of these functions are implemented in terms of an incomplete gamma function.
When the size of the trait (mantissa) is recognized (currently, for 53, 64 and 113-bit reals, plus 24-bit accuracy with an accuracy of 24 bits processed by promotion, it doubles), then a series of rational approximations developed by JM is used.
For z <= 0.5, then a rational approximation to erf is used, based on the observation that erf is an odd function, and therefore erf is calculated using:
erf(z) = z * (C + R(z*z));
where the rational approximation R (z * z) is optimized for the absolute error: if its absolute error is sufficiently small compared to the constant C, then any rounding error that occurs when calculating R (z * z) will effectively disappear from the result. As a result, the error for erf and erfc in this area is very low: the last bit is incorrect only in a very small number of cases.
For z> 0.5, we note that on a small interval [a, b):
erfc(z) * exp(z*z) * z ~ c
for some constant c.
Therefore, for z> 0.5, we calculate erfc using:
erfc(z) = exp(-z*z) * (C + R(z - B)) / z;
Again, R (z - B) is optimized for absolute error, and the constant C is the average value erfc (z) * exp (z * z) * z taken at the ends of the range. Once again, while the absolute error in R (z - B) is small compared to c, c + R (z - B) will be correctly rounded, and the error of the result will depend only on the accuracy of the exp function. In practice, in all cases, except for a very small number of cases, the error is limited to the last bit of the result. The constant B is chosen so that the left end of the range of rational approximation is 0.
For large z in the range [a, + ∞], the above approximation changes to:
erfc(z) = exp(-z*z) * (C + R(1 / z)) / z;
Rational approximations are explained in a painful detail . If you need more information, you can always look at the source code .