(35S) Erf(x)

05092020, 10:56 AM
(This post was last modified: 05092020 06:23 PM by lipoff.)
Post: #1




(35S) Erf(x)
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 builtin function. I could use the builtin 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 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. 

05232021, 08:52 PM
(This post was last modified: 05232021 11:18 PM by Albert Chan.)
Post: #2




RE: (35S) Erf(x)
(05092020 10:56 AM)lipoff Wrote: Erf(x) ≃ 1  (a₁t + a₂t + a₃t³) * Exp(x²) 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) (05212021 07:32 PM)Albert Chan Wrote: sqrt(pi)/2 * erf(w/2)/w Doing this with erf approximation formula (for real x, ε(x) ≤ 25e6) >>> 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 ≈ 56e6 + 132e6j, not as tight as for real input, but not bad ! Update: for comparison, with better Hastings formula (for real x, ε(x) ≤ 0.15e6) >>> 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.20e6 + 0.16e6j 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)