(35S) Erf(x)

[Albert Chan]

(05-09-2020 10:56 AM)lipoff Wrote:  Erf(x) ≃ 1 - (a₁t + a₂t + a₃t³) * Exp(-x²)
where t = 1/(1+px), a₁ = 0.3480242, a₂ = -0.0958798, a₃ = 0.7478556, and p = 0.47047

It seems formula also work with complex numbers, although not as good.
The bigger and more "complex" the input (phase angle close to ±pi/2), the worse the estimate.

I plotted |z| ≤ 3, 0 ≤ arg(z) < pi/2, estimate matched true erf(z) well.
Since formula is really ratio of 2 infinite length polynomials, we have erf_est(conj(z)) = conj(erf_est(z))
Since formula is designed for real x ≥ 0. For z with Re(z)<0, we compute via -erf(-z)
Thus, the plots really checked full circle.

Example, from my recent erf calculation: (|w| = 1, arg(w) = -pi/4)

(05-21-2021 07:32 PM)Albert Chan Wrote:  sqrt(pi)/2 * erf(w/2)/w
= 1/2                - (-i)/(3*1!*2^3)     + (-1)/(5*2!*2^5)    - (+i)/(7*3!*2^7)
+ 1/(9*4!*2^9) - (-i)/(11*5!*2^11) + (-1)/(13*6!*2^13) - (+i)/(15*7!*2^15) + ...
= (1/2 - 1/320 + 1/110592 - 1/76677120 + ...)
+ (1/24 - 1/5376 + 1/2703360 - 1/2477260800 + ...) * i
≈ 0.496884029215 + 0.0414810242685*i

Doing this with erf approximation formula (for real x, |ε(x)| ≤ 25e-6)

>>> from numpy import polyval, exp, sqrt, pi
>>> coef = [0.7478556, -0.0958798, 0.3480242, 0.]
>>> erf_est = lambda x: 1 - exp(-x*x) * polyval(coef, 1/(1+0.47047*x))
>>> w = exp(1j*-pi/4)

>>> sqrt(pi)/2 * erf_est(w/2)/w
(0.49682810049955073+0.041349472924459674j)

Error ≈ 56e-6 + 132e-6j, not as tight as for real input, but not bad !

Update: for comparison, with better Hastings formula (for real x, |ε(x)| ≤ 0.15e-6)

>>> coef2 = [1.061405429, -1.453152027, 1.421413741, -0.284496736, 0.254829592, 0.]
>>> erf2_est = lambda x: 1 - exp(-x*x) * polyval(coef2, 1/(1+0.3275911*x))

>>> sqrt(pi)/2 * erf2_est(w/2)/w
(0.49688622537023164+0.0414808658168212j)

Error ≈ -2.20e-6 + 0.16e-6j