python - 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()

对我来说，它给了

Integral vs exact: 6.81236e-16
math vs exact:     7.1137e-16
mpmath vs math:    4.90899e-14
mpmath vs integral:8.85422e-13

很明显，math 在它工作的地方给出了最好的精度，而 mpmath 在 math 工作的地方给出了几个数量级的误差，甚至更多更大的争论。

Answer 1

这是一个简单快速的实现，可在全局范围内提供 12-13 位数字的精度:

from scipy.special import exp, erfc

def erfcx(x):
    if x < 25:
        return erfc(x) * exp(x*x)
    else:
        y = 1. / x
        z = y * y
        s = y*(1.+z*(-0.5+z*(0.75+z*(-1.875+z*(6.5625-29.53125*z)))))
        return s * 0.564189583547756287