Analytic Error Function and Numeric Inverse Obtained by Geometric Means

1. Introduction

High-speed Monte Carlo simulations are used for a large spectrum of applications from mathematics to economy. As input for such simulations, the probability distribution are usually generated by pseudo-random number sambling, a method going back to a work of John von Neumann from 1951 [1]. In the era of “big data” such methods have to be fast and reliable, a sign of such neccessity being the release of the very first randomness Quside processing unit in 2023 [2]. Still, these samblings need to be cross-checked by exact methods, and for these the knowledge of analytical functions to describe the stochastic processes, among those the error function, are of tremendous importance.

By definition, a function is called analytic if it locally given by a converging Taylor series expansion. Even if a function itself turns out not to be analytic, its inverse can be analytic. The error function can be given analytically, one of these analytic expressions is the integral representation given by Craig in 1991 [3]. Craig mentioned this representation only in passing and did not give a derivation of it. In the following, there have been a couple of derivations of this formula [4,5,6]. In Sec. 2 we add a further one which is based on the same geometric considerations as employed in Ref. [7]. In Sec. 3 we give the series expansion for Craig’s integral representation and show the fast convergence of this series.

For the inverse error function, handbooks for special functions (cf. e.g., Ref. [8]) do not unveil such an analytic property. Instead, this function have to be approximated. Known approximations are dating back to the late 1960s and early 1970s [9,10]) and reach up to semi-analytical approximations by asymptotic expansion (cf., e.g., Refs. [11,12,13,14,15,16]. Using the same geometric considerations, in Sec. 4 we develop a couple of handy approximations which can easily be implemented in different computer languages, indicating the deviations from an exact treatment. In Sec. 5 we test the CPU time and give our conclusions.

2. Derivation of Craig’s integral representation

Ref. [7] provides an approximation for the Gaussian normal distribution obtained by geometric considerations. The same considerations apply to the error function $\mathrm{erf}\left(t\right)$ which is given by the normal distribution $P\left(t\right)$ via

$$\mathrm{erf}\left(t\right)=\frac{1}{\sqrt{\pi }}{\int }_{-t}^t{e}^{-x^2}dx=\frac{1}{\sqrt{2\pi }}{\int }_{-\sqrt{2}t}^{\sqrt{2}t}{e}^{-x^2/2}dx=P\left(\sqrt{2}t\right).$$

(1)

Translating the results of Ref. [7] to the error function, one obtains the approximation of order p to be

$${\mathrm{erf}}_p{\left(t\right)}^2=1-\frac{1}{N}\sum _{n=1}^N{e}^{-k_n^2{t}^2},$$

(2)

where the $N=2^p$ values $k_n$ ($n=1,2,\dots ,N$) are found in the intervals between $1/cos\left(\pi \right(n-1)/(4N\left)\right)$ and $1/cos(\pi n/(4N\left)\right)$ and can be optimised. In Ref. [7] it is shown that

$$|\mathrm{erf}\left(t\right)-\sqrt{1-e^{-k_{0,1}^2{t}^2}}|<0.0033$$

(3)

for $k_{0,1}=1.116$, and with $14\approx 0.0033/0.00024$ times larger precision

$$|\mathrm{erf}\left(t\right)-\sqrt{1-\frac{1}{2}(e^{-k_1^2{t}^2}+e^{-k_2^2{t}^2})}|<0.00024,$$

(4)

where $k_{1,1}=1.01$, $k_{1,2}=1.23345$. For the parameters $k_n$ taking the values $k_n=1/cos(\pi n/\left(4N\right))$ of the upper limits of those intervals, it can be shown that the deviation is given by

$$|\mathrm{erf}\left(t\right)-{\mathrm{erf}}_p\left(t\right)|<\frac{exp(-t^2)}{2N}\sqrt{1-exp(-t^2)}.$$

(5)

Given the values $k_n=1/cos\varphi \left(n\right)$ with $\varphi \left(n\right)=\pi n/\left(4N\right)$, in the limit $N\to \infty$ the sum over n in Equation (2) can be replaced by an integral with measure $dn=(4N/\pi )d\varphi \left(n\right)$ to obtain

$$\mathrm{erf}{\left(t\right)}^2=1-\frac{4}{\pi }{\int }_0^{\pi /4}exp\left(\frac{-t^2}{{cos}^2\varphi }\right)d\varphi .$$

(6)

3. Power series expansion

The integral in Equation (6) can be expanded into a power series in $t^2$,

$$\mathrm{erf}{\left(t\right)}^2=1-\frac{4}{\pi }\sum _{n=0}^{\infty }{c}_n\frac{{(-1)}^n}{n!}{\left(t^2\right)}^n$$

(7)

with

$$\begin{array}{ccc}\hfill c_n& =& {\int }_0^{\pi /4}\frac{d\varphi }{{cos}^{2n}\varphi }={\int }_0^{\pi /4}{(1+{tan}^2\varphi )}^{n}d\varphi ={\int }_0^1{(1+y^2)}^{n-1}dy\hfill \\ & =& \sum _{k=0}^{n-1}\left(\begin{array}{c}n-1\\ k\end{array}\right){\int }_0^1{y}^{2k}dy=\sum _{k=0}^{n-1}\frac{1}{2k+1}\left(\begin{array}{c}n-1\\ k\end{array}\right),\hfill \end{array}$$

(8)

where $y=tan\varphi$. The coefficients $c_n$ can be expressed by the hypergeometric function, $c_n={ }_2{F}_1(1/2,1-n;3/2;-1)$, also known as Barnes’ extended hypergeometric function. On the other hand, we can derive a constraint for the explicit finite series expression for $c_n$ that renders the series in Equation (7) to be convergent for all values of t. In order to be self-contained, intermediate steps to derive this constraint and to show the convergence are shown in the following. Necessary is Pascal’s rule

$$\begin{array}{c}{\displaystyle \left(\begin{array}{c}n\\ k\end{array}\right)+\left(\begin{array}{c}n\\ k-1\end{array}\right)=\frac{n!}{k!(n-k)!}+\frac{n!}{(k-1)!(n-k+1)!}}\\ =\frac{n!(n-k+1+k)}{k!(n-k+1)!}=\frac{(n+1)!}{k!(n+1-k)!}=\left(\begin{array}{c}n+1\\ k\end{array}\right)\end{array}$$

(9)

and the sum over the rows of Pascal’s triangle,

$$\sum _{k=0}^n\left(\begin{array}{c}n\\ k\end{array}\right)=2^n$$

(10)

which can be shown by mathematical induction. The base case $n=0$ is obvious, as $\left(\begin{array}{c}0\\ 0\end{array}\right)=1=2^0$. For the induction step from n to $n+1$ we write the first and last elements $\left(\begin{array}{c}n+1\\ 0\end{array}\right)=1$ and $\left(\begin{array}{c}n+1\\ n+1\end{array}\right)=1$ separately and use Pascal’s rule to obtain

$$\begin{array}{c}\hfill {\displaystyle \sum _{k=0}^{n+1}\left(\begin{array}{c}n+1\\ k\end{array}\right)=1+\sum _{k=1}^n\left(\begin{array}{c}n+1\\ k\end{array}\right)+1=}\\ & =& 1+\sum _{k=1}^n\left(\begin{array}{c}n\\ k\end{array}\right)+\sum _{k=1}^n\left(\begin{array}{c}n\\ k-1\end{array}\right)+1=2\sum _{k=0}^n\left(\begin{array}{c}n\\ k\end{array}\right)=2^{n+1}.\hfill \end{array}$$

(11)

This proves Equation (10). Returning to Equation (8), one has $0\le k\le n-1$ and, therefore,

$$\frac{1}{2n-1}\le \frac{1}{2k+1}\le 1.$$

(12)

For the result in Equation (8) this means that

$$\frac{1}{2n-1}\sum _{k=0}^{n-1}\left(\begin{array}{c}n-1\\ k\end{array}\right)\le c_n\le \sum _{k=0}^{n-1}\left(\begin{array}{c}n-1\\ k\end{array}\right)=2^{n-1},$$

(13)

i.e., the existence of a real number $c_n^{*}$ between $1/(2n-1)$ and 1 such that $c_n=c_n^{*}{2}^{n-1}$. One has

$${\mathrm{erf}}_p{\left(t\right)}^2=1-\frac{4}{\pi }\sum _{n=0}^N{c}_n\frac{{(-1)}^n}{n!}{\left(t^2\right)}^n=1-\frac{2}{\pi }\sum _{n=0}^N{c}_n^{*}\frac{{(-2t^2)}^n}{n!},$$

(14)

and because of $0\le c_n^{*}\le 1$ there is again a real number $c_N^{**}$ in the corresponding open interval so that

$$\frac{2}{\pi }\sum _{n=0}^N{c}_n^{*}\frac{{(-2t^2)}^n}{n!}=c_N^{**}\frac{2}{\pi }\sum _{n=0}^N\frac{{(-2t^2)}^n}{n!}<\frac{2}{\pi }\sum _{n=0}^N\frac{{(-2t^2)}^n}{n!}.$$

(15)

As the latter is the power series expansion of $(2/\pi )e^{-2t^2}$ which is convergent for all values of t, also the original series is convergent and, therefore, ${\mathrm{erf}}_p{\left(t\right)}^2$ with the limiting value shown in Equation (7). A more compact form of the power series expansion is given by

$$\mathrm{erf}{\left(t\right)}^2=\sum _{n=1}^{\infty }{c}_n\frac{{(-1)}^{n-1}}{n!}{\left(t^2\right)}^n,c_n=\sum _{k=0}^{n-1}\frac{1}{2k+1}\left(\begin{array}{c}n-1\\ k\end{array}\right).$$

(16)

4. Approximations for the inverse error function

Based on the geometric approach from Ref. [7], in the following we describe how to find simple, handy formulas that, guided by higher and higher orders of the approximation (2) for the error function lead to more and more advanced approximation of the inverse error function. Starting point is the degree $p=0$, i.e., the approximation in Equation (3). Inverting $E={\mathrm{erf}}_0\left(t\right)={(1-e^{-k_{0,1}^2{t}^2})}^{1/2}$ leads to $t^2=-ln(1-E^2)/k_{0,1}^2$, and using the parameter $k_{0,1}=1.116$ from Equation (3) gives

$$T_0=\sqrt{-ln(1-E^2)}/k_{0,1}^2.$$

(17)

For $0\le E\le 0.92$ the relative deviation $(T_{\left(0\right)}-t)/t$ from the exact value t is less than $1.11\%$, for $0\le E<1$ the deviation is less than $10\%$. Therefore, for $E>0.92$ a more precise formula has to be used. As such higher values for E appear only in $8\%$ of the cases, this will not essentially influence the CPU time.

Continuing with $p=1$, we insert $T_0=\sqrt{-ln(1-E^2)}/k_{0,1}^2$ into Equation (2) to obtain

$${\mathrm{erf}}_1\left(T_0\right)=\sqrt{1-\frac{1}{2}(e^{-k_{1,1}^2{T}_0^2}+e^{-k_{1,2}^2{T}_0^2})},$$

(18)

where $k_{1,1}=1.01$ and $k_{1,2}=1.23345$ are the same as for Equation (4). Taking the derivative of Equation (1) and approximating this by the difference quotient, one obtains

$$\frac{\mathrm{erf}\left(t\right)-\mathrm{erf}\left(T_0\right)}{t-T_0}=\frac{\Delta \mathrm{erf}\left(t\right)}{\Delta t}{|}_{t=T_0}\approx \frac{d\mathrm{erf}\left(t\right)}{dt}{|}_{t=T_0}=\frac{2}{\sqrt{\pi }}{e}^{-T_0^2},$$

(19)

leading to $t\approx T_1=T_0+\frac{1}{2}\sqrt{\pi }{e}^{T_0^2}(E-{\mathrm{erf}}_1\left(T_0\right))$. In this case, for in the larger interval $0\le E\le 0.995$ the relative deviation $(T_1-t)/t$ is less than $0.1\%$. Using ${\mathrm{erf}}_2\left(t\right)$ instead of ${\mathrm{erf}}_1\left(t\right)$ and inserting $T_1$ instead of $T_0$ one obtains $T_2$ with a relative deviation of maximally $0.01\%$ for the same interval. The results are shown in Figure 1.

The method to can be optimised by a method similar to the shooting method in boundary problems, giving dynamics to the calculation. Suppose that following one of the previous methods, for a particular argument E we have found an approximation $t_0$ for the value of the inverse error function at this argument. Using $t_1=1.01t_0$, one can adjust the improved result

$$t=t_0+A(E-\mathrm{erf}\left(t_0\right))$$

(20)

by inserting $E=\mathrm{erf}\left(t\right)$ and and calculating A for $t=t_1$. In general, this procedure gives a vanishing deviation close to $E=0$. In this case and for $t_0=T_1$, in the interval $0\le E\le 0.7$ the maximal deviation is slightly larger than ${10}^{-6}=0.0001\%$ while up to $E=0.92$ the deviation is restricted to ${10}^{-5}=0.001\%$. A more general ansatz

$$t=t_0+A(E-\mathrm{erf}\left(t_0\right))+B{(E-\mathrm{erf}\left(t_0\right))}^2$$

(21)

can be adjusted by inserting $E=\mathrm{erf}\left(t\right)$ for $t=1.01t_0$ and $t=1.02t_0$, and the system of equations

$$\Delta t=A\Delta E_1+B\Delta E_1^2,2\Delta t=A\Delta E_2+B\Delta E_2^2$$

(22)

with $\Delta t=0.01t_0$, $\Delta E_i=\mathrm{erf}\left(t_i\right)-\mathrm{erf}\left(t_0\right)$ can be solved for A and B to obtain

$$A=-\frac{(2\Delta E_1^2-\Delta E_2^2)\Delta t}{\Delta E_1\Delta E_2(\Delta E_1-\Delta E_2)},B=\frac{(-2\Delta E_1+\Delta E_2)\Delta t}{\Delta E_1\Delta E_2(\Delta E_1-\Delta E_2)}.$$

(23)

For $0\le E\le 0.70$ one obtains a relative deviation of $1.5·{10}^{-8}$, for $0\le E\le 0.92$ the maximal deviation is $5·{10}^{-7}$. Finally, the adjustment of

$$t=t_0+A(E-\mathrm{erf}\left(t_0\right))+B{(E-\mathrm{erf}\left(t_0\right))}^2+C{(E-\mathrm{erf}\left(t_0\right))}^3$$

(24)

leads to

$$\begin{array}{ccc}\hfill A& =& (3\Delta E_1^2\Delta E_2^2(\Delta E_1-\Delta E_2)-2\Delta E_1^2\Delta E_3^2(\Delta E_1-\Delta E_3)\phantom{(}\hfill \\ & & \phantom{(}+\Delta E_2^2\Delta E_3^2(\Delta E_2-\Delta E_3))\Delta t/D,\hfill \\ \hfill B& =& (-3\Delta E_1\Delta E_2(\Delta E_1^2-\Delta E_2^2)+2\Delta E_1\Delta E_3(\Delta E_1^2-\Delta E_3^2)\phantom{(}\hfill \\ & & \phantom{(}-\Delta E_2\Delta E_3(\Delta E_2^2-\Delta E_3^2))\Delta t/D,\hfill \\ \hfill C& =& (3\Delta E_1\Delta E_2(\Delta E_1-\Delta E_2)-2\Delta E_1\Delta E_3(\Delta E_1-\Delta E_3)\phantom{(}\hfill \\ & & \phantom{(}+\Delta E_2\Delta E_3(\Delta E_2-\Delta E_3))\Delta t/D,\hfill \end{array}$$

(25)

where $D=\Delta E_1\Delta E_2\Delta E_3(\Delta E_1-\Delta E_2)(\Delta E_1-\Delta E_3)(\Delta E_2-\Delta E_3)$. For $0\le E\le 0.70$ the relative deviation is restricted to $5·{10}^{-10}$ while up to $E=0.92$ the maximal relative deviation is $4·{10}^{-8}$. The results for the deviations of $T_{\left(n\right)}$ ($n=1,2,3$) for linear, quadratic and cubic dynamical approximation are shown in Figure 2.

5. Conclusions

In order to test the feasibility and speed, we have coded our algorithm in the computer language C under Slackware 15.0 on an ordinary hp laptop. The dependence of the CPU time for the calculation is estimated by calculating the value ${10}^6$ times in sequence. The speed of the calculation of course turns does not depend on the value for E, as the precision is not optimised. This would have to be neccessary for a practical application. For an arbitrary starting value $E=0.8$ we perform this test, and the results are given in Table 1. An analysis of the table shows that a further step in the degree p doubles the run time while the dynamics for increasing n adds a constant value of approximately $0.06$ seconds to the result. Despite the fact the increase of the dynamics needs the solution of a linear system of equations and the coding of the result, this endeavour is justified, as by using the dynamics one can increase the precision of the result without loosing calculational speed.