高精度误差函数计算

1 背景

有时需要进行高精度的误差函数计算,这时就需要误差函数的计算方式。

2 误差函数

2.1 性质

误差函数定义如下  
(1)
$$
\mathbf{erf}(x)=\dfrac{2}{\sqrt{\pi}}\int_0^x\exp(-t^2)dt
$$
而对于误差函数有如下性质  
(2)
$$
\begin{aligned}
\mathbf{erf}(x)&=-\mathbf{erf}(-x)
\\
\mathbf{erf}(0)&=0
\\
\mathbf{erf}(\infty)&=1
\\
-1 < \mathbf{erf}(x) &< 1
\\
\dfrac{\partial \mathbf{erf}(x)}{\partial x}&=\dfrac{2}{\sqrt{\pi}}\exp(-x^2)
\\
|\mathbf{erf}(20)-1| &< 10^{-175}
\\
|\mathbf{erf}(30)-1| &< 10^{-392}
\\
|\mathbf{erf}(40)-1| &< 10^{-696}
\\
|\mathbf{erf}(50)-1| &< 10^{-1087}
\end{aligned}
$$

2.2 函数展开

这里给出误差函数两种展开形式。  
(3)
$$
\begin{aligned}
\mathbf{erf}1(x)&=1-\dfrac{1}{e^{z^2}\sqrt{\pi}g(x)}
\\
&=1-\dfrac{1}{\sqrt{\pi}\exp(x^2+\ln(g(x)))}
\\
\\
g(x)&=x+\dfrac{1}{2(x+\dfrac{1}{x+\dfrac{3}{2(x+\dfrac{2}{x+\dfrac{5}{2(x+\dfrac{3}{x+\dfrac{7}{2(x+\dfrac{4}{x+\dfrac{9}{…}})}})}})}})}
\end{aligned}
$$

(4)
$$
\begin{aligned}
\mathbf{erf}2(x)&=\dfrac{2x}{\exp(x^2)\sqrt{\pi}(1-\dfrac{2x^2}{3+\dfrac{4x^2}{5-\dfrac{6x^2}{7+\dfrac{8x^2}{9-\dfrac{10x^2}{11+\dfrac{12x^2}{13-…}}}}}})}
\end{aligned}
$$

上面两种展开方式,计算误差函数时,收敛效率有差异。它们分别适用不同场景。

2.3 迭代计算

对于\(\mathbf{erf}1(x)\)可以构建如下的迭代算式  
(5)
$$
\begin{cases}
y_m=x+\dfrac{2m-1}{2x}
\\
\\
y_i=x+\dfrac{2i-1}{2(x+\dfrac{i}{y_{i+1}})} ,i=m-1,m-2,…,1
\\
\\
\mathbf{erf}1(x)=1-\dfrac{1}{y_1\sqrt{\pi}\exp(x^2)}=1-\dfrac{1}{\sqrt{\pi}\exp(x^2+\ln(y_1))}
\end{cases}
$$

对于\(\mathbf{erf}2(x)\)可以构建如下的迭代算式  
(6)
$$
\begin{cases}
y_m=2m-1
\\
\\
y_i=2i-1+(-1)^i\dfrac{2ix^2}{y_{i+1}},i=m-1,m-2,…,1
\\
\\
\mathbf{erf}2(x)=\dfrac{2x}{y_1\sqrt{\pi}\exp(x^2)}
\end{cases}
$$

2.4 截断与计算

由于\(\mathbf{erf}1(x)\)与\(\mathbf{erf}2(x)\)收敛效率有差异,因此这里根据经验,设有效位数为\(p\),则截断项数\(m\)根据不同的展开式进行获取  
(7)
$$
\mathbf{erf}(x)=\begin{cases}
x \leq 1 &,
\mathbf{erf}2(x),m=\lfloor 51.1267+0.35426p\rfloor
\\
1< x \leq 2 &,\mathbf{erf}2(x),m=\lfloor 83.1015+0.42509p\rfloor
\\
2< x \leq 3 &,\mathbf{erf}2(x),m=\lfloor 98.1029+0.49867p\rfloor
\\
3< x \leq 4 &,
\begin{cases}
q\leq 70&, \mathbf{erf}1(x) &, m=\lfloor 0.03074 + 0.02114p^2\rfloor
\\
q>70&, \mathbf{erf}2(x) &, m=\lfloor 130.000+0.55234p \rfloor
\end{cases}
\\
4< x \leq 5 &,
\begin{cases}
q\leq 110&, \mathbf{erf}1(x) &, m=\lfloor (0.119868p – 0.459333)^2 +2.409230\rfloor
\\
q>110&, \mathbf{erf}2(x) &, m=\lfloor 146.1909+ 0.607895p \rfloor
\end{cases}
\\
5< x \leq 10 &,
\begin{cases}
q\leq 450&, \mathbf{erf}1(x) &, m=\lfloor 0.003320773p^2-6.645615\rfloor
\\
q>450&, \mathbf{erf}2(x) &, m=\lfloor 493.4462 + 0.647483p \rfloor
\end{cases}
\\
10< x \leq 20 &,
\begin{cases}
q\leq 1790&, \mathbf{erf}1(x) &, m=\lfloor 0.0008286871p^2-25.99427\rfloor
\\
q>1790&, \mathbf{erf}2(x) &, m=\lfloor 1072.955+0.869305 p \rfloor
\end{cases}
\end{cases}
$$

因此,最终可以根据(3-7)结合性质(2)来进行误差函数计算。