Accurate Bernoulli numbers on the 41C, or "how close can you get"?

[Dieter]

(03-19-2014 06:47 AM)Ángel Martin Wrote:  Hi Dieter, glad to see you find the SandMath worth looking at. Which revision are you using? Are there two appendices 9a and 9b, or only one appendix 9?

The version I had was 4.44.44, now updated to 5.55.55 that indeed has additional information on the Normal quantile function.

Actually I do not use Sandmath. I heard about it quite often here, and when I tried a first look at MCODE I came across the manual. At the moment this is all what I got. ;-)

(03-19-2014 06:47 AM)Ángel Martin Wrote:  I say this because in the latest versions (Revision "M", see above link) there are a couple of functions directly related to the quantile, one (ICPF) using the inverse error function on a general Normal distribution (s,m), and another (QNTL) using Halley's iterative method on a standard Normal (0,1).

Halley's method is nice, but it can be done better and faster. ;-)

I now remember that some time ago I set up an algorithm specially taylored for the HP41 and other 10-digit machines. It first calculates an initial estimate with a very simple (1, 2) rational function, requiring just four constants with few digits, and then one single correction step is sufficient for a result with error less than 0,055 ULP, i.e. 5,5 units in the 12th significant digit. All it needs is a sufficiently accurate CDF (cumulative distribution function).

The following method assumes p < 0,5. The rest is handled by CDF(-x) = 1 - CDF(x).

\(u := \sqrt{-2 ln p}\)

\(x := max( 0 ,  u - \large \frac{2,374 + 0,39 u}{1 + 1,11 u + 0,07157 u^2} )\)

\(t := \large \frac{CDF(x) - p}{PDF(x)}\)

Now apply a third order (!) correction, cf. Abramowitz & Stegun, 10th ed., p. 954:

\(x := x + t + \frac{x}{2}t^2 + \frac{2x^2+1}{6}t^3\)

If evaluated exactly, the result has an error within approx. +0 and -5,5 units in the 12th significant digit throughout the 41's working range (even down to 1E-121).

Here is an implementation in Visual Basic for Excel:

Code:

Function icdf10(p)
  signflag = (p > 0.5)
  If signflag Then q = 1 - p Else q = p
  u = Sqr(-2 * Log(q))
  x = u - (2.374 + 0.39 * u) / ((0.07157 * u + 1.11) * u + 1)
  If x < 0 Then x = 0
  t = (normcdf(x) - q) / normpdf(x)
  x = t * t * t * (2 * x * x + 1) / 6 + t * t * x / 2 + t + x
  If signflag Then x = -x
  icdf10 = x
End Function

Here normcdf() and normpdf() refer to the Normal CDF resp. PDF. Please note that the PDF refers to the upper integral, i.e. from x to infinity. Otherwise two small adjustments are required.

The essential magic is in the combination of a dedicated initial guess and a very effective correction afterwards. If the term at t³ is omitted, the result quite exactly matches that of Halley's method and the error is less than 2 units in the 9th digit.

All this should work for p = 0,2 as well as p = 1E-12 or 1E-99 (that's why \(erf^{-1}(2p-1)\) is not a good idea). Cases close to the distribution center and and the lower working limit (1E-99) require some care at the calculation of CDF(x) - p since this difference may become numerically zero due to underflow or digit cancellation. But there are ways to handle this.

Dieter

Edit: corrected an error and tweaked one of the constants for even better accuracy.