In matlab there is a special function which is not available in any of the collections for the Python I know (numpy, scipy, mpmath, ...).
Probably there are other places where functions like this one may be found?
UPD For all who think that the question is trivial, please try to compute this function for argument ~30 first.
UPD2 Arbitrary precision is a nice workaround, but if possible I would prefer to avoid it. I need a "standard" machine precision (no more no less) and maximum speed possible.
UPD3 It turns out, mpmath gives surprisingly inaccurate result. Even where standard python math works, mpmath results are worse. Which makes it absolutely worthless.
UPD4 The code to compare different ways to compute 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()
For me, it gives
Integral vs exact: 6.81236e-16
math vs exact: 7.1137e-16
mpmath vs math: 4.90899e-14
mpmath vs integral:8.85422e-13
Obviously, math gives best possible precision where it works while mpmath gives error couple orders of magnitude larger where math works and even more for larger arguments.
Here is a simple and fast implementation giving 12-13 digit accuracy globally:
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