(35S) Erf(x)

[lipoff]

I needed to be able to compute the error function, Erf(x) on my HP 35s. Since it's not a WP 31s or a WP 34s, Erf(x) is not a built-in function. I could use the built-in numerical integrator to calculate the definite integral:

Erf(x) = ∫ 2/√π * (e^(-T^2)) dT, from T = 0 to T = x

But it is quite slow; it takes about 6 seconds on FIX 4 and almost 25 seconds on FIX 9.

So I got out my copy of the Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables by Abramowitz and Stegun --- better known as Abramowitz and Stegun --- and found a polynomial approximation to Erf(x) with enough precision for my purposes:

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

Let's implement that now as an RPN program:

Code:

E001     LBL E
E002     STO X
E003     0.3480242
E004     0.47047
E005     RCL X
E006     *
E007     1
E008     +
E009     1/x
E010     *
E011     0.0958798
E012     0.47047
E013     RCL X
E014     *
E015     1
E016     +
E017     1/x
E018     x^2
E019     *
E020     -
E021     0.7478556
E022     0.47047
E023     RCL X
E024     *
E025     1
E026     +
E027     1/x
E028     3
E029     y^x
E030     *
E031     +
E032     RCL x
E033     x^2
E034     +/-
E035     e^x
E036     *
E037     1
E038     x<>y
E039     -
E040     RTN

Of course, on the HP 35s, you can also implement it as an equation in algebraic mode:

Code:

1-(0.348024*(1÷(1+0.47047*X))-0.0958798*SQ(1÷(1+0.47047*X))+0.7478556*(1÷(1+0.47047*X))^3*EXP(⁻SQ(X)))

It's just slightly faster in the RPN implementation; just over 1 second in RPN and just under 2 seconds in algebraic. Both are much faster than evaluating the integral!

This is one thing I really enjoy about the HP 35s --- there are so many ways to accomplish the same thing . . . but almost none of them are equally fast.

Anyway, I hope this is helpful if anyone else needs to calculate Erf(x), and note that Abramowitz and Stegun does have more accurate (but more complicated) polynomial approximations to Erf(x) and many other functions should you need to implement them. You can find an online copy as well.

[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