(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.