Evaluation of Generalized Error Function via Fast-Converging Power Series

Serdar Beji

Faculty of Naval Architecture and Ocean Engineering, Istanbul Technical University, Istanbul, Türkiye.

DOI: 10.4236/apm.2024.146028                

Abstract

A generalized form of the error function, $G_p\left(x\right)=\frac{p}{\Gamma \left(1/p\right)}{\displaystyle {\int }_0^x{\text{e}}^{-t^p}\text{d}t}$ , which is directly associated with the gamma function, is evaluated for arbitrary real values of $p>1$ and $0<x\le +\infty$ by employing a fast-converging power series expansion developed in resolving the so-called Grandi’s paradox. Comparisons with accurate tabulated values for well-known cases such as the error function are presented using the expansions truncated at various orders.

Keywords

Generalized Error Function, Gamma Function, Grandi’s Paradox, Fast-Converging Power Series

1. Introduction

In the fourth volume of the Mécanique Céleste, Laplace (1749-1827), treating the refraction of light on the hypothesis of uniform temperature, makes use of the integral ${\displaystyle {\int }_T^{\infty }}{\text{e}}^{-t^2}\text{d}t$ and, besides stating that ${\displaystyle {\int }_0^{\infty }}{\text{e}}^{-t^2}\text{d}t=\sqrt{\pi }/2$ , develops a continued fraction solution to it ([1], pp. 485-493):

${\displaystyle {\int }_T^{\infty }}{\text{e}}^{-t^2}\text{d}t=\frac{{\text{e}}^{-T^2}/2T}{1+\frac{q}{1+\frac{2q}{1+\frac{3q}{1+\ddots }}}}$ (1)

where $q=1/2T^2$ and q should not exceed 1/4 or equivalently T should be greater than $\sqrt{2}$ to get good results; otherwise, alternative series solutions are suggested. The most commonly used form of this particular integral is known as the error function:

$\text{erf}\left(x\right)=\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^x}{\text{e}}^{-t^2}\text{d}t$ (2)

for which various series solutions and approximations are available. Two generalized forms of the error function are defined as ([2], p. 297):

$y={\text{e}}^{x^p}{\displaystyle {\int }_x^{\infty }}{\text{e}}^{-t^p}\text{d}t,y={\text{e}}^{x^{-p}}{\displaystyle {\int }_0^x}{\text{e}}^{t^p}\text{d}t.$ (3)

where $p=\mathrm{2,3,4,}\cdots$ . Euler’s integral for the gamma function ([3], p. 500) is

$\Gamma \left(z\right)={\displaystyle {\int }_0^{+\infty }}{\text{e}}^{-s}{s}^{z-1}\text{d}s,\Re \left(z\right)>0.$ (4)

Let $z=1/p$ where, like z, p may be complex although treated here as real, and introduce a change of variable such that $s=t^p$ so that $\text{d}s=p{t}^{p-1}\text{d}t$ . Noting that $s^{z-1}=s^{1/p-1}=t^{p\left(1/p-1\right)}=t^{1-p}$ and substituting the new variables into (4) give

$\left(1/p\right)\Gamma \left(1/p\right)={\displaystyle {\int }_0^{+\infty }}{\text{e}}^{-t^p}\text{d}t$ (5)

which is introduced here as a new definition of the gamma function. In accord with Equation (5), the present work goes a step further and considers the normalized integral function

$G_p\left(x\right)=\frac{1}{G_p\left(+\infty \right)}{\displaystyle {\int }_0^x}{\text{e}}^{-t^p}\text{d}t=\frac{p}{\Gamma \left(1/p\right)}{\displaystyle {\int }_0^x}{\text{e}}^{-t^p}\text{d}t$ (6)

with real $p>1$ and $0<x\le +\infty$ . The normalized integral function is designated by $G_p\left(x\right)$ due to its close affinity with the gamma (factorial) function $\Gamma \left(1/p\right)=p\cdot G_p\left(+\infty \right)$ as can be seen from (5) and (6). Furthermore, for $p=2$ Equation (6) becomes identical with the error function given in (2). Here, a fast-converging power series solution valid for the entire range of $0<x\le +\infty$ is developed for $G_p\left(x\right)$ by making use of a consistent and convergence improved series expansion introduced in the resolution of Grandi’s paradox [4]. The same approach was previously used for obtaining a series solution to the exponential integral $E_1\left(x\right)={\displaystyle {\int }_x^{\infty }}{t}^{-1}{\text{e}}^{-t}\text{d}t$ with very satisfactory results [5]. This particular integral, which is also related to the gamma function, was transformed to a form that required a fast-converging and accurate expansion of $1/\left(1+u\right)$ for its evaluation. Equation (6) is likewise elaborated into integrals containing $1/{\left(1+u\right)}^{\beta }$ with $\beta =1-1/p$ hence is evaluated with the aid of expansions introduced in Beji [4] for the resolution of Grandi’s paradox. In order to clarify the underlying methodology, a brief account of the Grandi’s paradox and its resolution is addressed next.

2. Grandi’s Paradox and Its Resolution

Grandi (1671-1742) is known by his book entitled briefly as Quadratura, which essentially lacks originality except for two particular items: construction of a curve later named as the Witch of Agnesi and statement of a paradox associated with the series expansion of $1/\left(1+u\right)$ within $0\le u\le 1$ . The paradox has become well-known and named after Grandi [6]. The series expansion of $1/\left(1+u\right)$ , either by simple division or Maclaurin series, can be easily established for $0\le u\le 1$ and a straightforward manipulation described in [4] (multiply both sides by u and redefine $u\to 1/u$ ) gives the corresponding convergent series for $1\le u\le +\infty$ :

$\frac{1}{1+u}=1-u+u^2-u^3+u^4-u^5+\cdots \text{for}0\le u\le 1$ (7a)

$\frac{1}{1+u}=\frac{1}{u}-\frac{1}{u^2}+\frac{1}{u^3}-\frac{1}{u^4}+\frac{1}{u^5}-\cdots \text{for}1\le u\le +\infty$ (7b)

The paradox arises when u is set to unity, $u=1$ , which renders the left side 1/2 while the right side becomes either 0 or 1, depending on the number of terms kept in the expansion. Beji [4] introduced the consistent truncation and convergence improvement technique to resolve the paradox. Thus, instead of equations given in (7), the following series expansions were developed.

$\frac{1}{1+u}=1+\frac{1}{2^N}{\displaystyle \sum _{m=1}^N}\left[{\displaystyle \sum _{k=m}^N}\left(\begin{array}{c}N\\ k\end{array}\right)\right]{\left(-1\right)}^m{u}^m\text{for}0\le u\le 1$ (8a)

$\frac{1}{1+u}=\frac{1}{u}+\frac{1}{2^N}{\displaystyle \sum _{m=1}^N}\left[{\displaystyle \sum _{k=m}^N}\left(\begin{array}{c}N\\ k\end{array}\right)\right]{\left(-1\right)}^m\frac{1}{u^{m+1}}\text{for}1\le u\le +\infty$ (8b)

where N is the order of truncation. The new forms (8) give exactly 1/2 when u is set to unity for all the orders of expansion. A remarkable feature of the new series expansions is the change of coefficients depending on the order of truncation. Specifically, the coefficients are re-adjusted for each truncation order to give a highly accurate expansion of the generating function. Note also that the convergence of series (8a) (8b) are ensured since the positive powers for $u\le 1$ and negative powers for $u\ge 1$ are used. As to the case $u=1$ , both expansions yield exactly 1/2 for any order of truncation N and terminate the inconsistency, which is the origin of the paradox. Applications of this technique to various Grandi-type functions with real and complex arguments can be found in [4] [7].

Figure 1 compares the third- and fourth-order series expansions of (7) and (8) against the exact expression $1/\left(1+u\right)$ . Smooth and accurate approximations obtained from (8), which are virtually inseparable from $1/\left(1+u\right)$ , are quite satisfactory when compared with (7).

3. Development of Series Solution

Integral stated in (4) is first split into two parts as

$G_p\left(+\infty \right)\cdot G_p\left(x\right)={\displaystyle {\int }_0^x}{\text{e}}^{-t^p}\text{d}t={\displaystyle {\int }_0^1}{\text{e}}^{-t^p}\text{d}t+{\displaystyle {\int }_1^x}{\text{e}}^{-t^p}\text{d}t=I_0+I_1$ (9)

where, by definition, $I_0={\displaystyle {\int }_0^1}{\text{e}}^{-t^p}\text{d}t$ and $I_1={\displaystyle {\int }_1^x}{\text{e}}^{-t^p}\text{d}t$ . Integral $I_0$ , being in the range of 0 to 1, is fairly easily evaluated by employing a straightforward series expansion of ${\text{e}}^{-t^p}$ . On the other hand, integral $I_1$ , which implicitly assumes $x>1$ , requires some novelty for its evaluation.

3.1. Evaluation of $I_0={\displaystyle {\int }_0^1}{e}^{-t^p}dt$

Using the standard series expansion of exponential function gives

Figure 1. Function $1/\left(1+u\right)$ (black), third- and fourth-order series expansions according to (0) (blue, green) and (0) (red, orange).

$\begin{array}{c}{I}_0={\displaystyle {\int }_0^1}{\text{e}}^{-t^p}\text{d}t={\displaystyle {\int }_0^1}\left(1-\frac{t^p}{1!}+\frac{t^{2p}}{2!}-\frac{t^{3p}}{3!}+\frac{t^{4p}}{4!}-\frac{t^{5p}}{5!}+\cdots \right)\text{d}t\\ ={t-\frac{t^{p+1}}{\left(p+1\right)\cdot 1!}+\frac{t^{2p+1}}{\left(2p+1\right)\cdot 2!}-\frac{t^{3p+1}}{\left(3p+1\right)\cdot 3!}+\cdots |}_0^1\\ =1-\frac{1}{\left(p+1\right)\cdot 1!}+\frac{1}{\left(2p+1\right)\cdot 2!}-\frac{1}{\left(3p+1\right)\cdot 3!}+\cdots \end{array}$ (10)

which is absolutely convergent as long as the upper limit of integration is equal to or less than one, as is the case here. A quite similar expression corresponding to the error function, $p=2$ , for an arbitrary upper limit can be found in ([2], p. 297).

3.2. Evaluation of $I_1={\displaystyle {\int }_1^x}{e}^{-t^p}dt$

In order to evaluate $I_1$ some preliminary arrangements and splitting of it into separate integrals are necessary. First, introduce a change of variable $t^p=v$ hence $\text{d}t=\text{d}v/\left[p{v}^{1-1/p}\right]$ and define $\beta =1-1/p$ for ease of notation so that

$I_1={\displaystyle {\int }_1^x}{\text{e}}^{-t^p}\text{d}t=\frac{1}{p}{\displaystyle {\int }_1^{x^p}}\frac{{\text{e}}^{-v}}{v^{\beta }}\text{d}v$ (11)

Note that $0<\beta <1$ since $p>1$ and for the special case of the error function $p=2$ hence $\beta =1/2$ . Let $v=u+1$ ,

$I_1=\frac{1}{p}{\displaystyle {\int }_0^{x^p-1}}\frac{{\text{e}}^{-\left(1+u\right)}\text{d}u}{{\left(1+u\right)}^{\beta }}=\frac{1}{p\cdot \text{e}}\left[{\displaystyle {\int }_0^1}\frac{{\text{e}}^{-u}\text{d}u}{{\left(1+u\right)}^{\beta }}+{\displaystyle {\int }_1^{x^p-1}}\frac{{\text{e}}^{-u}\text{d}u}{{\left(1+u\right)}^{\beta }}\right]$ (12)

Define the first integral within the brackets as

$I_{10}={\displaystyle {\int }_0^1}\frac{{\text{e}}^{-u}\text{d}u}{{\left(1+u\right)}^{\beta }}$ (13)

Since $0\le u\le 1$ , $I_{10}$ is going to be treated by employing series expansion (8a). The second integral within the brackets is split into two parts as

${\displaystyle {\int }_1^{x^p-1}}\frac{{\text{e}}^{-u}\text{d}u}{{\left(1+u\right)}^{\beta }}={\displaystyle {\int }_1^{x^p}}\frac{{\text{e}}^{-u}\text{d}u}{{\left(1+u\right)}^{\beta }}-{\displaystyle {\int }_{x^p-1}^{x^p}}\frac{{\text{e}}^{-u}\text{d}u}{{\left(1+u\right)}^{\beta }}=I_{11}-I_{12}$ (14)

where, by definition,

$I_{11}={\displaystyle {\int }_1^{x^p}}\frac{{\text{e}}^{-u}\text{d}u}{{\left(1+u\right)}^{\beta }},I_{12}={\displaystyle {\int }_{x^p-1}^{x^p}}\frac{{\text{e}}^{-u}\text{d}u}{{\left(1+u\right)}^{\beta }}$ (15)

$I_{12}$ is elaborated further with two successive changes of variables; first $u=y+x^p-1$ then $y=x^p\overline{u}$ , and finally $\overline{u}$ is renamed as u so that

$I_{12}={\displaystyle {\int }_{x^p-1}^{x^p}}\frac{{\text{e}}^{-u}\text{d}u}{{\left(1+u\right)}^{\beta }}=x\cdot {\text{e}}^{1-x^p}{\displaystyle {\int }_0^{1/x^p}}\frac{{\text{e}}^{-x^{p}u}\text{d}u}{{\left(1+u\right)}^{\beta }}$ (16)

Considering the entire $I_1$ ,

$\begin{array}{c}\left(p\cdot \text{e}\right)I_1={\displaystyle {\int }_0^1}\frac{{\text{e}}^{-u}\text{d}u}{{\left(1+u\right)}^{\beta }}+{\displaystyle {\int }_1^{x^p}}\frac{{\text{e}}^{-u}\text{d}u}{{\left(1+u\right)}^{\beta }}-x\cdot {\text{e}}^{1-x^p}{\displaystyle {\int }_0^{1/x^p}}\frac{{\text{e}}^{-x^{p}u}\text{d}u}{{\left(1+u\right)}^{\beta }}\\ =I_{10}+I_{11}-I_{12}\end{array}$ (17)

Integrals $I_{10}$ , $I_{11}$ , and $I_{12}$ are now evaluated in succession by employing the improved series expansions defined in (8) for $1/\left(1+u\right)$ .

3.2.1. Evaluation of $I_{10}={\displaystyle {\int }_0^1}\left[e^{-u}/{\left(1+u\right)}^{\beta }\right]du$

In order to evaluate $I_{10}$ it is necessary to express $1/{\left(1+u\right)}^{\beta }$ as a power series and to do so we first employ the improved series expansion of $1/\left(1+u\right)$ as given in Equation (8a) for the range $0\le u\le 1$ . Let $g\left(u\right)=1/\left(1+u\right)$ so that

$g\left(u\right)=\frac{1}{1+u}=1+\frac{1}{2^N}{\displaystyle \sum _{m=1}^N}\left[{\displaystyle \sum _{k=m}^N}\left(\begin{array}{c}N\\ k\end{array}\right)\right]{\left(-1\right)}^m{u}^m=C_0+{\displaystyle \sum _{m=1}^N}{C}_m{u}^m$ (18)

where, by definition,

$C_0=1,C_m=\frac{{\left(-1\right)}^m}{2^N}{\displaystyle \sum _{k=m}^N}\left(\begin{array}{c}N\\ k\end{array}\right)\text{for}m=\mathrm{1,2,}\cdots \mathrm{,}N$ (19)

We now establish $f\left(u\right)=g^{\beta }\left(u\right)=1/{\left(1+u\right)}^{\beta }={\left[C_0+{\displaystyle {\sum }_{m=1}^N}{C}_m{u}^m\right]}^{\beta }$ as a Maclaurin series truncated at a definite order N:

$f\left(u\right)=f\left(0\right)+f^{\prime }\left(0\right)u+\frac{{f^{\prime }}^{\prime }\left(0\right)}{2!}{u}^2+\frac{{{f^{\prime }}^{\prime }}^{\prime }\left(0\right)}{3!}{u}^3+\cdots +\frac{f^{\left(N\right)}\left(0\right)}{N!}{u}^N$ (20)

As a demonstration, taking $N=3$ gives

$\begin{array}{l}f=g^{\beta }=1,f^{\prime }=\beta g^{\beta -1}{g}^{\prime },{f^{\prime }}^{\prime }=\beta \left(\beta -1\right)g^{\beta -2}{g^{\prime }}^2+\beta g^{\beta -1}{g^{\prime }}^{\prime },\\ {{f^{\prime }}^{\prime }}^{\prime }=\beta \left(\beta -1\right)\left(\beta -2\right)g^{\beta -3}{g^{\prime }}^3+3\beta \left(\beta -1\right)g^{\beta -2}{g}^{\prime }{g^{\prime }}^{\prime }+\beta g^{\beta -1}{{g^{\prime }}^{\prime }}^{\prime }\end{array}$ (21)

in which the functions are to be evaluated at $u=0$ . Use of $g\left(u\right)=C_0+{\displaystyle {\sum }_{m=1}^N}{C}_m{u}^m$ results in

$\begin{array}{l}f\left(0\right)={\widehat{C}}_0=\mathrm{0<mo>!</mo>}\cdot C_0=1,f^{\prime }\left(0\right)=\beta {\widehat{C}}_0{\widehat{C}}_1,{f^{\prime }}^{\prime }\left(0\right)=\beta \left(\beta -1\right){\widehat{C}}_0{\widehat{C}}_1^2+\beta {\widehat{C}}_0{\widehat{C}}_2,\\ {{f^{\prime }}^{\prime }}^{\prime }\left(0\right)=\beta \left(\beta -1\right)\left(\beta -2\right){\widehat{C}}_0{\widehat{C}}_1^3+3\beta \left(\beta -1\right){\widehat{C}}_0{\widehat{C}}_1{\widehat{C}}_2+\beta {\widehat{C}}_0{\widehat{C}}_3\end{array}$ (22)

where ${\widehat{C}}_m=m\mathrm{<mo>!</mo>}\cdot C_m$ and $C_m$ ’s are as defined in Equation (19). $m\mathrm{<mo>!</mo>}$ multiplying $C_m$ arises from successive differentiations of the polynomial function $g\left(u\right)$ given in (18). A practical method of carrying out higher derivatives of $f=g^{\beta }$ evaluated at $u=0$ is developed in Appendix A. For ease of notation we define a set of new coefficients:

$\begin{array}{l}{D}_0=f\left(0\right)/\mathrm{0<mo>!</mo>}={\widehat{C}}_0=1\\ D_1=f^{\prime }\left(0\right)/\mathrm{1<mo>!</mo>}=\beta {\widehat{C}}_0{\widehat{C}}_1/\mathrm{1<mo>!</mo>}\\ D_2={f^{\prime }}^{\prime }\left(0\right)/\mathrm{2<mo>!</mo>}=\left[\beta \left(\beta -1\right){\widehat{C}}_0{\widehat{C}}_1^2+\beta {\widehat{C}}_0{\widehat{C}}_2\right]/\mathrm{2<mo>!</mo>}\mathrm{,}\\ D_3={{f^{\prime }}^{\prime }}^{\prime }\left(0\right)/\mathrm{3<mo>!</mo>}=\left[\beta \left(\beta -1\right)\left(\beta -2\right){\widehat{C}}_0{\widehat{C}}_1^3+3\beta \left(\beta -1\right){\widehat{C}}_0{\widehat{C}}_1{\widehat{C}}_2+\beta {\widehat{C}}_0{\widehat{C}}_3\right]/\mathrm{3<mo>!</mo>}\\ ⋮\\ D_N=f^{\left(N\right)}\left(0\right)/N!=\cdots \end{array}$ (23)

With the use of above coefficients in Equation (20) a power series for $1/{\left(1+u\right)}^{\beta }$ for the range $0\le u\le 1$ is obtained. A simple manipulation, described for Equations (7b) and (8b), gives the corresponding expression for the range $1\le u\le +\infty$ . First, changing the variable by setting $u=1/\overline{u}$ , then dividing both sides by ${\overline{u}}^{\beta }$ and finally renaming $\overline{u}$ as u again, renders the series valid for $1\le u\le +\infty$ . Thus,

$\frac{1}{{\left(1+u\right)}^{\beta }}=D_0+D_{1}u+D_2{u}^2+D_3{u}^3+\cdots +D_N{u}^N\text{for}0\le u\le 1$ (24a)

$\frac{1}{{\left(1+u\right)}^{\beta }}=\frac{D_0}{u^{\beta }}+\frac{D_1}{u^{\beta +1}}+\frac{D_2}{u^{\beta +2}}+\frac{D_3}{u^{\beta +3}}+\cdots +\frac{D_N}{u^{\beta +N}}\text{for}1\le u\le +\infty$ (24b)

The accuracy of equations given in (24) to represent $1/{\left(1+u\right)}^{\beta }$ for various $\beta =1-1/p$ values may be seen in Figure 2, where the exact expression $1/{\left(1+u\right)}^{\beta }$ and Equation (24) are depicted for $p=3/2,2,3$ , and 10 which correspond to $\beta =1/3,1/2,2/3$ , and 9/10, respectively. Figure 2 is drawn by the use of (24a) for $0\le u\le 1$ and (24b) for $1\le u\le 5$ and both series are truncated at $N=4$ .

Figure 2. Function $1/{\left(1+u\right)}^{\beta }$ (red) and Equation (24) for $N=4$ terms (blue dots) compared for $p=3/2,\beta =1/3$ ; $p=2,\beta =1/2$ ; $p=3,\beta =2/3$ ; $p=10,\beta =9/10$ values.

Equation (24a) is now used for the evaluation of $I_{10}$ for the range $0\le u\le 1$ while Equation (24b) is employed for the evaluation of $I_{11}$ in §3.2.2 for the range $1\le u\le x^p$ where $x^p>1$ and may be set to $+\infty$ .

$I_{10}={\displaystyle {\int }_0^1}\frac{{\text{e}}^{-u}\text{d}u}{{\left(1+u\right)}^{\beta }}={\displaystyle {\int }_0^1}\left(D_0+D_{1}u+D_2{u}^2+D_3{u}^3+\cdots +D_N{u}^N\right){\text{e}}^{-u}\text{d}u$ (25)

Successively integrating by parts yields

$\begin{array}{c}{I}_{10}=-{D_0{\text{e}}^{-u}|}_0^1-{D_1\left(1+u\right){\text{e}}^{-u}|}_0^1-{D_2\left(2+2u+u^2\right){\text{e}}^{-u}|}_0^1\\ -{D_3\left(6+6u+3u^2+u^3\right){\text{e}}^{-u}|}_0^1-{D_4\left(24+24u+12u^2+4u^3+u^4\right){\text{e}}^{-u}|}_0^1-\cdots \end{array}$ (26)

which may be easily generalized as

$I_{10}={\displaystyle \sum _{m=0}^N}\left(m!D_m-{\text{e}}^{-1}{\displaystyle \sum _{k=m}^N}\frac{k!}{m!}{D}_k\right)$ (27)

We have thus completed the evaluation of $I_{10}$ in terms of the coefficients of the power series given in Equation (24).

3.2.2. Evaluation of $I_{11}={\displaystyle {\int }_1^{x^p}}\left[{\text{e}}^{-u}/{\left(1+u\right)}^{\beta }\right]\text{d}u$

For the evaluation of $I_{11}$ we make use of Equation (24b) as $1\le u\le x^p$ in this case.

$I_{11}={\displaystyle {\int }_1^{x^p}}\frac{{\text{e}}^{-u}\text{d}u}{{\left(1+u\right)}^{\beta }}={\displaystyle {\int }_1^{x^p}}\left(\frac{D_0}{u^{\beta }}+\frac{D_1}{u^{\beta +1}}+\frac{D_2}{u^{\beta +2}}+\cdots +\frac{D_N}{u^{\beta +N}}\right){\text{e}}^{-u}\text{d}u$ (28)

The first integral ${\displaystyle {\int }_1^{x^p}}\left({\text{e}}^{-u}/u^{\beta }\right)\text{d}u$ is equal to $p\cdot I_1$ as can bee seen from Equation (13) and cannot be evaluated directly; therefore, we keep it as $p\cdot I_1$ . On the other hand, the remaining integrals are worked out by partial integration so that

${\displaystyle {\int }_1^{x^p}}\frac{{\text{e}}^{-u}}{u^{\beta +1}}={\text{e}}^{-1}\left[\frac{1}{\beta }\right]-x{\text{e}}^{-x^p}\left[\frac{x^{-p}}{\beta }\right]-\frac{p\cdot I_1}{\beta }$ (29a)

$\begin{array}{c}{\displaystyle {\int }_1^{x^p}}\frac{{\text{e}}^{-u}}{u^{\beta +2}}={\text{e}}^{-1}\left[\frac{1}{\beta +1}-\frac{1}{\beta \left(\beta +1\right)}\right]-x{\text{e}}^{-x^p}\left[\frac{x^{-2p}}{\beta +1}-\frac{x^{-p}}{\beta \left(\beta +1\right)}\right]\\ +\frac{p\cdot I_1}{\beta \left(\beta +1\right)}\end{array}$ (29b)

$\begin{array}{c}{\displaystyle {\int }_1^{x^p}}\frac{{\text{e}}^{-u}}{u^{\beta +3}}={\text{e}}^{-1}\left[\frac{1}{\beta +2}-\frac{1}{\left(\beta +1\right)\left(\beta +2\right)}+\frac{1}{\beta \left(\beta +1\right)\left(\beta +2\right)}\right]\\ -x{\text{e}}^{-x^p}\left[\frac{x^{-3p}}{\beta +2}-\frac{x^{-2p}}{\left(\beta +1\right)\left(\beta +2\right)}-\frac{x^{-p}}{\beta \left(\beta +1\right)\left(\beta +2\right)}\right]\\ -\frac{p\cdot I_1}{\beta \left(\beta +1\right)\left(\beta +2\right)}\end{array}$ (29c)

Continuing in this manner and adding all the results to establish $I_{11}$ as expressed in (28) give

$\begin{array}{c}{I}_{11}=p\cdot I_1\left[D_0-\frac{D_1}{\beta }+\frac{D_2}{\beta \left(\beta +1\right)}-\frac{D_3}{\beta \left(\beta +1\right)\left(\beta +2\right)}+\cdots \right]\\ +{\text{e}}^{-1}{D}_1\left[\frac{1}{\beta }\right]+{\text{e}}^{-1}{D}_2\left[\frac{1}{\beta +1}-\frac{1}{\beta \left(\beta +1\right)}\right]\\ +{\text{e}}^{-1}{D}_3\left[\frac{1}{\beta +2}-\frac{1}{\left(\beta +1\right)\left(\beta +2\right)}+\frac{1}{\beta \left(\beta +1\right)\left(\beta +2\right)}\right]+\cdots \\ -x{\text{e}}^{-x^p}\left[\frac{x^{-p}}{\beta }\right]-x{\text{e}}^{-x^p}\left[\frac{x^{-2p}}{\beta +1}-\frac{x^{-p}}{\beta \left(\beta +1\right)}\right]\\ -x{\text{e}}^{-x^p}\left[\frac{x^{-3p}}{\beta +2}-\frac{x^{-2p}}{\left(\beta +1\right)\left(\beta +2\right)}+\frac{x^{-p}}{\beta \left(\beta +1\right)\left(\beta +2\right)}\right]-\cdots \end{array}$ (30)

The general form of (30) may be expressed as

$\begin{array}{c}{I}_{11}=p\cdot I_1\left[D_0+{\displaystyle \sum _{m=1}^N}{\left(-1\right)}^m\left({\displaystyle \prod _{k=1}^m}{\left(\beta +m-k\right)}^{-1}\right)D_m\right]\\ +{\text{e}}^{-1}\left\{{\displaystyle \sum _{m=1}^N}\left[{\displaystyle \sum _{k=1}^m}{\left(-1\right)}^{k+1}\left({\displaystyle \prod _{l=1}^k}{\left(\beta +m-l\right)}^{-1}\right)\right]D_m\right\}\\ -x{\text{e}}^{-x^p}\left\{{\displaystyle \sum _{m=1}^N}\left[{\displaystyle \sum _{k=1}^m}{\left(-1\right)}^{k+1}{x}^{-p\left(m-k+1\right)}\left({\displaystyle \prod _{l=1}^k}{\left(\beta +m-l\right)}^{-1}\right)\right]D_m\right\}\end{array}$ (31)

3.2.3. Evaluation of $I_{12}=x\cdot {\text{e}}^{1-x^p}{\displaystyle {\int }_0^{1/x^p}}\left[{\text{e}}^{-x^{p}u}/{\left(1+u\right)}^{\beta }\right]\text{d}u$

To complete the evaluation of $I_1$ the only remaining integral is

$I_{12}=x\cdot {\text{e}}^{1-x^p}{\displaystyle {\int }_0^{1/x^p}}\left[{\text{e}}^{-x^{p}u}/{\left(1+u\right)}^{\beta }\right]\text{d}u$ as defined by (16). Considering only the integral itself and noting that $1/x^p<1$ , the series expansion (24a) valid for $0\le u\le 1$ can be employed.

${\displaystyle {\int }_0^{1/x^p}}\frac{{\text{e}}^{-x^{p}u}\text{d}u}{{\left(1+u\right)}^{\beta }}={\displaystyle {\int }_0^{1/x^p}}\left(D_0+D_{1}u+D_2{u}^2+D_3{u}^3+\cdots +D_N{u}^N\right){\text{e}}^{-x^{p}u}\text{d}u$ (32)

Proceeding very similar to the evaluation of $I_{10}$ , we get

$\begin{array}{c}{\displaystyle {\int }_0^{1/x^p}}\frac{{\text{e}}^{-x^{p}u}\text{d}u}{{\left(1+u\right)}^{\beta }}=-{\frac{D_0}{x^p}{\text{e}}^{-x^{p}u}|}_0^{x^{-p}}-{\frac{D_1}{x^p}\left(x^{-p}+u\right){\text{e}}^{-x^{p}u}|}_0^{x^{-p}}\\ -{\frac{D_2}{x^p}\left(2x^{-2p}+2x^{-p}u+u^2\right){\text{e}}^{-x^{p}u}|}_0^{x^{-p}}\\ -{\frac{D_3}{x^p}\left(6x^{-3p}+6x^{-2p}u+3x^{-p}{u}^2+u^3\right){\text{e}}^{-x^{p}u}|}_0^{x^{-p}}-\cdots \end{array}$ (33)

Generalizing the above expression and multiplying by $x\cdot {\text{e}}^{1-x^p}$ give $I_{12}$ as

$I_{12}=x\cdot {\text{e}}^{1-x^p}{\displaystyle \sum _{m=0}^N}\left(x^{-p\left(m+1\right)}m!D_m-{\text{e}}^{-1}{\displaystyle \sum _{k=m}^N}{x}^{-p\left(k+1\right)}\frac{k!}{m!}{D}_k\right)$ (34)

Finally, $I_1$ can be obtained by the use of Equation (17) $\left(p\cdot e\right)I_1=I_{10}+I_{11}-I_{12}$ ,

$\begin{array}{c}\left(p\cdot \text{e}\right)I_1={\displaystyle \sum _{m=0}^N}\left(m!D_m-{\text{e}}^{-1}{\displaystyle \sum _{k=m}^N}\frac{k!}{m!}{D}_k\right)\\ +p\cdot I_1\left[D_0+{\displaystyle \sum _{m=1}^N}{\left(-1\right)}^m\left({\displaystyle \prod _{k=1}^m}{\left(\beta +m-k\right)}^{-1}\right)D_m\right]\\ +{\text{e}}^{-1}\left\{{\displaystyle \sum _{m=1}^N}\left[{\displaystyle \sum _{k=1}^m}{\left(-1\right)}^{k+1}\left({\displaystyle \prod _{l=1}^k}{\left(\beta +m-l\right)}^{-1}\right)\right]D_m\right\}\\ -x{\text{e}}^{-x^p}\left\{{\displaystyle \sum _{m=1}^N}\left[{\displaystyle \sum _{k=1}^m}{\left(-1\right)}^{k+1}{x}^{-p\left(m-k+1\right)}\left({\displaystyle \prod _{l=1}^k}{\left(\beta +m-l\right)}^{-1}\right)\right]D_m\right\}\\ -x\cdot {\text{e}}^{1-x^p}{\displaystyle \sum _{m=0}^N}\left(x^{-p\left(m+1\right)}m!D_m-{\text{e}}^{-1}{\displaystyle \sum _{k=m}^N}{x}^{-p\left(k+1\right)}\frac{k!}{m!}{D}_k\right)\end{array}$ (35)

which can be re-arranged for $I_1$ as

$\begin{array}{l}p\cdot I_1\left[\text{e}-D_0-{\displaystyle \sum _{m=1}^N}{\left(-1\right)}^m\left({\displaystyle \prod _{k=1}^m}{\left(\beta +m-k\right)}^{-1}\right)D_m\right]\\ ={\displaystyle \sum _{m=0}^N}\left(m!D_m-{\text{e}}^{-1}{\displaystyle \sum _{k=m}^N}\frac{k!}{m!}{D}_k\right)\\ +{\text{e}}^{-1}\left\{{\displaystyle \sum _{m=1}^N}\left[{\displaystyle \sum _{k=1}^m}{\left(-1\right)}^{k+1}\left({\displaystyle \prod _{l=1}^k}{\left(\beta +m-l\right)}^{-1}\right)\right]D_m\right\}\\ -x{\text{e}}^{-x^p}\left\{{\displaystyle \sum _{m=1}^N}\left[{\displaystyle \sum _{k=1}^m}{\left(-1\right)}^{k+1}{x}^{-p\left(m-k+1\right)}\left({\displaystyle \prod _{l=1}^k}{\left(\beta +m-l\right)}^{-1}\right)\right]D_m\right\}\\ -x\cdot {\text{e}}^{1-x^p}{\displaystyle \sum _{m=0}^N}\left(x^{-p\left(m+1\right)}m!D_m-{\text{e}}^{-1}{\displaystyle \sum _{k=m}^N}{x}^{-p\left(k+1\right)}\frac{k!}{m!}{D}_k\right)\end{array}$ (36)

Computing $I_0$ from (10) and $I_1$ from (36) enable us to compute the integral sought $G_p\left(+\infty \right)\cdot G_p\left(x\right)={\displaystyle {\int }_0^x{\text{e}}^{-t^p}\text{d}t}=I_0+I_1$ as stated in (9).

4. Numerical Results

Formulations developed for computing $G_p\left(x\right)=\frac{p}{\Gamma \left(1/p\right)}{\displaystyle {\int }_0^x{\text{e}}^{-t^p}\text{d}t}$ are programmed in FORTRAN, which is included in Appendix B. In this section, various comparisons with available numerical values are presented. First, the tables given in ([2], pp. 310-311) for the error function $\text{erf}\left(x\right)=\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^x{\text{e}}^{-t^2}\text{d}t}$ are used for comparisons with the computed results obtained for $N=4$ . Then, quite similar comparisons, again using the tables of Abramowitz and Stegun ([2], p. 320) are presented for $G_3\left(x\right)=\frac{3}{\Gamma \left(1/3\right)}{\displaystyle {\int }_0^x{\text{e}}^{-t^3}\text{d}t}$ . For different truncation orders N and a range of p values, §4 presents a table listing computed $G_p\left(+\infty \right)$ values, all of which must be exactly unity. Two graphical representations of the normalized integral function $G_p\left(x\right)$ for $1\le p\le 10$ and $0\le x\le 5$ are also included.

4.1. Comparisons for $\text{erf}\left(x\right)=\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^x{\text{e}}^{-t^2}\text{d}t}$ or $G_2\left(x\right)=\frac{2}{\Gamma \left(1/2\right)}{\displaystyle {\int }_0^x{\text{e}}^{-t^2}\text{d}t}$

The well-known error function corresponds to $p=2$ case hence it is considered first. Table 1 compares the tabulated values of ([2], pp. 310-11), labeled “exact”, with computations from present formulation for the truncation order $N=4$ . Integral value at infinity can be obtained either by simply setting the upper limit contribution of $I_{11}$ to zero or taking a relatively large x value such as $x=100$ ; but even $x=5$ is sufficient to obtain a perfectly acceptable result for $x=+\infty$ . For the present case of $p=2$ this value corresponds to $\Gamma \left(1/2\right)/2=\sqrt{\pi }/2$ and the values in Table 1 are normalized by multiplying with $2/\Gamma \left(1/2\right)=2/\sqrt{\pi }$ . Finally, the relative error percentages are computed from 100∙(1 – Approximate/Exact).

Table 1. Comparisons for $G_2\left(x\right)=\text{erf}\left(x\right)=\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^x{\text{e}}^{-t^2}\text{d}t}$ .

	Abramowitz and Stegun [2]	Present work	Relative error
x	(Exact)	(N = 4)	(%)
0.00	0.00000000	0.00000000	−0.0000
0.25	0.27632639	0.27632639	−0.0000
0.50	0.52049988	0.52049988	−0.0000
0.75	0.71115563	0.71115563	−0.0000
1.00	0.84270079	0.84270079	−0.0000
1.25	0.92290013	0.92227506	−0.0677
1.50	0.96610514	0.96578097	−0.0336
1.75	0.98667167	0.98657212	−0.0101
2.00	0.99532227	0.99532944	−0.0072
$+\infty$	1.00000000	1.00006365	−0.0064

For $0.00\le x\le 1.00$ the computed results agree perfectly well with the tabulated values because both are obtained from the series expansion in (10) with upper limit of integration set to the given x value. In computations the number of terms was restricted to 150 due to the limited capability of the machine to compute large factorials. Tables of Abramowitz and Stegun [2] terminate at $x=2$ ; therefore, the comparisons accord with this upper limit.

For $x\ge 1$ the present formulation uses (36) with $N=4$ as no improvement in accuracy is observed for higher truncation orders $N=\mathrm{5,6,7,8}$ . However, as Table 3 in §4 reveals, higher truncation orders give improved results for higher p values, especially for $x=+\infty$ . The present comparisons have the same tendency concerning x; namely, increasing the upper limit of integration decreases the relative error percentage appreciably. For instance, the error is approximately 0.07% for $x=1.25$ ; but it reduces tenfold to −0.007% for $x=2.00$ . Despite these moderately satisfactory results, even for a truncation order of as small as $N=4$ , the lack of improvement in accuracy for higher truncation orders indicates an inherently difficult problem, which is also observed for numerical integrations performed by Simpson’s rule of integration for purposes of extra checking. High precision computations in the range $0<x\le 1$ did not yield results as good as the series solutions despite taking an extremely small step of $\Delta x=1/1000000$ .

4.2. Comparisons for $G_3\left(x\right)=\frac{3}{\Gamma \left(1/3\right)}{\displaystyle {\int }_0^x{\text{e}}^{-t^3}\text{d}t}$

Another group of tables given in ([2], p. 320] is for the integral

$G_3\left(x\right)=\frac{3}{\Gamma \left(1/3\right)}{\displaystyle {\int }_0^x{\text{e}}^{-t^3}\text{d}t}$ , which corresponds to $p=3$ in the present notation. Table 2 shows the comparisons up to $x=2.30$ , which is the highest value included in [2]. Error percentages show similar characteristics with those of $p=2$ case; yet, overall the errors are smaller. The listed highest error is 0.056% for $x=1.20$ and it reduces to 0.004% for $x=2.30$ . Again, values for $0.00\le x\le 1.00$ agree perfectly well because they are computed from the series expansion in complete agreement with the one used in [2].

Table 2. Comparisons for normalized integral $G_3\left(x\right)=\frac{3}{\Gamma \left(1/3\right)}{\displaystyle {\int }_0^x{\text{e}}^{-t^3}\text{d}t}$ .

	Abramowitz and Stegun [2]	Present work	Relative error
x	(Exact)	(N = 4)	(%)
0.00	0.0000000	0.0000000	0.0000
0.30	0.3337037	0.3337037	0.0000
0.70	0.7227669	0.7227669	0.0000
1.00	0.9042886	0.9042886	0.0000
1.20	0.9641064	0.9635709	0.0555
1.50	0.9951149	0.9949880	0.0128
1.70	0.9991499	0.9990934	0.0336
2.10	0.9999925	0.9999496	0.0057
2.30	0.9999997	0.9999569	0.0043
$+\infty$	1.0000000	0.9999573	0.0043

4.3. Computations for Various p and x Values

For different truncation orders $N=\mathrm{4,}\cdots \mathrm{,8}$ , Table 3 presents computational results for the normalized integral $G_p\left(+\infty \right)$ , which is exactly equal to unity for any p. Deviation from the unity corresponds to the absolute error in computational result. It is clear from Table 3 that, regardless of the truncation number,

Table 3. Normalized integral $G_p\left(+\infty \right)=\frac{p}{\Gamma \left(1/p\right)}{\displaystyle {\int }_0^{+\infty }{\text{e}}^{-t^p}\text{d}t}$ for various p values and truncation orders N.

p	N = 4	N = 5	N = 6	N = 7	N = 8
1.20	1.000544	0.999778	1.001037	0.998898	1.001525
1.50	1.000405	1.000205	1.001224	0.998582	1.002317
2.00	1.000064	1.000844	1.000719	0.999117	1.001955
3.00	0.999957	1.001095	1.000106	0.999937	1.000867
4.25	1.000047	1.000912	0.999939	1.000216	1.000303
5.75	1.000122	1.000685	0.999931	1.000251	1.000090
7.80	1.000161	1.000485	0.999961	1.000212	1.000013
10.00	1.000168	1.000358	0.999985	1.000167	0.999995

computational accuracy increases with increasing p values. In particular, $N=8$ case with four/five decimal places accuracy is remarkably good for $p\ge 5.75$ . On the other hand, lower truncation orders $N=\mathrm{4,5}$ appear to perform better for relatively smaller $p\le 3$ values. Nevertheless, there appears no definite trend of truncation orders. Figure 3 depicts $G_p\left(x\right)$ normalized integral function within the range $0\le x\le 5$ for five different p values. Computations for $p=1$ may be performed by setting p to a slightly higher value such as $p=1+{10}^{-6}$ ; otherwise, $\beta =1-1/p$ becomes exactly zero and the scheme fails. As expected, all the curves tend to unity with increasing x, even as low as x = 5.

A perspective view of $G_p\left(x\right)$ function for x and p pairs is presented in Figure 4. The most notable characteristic of both Figure 3 and Figure 4 is the sharp increase observed in $G_p\left(x\right)$ in the neighborhood of $x=1$ as p values increase.

Figure 3. Normalized integral function $G_p\left(x\right)=\frac{p}{\Gamma \left(1/p\right)}{\displaystyle {\int }_0^x{\text{e}}^{-t^p}\text{d}t}$ for a number of p values within $0\le x\le 5$ .

Figure 4. A perspective view of normalized integral function $G_p\left(x\right)=\frac{p}{\Gamma \left(1/p\right)}{\displaystyle {\int }_0^x{\text{e}}^{-t^p}\text{d}t}$ for $1\le p\le 10$ and $0\le x\le 5$ values.

5. Concluding Remarks

A generalized form of the error function, defined as $G_p\left(x\right)=\frac{p}{\Gamma \left(1/p\right)}{\displaystyle {\int }_0^x{\text{e}}^{-t^p}\text{d}t}$ , is

evaluated by means of fast-converging novel series expansions developed in resolving the so-called Grandi’s paradox. Computations indicate that for as low truncation order as $N=4$ , the technique introduced here gives quite acceptable results when compared with available tables; however, the lack of substantial improvement in accuracy for higher truncation orders $N=\mathrm{5,}\cdots \mathrm{,8}$ poses a drawback. On the other hand, visibly better results for higher p and N values may indicate a potential for further improvement with even higher truncation orders $N>8$ , which are not pursued here since the expansions become quite cumbersome as can be seen in the FORTRAN code for $N=8$ . A final criticism may be raised against the idea of pursuing such analytical solutions at the present advanced age of computers. It is perfectly possible to do numerical integrations with very good accuracy by various methods; but restricting our interests only to available tools would unquestionably limit the mathematical progress which may otherwise produce fruitful extensions to new areas.

Appendix A: Derivatives of $f=g^{\beta }$ for $g=C_0+{\displaystyle {\sum }_{m=0}^N}{C}_m{u}^m$

A practical methodology is developed here for obtaining higher derivatives of $f=g^{\beta }$ at $u=0$ with $g=C_0+{\displaystyle {\sum }_{m=0}^N}{C}_m{u}^m$ . The approach is essentially based on a notational procedure and considerably simpler than the direct way. Let us begin with $f\left(0\right)=g^{\beta }\left(0\right)={\widehat{C}}_0=\mathrm{0<mo>!</mo>}\cdot C_0$ , recall that in general ${\widehat{C}}_m=m\mathrm{<mo>!</mo>}\cdot C_m$ as defined in §3.2.1, and introduce a differentiation scheme as follows

${\widehat{f}}^{\prime }\left(0\right)={{\widehat{C}}^{\prime }}_0={\widehat{C}}_0{\widehat{C}}_1$ (37)

where ${\widehat{C}}_1=\mathrm{1<mo>!</mo>}\cdot C_1$ indicates the first derivative of ${\widehat{C}}_0=\mathrm{0<mo>!</mo>}\cdot C_0$ , which is no longer considered as a constant but the value of $g^{\beta }\left(0\right)$ . For the same reason, the derivative of ${\widehat{C}}_0$ is always evaluated as ${\widehat{C}}_0{\widehat{C}}_1$ which is equal to $g^{\beta }\left(0\right)g^{\prime }\left(0\right)$ . Equation (37) is not complete however; we must have the multiplier $\beta$ associated with the first derivative ${\widehat{C}}_1$ so that the correct expression becomes $f^{\prime }\left(0\right)=\beta {\widehat{C}}_0{\widehat{C}}_1$ . Nevertheless, presently we shall ignore these multipliers of $\beta$ , $\left(\beta -1\right)$ , etc. and consider them in a systematic manner afterwards. Thus, proceeding from (37),

${{\widehat{f}}^{\prime }}^{\prime }\left(0\right)={\left({\widehat{C}}_0{\widehat{C}}_1\right)}^{\prime }\mathrm{<mo>=</mo>}{{\widehat{C}}^{\prime }}_0{\widehat{C}}_1+{\widehat{C}}_0{{\widehat{C}}^{\prime }}_1={\widehat{C}}_0{\widehat{C}}_1{\widehat{C}}_1+{\widehat{C}}_0{\widehat{C}}_2={\widehat{C}}_0{\widehat{C}}_1^2+{\widehat{C}}_0{\widehat{C}}_2$ (38)

The second power appearing in ${\widehat{C}}_1^2$ indicates a double differentiation hence requires the multiplier $\beta \left(\beta -1\right)$ while the first power in ${\widehat{C}}_2$ corresponds to the multiplier $\beta$ only. The complete expression is then ${f^{\prime }}^{\prime }\left(0\right)=\beta \left(\beta -1\right){\widehat{C}}_0{\widehat{C}}_1^2+\beta {\widehat{C}}_0{\widehat{C}}_2$ as given in Equation (22). Differentiation of (38) yields

${{{\widehat{f}}^{\prime }}^{\prime }}^{\prime }\left(0\right)={\left({\widehat{C}}_0{\widehat{C}}_1^2+{\widehat{C}}_0{\widehat{C}}_2\right)}^{\prime }={\widehat{C}}_0{\widehat{C}}_1^3+3{\widehat{C}}_0{\widehat{C}}_1{\widehat{C}}_2+{\widehat{C}}_0{\widehat{C}}_3$ (39)

which in turn corresponds to

${{f^{\prime }}^{\prime }}^{\prime }\left(0\right)=\beta \left(\beta -1\right)\left(\beta -2\right){\widehat{C}}_1^3+3\beta \left(\beta -1\right){\widehat{C}}_1{\widehat{C}}_2+\beta {\widehat{C}}_3$ (40)

Note that total powers of the coefficients determines the number of derivatives; for example, ${\widehat{C}}_1^3$ has a power of 3 hence is multiplied by $\beta \left(\beta -1\right)\left(\beta -2\right)$ . Likewise, ${\widehat{C}}_1^1{\widehat{C}}_2^1$ has total $1+1=2$ powers and is multiplied by $\beta \left(\beta -1\right)$ . It is also relatively straightforward to establish the corresponding functional form from Equation (40) as

${{f^{\prime }}^{\prime }}^{\prime }\left(0\right)=\beta \left(\beta -1\right)\left(\beta -2\right)g^{\beta -3}{g^{\prime }}^3+3\beta \left(\beta -1\right)g^{\beta -2}{g}^{\prime }{g^{\prime }}^{\prime }+\beta g^{\beta -1}{{g^{\prime }}^{\prime }}^{\prime }$ (41)

Demonstration is concluded by giving the resulting expressions for the fourth-derivative.

${\widehat{f}}^{iv}\left(0\right)={\widehat{C}}_0{\widehat{C}}_1^4+6{\widehat{C}}_0{\widehat{C}}_1^2{\widehat{C}}_2+4{\widehat{C}}_0{\widehat{C}}_1{\widehat{C}}_3+3{\widehat{C}}_0{\widehat{C}}_2^2+{\widehat{C}}_0{\widehat{C}}_4$ (42)

which in turn corresponds to

$\begin{array}{c}{f}^{iv}\left(0\right)=\beta \left(\beta -1\right)\left(\beta -2\right)\left(\beta -3\right){\widehat{C}}_1^4+6\beta \left(\beta -1\right)\left(\beta -2\right){\widehat{C}}_1^2{\widehat{C}}_2\\ +4\beta \left(\beta -1\right){\widehat{C}}_1{\widehat{C}}_3+3\beta \left(\beta -1\right){\widehat{C}}_2^2+\beta {\widehat{C}}_4\end{array}$ (43)

and in terms of function g reads

$\begin{array}{c}{f}^{iv}\left(0\right)=\beta \left(\beta -1\right)\left(\beta -2\right)\left(\beta -3\right)g^{\beta -4}{g^{\prime }}^4+6\beta \left(\beta -1\right)\left(\beta -2\right)g^{\beta -3}{g^{\prime }}^2{g^{\prime }}^{\prime }\\ +4\beta \left(\beta -1\right)g^{\beta -2}{g}^{\prime }{{g^{\prime }}^{\prime }}^{\prime }+3\beta \left(\beta -1\right)g^{\beta -2}{{g^{\prime }}^{\prime }}^2+\beta g^{\beta -1}{g}^{iv}\end{array}$ (44)

Higher-order derivatives should be obtained fairly easily by following the same methodology.

Appendix B: FORTRAN Code for Computation of $G_p\left(x\right)=\frac{p}{\Gamma \left(1/p\right)}{\displaystyle {\int }_0^x{\text{e}}^{-t^p}\text{d}t}$

use msimsl

parameter(n=4)

double precision c(0:50),d(0:50)

double precision bt,cfn,ep,em,one,p,pi,x

double precision xmlt,sgn,sum,sumlt,twopwn

double precision XI,XI1,XI2,XI0,XI10,XI1a,XI1b,XI1c

c

c Pi & Exponential values

pi=4*datan(one); one=1; ep=dexp(one); em=dexp(-one)

c

c Power p and beta

p=2.0; bt=1-1/p

c

c Normalization coefficient p/Gamma(1/p)

cfn=p/dgamma(1/p)

c

c Upper limit of integration

x=100

c

c Regularized coefficients

c(0)=1; twopwn=2**dfloat(n)

c

do m=1,n; sum=0

c

do k=m,n

sum=sum+dbinom(n,k)

enddo

c

sgn=1-2*mod(m,2)

c(m)=sgn*sum/twopwn

enddo

c

c Coefficients for derivatives

do m=0,n; c(m)=dfac(m)*c(m); enddo

c

c Call Maclaurin coefficients

call mac(n,bt,c,d)

c

c Integral Parts

c

c Computation of I_0

c

if(x.le.one) then; XI0=x

do m=1,150; sgn=1-2*mod(m,2)

XI0=XI0+sgn*(x**(1+p*dfloat(m)))/(1+p*dfloat(m))/dfac(m)

enddo; goto 50

else; XI0=1.0

do m=1,150; sgn=1-2*mod(m,2)

XI0=XI0+sgn/(1+p*dfloat(m))/dfac(m)

enddo; endif

c c Computation of I_10

XI10=0.0

do m=0,n; XI10=XI10+dfac(m)*d(m)

do k=m,n; XI10=XI10-em*(dfac(k)/dfac(m))*d(k)

enddo; enddo

c

c Computation of I_11

c

c I1a First part of I1

XI1a=0.0 do m=1,n; sumlt=0.0

do k=1,m; sgn=1-2*mod(k+1,2); xmlt=1.0

do l=1,k; xmlt=xmlt/(bt+dfloat(m)-dfloat(l))

enddo; sumlt=sumlt+sgn*xmlt; enddo

XI1a=XI1a+sumlt*d(m); enddo

XI1a=em*XI1a

c

c I1b Second part of I1

XI1b=0.0; do m=1,n; sumlt=0.0

do k=1,m; sgn=1-2*mod(k+1,2); xmlt=1.0

do l=1,k; xmlt=xmlt/(bt+dfloat(m)-dfloat(l)); enddo

sumlt=sumlt+sgn*xmlt*(x**(-p*(1+dfloat(m)-dfloat(k))))

enddo; XI1b=XI1b+sumlt*d(m); enddo

XI1b=x*dexp(-x**p)*XI1b

c

c I1c Coefficient of I1

XI1c=d(0); do m=1,n; sgn=1-2*mod(m,2); xmlt=1.0

do l=1,m; xmlt=xmlt/(bt+dfloat(m)-dfloat(l)); enddo

XI1c=XI1c+sgn*xmlt*d(m); enddo

c

c Computation of I_2

XI2=0.0; do m=0,n; XI2=XI2+(x**(-p*(dfloat(m)+1.)))*dfac(m)*d(m)

do k=m,n

XI2=XI2-em*(x**(-p*(dfloat(k)+1.)))*(dfac(k)/dfac(m))*d(k)

enddo; enddo

XI2=x*dexp(1-x**p)*XI2

c

c Final Value

XI1=XI10+XI1a-XI1b-XI2

XI1=XI1/(ep-XI1c)/p

XI=XI0+XI1; goto 60

c

50 XI=XI0

c

60 write(*,40)p,x,cfn*XI

c

40 format(8(2x,f14.6))

c

stop; end

c subroutine mac(n,bt,c,d)

use msimsl

double precision bt,c(0:50),d(0:50)

c

d(0)=c(0)

c d(1)=c(1)

c

d(2)=(bt-1)*c(1)*c(1)+c(2)

c

d(3)=(bt-1)*(bt-2)*c(1)*c(1)*c(1)+3*(bt-1)*c(1)*c(2)+c(3)

c

d(4)=(bt-1)*(bt-2)*(bt-3)*c(1)*c(1)*c(1)*c(1)

d(4)=d(4)+6*(bt-1)*(bt-2)*c(1)*c(1)*c(2)

d(4)=d(4)+3*(bt-1)*c(2)*c(2)+4*(bt-1)*c(1)*c(3)+c(4)

c

d(5)=(bt-1)*(bt-2)*(bt-3)*(bt-4)*c(1)*c(1)*c(1)*c(1)*c(1)

d(5)=d(5)+10*(bt-1)*(bt-2)*(bt-3)*c(1)*c(1)*c(1)*c(2)

d(5)=d(5)+10*(bt-1)*(bt-2)*c(1)*c(1)*c(3)

d(5)=d(5)+15*(bt-1)*(bt-2)*c(1)*c(2)*c(2)

d(5)=d(5)+5*(bt-1)*c(1)*c(4)+10*(bt-1)*c(2)*c(3)+c(5)

c

d(6)=(bt-1)*(bt-2)*(bt-3)*(bt-4)*(bt-5)

d(6)=d(6)*c(1)*c(1)*c(1)*c(1)*c(1)*c(1)

d(6)=d(6)+15*(bt-1)*(bt-2)*(bt-3)*(bt-4)*c(1)*c(1)*c(1)*c(1)*c(2)

d(6)=d(6)+45*(bt-1)*(bt-2)*(bt-3)*c(1)*c(1)*c(2)*c(2)

d(6)=d(6)+20*(bt-1)*(bt-2)*(bt-3)*c(1)*c(1)*c(1)*c(3)

d(6)=d(6)+15*(bt-1)*(bt-2)*c(2)*c(2)*c(2)

d(6)=d(6)+60*(bt-1)*(bt-2)*c(1)*c(2)*c(3)

d(6)=d(6)+15*(bt-1)*(bt-2)*c(1)*c(1)*c(4)

d(6)=d(6)+10*(bt-1)*c(3)*c(3)+15*(bt-1)*c(2)*c(4)

d(6)=d(6)+6*(bt-1)*c(1)*c(5)+c(6)

c

d(7)=(bt-1)*(bt-2)*(bt-3)*(bt-4)*(bt-5)*(bt-6)

d(7)=d(7)*c(1)*c(1)*c(1)*c(1)*c(1)*c(1)*c(1)

d(7)=d(7)+21*(bt-1)*(bt-2)*(bt-3)*(bt-4)*(bt-5)*(c(1)**5)*c(2)

d(7)=d(7)+35*(bt-1)*(bt-2)*(bt-3)*(bt-4)*(c(1)**4)*c(3)

d(7)=d(7)+105*(bt-1)*(bt-2)*(bt-3)*(bt-4)*(c(1)**3)*c(2)**2

d(7)=d(7)+35*(bt-1)*(bt-2)*(bt-3)*(c(1)**3)*c(4)

d(7)=d(7)+210*(bt-1)*(bt-2)*(bt-3)*c(1)*c(1)*c(2)*c(3)

d(7)=d(7)+21*(bt-1)*(bt-2)*c(1)*c(1)*c(5)

d(7)=d(7)+105*(bt-1)*(bt-2)*(bt-3)*c(1)*c(2)**3

d(7)=d(7)+105*(bt-1)*(bt-2)*c(1)*c(2)*c(4)

d(7)=d(7)+70*(bt-1)*(bt-2)*c(1)*c(3)*c(3)+7*(bt-1)*c(1)*c(6)

d(7)=d(7)+105*(bt-1)*(bt-2)*c(2)*c(2)*c(3)+21*(bt-1)*c(2)*c(5)

d(7)=d(7)+35*(bt-1)*c(3)*c(4)+c(7)

c

d(8)=(bt-1)*(bt-2)*(bt-3)*(bt-4)*(bt-5)*(bt-6)*(bt-7)

d(8)=d(8)*c(1)*c(1)*c(1)*c(1)*c(1)*c(1)*c(1)*c(1)

&+28*(bt-1)*(bt-2)*(bt-3)*(bt-4)*(bt-5)*(bt-6)*(c(1)**6)*c(2)

d(8)=d(8)+56*(bt-1)*(bt-2)*(bt-3)*(bt-4)*(bt-5)*(c(1)**5)*c(3)

d(8)=d(8)+210*(bt-1)*(bt-2)*(bt-3)*(bt-4)*(bt-5)*(c(1)**4)*c(2)**2

d(8)=d(8)+70*(bt-1)*(bt-2)*(bt-3)*(bt-4)*(c(1)**4)*c(4)

d(8)=d(8)+560*(bt-1)*(bt-2)*(bt-3)*(bt-4)*c(1)*c(1)*c(1)*c(2)*c(3)

d(8)=d(8)+56*(bt-1)*(bt-2)*(bt-3)*c(1)*c(1)*c(1)*c(5)

d(8)=d(8)+420*(bt-1)*(bt-2)*(bt-3)*(bt-4)*c(1)*c(1)*c(2)**3

d(8)=d(8)+420*(bt-1)*(bt-2)*(bt-3)*c(1)*c(1)*c(2)*c(4)

d(8)=d(8)+280*(bt-1)*(bt-2)*(bt-3)*c(1)*c(1)*c(3)*c(3)

d(8)=d(8)+28*(bt-1)*(bt-2)*c(1)*c(1)*c(6)

d(8)=d(8)+840*(bt-1)*(bt-2)*(bt-3)*c(1)*c(2)*c(2)*c(3)

d(8)=d(8)+168*(bt-1)*(bt-2)*c(1)*c(2)*c(5)

d(8)=d(8)+280*(bt-1)*(bt-2)*c(1)*c(3)*c(4)+8*(bt-1)*c(1)*c(7)

&+105*(bt-1)*(bt-2)*(bt-3)*c(2)**4+210*(bt-1)*(bt-2)*c(2)*c(2)*c(4)

d(8)=d(8)+280*(bt-1)*(bt-2)*c(2)*c(3)**2+28*(bt-1)*c(2)*c(6)

&+56*(bt-1)*c(3)*c(5)+35*(bt-1)*c(4)*c(4)+c(8)

c

do m=1,n; d(m)=bt*d(m)/dfac(m); enddo

c

return; end

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Laplace, P.S. (1839) Méchanique Céleste, Vol. IV. Little, Brown and Company.
[2]	Abramowitz, M. and Stegun, I.A. (1972) Handbook of Mathematical Functions. Dover Publications.
[3]	Arfken, G.B. and Weber, H.J. (2005) Mathematical Methods for Physicists. Elsevier Academic Press.
[4]	Beji, S. (2020) Resolution of Grandi’s Paradox and Investigations on Related Series. Applied Mathematics E-Notes, 20, 265-277.
[5]	Beji, S. (2021) Evaluation of Exponential Integral by Means of Fast-Converging Power Series. Advances in Pure Mathematics, 11, 101-108. https://doi.org/10.4236/apm.2021.111006
[6]	O’Conner, J.J. and Robertson, E.F. (2024) McTutor History of Mathematics Archive.
[7]	Beji, S. (2020) Resolution of Grandi’s Paradox as Extended to Complex Valued Functions. Advances in Pure Mathematics, 10, 447-463. https://doi.org/10.4236/apm.2020.108027.