高斯误差函数erf的数值计算方法（C++实现）

由于毕设的数学推导中涉及了 e r f mathrm{erf} erf 函数，关于其他函数的渐近计算推导见链接类指数级数（指数积分函数的变体）数值计算算法的C++实现。

反正闲得无聊，虽然知道这种函数肯定有现成的轮子了，然而我是情报弱者。再加上最后我的算法是要在 C++ 平台上进行实现的，不如自己造一手轮子。

注意：因为我的场景只涉及 x ⩾ 0 xgeqslant0 x⩾0 的情形，所以只针对这种情况进行了考虑。事实上，根据对称性 x < 0 x<0 x<0 ，直接用 e r f ( x ) = 1 − e r f ( − x ) mathrm{erf}left(xright)=1-mathrm{erf}left(-xright) erf(x)=1−erf(−x) 即可

注：为了与已有函数区分，函数名前加了Hsk，Hsk是我的ID的缩写(羽咲->Hasaki->Hsk)

话不多说，先给代码再解释。这里的计算精度为 ε < 1 0 − 10 varepsilon<10^{-10} ε<10−10 ，如果想要修改计算精度，应当调节 HSK_EPSILON 和 HSK_INF_SERIES 这两个宏。

#include

#define HSK_EPSILON 1e-10
#define HSK_INF_SERIES 3.5  // HSK_INF_SERIES = log10(1/HSK_EPSILON) / 4 + 1
#define HSK_PI 3.14159265358979323846

double HskErf (double x) {
    double ans = 0;
    static double k = 2 / sqrt(HSK_PI);

    // 使用无穷展开式
    if (x >= HSK_INF_SERIES) {
        double x2 = - 1 / x / x;
        double an = k * exp(-x * x) / x / 2;    // a1

        double n;
        for (n = 1; n < 5; ++ n) {
            ans += an;
            an *= (n - 0.5) * x2;  // a{n}(x) <- a{n-1}(x)
        }
        return 1 - ans;
    }
    // 使用原点展开式
    else {
        double x2 = - x * x;
        double an = k * x;

        double n;
        for (n = 1; 2 * abs(an) > HSK_EPSILON; ++ n) {
            ans += an;
            an *= (1 - 1 / (n + 0.5)) * x2 / n;  // a{n+1}(x) <- a{n}(x)
        }
        return ans;
    }
}

e r f ( x ) = 2 π ∫ 0 x exp ⁡ ( − t 2 ) d t (4.1) mathrm{erf}left(xright)=frac{2}{sqrt{pi}}int_0^x{expleft(-t^2right)mathrm{d}t}tag{4.1} erf(x)=π ​2​∫0x​exp(−t2)dt(4.1)

令
a n ( x ) = ( − 1 ) n x 2 n + 1 n ! ( 2 n + 1 ) (4.2) a_nleft( x right) =frac{left( -1 right) ^nx^{2n+1}}{n!left(2n+1right)}tag{4.2} an​(x)=n!(2n+1)(−1)nx2n+1​(4.2)
有递推式
a n ( x ) = ( − 1 ) n x 2 n + 1 n ! ( 2 n + 1 ) = − ( 2 n − 1 ) x 2 ( 2 n + 1 ) n ( − 1 ) n − 1 x 2 n − 1 ( n − 1 ) ! ( 2 n − 1 ) = ( 1 − 2 2 n + 1 ) − x 2 n a n ( x ) (4.3) begin{aligned} a_nleft( x right) &=frac{left( -1 right) ^nx^{2n+1}}{n!left( 2n+1 right)}\ &=-frac{left( 2n-1 right) x^2}{left( 2n+1 right) n}frac{left( -1 right) ^{n-1}x^{2n-1}}{left( n-1 right) !left( 2n-1 right)}\ &=left( 1-frac{2}{2n+1} right) frac{-x^2}{n}a_nleft( x right)\ end{aligned}tag{4.3} an​(x)​=n!(2n+1)(−1)nx2n+1​=−(2n+1)n(2n−1)x2​(n−1)!(2n−1)(−1)n−1x2n−1​=(1−2n+12​)n−x2​an​(x)​(4.3)
原点展开式为
e r f ( x ) = 2 π ∫ 0 x exp ⁡ ( − t 2 ) d t = 2 π ∑ n = 0 ∞ ∫ 0 x ( − 1 ) n t 2 n n ! d t = 2 π ∑ n = 0 ∞ ( − 1 ) n x 2 n + 1 n ! ( 2 n + 1 ) = 2 π ∑ n = 0 ∞ a n ( x ) (4.4) begin{aligned} mathrm{erf}left( x right) &=frac{2}{sqrt{pi}}int_0^x{exp left( -t^2 right) mathrm{d}t}\ &=frac{2}{sqrt{pi}}sum_{n=0}^{infty}{int_0^x{frac{left( -1 right) ^nt^{2n}}{n!}mathrm{d}t}}\ &=frac{2}{sqrt{pi}}sum_{n=0}^{infty}{frac{left( -1 right) ^nx^{2n+1}}{n!left(2n+1right)}}\ &=frac{2}{sqrt{pi}}sum_{n=0}^{infty}{a_nleft( x right)}\ end{aligned} tag{4.4} erf(x)​=π ​2​∫0x​exp(−t2)dt=π ​2​n=0∑∞​∫0x​n!(−1)nt2n​dt=π ​2​n=0∑∞​n!(2n+1)(−1)nx2n+1​=π ​2​n=0∑∞​an​(x)​(4.4)
余项估计（假设 N + 2 > 2 x 2 N+2> 2x^2 N+2>2x2 ）(详见类指数级数（指数积分函数的变体）数值计算算法的C++实现的第三章）
∣ a N + 1 + p ( x ) ∣ = ∣ x ∣ 2 ( N + 1 + p ) + 1 ( N + 1 + p ) ! = ∣ x ∣ 2 N + 1 ( N + 1 ) ! ( x 2 ) p ( N + 2 ) ⋯ ( N + 1 + p ) ⩽ ∣ x ∣ 2 N + 1 ( N + 1 ) ! ( x 2 N + 2 ) p < ∣ a N + 1 ( x ) ∣ × ( 1 2 ) p (4.5) begin{aligned} left| a_{N+1+p}left( x right) right| &=frac{left| xright| ^{2left(N+1+pright)+1}}{left( N+1+p right) !}\ &=frac{left| x right|^{2N+1}}{left( N+1 right) !}frac{left( x^2 right)^p}{left( N+2 right) cdots left( N+1+p right)}\ &leqslant frac{left| x right|^{2N+1}}{left( N+1 right) !}left( frac{x^2}{N+2} right) ^p\ &

2 ∣ a N + 1 ( x ) ∣ ⩾ 2 ∣ a N + 1 ( N + 2 2 ) ∣ = 2 ( N + 2 2 ) 2 N + 3 ( N + 1 ) ! > 2 (4.6) begin{aligned} 2left| a_{N+1}left( x right) right|&geqslant 2left| a_{N+1}left( sqrt{frac{N+2}{2}} right) right|\ &=frac{2left( sqrt{frac{N+2}{2}} right) ^{2N+3}}{left( N+1 right) !}\ &>2\ end{aligned} tag{4.6} 2∣aN+1​(x)∣​⩾2∣∣∣∣∣​aN+1​(2N+2​ ​)∣∣∣∣∣​=(N+1)!2(2N+2​ ​)2N+3​>2​(4.6)

渐近展开
e r f ( x ) = 2 π ∫ 0 x exp ⁡ ( − t 2 ) d t = 2 π ( π 2 − exp ⁡ ( − x 2 ) ∑ n = 0 ∞ ( − 1 ) n ( 2 n − 1 ) ! ! 2 n + 1 1 x 2 n + 1 ) = 1 − 2 π exp ⁡ ( − x 2 ) ∑ n = 0 ∞ ( − 1 ) n ( 2 n − 1 ) ! ! 2 n + 1 1 x 2 n + 1 (4.8) begin{aligned} mathrm{erf}left( x right) &=frac{2}{sqrt{pi}}int_0^x{exp left( -t^2 right) mathrm{d}t}\ &=frac{2}{sqrt{pi}}left( frac{sqrt{pi}}{2}-exp left( -x^2 right) sum_{n=0}^{infty}{left( -1 right) ^nfrac{left( 2n-1 right) !!}{2^{n+1}}frac{1}{x^{2n+1}}} right)\ &=1-frac{2}{sqrt{pi}}exp left( -x^2 right) sum_{n=0}^{infty}{left( -1 right) ^nfrac{left( 2n-1 right) !!}{2^{n+1}}frac{1}{x^{2n+1}}}\ end{aligned}tag{4.8} erf(x)​=π ​2​∫0x​exp(−t2)dt=π ​2​(2π ​​−exp(−x2)n=0∑∞​(−1)n2n+1(2n−1)!!​x2n+11​)=1−π ​2​exp(−x2)n=0∑∞​(−1)n2n+1(2n−1)!!​x2n+11​​(4.8)
注意到无穷展开式并非展开阶数越高越好，当阶数较大时，系数会以比指数更快的速度发散。经过 mathematica 测试，发现展开到 n = 4 n=4 n=4 （即 x − 9 x^{-9} x−9 ）项时效果最好。因此实际渐近展开式为
e r f ( x ) ≈ 1 − 2 π exp ⁡ ( − x 2 ) [ 1 2 x − 1 4 x 3 + 3 8 x 5 − 15 16 x 7 + 105 32 x 9 ] (4.9) mathrm{erf}left( x right) approx 1-frac{2}{sqrt{pi}}exp left( -x^2 right) left[ frac{1}{2x}-frac{1}{4x^3}+frac{3}{8x^5}-frac{15}{16x^7}+frac{105}{32x^9} right] tag{4.9} erf(x)≈1−π ​2​exp(−x2)[2x1​−4x31​+8x53​−16x715​+32x9105​](4.9)
其误差近似小于 1 0 − 4 ( x − 1 ) 10^{-4left(x-1right)} 10−4(x−1) ，我们可以通过 mathematica 画图证明这一近似误差估计式，其拟合效果非常好

Plot[{Log10[
   Abs[Erf[x] - (1 - 
       2/Sqrt[Pi]
         Exp[-x^2] (1/(2 x) - 1/(4 x^3) + 3/(8 x^5) - 15/(16 x^7) + 
          105/(32 x^9)))]], -4 x + 4}, {x, 1, 5}]

也就是说，当 x ⩾ log ⁡ 10 ( 1 / ε ) 4 + 1 xgeqslant frac{log_{10}left(1/varepsilonright)}{4}+1 x⩾4log10​(1/ε)​+1 时使用无穷渐近展开，反之使用原点渐近展开。

本节采用泰勒展开的方式推导了 erf 函数的展开式，这一展开式在 x x x 较小和较大时计算速度都较快，但 2 < x < 3.5 2