互补误差函数erfc是与标准正态分布密切相关的特殊函数.它经常用于统计学和自然科学(例如扩散问题),其中需要考虑该分布的"尾部",因此使用误差函数erf是不合适的.
互补误差函数是在ISO C99标准数学库提供作为功能erfcf,erfc和erfcl; 这些随后也被采用到ISO C++中.因此,源代码可以很容易地在该库的开源实现中找到,例如在glibc中.
然而,许多现有的实现本质上是标量的,而现代处理器硬件是面向SIMD的(显式地,如在x86 CPU中,或隐式地,如在GPU中).出于性能原因,因此非常需要可矢量化的实现.这意味着需要避免分支,除非作为选择分配的一部分.同样,没有指出表的广泛使用,因为并行查找通常是低效的.
如何构建单精度函数的高效矢量化实现erfcf()？测量的准确度ulp应与glibc的标量实现大致相同,其最大误差为3.12575 ulps(通过详尽测试确定).可以假设融合乘法加法(FMA)的可用性,因为此时所有主要处理器架构(CPU和GPU)都提供它.处理浮点状态标志时errno可以忽略,非正规,无穷大和NaN应根据ISO 75的IEEE 754绑定进行处理.
在matlab中有一个特殊的函数,在我所知道的Python(numpy,scipy,mpmath,...)的任何集合中都没有.
可能还有其他地方可以找到像这样的功能吗？
UPD对于所有认为这个问题很简单的人,请首先尝试为参数~30计算此函数.
UPD2任意精度是一个很好的解决方法,但如果可能的话我宁愿避免它.我需要一个"标准"的机器精度(不多也不少)和最高速度.
UPD3事实证明,mpmath结果令人惊讶地不准确.即使标准python math工作,mpmath结果也会更糟.这使它绝对毫无价值.
UPD4用于比较计算erfcx的不同方法的代码.
import numpy as np
def int_erfcx(x):
"Integral which gives erfcx"
from scipy import integrate
def f(xi):
return np.exp(-x*xi)*np.exp(-0.5*xi*xi)
return 0.79788456080286535595*integrate.quad(f,
0.0,min(2.0,50.0/(1.0+x))+100.0,limit=500)[0]
def my_erfcx(x):
"""M. M. Shepherd and J. G. Laframboise,
MATHEMATICS OF COMPUTATION 36, 249 (1981)
Note that it is reasonable to compute it in long double
(or whatever python has)
"""
ch_coef=[np.float128(0.1177578934567401754080e+01),
np.float128( -0.4590054580646477331e-02),
np.float128( -0.84249133366517915584e-01),
np.float128( 0.59209939998191890498e-01),
np.float128( -0.26658668435305752277e-01),
np.float128( …