# efficient and accurate approximation of error function

I am looking for the numerical approximation of error function, which must be efficient and accurate. Thanks in advance $$\mathrm{erf}(z)=\frac2{\sqrt\pi}\int_0^z e^{-t^2} \,\mathrm dt$$

"Efficient and accurate" is probably contradictory... Have you tried the one listed in http://en.wikipedia.org/wiki/Error_function#Approximation_with_elementary_functions ?

• yes, I have tried this. Its accuracy is up to 2 decimal places. Do we have more than this? – shaikh Jun 3 '11 at 2:40
• @shaikh, C99 has an erf function, which should be quite accurate. – lhf Jun 3 '11 at 2:42
• where can I find the derivation of C99 erf...? – shaikh Jun 3 '11 at 2:45
• @shaikh, oh, you have to read the code. See for instance the cephes library. – lhf Jun 3 '11 at 2:48
• @shaikh: Or boost's implementation – ziyuang Jun 3 '11 at 2:51
• I will try this. thx – shaikh Jun 3 '11 at 2:42

The oldest approximations to the error function $\mathrm{erf}$ encountered in existing software come from:

Cecil Hastings, Jr., Approximations for Digital Computers, Princeton, NJ: Princeton University Press 1955, pp. 167-169, pp. 185-187

Hastings gives six minimax approximations designed for absolute error bounds, the most accurate of which bounds the error by $1.5 \times {10}^{-7}$. Four of the approximations were later incorporated into:

Milton Abramowitz & Irene A. Stegun (eds.), Handbook of Mathematical Functions, Washington, D.C.: National Bureau of Standards 1964, sections 7.1.25 - 7.1.28 (online)

The first approximations designed with a relative error bound were published by William J. Cody:

William J. Cody, "Rational Chebyshev approximations for the error function." Mathematics of Computation, Vol. 23, No. 107, July 1969, pp. 631-637 (online)

Cody used these rational minimax approximations as the basis for a Fortran implementation of the error function $\mathrm{erf}$, complementary error function $\mathrm{erfc}$ and the scaled complementary error function $\mathrm{erfcx}$ which he contributed to SPECFUN (online) during the 1980s, and improved until 1990. This work is described in:

William J. Cody, "Performance Evaluation of Programs for the Error and Complementary Error Functions." ACM Transactions on Mathematical Software, Vol. 16, No. 1, March 1990, pp. 29-37

In testing with 500 million random test vectors I found no errors larger than 6.03 ulps in Cody's SPECFUN implementation of $\mathrm{erf}$.

Around the same time Stephen Moshier released the Cephes library (online), which he maintained well into the 2000s. This library is written in C and contains an implementation of $\mathrm{erf}$ that is based on rational approximations as well. The library was documented in:

Stephen Lloyd Baluk Moshier, Methods and Programs for Mathematical Functions, Ellis Horwood Limited, Chichester, 1989

Testing with 500 million random test vectors I found no errors in excess of 2.95 ulps in Moshier's Cephes implementation of $\mathrm{erf}$.

In 1993 Sun Microsystems, Inc released the Freely Distributable Math Library, or FDLIBM, which is based on the work of K. C. Ng. This library, written in C, aimed to provide a faithfully-rounded implementation of $\mathrm{erf}$ (online), based on rational minimax approximations. A faithfully-rounded function has an error bound of 1 ulp or less. The result returned is either the floating-point number closest to the mathematical result, or its immediate neighbor.

Based on my tests with 500 million random test vectors, the FDLIBM implementation falls slightly short of its accuracy goal, in that errors in excess of 1 ulp do occur but no errors larger than 1.01 ulp were observed.

Using modern tools and processor architectures, it is not too hard to create one's own competitive implementation of $\mathrm{erf}$. Rational approximations have largely fallen out of favor because of the high computational cost of division compared with additions and multiplications. Polynomial minimax approximations are desired, and can readily be generated with tools such as Maplesoft's Maple and Wolfram's Mathematica, both of which have dedicated functions for this: the numapprox command in Maple and the MiniMaxApproximation command in Mathematica.

It should be noted that these generate minimax approximations without taking into account that the coefficients will be rounded to target precision, e.g. IEEE-754 double precision. The Sollya tool with its fpminimax command is unique in offering that functionality, it typically generates superior approximations.

To my knowledge, all these tools use variants of the Remez algorithm which may be augmented by additional heuristics. It is quite possible to write one's own implementation of this algorithm, this basically requires a solver for systems of linear equations and an arbitrary-precision library. I use Richard Brent's MP library for the latter, but would not recommended it for new software since the design is outdated.

Looking at the graph of $\mathrm{erf}$ we find that it is symmetric about the origin, so approximations can be restricted to the positive half-plane. The graph further suggest two basic approximation intervals. For smaller arguments, the function closely follows its tangent through the origin, while for larger arguments it approaches $y = 1$ in the shape of exponential decay. The crossover between the two regions is roughly around $x = 1$. For simplicity, we might choose that as the actual crossover point.

The series expansion for $\mathrm{erf}$ suggests a polynomial approximation of the form $x p({x}^{2})$ for $x \in [0, 1]$, while a quick numerical check shows that for IEEE-754 double precision, erf(x) == 1.0 holds for $x > 5.921875$, so we can use an approximation $1 - \exp (x q(x))$ for $x \in (1, 5.921875]$; $p$ and $q$ are polynomials.

The fused-multiply add operation (FMA) enables the evaluation of the polynomials with improved accuracy, by reducing rounding error and providing protection against subtractive cancellation, as the intermediate product is not rounded. An FMA operation is provided by most modern processor architectures, and is exposed in ISO C/C++ via the standard math function fma().

Below is exemplary C code for a IEEE-754 double-precision implementation that follows the design guidelines sketched above. The accuracy of this implementation depends somewhat on the accuracy of the exponential function exp(). With a faithfully-rounded implementation of exp() I observed no errors in excess of 1.07 ulp using 1 billion random test vectors. With some tweaking, a faithfully-rounded approximation should be possible, as I demonstrated for equivalent single-precision code.

double my_erf (double a)
{
double r, s, t, u;

t = fabs (a);
s = a * a;
if (t >= 1.0) {
// -0x1.9f363bac00b3ap-58, 0x1.17f6b1da05762p-51
r = fma (-5.6271698458222802e-018, t, 4.8565951833159269e-016);
// -0x1.66b85b7ea724ep-46, 0x1.22907eeb53bbdp-41
u = fma (-1.9912968279795284e-014, t, 5.1614612430130285e-013);
r = fma (r, s, u);
r = fma (r, t, -9.4934693735334407e-012);  // -0x1.4e0591fcf7297p-37
r = fma (r, t,  1.3183034417266867e-010);  //  0x1.21e5e2d8345c4p-33
r = fma (r, t, -1.4354030030124722e-009);  // -0x1.8a8f81b7cd224p-30
r = fma (r, t,  1.2558925114367386e-008);  //  0x1.af85793b8cf2fp-27
r = fma (r, t, -8.9719702096026844e-008);  // -0x1.8157db0ed6deap-24
r = fma (r, t,  5.2832013824236141e-007);  //  0x1.1ba3c45370f62p-21
r = fma (r, t, -2.5730580226095829e-006);  // -0x1.595999b7e9e12p-19
r = fma (r, t,  1.0322052949682532e-005);  //  0x1.5a59c27b3c70ep-17
r = fma (r, t, -3.3555264836704290e-005);  // -0x1.197b61ee3732bp-15
r = fma (r, t,  8.4667486930270974e-005);  //  0x1.631f0597f6424p-14
r = fma (r, t, -1.4570926486272249e-004);  // -0x1.319310dfb85f3p-13
r = fma (r, t,  7.1877160107951816e-005);  //  0x1.2d798353da7c3p-14
r = fma (r, t,  4.9486959714660115e-004);  //  0x1.037445e25d35dp-11
r = fma (r, t, -1.6221099717135142e-003);  // -0x1.a939f51db8bcbp-10
r = fma (r, t,  1.6425707149019371e-004);  //  0x1.5878d80188692p-13
r = fma (r, t,  1.9148914196620626e-002);  //  0x1.39bc5e0e9e090p-6
r = fma (r, t, -1.0277918343487556e-001);  // -0x1.a4fbc8f8ff7d7p-4
r = fma (r, t, -6.3661844223699315e-001);  // -0x1.45f2da3ae06f8p-1
r = fma (r, t, -1.2837929411398119e-001);  // -0x1.06ebb92d9ffa8p-3
r = fma (r, t, -t);
r = 1.0 - exp (r);
r = copysign (r, a);
} else {
r =            -7.7794684889579207e-10; // -0x1.abae491c43c5cp-31
r = fma (r, s,  1.3710980398029044e-8); //  0x1.d71b0f1b10b88p-27
r = fma (r, s, -1.6206313758493333e-7); // -0x1.5c0726f04dd28p-23
r = fma (r, s,  1.6447131571278885e-6); //  0x1.b97fd3d992901p-20
r = fma (r, s, -1.4924712302009674e-5); // -0x1.f4ca4d6f3e2a0p-17
r = fma (r, s,  1.2055293576900647e-4); //  0x1.f9a2baa8fede1p-14
r = fma (r, s, -8.5483259293145256e-4); // -0x1.c02db03dd71f5p-11
r = fma (r, s,  5.2239776061184708e-3); //  0x1.565bccf92b2f6p-8
r = fma (r, s, -2.6866170643111489e-2); // -0x1.b82ce311fa944p-6
r = fma (r, s,  1.1283791670944184e-1); //  0x1.ce2f21a040d15p-4
r = fma (r, s, -3.7612638903183515e-1); // -0x1.812746b0379bcp-2
r = fma (r, s,  1.2837916709551256e-1); //  0x1.06eba8214db68p-3
r = fma (r, a, a);
}
return r;
}


I know its an old post, but others stumbling upon this post might find it helpful. You can try the accurate approximate analytical expression for faster numerical evaluation such as this or answers in this post especially by Ron Gordon.

It can be computed by using the complex error function (aka the Faddeeva function):

$$1-{\rm{erf}}(z)=e^{-z^2}w(iz)$$

Matlab and C packages for the Faddeeva function are available in the Matlab Central.