|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 upper-tail 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 J-M 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 10-digit HPs I get the following:
47 ENTER 47.0001 / 10 * gives 9.999978723
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.
P.S. The analogous discrepancy on my 12 digit calcs is similar but not as extreme--9.99997872345 with normal division, 9.99997872344 with the reciprocal and multiply technique. Normal division is again the correct approximation to the actual result.