Analytic error function and numeric inverse obtained by geometric means

Using geometric considerations, we provide a clear derivation of the integral representation for the error function, known as the Craig formula. We calculate the corresponding power series expansion and prove the convergence. The same geometric means finally help to systematically derive handy formulas that approximate the inverse error function. Our approach can be used for applications in e.g.\ high-speed Monte Carlo simulations where this function is used extensively.


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 pseudorandom 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 eestidima@gmail.com, stefan.groote@ut.ee. 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 semianalytical 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 discuss our result and test the CPU time. Sec. 6 contains our conclusions.

Derivation of Craig's integral representation
Ref. [7] provides an approximation for the integral over the Gaussian standard normal distribution obtained by geometric considerations that is related to the cumulative distribution function via is the Laplace function. The same considerations apply to the error function erf(t) which is related to P (t) via Translating the results of Ref. [7] to the error function, one obtains the approximation of order p to be where the N = 2 p values k p,n (n = 1, 2, . . . , N) are found in the intervals between 1/ cos(π(n − 1)/(4N)) and 1/ cos(πn/(4N)). The way of selecting those values is extensively described in Ref. [7] where it is shown that for k 0,1 = 1.116, and with 14 ≈ 0.0033/0.00024 times larger precision for k 1,1 = 1.01, k 1,2 = 1.23345. For the parameters taking the values k p,n = 1/ cos(πn/(4N)) of the upper limits of those intervals, it can be shown that the deviation is given by Given the values k p,n = 1/ cos φ(n) with φ(n) = πn/(4N), in the limit N → ∞ the sum over n in Eq.

Power series expansion
The integral in Eq. (6) can be expanded into a power series in t 2 , where y = tan φ. 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 Eq. (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 and the sum over the rows of Pascal's triangle, n + 1 k This proves Eq. (10). Returning to Eq. (8), one has 0 ≤ k ≤ n − 1 and, therefore, For the result in Eq. (8) this means that 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 and because of 0 ≤ c * n ≤ 1 there is again a real number c * * N in the corresponding open interval so that (15) 2 As the latter is the power series expansion of (2/π)e −2t 2 which is convergent for all values of t, also the original series is convergent and, therefore, erf p (t) 2 with the limiting value shown in Eq. (7). A more compact form of the power series expansion is given by

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 Eq. (3).
(2) to obtain where k 1,1 = 1.01 and k 1,2 = 1.23345 are the same as for Eq. (4). Taking the derivative of Eq. (1) and approximating this by the difference quotient, one obtains ). In this case, for in the larger interval 0 ≤ E ≤ 0.995 the relative deviation (T 1 − t)/t is less than 0.1%. Using erf 2 (t) instead of erf 1 (t) 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 Fig. 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 by inserting E = erf(t) 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 ≤ E ≤ 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 can be adjusted by inserting E = erf(t) for t = 1.01t 0 and t = 1.02t 0 , and the system of equations can be solved for A and B to obtain For 0 ≤ E ≤ 0.70 one obtains a relative deviation of 1.5 · 10 −8 , for 0 ≤ E ≤ 0.92 the maximal deviation is 5·10 −7 . Finally, the adjustment of . For 0 ≤ E ≤ 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 (n) (n = 1, 2, 3) for linear, quadratic and cubic dynamical approximation are shown in Fig. 2.

Discussion
In order to test the feasibility and speed, we have coded our algorithm in the computer language C under Slackware 15.0 (linux 5.15.19) on an ordinary hp laptop with Intel® Core™2 Duo CPU P8600 @ 2.4GHz with 3MiB memory used. 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 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 this 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 that 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.  (1) The results for the deviations in Figs. 1 and 2 are multiplied by increasing decimal powers in order to make the result comparable. This fact indicates that the convergence is improved in each of the steps for p or n at least by the corresponding inverse power. While the approximations in the statical approximations n = 0 in Fig. 1 show both deviations close to E = 0 and for higher values of E, the dynamical approximations in Fig. 2 show no deviation at E = 0 and moderate deviations for higher values. On the other hand, the costs for an improvement step in either p or n is at most a factor of two higher CPU time. This means that the calculation and coding of expressions like Eqs. (25) is justified by the increase of precision gained. Given the goals for the precision, the user can decide to which degrees p and n the algorithm should be developed. In order to prove the precision, in Table 2 we show the convergence of our procedure for p = 2 fixed and increasing values of n. The last column shows the CPU times for 10 6 runs of the algorithm proposed in Ref. [12] with N given in the last column of the table in Ref. [12], as coded in C.

Conclusions
In this paper we developed and described an approximative algorithm for the determination of the error function which is based on geometric considerations. Along the lines explained in this paper, the algorithm can be easily implemented and extended. We have shown that each improvement step gains an improvement of the precision of at least a factor of ten, at the cost of at most a factor of two more CPU time. As a bonus, we have given a geometric derivation of Craig's integral representation of the error function and a converging power series expansion for this formula. Table 2. Results for p = 2 and increasing values of n for values of E approaching E = 1. The last column shows the CPU time for 10 6 runs according to the algorithm proposed in Ref. [12] for the values of N given in the last column of the table displayed in Ref. [12]. E = n = 0 n = 1 n = 2 n = 3 n = 4 n = 5 [12]