HP Forum Archive 17

[ Return to Index | Top of Index ]

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 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 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.

Re: a/b != a*(1/b) on HP41CV, 33C, 15C, 67.... Message #2 Posted by Les Wright on 11 Apr 2007, 8:03 a.m., in response to message #1 by Les Wright Before someone points out the "why" I think I have an idea what is going on. The "reciprocal then multiply" approach is a two step process. The intermediate step of taking the reciprocal of 47.0001 basically rounds that result to 2.127655048e-2, since every time a result is placed on the the stack the three internal guard digits are lost. Multiply that by 470, round to 10 digits, place on stack--voila! Inherited rounding error. 9.999978726 is not the properly rounded 10 digit result of 470/47.0001, but it is the properly round result of 470 * 2.127655048e-2. The division is a single step with all of the intermediate computations carried on internally to 13 digits before the final 10 digit rounded returned to the stack. So it gets the result right. The lesson seems to be that if maximum possible accuracy is desired in division, I should just do a division, since adding extra steps just eliminates guard digits that could prove very important in certain situations. Les

[ Return to Index | Top of Index ]

Go back to the main exhibit hall