{'fraction_non_alphanumeric': 0.05268697374938088, 'fraction_numerical': 0.030151064883605746, 'mean_word_length': 3.858045112781955, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Evaluation of the Gauss integral

The normal or Gaussian distribution plays a prominent role in almost all fields of science. However, it is well known that the Gauss (or Euler--Poisson) integral over a finite boundary, as it is necessary for instance for the error function or the cumulative distribution of the normal distribution, cannot be expressed by analytic functions. This is proven by the Risch algorithm. Still, there are proposals for approximate solutions. In this paper, we give a new solution in terms of normal distributions by applying a geometric procedure iteratively to the problem.


Introduction
The normal or Gaussian distribution plays a prominent role in almost all fields of science, as the sum of random variables tends to the normal distribution if the quite general conditions of the central limit theorem [1] are satisfied. Besides the unbounded normal integral, the bounded integral or error function is crucial for the determination of probabilities. However, there is no analytic expression found for this function, a fact that can be tested by using the Risch algorithm [2,3]. Powerful modern computer facilities but also simple personal computers allow for a numerical calculation of the error function with any needed precision. However, if a multitude of such calculations has to be performed in a limited time, for instance, in Monte Carlo simulations, the processing time becomes essential. In order to speed up these calculations, simple and more educated approximations have been proposed in the literature. The spectrum of approximations contains, for instance, the Gaussian exponential function, including either numerical constants [5] or powers and square roots [4], approximations using ordinary and hyperbolic tangent functions [6], a rational function of an exponential function with the exponent given by a power series [7], eestidima@gmail.com. or an approximation by Jacobi theta functions [8]. Without knowing the error function explicitly, expectation values can be calculated by an approximation of the normal distribution by a series in ordinary exponential functions [9].
The present paper contains a continuation of this topic. Employing a geometric approach, we provide an approximation of the squared error function by a finite sum of N Gaussian exponential functions with different widths, where the values of which are constrained to fixed intervals. We show that, by fine-tuning these width parameters, one can optimise the precision, which, even for the leading order N = 1, is better than the error estimates given by the constraints in Ref. [4]. In addition, by choosing N as appropriately large, one can afford an arbitrary precision. On the other hand, even on a personal computer, the calculation with our leading order approximation is obtained 34 times faster than an exact numerical calculation, the processing time for higher orders being multiplied by N.
Our paper is organised as follows. In Sec. 2 we introduce the basic concepts for the calculation of the Gaussian integral that are necessary for the understanding of our geometric approach. The precisions of the leading order approximation obtained here and simple, straightforward extensions of this approximation are discussed in Sec. 3. In Sec. 4 we explain the geometric background for our approximation and provide a systematic way to create higher order approximations. The iterative construction of higher order approximations is explained in general in Sec. 5 in terms of partitions before we turn to the partition into N = 2 p intervals for increasing values of p. In Sec. 6 we explain a similar ternary construction. In Sec. 7 we provide our conclusions and an outlook on possible extensions. The convergence of our iterative procedure is discussed in more detail in Appendix A. In addition, we discuss the continuum limit, which, of course, cannot be part of the algorithm but allows, as a bonus, for a different representation of the error function.

Basic concepts
The error function is based on the standard normal density distribution which does not have a direct practical meaning, while it is desirable to evaluate the integral of this function over a bounded interval [−t, t], leading to the probability P (t) to find the result within this interval, From Eqs. (1) and (2), the square of probability is given by (3) where the integration area is a square in Fig. 1(A). Introducing polar coordinates x = r cos ϕ and y = r sin ϕ, one obtains The integral in Eq. (4) is analytically calculable if the integration is performed over the interior of some circle with radius R. Indeed, Here the function I(R) increases monotonically with R as the integral in Eq. (4) is taken over a positive function. This is why I(R = m) < P < I(R = M), with m = t and M = t √ 2, see Fig. 1(A). Therefore, Using a PC for analyzing the set of Eqs. (1), (2), and (6), one concludes that k(t) is even more constrained by 1 < k(t) < 4/π. Hence, and P M (t) − P m (t) has a maximum of 0.0592 at t = t 0 = 1.0668. The inequality (7) proves to be incomparably elegant, easy to remember, and much more accurate than the best result of Ref. [4], which, if transformed into the present formalism, will be The largest range P M (t) − P m (t) ≈ 0.330 for these constraints occurs for t = 0. Compared to this, even at the leading order observed so far, our value for P M (t) − P m (t) < 0.0592 ≪ 0.330 is more restrictive. In more detail, if, in Eq. (6), we choose k = 4/π, the error will be below 0.006, but if we take k = 1.116, the maximal error is only 0.0033. Even modern reviews on this subject do not have better results [10].

Simple extensions
By adding additional terms to the leading order approximation, one can increase the precision further. For the normal integral and the leading order approximation, for k 1 = 1.116, one has However, the precision increases by a factor of 14 ≈ 0.0033/0.00024 by using where k 1 = 1.01, k 2 = 1.23345. The next order of precision has where k 1 = 1.02335, k 2 = 1.05674, and k 3 = 1.28633. Therefore, this formula with three exponentials is at least 8 times more precise than the one with two exponentials, and it is at least 110 times more precise than Eq. (11) with one exponential only. Finally, it is 11 000 ≈ 0.33/0.00003 times more precise than the approximation in Ref. [4]. As it turns out, the values for the parameters k i for i running from 1 to N take values between 1 and √ 2, while the sum of the exponential factors is divided by N. Still, there is a degree of arbitrariness in the determination of these parameters. In order to remove this arbitrariness, in the following, we develop an iterative method based on geometry.

Geometric background of our procedure
In order to understand our method, we refer to Fig. 1(B) for the first step. Starting with the square with side length 2t representing the square P 2 (t) of the probability, we turn this square by an angle of π/4 to obtain P 2 (t) again. In union and intersection, these two overlayed squares construct two 8-angle figures. In order to obtain the area of the larger figure, one has to subtract the area of the smaller figure from the twofold square area, as this smaller figure is covered twice by the two squares. Accordingly, for the integrals over the probability density, one obtains the relation between the probabilities. Now, and the angle θ = π/8 is enclosed between the x axis and the vector u shown in Fig. 1(B). We study ∆(k 1 , k 2 , t) = P (k 1 , k 2 , t) − P (t) with (17) P (k 1 , k 2 , t) = 1 − 1 2 e −k 2 1 t 2 /2 + e −k 2 2 t 2 /2 . Drawing three-dimensional graphics and looking for a minimum of |∆(k 1 , k 2 , t)|, one obtains where k 1 = 1.01 and k 2 = 1.23345. This is the starting point.

Basic construction of the procedure
In order to construct the iteration, one performs a partition of the figure describing ω(t), Ω(t), or both of these, by repeating the geometric construction shown before. For instance, taking only the larger 8-angle figure describing Ω(t), one can turn this figure by an angle θ = π/16 and overlay the new figure with the old one. In doing so, one can separate a new larger and smaller 16-angle figure in the same way as was carried out before for the 8-angle figures. Accordingly, by geometric means, one obtains new constraints. In order to describe the procedure in a unique way, in each iterative step, we rename k n by k 2n , and, if this new k 2n is subject to a partition, the smaller and larger figure of this partition are related to the values k 2n−1 and k 2n , respectively.
Using the case in the previous section as an illustrative example for the procedure, we might keep the smaller 8-angle figure related to ω(t) but apply a partition to the larger 8-angle figure related to Ω(t). Accordingly, k 1 is replaced by k 2 and k 2 is replaced by k 4 , but this new k 4 is again split up into k 3 and k 4 . The constraint for the lowest parameter k 2 (the former k 1 ) remains the same, The intervals are consecutive, but π/8 is replaced by 2π/16 in order to indicate the new partition. Finally, the upper limit stays at 1/ cos(4π/16) = 1/ cos(π/4) = √ 2. For these values k 2 , k 3 , and k 4 , one obtains the approximation because of the geometric transformations of Fig. 1 and the corresponding double use of Eq. (14). Note that the set of parameters k 2 , k 3 , k 4 is different from the set k 1 , k 2 and k 3 in Eq. (13). Indeed, if, for Eq. (23), one uses k 2 = 1.025187, k 3 = 1.1249, and k 4 = 1.31336, the precision improves to 0.000015. Still, it is obvious that this example is only half of an iteration step, and one could do definitely achieve a greater result by also performing the partition for the smaller 8-angle figure, leading to four parameters separated uniformly, Therefore, a full iteration step is increasing the number of parameters k n by a factor of two, and, after p full iteration steps, one has N = 2 p parameters. Each iteration step is finalized by optimizing the N (or less) parameters k n . For any finite (or even very large) N, the constraints (25) k min n (N) ≤ k n ≤ k max n (N) with k min n (N) = k max n−1 (N) can be calculated from geometry observations in a similar fashion. In practice, for a small set of parameters, we use a graphical method. For instance, the method applied to obtain the three values k 2 , k 3 , and k 4 in Eq. (23) was to look for the solution of the system of three equations The values t = 1, √ 2, and 2 are used as nodes for this approximation. Their choice depends on the application of the approximation and has to be adjusted to the number of width parameters to be determined. Each equation in (26) can be treated individually. Therefore, the solution is very easy to find. From Q(t = 1) = 0, one extracts the function k 2 = k 2 (k 3 , k 4 ). Inserting these solutions into Q(t = √ 2) = 0, one extracts the two positive functions k 3 = k 3 (k 4 ) and k 3 =k 3 (k 4 ). Inserting these solutions into Q(t = 2) = 0 and plotting the function Q(t = 2, k 4 ), one finds the position of the zero, which proves to be k 4 = 1.31336. Using this knowledge, one obtains k 3 = k 3 (k 4 ) and k 2 = k 2 (k 3 , k 4 ) as well. However, as k 3 =k 3 (k 4 ) is given for k 4 < 1 only, this is not a valid solution, as k n ≥ 1 for all n. Note that the graphical method cannot be applied any more for N ≥ 4. Instead, we used a random number generator to create values for the parameters k n in the respective intervals in Eq. (24). Proceeding in this way, for N = 4 (p = 2), we obtain the values k 1 = 1.00725, k 2 = 1.04665, k 3 = 1.12192, and k 4 = 1.3129, and a precision of 0.00001, which is, again, the lowest precision for a given N. As becomes obvious, the lowest precisions are obtained for uniform partitions. This is not only the case for N being a power of 2 but also for N being a power of 3, as discussed in the next section. Figure 2. The case of three "boxes" and accordingly three approximating exponents. The picture becomes more and more rotationally symmetric as the number of boxes grows.

A similar ternary procedure
As the approximation (13) gained high precision, we tried and succeeded in finding a geometric interpretation for this, as is shown in Fig. 2. In this ternary approach, the initial step is to rotate the square not by an angle of π/4 as in the previous approach but by an angle of π/6. The overlapping squares in Fig. 2 can be split up into three 12-angle figures that, at the same time, determine the constraints for the parameters k i , (28) 1 ≤ k 1 ≤ 1/ cos(π/12) ≤ k 2 ≤ 1/ cos(π/6) ≤ k 3 ≤ 1/ cos(π/4).
Note that the values k 1 = 1.02335, k 2 = 1.05674 and k 3 = 1.28633 chosen in Eq. (13) fit into these intervals. This procedure can be continued iteratively in a ternary way, i.e., turning the 12-angle figures by an angle of π/18, and generally by the angle α = π/(2 · 3 p ). In the next section, we deal with the convergence of this and the previous procedure for increasing values of p.

Conclusions and Outlook
In this paper, we have given an approximation for the Gauss integral with a finite boundary in terms of the square root of a normalized sum of normal distributions plus one, each of those distributions depending on the (symmetric) boundary [−t, t] of the integral and a set of maximally N parameters k n . By simple geometrical means, it is shown that these parameters are constrained to intervals given by the inverse cosine with equally distributed angles. We performed this approximation procedure in both a binary (N = 2 p ) and a ternary way (N = 3 p ) and showed that the procedure converges for an increasing degree p. The continuum limit leads to a further approximation. ropean Regional Development Fund under Grant No. TK133. For geometry reasons, in the limit N → ∞, the largest parameter in the infinite set {k 1 , k 2 , . . . , k ∞ } must be √ 2, whereas the lowest one must be 1. The reason is that, in using the technique as in Fig. 1(B) over and over again, the final areas of integration turn to perfect circles between the radii t and t √ 2. The convergence of this method becomes obvious by considering the backstep iteration. Suppose we start with an approximation for a given set of parameters k n with a given precision. The degeneration of two adjacent parameters means that a partition is skipped, leading to a more imprecise approximation as the degree of freedom in choosing different parameter values is lost.
For the general analysis we calculate the convergence by fixing the parameters in Eq. (29) to the upper boundary, k n = 1/ cos(πn/(4 · 2 p )), and analyse (30) ∆ N (N, t) = ∆(k 1 , . . . , k N , t) = P (k 1 , . . . , k N , t) − P (t) for N = 2 p and a fixed value of t, e.g., t 0 = 1.0668, at which, the uncertainty range of Eq. (7) turns out to be maximal. One obtains the values in Table 1 demonstrating the convergence of the approximations. For the ternary procedure we again fix the parameters to the upper boundary, k n = 1/ cos(πn/(4 · 3 p )), and analyse ∆ N (N, t) for N = 3 p and for the same fixed value t 0 = 1.0668. One obtains the values in Table 2. Table 2. Deviations in the uniform N = 3 p approximations for increasing p. Note that values higher than p = 10 could not be checked with the PC at hand.  Tables 1 and 2 can be approximated by the formula |∆ p (N, 1)| < 0.09/N, i.e., the deviation is inversely proportional to N. This can be seen as follows. The worst error of the squared Gauss integral P 2 (t) is given by using the approximations where the k n takes the maximal or minimal values, respectively. The difference between these squares of extremal values is given by Finally, one can use P (t) > P min (t) = √ 1 − e −t 2 /2 to obtain where H(t) √ 1 − e −t 2 /2 ≤ 2 2 3 3 /5 5 is used. The error for P (t) itself is, at most, the difference between the two extremal values. Eq. (29) can be considered as the discretised form of the Gauss integral. Applying the continuum limit i f i (z i )∆z i → f (z)dz, one obtains (34) P 2 (t) = 1 − 4 π π/4 0 exp −t 2 2 cos 2 φ dφ.
The exponential function can be expanded into a series of finite degree N. Again, we obtain an approximation, as