Accurate Normal Distribution for the HP67/97

[Albert Chan]

Discovered a simpler and faster 1 Exp method.

Instead of correcting Z(x) to Z(x+h) via Taylor expansion of Z(x+h), do this:

Z(x+h) / Z(x) = exp(-½(x+h)² + ½(x²)) = exp(-x h - h²/2) = exp(y)

exp(y) = 1 + y + y²/2 + y³/6 + ...

|y| ~ |x h|, a very small number, above ... can be ignored.

PHP Code:

double pdf(double z)             // 1 Exp method, revised
{
    double x = (float) z;        // clear low bits
    double y = 0.5*(x-z)*(x+z);  // PDF(z) = exp(y) PDF(x)
    double y2 = y*y;
    z = 0.3989422804014327 * exp(-0.5*x*x);
    y += (y*(1./6) + 0.5) * y2;  // ≈ exp(y) - 1
    return z + y * z;
} 


Note: for calculator program, return z * (1. + y) instead
Calculator might underflow the correction term.