How to accurately calculate the error function $\operatorname{erf}(x)$ with a computer?

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

$$\operatorname{erf}(x)=\frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2}\ dt$$

I have tried using the following formula,

(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. Jul 20, 2010 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.
Jul 20, 2010 at 20:49
• Have you tried numerical integration? Gaussian Quadrature is an accurate technique Aug 28, 2010 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. Aug 29, 2010 at 23:07
• @Evgeny this was probably python but I have honestly no idea what I was up to seven years ago
Nov 9, 2017 at 8:54

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. :)
Jul 30, 2010 at 20:02
• +1 for the Winitzki reference; I've seen that 2nd paper before + it's a nice one. Aug 6, 2010 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. Aug 6, 2010 at 13:10
• Web Archive of Winitzki Paper: web.archive.org/web/20100601094000/https://…, it´s Ansätze
– user491315
Jan 25, 2019 at 13:45

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. Sep 1, 2011 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. Sep 1, 2011 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. Sep 1, 2011 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 Sep 1, 2011 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

$$\text{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)

I'm surprised that nobody proposed the obvious option: Numerical quadrature. One of the simplest methods would be the composed trapezoidal rule, where you take equally spaced points $$0=x_0, x_1, \cdots, x_n = x$$, with spacing $$h$$, and obtain the estimate $$\mathrm{erf} (x) \approx \frac{h}{\sqrt{\pi}}\left(1+e^{-x^2} + 2 \sum_{i=1}^{n-1} e^{-x_i^2}\right)$$

You can use the error formula to decide which $$n$$ to use for a particular target precision.

• Could you also explore the computational complexity to reach the same accuracy as the other methods mentioned, starting with the error $10^{-7}$ method given in the question? May 2, 2019 at 8:48
• For this specific method, the error function would be given by $$|\varepsilon(x)| = \left| \frac{x^3}{12n^2} f''(\xi) \right| \leq \frac{x^3}{6n^2}$$ So, for a given $x$, you can choose $n$ so that any prescribed accuracy is reached. Performing Richardson extrapolation could accelerate the process (Romberg's method). I'm not sure this was your question though.@Lutz May 2, 2019 at 10:28
• So to get an error of $10^{-7}$ at $x=10$ you need $n\simeq 4\cdot 10^4$ for the trapzoidal rule or possibly $n\sim 10^2$ for the Simpson rule. Going down further will probably increase the constants too much due to insufficient sampling density. This is still a lot more effort than the dedicated approximation formulas, more so when $x$ is large. Any true improvement, variable step size, re-parametrization of the integral to get a more linear integrand etc. will lead to methods mentioned in the other answers. So your "obvious option" is not so obvious after all. May 2, 2019 at 10:49
• @LutzL I see your point and kind of agree, but I still think it is obvious, in the sense that is a very natural first approach that will work for a wide range of problems. Considering that for an error of $10^{-7}$ you don't need to compute for any $x>4$ (you can set the value to $1$) and using Romberg's method I expect this solution to be competitive, while having the advantage of requiring no extra work should you decide to make any changes in the integrand. May 2, 2019 at 11:17
• Can you do the same thing for the inverse of the error function? Oct 24, 2019 at 23:12

You could get good theoretical approximations of the error function using $$\text{erf}\left(x\right)\sim \text{sgn}(x)\sqrt{1-\exp\Big(-\frac 4 {\pi}\,\frac{1+P_n(x^2)}{1+Q_n(x^2)}\,x^2 \Big)}$$ For example $$P_1(x^2)=\frac{10-\pi ^2}{5 (\pi -3) \pi }x^2\qquad Q_1(x^2)=\frac{120-60 \pi +7 \pi ^2}{15 (\pi -3) \pi }$$ while $$P_2(x^2)=\frac{\left(105840-110880 \pi +37800 \pi ^2-4260 \pi ^3+69 \pi ^4-17 \pi ^5\right) }{3 \pi \left(-12600+12600 \pi -3360 \pi ^2-30 \pi ^3+73 \pi ^4\right)}x^2+$$ $$\frac{-2116800+1270080 \pi ^2-504000 \pi ^3+48510 \pi ^4+503 \pi ^6}{315 \pi ^2 \left(-12600+12600 \pi -3360 \pi ^2-30 \pi ^3+73 \pi ^4\right)}x^4$$

$$Q_2(x^2)=\frac{60480-70560 \pi +27720 \pi ^2-3600 \pi ^3-143 \pi ^4+43 \pi ^5}{\pi \left(-12600+12600 \pi -3360 \pi ^2-30 \pi ^3+73 \pi ^4\right)}x^2+$$ $$\frac{-6350400+8467200 \pi -4515840 \pi ^2+1192800 \pi ^3-145320 \pi ^4+2380 \pi ^5+793 \pi ^6}{105 \pi ^2 \left(-12600+12600 \pi -3360 \pi ^2-30 \pi ^3+73 \pi ^4\right)}x^4$$ and so on.

The first one gives a maximum absolute error of $$0.000150$$ while the second gives a maximum absolute error of $$0.000012$$.

With regard to the infinite norms, they are respectively $$3.04\times 10^{-8}$$ and $$1.20\times 10^{-10}$$.