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

03192014, 01:37 PM
(This post was last modified: 03192014 06:54 PM by Dieter.)
Post: #18




RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03192014 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. ;) (03192014 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 10digit 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 1E121). Here is an implementation in Visual Basic for Excel: Code: Function icdf10(p) 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 = 1E12 or 1E99 (that's why \(erf^{1}(2p1)\) is not a good idea). Cases close to the distribution center and and the lower working limit (1E99) 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. 

« Next Oldest  Next Newest »

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