Accurate Normal Distribution for the HP67/97

[Albert Chan]

(12-12-2018 10:49 PM)John Keith Wrote:  My implementation is straightforward except for the second line,
double x = (float) z;
which turned into -6 RND i.e. round to 6 digits.

Round to 6 digits does not work. Example: 4.56789²/2 = 10.4328 095260 5
Try round z to FIX-4, we get:

|y| < 21.4 * 5e-5 = 0.00107, y⁴/24 < 5.5e-14

3rd order correction also enough for C double:
From header file float.h: DBL_EPSILON = 2^-52 ~ 2.22e-16

|y| < 38.6 * 2^-19 ~ 0.000074, y⁴/24 < 1.2e-18

Quote:I did not experience any underflow, just a less accurate result.

An example of underflow, try Z(20.5000 00002):

z = Z(20.5) = 2.21198 438021 e-92
y = 0.5*(20.5 - 20.5000 00002)*(20.5 + 20.5000 00002) ~ -4.1e-8
Other correction terms y²/2, y³/6 are too small, and can be ignored.

-> correction = y z ~ -9.069e-100, underflow to 0

That was why I mentioned do z*(1. + y) for calculator program.

-> Z(20.5000 00002) = Z(20.5) * 0.99999 9959 = 2.21198 428952 e-92, matched true value