Accurate Normal Distribution for the HP67/97

[Albert Chan]

We can use Taylor expansion to calculate Z(x+h) with 1 Exp call

Z(x+h) = Z(x) + Z'(x) h + Z''(x)/2 h² + Z'''(x)/6 h³ + ...
Z(x+h) / Z(x) = 1 - x h + (x² - 1)/2 h² - x (x² - 3)/6 h³ + ...

Example: calculate Z(20.3333 333333)

x = z-score rounded to FIX-4, thus exp(-x²/2) can be calculate accurately.
h = z-score - x, thus |h| ≤ 5e-5

-> above example, x=20.3333, h=3.33333e-5
-> Z(x) = exp(-x²/2) / √(2 Pi) = 6.65095 520534 e-91

Correction to get to Z(x+h):

Code:

c1 = -x h =          -677775.98889 e-9
c2 = (c1+h)(c1-h) / 2 = +229.13459 e-9
c3 = c1 (c2 - h²) / 3 =   -0.05152 e-9
c  = c3 + c2 + c1 =  -677546.90582 e-9

Z(x+h) ~ Z(x) + c Z(x) = 6.64644 887122 e-91, matched true value