在python中是否有可伸缩的互补误差函数？

Question

在matlab中有一个special function，这在我所知道的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(   0.9074997670705265094e-02),
             np.float128(  -0.2413163540417608191e-02),
             np.float128(    0.490775836525808632e-03),
             np.float128(    -0.69169733025012064e-04),
             np.float128(      0.4139027986073010e-05),
             np.float128(       0.774038306619849e-06),
             np.float128(      -0.218864010492344e-06),
             np.float128(        0.10764999465671e-07),
             np.float128(         0.4521959811218e-08),
             np.float128(         -0.775440020883e-09),
             np.float128(          -0.63180883409e-10),
             np.float128(           0.28687950109e-10),
             np.float128(             0.194558685e-12),
             np.float128(            -0.965469675e-12),
             np.float128(              0.32525481e-13),
             np.float128(              0.33478119e-13),
             np.float128(              -0.1864563e-14),
             np.float128(              -0.1250795e-14),
             np.float128(                 0.74182e-16),
             np.float128(                 0.50681e-16),
             np.float128(                 -0.2237e-17),
             np.float128(                 -0.2187e-17),
             np.float128(                    0.27e-19),
             np.float128(                    0.97e-19),
             np.float128(                     0.3e-20),
             np.float128(                    -0.4e-20)]
    K=np.float128(3.75)
    y = (x-K) / (x+K)
    y2 = np.float128(2.0)*y
    (d, dd) = (ch_coef[-1], np.float128(0.0))
    for cj in ch_coef[-2:0:-1]:             
        (d, dd) = (y2 * d - dd + cj, d)
    d = y * d - dd + ch_coef[0]
    return d/(np.float128(1)+np.float128(2)*x)

def math_erfcx(x):
    import scipy.special as spec
    return spec.erfc(x) * np.exp(x*x)

def mpmath_erfcx(x):
    import mpmath
    return mpmath.exp(x**2) * mpmath.erfc(x)

if __name__ == "__main__":
    x=np.linspace(1.0,26.0,200)
    X=np.linspace(1.0,100.0,200)

    intY  = np.array([int_erfcx(xx*np.sqrt(2)) for xx in X])
    myY   = np.array([my_erfcx(xx) for xx in X])
    myy   = np.array([my_erfcx(xx) for xx in x])
    mathy = np.array([math_erfcx(xx) for xx in x])
    mpmathy = np.array([mpmath_erfcx(xx) for xx in x])
    mpmathY = np.array([mpmath_erfcx(xx) for xx in X])

    print ("Integral vs exact: %g"%max(np.abs(intY-myY)/myY))
    print ("math vs exact:     %g"%max(np.abs(mathy-myy)/myy))
    print ("mpmath vs math:    %g"%max(np.abs(mpmathy-mathy)/mathy))
    print ("mpmath vs integral:%g"%max(np.abs(mpmathY-intY)/intY))

exit()

对我来说，它给了我

^{pr2}$

显然，math在它工作的地方提供了最好的精度，而mpmath在math起作用的地方给出了更大数量级的误差偶，对于更大的参数，甚至更高。在

Answer 1

一个高度优化的C++实现Erfcx（对于真实和复杂的参数）最近是{A1}，应该在SciPy版本0.12中。在

Answer 2

下面是一个简单而快速的实现，在全球范围内提供12-13位精度：

^{1}$

Answer 3

gmpy2库提供对MPFR多精度库的访问。对于正常精度，它几乎比mpmath快5倍。在

^{1}$

两个库返回30的相同结果。在

^{pr2}$

免责声明：我维护gmpy2。我正在积极努力开发一个新的版本，但是对于这个计算来说，当前版本不应该有任何问题。在

编辑：我很好奇调用两个函数而不是只调用一个函数的开销，所以我完全用C实现了gmpy2.erfcx（），但仍然使用MPFR来执行计算。改进的程度低于我的预期。如果您认为erfcx（）有用，我可以将其添加到下一个版本中。在

$ py27 -m timeit -s "import gmpy2" "gmpy2.erfcx(30)"
100000 loops, best of 3: 9.45 usec per loop