StackOverflow2013

Note that there are some explanatory texts on larger screens.

plurals

POIs there a scaled complementary error function in python available?
text
Body
copied!<p>In matlab there is a <a href="http://www.mathworks.com/help/techdoc/ref/erfcx.html" rel="nofollow">special function</a> which is not available in any of the collections for the Python I know (numpy, scipy, mpmath, ...). </p> <p>Probably there are other places where functions like this one may be found?</p> <p><strong>UPD</strong> For all who think that the question is trivial, please try to compute this function for argument ~30 first. </p> <p><strong>UPD2</strong> 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. </p> <p><strong>UPD3</strong> It turns out, <code>mpmath</code> gives surprisingly inaccurate result. Even where standard python <code>math</code> works, <code>mpmath</code> results are worse. Which makes it absolutely worthless. </p> <p><strong>UPD4</strong> The code to compare different ways to compute erfcx.</p> <pre><code>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() </code></pre> <p>For me, it gives </p> <pre><code>Integral vs exact: 6.81236e-16 math vs exact: 7.1137e-16 mpmath vs math: 4.90899e-14 mpmath vs integral:8.85422e-13 </code></pre> <p>Obviously, <code>math</code> gives best possible precision where it works while <code>mpmath</code> gives error couple orders of magnitude larger where <code>math</code> works and even more for larger arguments. </p>

Querying!

Guidance

An individual column

Larger individual text columns get their own page to allow for proper reading.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload