a/b != a*(1/b) on HP41CV, 33C, 15C, 67.... Message #1 Posted by Les Wright on 11 Apr 2007, 7:42 a.m.
I am sure this is in the area of Numeric Programming Skills 101, but I thought I would share an observation.
I am working on an HP41 program that, among other things, computes the uppertail probability of the Student t distribution given degrees of freedom and t statistic value. (The t distribution is a specific case of the incomplete beta function, and the program has code to compute that. I like JM Baillard's programs in the Software Library, but I thought I would write my own routine using the Modified Lentz procedure in Numerical Recipes that didn't have dependency on other routines save a reliable way to compute the Gamma function.) In certain cases, such as when degrees of freedom are relatively large and the other argument is very small, I may need to compute a quotient that turns out to be pretty close to unity. For such arguments I was getting more digit loss in the final result than I expected, and found that the situation improved greatly when I computed such quotients as divisions a/b rather than multiplications (1/b)*a. I thought the calculator would treat these situations as identical, but I was wrong.
For example, on all of my 10digit HPs I get the following:
47 ENTER 47.0001 / 10 * gives 9.999978723
BUT
47.0001 1/x 47 * 10 * gives 9.999978726
In Mathematica, the quotient to 20 digits is 9.9999787234495245755, so the straight division is the "more correct" result here.
It turns out that this difference of 3 ULP was having a profound effect on later computations such that 2 or even three digits were lost. When I programmed the division as, well, a division, the results were much better.
I am sure this has been discussed before, but I share it in case it has not come up in a while. It is a common step saving device in RPN programs to compute quotients as (1/b)*a to avoid excessive stack lift and register swapping. I have used it a lot myself. But I will be more careful now, especially when I am interested in preserving as much digit accuracy as possible.
Les
P.S. The analogous discrepancy on my 12 digit calcs is similar but not as extreme9.99997872345 with normal division, 9.99997872344 with the reciprocal and multiply technique. Normal division is again the correct approximation to the actual result.
