Accurate Normal Distribution for the HP67/97

[Massimo Gnerucci]

(12-15-2018 12:06 PM)pier4r Wrote:  It is a bit like paper dictionaries. You search a word and you may stumble in others while you search it. Or you can do a bit of random looking.
Of course digital dictionaries can simulate this, but not as good as paper ones yet.

The same with the mathematical tables. Sometimes curiosity arises from what one experiences.

Serendipity.

[pier4r]

(12-15-2018 08:22 PM)Massimo Gnerucci Wrote:  Serendipity.

Almost.

Here (see n1) the first meaning given is An unsought, unintended, and/or unexpected, but fortunate, discovery and/or learning experience that happens by accident..
What I mean is something different.

Example: you know that something in the "working tools" section of the department store will catch your attention, you just don't know which tool. You go there and look around until you find something (if something is there). You don't go looking in the sport section though.

Or: you know that if you go on ebay searching for calculators you may find something, but you don't know precisely if it is a 12c, a 50g, a 48g, a 35, a 35s, a 42s a dm42, a casio (!), a ti (!!), and so on. You search for calculators, not for a random object.

So it is not really an unintended search. The same happens when I look for words in a dictionary or numbers in numerical tables. I am searching for something of that genre, only I cannot define it until I see it.

n1: https://en.wiktionary.org/wiki/serendipity (wiki wiki wiki!)

[John Keith]

(12-15-2018 06:38 PM)Dieter Wrote:  Then you have found your best way to round the input. And indeed rounding to 6 significant digits may be the best choice here. The square of such a number is exact, and the division by 2 will still be exact if the mantissa of x² does not exceed 2. This is true if the mantissa of x does not exceed √20, i.e. up to 4,47213. Which accounts for 44,72% of all possible x. The remaining 55,28% are also exact if the final digit of x² is even. Otherwise the value of x²/2 may be 1/2 ULP high. So an error only occurs in only 27,64% of all cases, and even if this happens the result if only half an ULP off. I have not analyzed what this means for the final Z(z) result, but a slight error in the last digit may always remain in such calculations.

After all there is no guarantee that the e^x function is accurate to half an ULP. Take a look at the 15C Advanced Functions Handbook where the possible errors of different functions are explained in detail. Here an error of slightly more than 1 ULP (but less than 3 ULP) is possible.

Please correct me if I'm wrong here. ;-)

Dieter

I certainly did not make that detailed of an analysis but it seemed to me intuitively that Albert's method double x = (float) z discards exactly half of the significant bits. Similarly -6 RND discards half of the significant (BCD) digits in a 12-digit calculator.

Nonetheless, the 1-Exp method implemented as above had errors as large as 34 ULPs with some inputs, whereas the 2-Exp method never had more than 1 ULP of error in all of the inputs I tried. I certainly may have made some errors in my implementation but the 2-Exp method works better for me even though the program is larger and slower.

I also noticed that Albert made a couple of new posts in the HP 50 thread which I will have to check out.

[Dieter]

(12-16-2018 07:53 PM)John Keith Wrote:  I certainly did not make that detailed of an analysis but it seemed to me intuitively that Albert's method double x = (float) z discards exactly half of the significant bits. Similarly -6 RND discards half of the significant (BCD) digits in a 12-digit calculator.

IEEE 754 single precision numbers are stored with a 23 bit mantissa, compared to 52 bits in double precision. So actually more than half of the mantissa bits are discarded. That's about 7 vs. 15-16 decimal digits.

Rounding to a certain number of significant digits (RND with negative argument) is the more flexible method. So stay with –6 RND, this seems to produce results close to the theoretical optimum:

(12-16-2018 07:53 PM)John Keith Wrote:  Nonetheless, the 1-Exp method implemented as above had errors as large as 34 ULPs with some inputs, whereas the 2-Exp method never had more than 1 ULP of error in all of the inputs I tried.

That's as good as it gets.

Dieter

[Albert Chan]

Hi, John Keith

Regarding how to split z for 1 Exp Method Revision 1, note the reason for the split:

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

The split have to ensure *BOTH* good Z(x) and exp(y).
If |z| is limited to below √2000 ~ 44.7, FIX-4 work well.

Above this z range, x = -5 RND of z, (to evaluate Z(x) correctly).
|y| is now too big, exp(y) need more terms to get "good enough".

Instead of updating revision 1 with above, I choose an easier way.
1 Exp Method Revision 2 only need a linear correction, and work for even bigger |z|

Example: Z(99.1234567), assuming calculator can handle bigger exponent

B = z²/2 = 4912.72 983408
D = exp(-B) / √(2 Pi) = 1.07016 820397 e-2134
x = z to 5 digits = 99.123
h = z - x = 0.0004567
y = B - x²/2 - x h - h²/2 = 161.2555 e-11

Z(z) = D + y D = D + 173 ULP = 1.07016 820570 e-2134, error = +1 ULP