# How to accurately calculate the error function erf(x) with a computer?

I am looking for an accurate algorithm to calculate the error function

I have tried using [this formula]

(http://stackoverflow.com/a/457805) (Handbook of Mathematical Functions, formula 7.1.26), but the results are not accurate enough for the application.

-
You may want to take a look at python's code.google.com/p/mpmath or other libraries that advertise a "multiple precision" feature. Also, this may be a better question for stack overflow instead, since it's more of a computer science thing. – Jon Bringhurst Jul 20 '10 at 20:26
@Jon: Nope, I'm not interested in a library, there is no such library for the language I'm writing in (yet). I need the mathematical algorithm. – badp Jul 20 '10 at 20:49
Have you tried numerical integration? Gaussian Quadrature is an accurate technique – Digital Gal Aug 28 '10 at 1:25
GQ is nice, but with (a number of) efficient methods for computing $\mathrm{erf}$ already known, I don't see the point. – J. M. Aug 29 '10 at 23:07

I am assuming that you need the error function only for real values. For complex arguments there are other approaches, more complicated than what I will be suggesting.

If you're going the Taylor series route, the best series to use is formula 7.1.6 in Abramowitz and Stegun. It is not as prone to subtractive cancellation as the series derived from integrating the power series for $\exp(-x^2)$. This is good only for "small" arguments. For large arguments, you can use either the asymptotic series or the continued fraction representations.

Otherwise, may I direct you to these papers by S. Winitzki that give nice approximations to the error function.

I wrote about the computation of the (complementary) error function (couched in different notation) in this answer to a CV question.

-
Assumption correct. :) – badp Jul 30 '10 at 20:02
+1 for the Winitzki reference; I've seen that 2nd paper before + it's a nice one. – Jason S Aug 6 '10 at 13:05
The only thing infuriating for me about the Winitzki paper is that apart from the Hermite-Pade ansatzes (or however you pluralize these things :P), he gives no indication on how he arrived at these simple approximants. – J. M. Aug 6 '10 at 13:10
@Ben: Thanks a bunch! – J. M. May 4 '11 at 5:02

You can use a Taylor polynomial of sufficient degree to guarantee the accuracy that you need. (The Taylor series for erf(x) is given on the Wikipedia page to which you linked.) The Lagrange Remainder term can be used to bound the error in the Taylor series approximation.

-

Here's a link to the boost c++ math library documentation. They use their implementation of the incomplete gamma function, which in turn uses a mixed approach depending on the argument. It's all fairly well documented should you care to duplicate their method. And it looks like their error is within a few multiples of the machine epsilon.

Other than that, I would try the Taylor series. Numerical approximation might lead to a larger error term than the analytic one though, and it will only be valid in a neighborhood of 0. For larger values you could use the asymptotic series.

Another idea would be to restrict the domain to a closed interval. If you size it properly, then the function will appear constant with respect to your machine precision outside of this interval. Once you have a compact domain, you can know exactly how many Taylor terms you need, or you can use other types of spline interpolation. Chebyshev polynomials come to mind.

As for the problem that the language your writing in has no such library already: for me that is probably not as big of a deal as you think. Most languages seem to have a way to link in C functions, and if that is the case, then there is an open source implementation somewhere out there.

-
The naïve (alternating) Maclaurin series is not really that numerically sound; I had already mentioned in my answer the modified series that has much better properties for computing $\mathrm{erf}(x)$ near the origin. The (Laplace) continued fraction tends to be slightly easier to handle than the asymptotic series for medium-to-large arguments. – J. M. Sep 1 '11 at 10:34
If you're going for approximations of fixed degree near the origin, constructing a Padé approximant is slightly better than using a truncated Maclaurin series. – J. M. Sep 1 '11 at 10:35
I'll agree with that assessment. I thought about mentioning the numerical instability, but the post was already long. I think the best bet is to use a hybrid approach depending on the size of the argument. That way you can make an appropriate trade off of precision versus speed. Which method you use for which intervals is down to experimentation. – Tim Seguine Sep 1 '11 at 10:51
A lot of this comes down to the desired accuracy and how fast it must be... I think Chebyshev interpolation is worth looking into in any case – Tim Seguine Sep 1 '11 at 10:56

A simple way of computing error function is to use Kummer's equation of the form of,

$M(a,b,z)=\sum_{s=0}^{\infty} \frac{(a)_s}{(b)_s s!}z^s=1+\frac{a}{b}z+\frac{a(a+1)}{b(b+1)2!}z^2+...,\quad z \in \mathbb{C}$

and

$M(a,b,z)=\sum_{s=0}^{\infty}\frac{(a)_s}{\Gamma{(b+s)}s!}z^s$

and

$erf(z)=\frac{2z}{\sqrt{\pi}}M(.5,1.5,-z^2)=\frac{2z}{\sqrt{\pi}}e^{-z^2}M(1,1.5,z^2)$

Here is the R code,

f<-function(a,b,z,maxt=5){
s=1:maxt
a=a+s
b=b+s
ss=1
for(i in s){
mt=prod(a[1:i]/b[1:i])
ss=ss+mt *z^(2*i)/factorial(i)
}
ss=2*z/sqrt(pi)*exp(-z^2)*ss
return(ss)
}


here is the plot (for real z)

-