are there any details available for the algorithm behind the erf-function of boost? The documentation of the module is not very precise. All I found out is that several methods are mixed. For me it looks like variations of Abramowitz and Stegun.
Sebastian
The docs for Boost Math Toolkit has a long list of references, among which Abramowitz and Stegun. The erf-function interface contains a policy template parameter that can be used to control the numerical precision (and hence its run-time complexity).
#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 a verbatim copy of the section "Implementation" of the earlier provided reference 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
The generic versions of these functions are implemented in terms of the incomplete gamma function.
When the significand (mantissa) size is recognised (currently for 53, 64 and 113-bit reals, plus single-precision 24-bit handled via promotion to double) then a series of rational approximations devised by JM are 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 optimised for absolute error: as long as its absolute error is small enough compared to the constant C, then any round-off error incurred during the computation of R(z*z) will effectively disappear from the result. As a result the error for erf and erfc in this region is very low: the last bit is incorrect in only a very small number of cases.
For z > 0.5 we observe that over a small interval [a, b) then:
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 optimised for absolute error, and the constant C is the average of erfc(z) * exp(z*z) * z taken at the endpoints of the range. Once again, as long as the absolute error in R(z - B) is small compared to c then c + R(z - B) will be correctly rounded, and the error in the result will depend only on the accuracy of the exp function. In practice, in all but a very small number of cases, the error is confined to the last bit of the result. The constant B is chosen so that the left hand end of the range of the rational approximation is 0.
For large z over a range [a, +∞] the above approximation is modified to:
erfc(z) = exp(-z*z) * (C + R(1 / z)) / z;
The rational approximations are explained in excruciating detail. Tf you need more details, you can always look at the source code.