Algorithm of boost::math::erf

Question

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.

• Which methods are mixed?
• How are the methods mixed?
• What is the complexity of the erf-function (constant time)?

Sebastian

Answer 1

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.