a/b != a*(1/b) on HP41CV, 33C, 15C, 67....


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.


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.


Possibly Related Threads…
Thread Author Replies Views Last Post
  hp41cv memory modules tim peterson 10 3,658 11-07-2013, 11:09 PM
Last Post: Garth Wilson
  HP41CV Mike Powell 12 3,506 10-30-2013, 08:17 AM
Last Post: Mike Powell
  HP41CV Bry Sun 11 3,304 10-22-2013, 12:52 PM
Last Post: Geoff Quickfall
  Battery holder for my 33C Palmer O. Hanson, Jr. 1 1,095 07-08-2013, 04:03 PM
Last Post: Gerson W. Barbosa
  HP41cv expiriences of a new user Wolfgang 28 7,123 02-13-2013, 08:02 PM
Last Post: Mike Morrow
  A simple question about the HP41C or HP41CV Antoine M. Couëtte 6 2,389 12-16-2012, 04:06 AM
Last Post: Antoine M. Couëtte
  Problem with HP41CV Bruce Larrabee 19 4,978 09-18-2012, 09:29 AM
Last Post: Luiz C. Vieira (Brazil)
  HP41CV drains battery. Schematics? Harald 15 4,265 03-18-2012, 09:22 PM
Last Post: Luiz C. Vieira (Brazil)
  HP41CV batteries John Abbott (S. Africa) 3 1,497 03-02-2012, 06:15 AM
Last Post: John Abbott (S. Africa)
  HP41CV faulty Steve Hunt 6 2,010 12-27-2011, 12:45 AM
Last Post: Kerem Kapkin (Silicon Valley, CA)

Forum Jump: