What should you get?



#24

There has been some recent discussion of the HP33S's errors when computing the TAN function.

The calculator forensics page has a result for the HP33S:
http://www.rskey.org/~mwsebastian/miscprj/forensics.htm

The philosophy behind the math routines in the HP calculators after the influence of Prof. Kahan began to be felt is expressed like this:

"Each calculation is to return a result as though the calculation were done to infinite precision and the result properly rounded (round to even on most Saturn calcs) to an n-digit result."

The HP15C Advanced Functions Handbook describes which functions of that calculator met this goal. Some of the mathematical functions such as the transcendentals and trigonometrics don't always meet the goal.

An interesting question for any calculator is, "How closely does the calculator's math approach this goal"?

The so-called Forensics Algorithm at the site mentioned above is this calculation:

arcsin(arccos(arctan(tan(cos(sin(9))))))

Now for a mini challenge. What result should you expect for an n-digit calculator, 8 <= n <= 16, if each calculation is done according to the goal described above?

Hint: the answer is not 9.

Look through the results at the Forensics web site, and see if any calculators meet the goal. Which ones are they?


#25

Hi, Rodger;

I never took a closer look at the forensics results for

arcsin(arccos(arctan(tan(cos(sin(9))))))
before, and after that 'search', I am confused... Indeed, the many 'out-of-range' results may be espected, depending on how algorithms were writen and how do they deal with exceptions, so the few 'error' and '0.0' may have their excuse. Some interesting results in the range '10.xxx' to '15.xxx' are worth knowing about. The '71.252182'
found by Katie Wasserman in a 'National Semiconductor Scientist' is actually intriguing, mainly when we find the exact same result in both Novus 4520 Scientist, Novus 4525 Scientist PR and Privileg PR-55NC. Same algorithm?

One particular result called my attention: 9.999999999833 found in one Psion Diamond Mako/Revo Plus. The other three models, namely Psion Series 3a, Psion Series 5 (default application) and Psion Series 5 dCalc V. 2.02 (PocketIQ), are all in the range 8.99999999983xxx. any chance that the we have a typo in the Psion Diamond Mako/Revo Plus result?

Interesting collection of forensics. I guess that this sort of comparison table might be a good reference for those dealing with precision computation. About the 'new breed' HP models with a 9.0000000... result. One of them - HP30S - had already been thge subject of an earlier thread (some time ago, IIRC was Karl the one who pointed that out) where its 32-bit internal precision was the issue, so I was not so taken by surprise in this case. But the HP9G... I didn´t see that comming.

Thanks, Rodger. And thaks to all other contributors who share their particular findings with us all here.

Cheers.

Luiz (Brazil)


#26

Hi, Luiz --

Quote:
arcsin(arccos(arctan(tan(cos(sin(9))))))

About the 'new breed' HP models with a 9.0000000... result. One of them - HP30S - had already been thge subject of an earlier thread (some time ago, IIRC was Karl the one who pointed that out) where its 32-bit internal precision was the issue, so I was not so taken by surprise in this case.


(NOTE: Corrected and updated with information from Rodger Rosembaum's post in reply.)

Several threads are found in Archive #15, including Rodger's original post on this topic.

http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv015.cgi?read=85973#85973

http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv015.cgi?read=88282#88282

http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv015.cgi?read=83934#83934

The "forensics" expression, while straightforward, always seemed a bit arbitrary to me.

The HP-30S, with only two exponent digits, apparently uses an 80-bit extended-double floating-point representation. It allows extra precision -- up to 24 digits -- for the mantissa. The Saturn-processor 64-bit binary-coded decimal (BCD) representation with three exponent digits allows only 12 digits for the mantissa.

Here are several useful links:

http://en.wikipedia.org/wiki/IEEE_Floating_Point_Standard

http://en.wikipedia.org/wiki/Binary_coded_decimal


Part of the reason that the HP-30S returns the exact answer of 9 for the forensics test is that the HP-30S takes the liberty of rounding results to integers that are in very close proximity, which provides answers that are reassuring to novice users. This was noted before by others, but here's my favorite example:

On a Saturn-processor calculator (or their KinHPo descendants), take "sin 3.14159265358" in radians mode. You will get the numerically-correct result of 9.79323846264 x 10-12 -- which not-coincidentally are the next 12 significant digits of pi, given that the input was not, in fact, exactly pi.

Now, try the same on the HP-30S. Start with "sin 3.1415926535". The displayed result is 8.979323846 x 10-11 -- correct to its 10-digit display. Then try "sin 3.14159265358". The answer returned is exactly zero! What happened? Rounding, for the sake of reassurance -- "This answer must actually be zero, so let's return that result to gratify the user."

"EXACTITUDE OF REPRESENTATION" -- BCD VS. FLOATING POINT

A nice thing about BCD is its exactitude and direct correspondence to the encoded value. This makes it well-suited for calculators, which usually display a result immediately after each calculation. In programs, a "test for zero" can't be bollixed by a tiny floating-point argument that wasn't quite represented as zero. In HP's 64-bit BCD, the mantissa of a zero is simply a string of 48 cleared bits.

Retrieve pi on a BCD calculator, then subtract the integer part and multiply the result by 10. Keep repeating, and you will see no additional digits: Pi to 12 digits is given as 3.14159265359, and nothing more.

Now, try the HP-30S. Retrieve pi, and keep revealing the extra digits as above. You'll obtain

3.14159 26535 89793 23846 264 8727...

The first 24 digits are correct, but the next four that follow should have been "3383". The limits of resolution have been exceeded, but the HP-30S will keep giving more incorrect digits if the procedure is continued.

HP-30S -- REPRESENTATION VS. COMPUTATION

The HP-30S' representational accuracy of 24 significant digits does not extend to the same computational accuracy.

The digit-revealing procedure after computing e^(1) on the HP-30S produces

2.7182 8182 8459 0452 018...

which strays from the correct value in the 18th significant digit:

2.7182 8182 8459 0452 353...

But, of course, there's a bottomless well of incorrect digits to be retrieved.


Going back to the forensics example, one can see on the HP-30S that

  • arcsin(arccos(arctan(tan(cos(sin(8.999999999)))))) = 8.9999999990528...

  • arcsin(arccos(arctan(tan(cos(sin(9.000000001)))))) = 9.0000000010528...

    with error appearing in the 12th significant digit.


    In summary, if there's a sound argument for using floating-point representation instead of BCD on a calculator (and particularly the slow, non-programmable, non-graphing HP-30S), I'd like to see it.

    -- KS


    Edited: 2 Dec 2006, 9:40 p.m. after one or more responses were posted


  • #27

    Quote:
    I can't seem to dig up the thread, but I do vaguely remember it. The "forensics" expression, while straightforward, always seemed a bit arbitrary to me (Why 9? Why not 7, or 13? Why not reverse the positions of SIN and TAN?)

    I think the thread is:
    http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv015.cgi?read=85973#85973

    On the Forensics site, the author explains "Why 9?". It was because that was the largest integer he could enter with a single key press.

    As for "Why not reverse the positions of SIN and TAN?", I guess he just liked to press the keys in a left to right and then back again order (or perhaps the reverse on some calcs).

    Quote:
    I believe that the HP-30S, with only two exponent digits, uses a 64-bit floating-point representation. It is not the IEEE standard, but it does allow extra precision -- up to 24 digits -- for the mantissa.

    It was established shortly after the HP-30S came out that it uses an 80 bit internal (mantissa) arithmetic. This gives 80 / (LOG(10)/LOG(2)) ~= 24 digits decimal.

    Quote:
    Now, try the HP-30S. Retrieve pi, and keep revealing the extra digits as above. You'll obtain 3.14159 26535 89793 23846 264 87271...

    The explanation for this is the following:

    If you convert PI to its binary representation, round to 77 bits, and then convert back to decimal you get
    3.14159 26535 89793 23846 26487 27060 97060 62285 40447 09129 08436 73050 40359 49707 03125 00000 00000 00000...

    The HP30S is just giving back the correct decimal digits for a 77 bit internal binary representation of PI.

    However, the strange thing is, if I follow your procedure, I get 3.14159 26535 89793 23846 19869 82570 92818 40886 43320 13763... which corresponds to PI converted to binary, rounded to 70 bits, and converted back to decimal. I suspect that if I kept going, I would continue to get the correct digits for the 70 bit PI, but there are 18 more and I'm tired of pressing buttons!

    I wonder why you got the decimal digits corresponding to a 77 bit approximation to PI, and I got those for a 70 bit approximation. What I did to "retrieve pi" was to press the PI button just to the left of the "enter" button. That was my starting number for the process. Did you do something different?

    Quote:
    The digit-revealing procedure after computing e^(1) on the HP-30S produces 2.7182 8182 8459 0452 018...

    If I carry the process further, I get:
    2.7182 8182 8459 0452 0181 7900 7609... which is the correct decimal re-conversion of a properly rounded 55-bit representation of e. Don't ask me why 55 bits instead of 80 bits. Or why the 77 bits and 70 bits for PI above. I've noticed this behavior before when investigating the internals of the HP30S.


    #28

    Rodger --

    Thank for the clarifications and information. I have incorporated corrections and enhancements to my post.

    I'm not sure what accounted for the differences in the low-order digits of pi, but we have the same procedure in mind.

    I'm also not sure about internal rounding in the HP-30S, or how many bits, but I think we both agree that the HP-30S' calculations based on its 80-bit internal representation are not as robust as they ought to be. Frankly, I'm not sure what purpose was served by abandoning the 64-bit BCD and all the Saturn-processor mathematical algorithms developed for it. The HP-30S is no computer -- it's a slow, non-graphing, non-programmable calculator. Whatever benefits an 80-bit floating-point representation would offer were just not applicable here.

    -- KS


    #29

    Do you think Jean-Yves would in fact know the answer to the question "why?"

    #30

    Well, as I said in another thread:

    "Rather than do proper arithmetic, just overwhelm the problem by carrying many extra digits, and throw a lot of them away in the final results.

    This is easier to do nowadays given the state of the art in microelectronics. But in the early days of calculators, carrying 25 digits to get a final result of 12 digits was a (not very elegant) shortcut the designers of the day couldn't afford."


    #31

    This compounds the difficulty of teaching kids these days the concept of significant digits, the use of estimation in quickly carrying out simple mental calculations (what should the tip be on tonight's check?), doubtful digits, error margins, propagation of errors, and such.

    I hope hatches on a spaceship or submarine can keep tight these days... or for that matter, a car door!

    Have any of you seen calculations by college students these days? It is exactly "... overwhelm(ing) the problem by carrying many extra digits, and throw a lot of them away in the final results... " (parenthetical addition mine), and it is upsetting and disgusting and depressing that they neither can nor really want to tell what figures are... well, real, and which are not, and let's not even get into real and complex. They just either copy down all ten or twelve digits spat out by their TIs or Casios or Sharps or blindly round (truncate, really, in the case of many of them) off to give the "final answer". Aaagghh!

    #32

    Hi, Karl;

    thanks for pointing that out... 32-bit internal precision! Where was my head placed at that time... I wonder where did I take that figure from. When writing, I was sure I was refering to a 64-bit precision, that is even wrong when compared to the actual 80-bit.

    Thanks again.

    Very good post this one of yours, as always.

    Best regards.

    Luiz (Brazil)

    #33

    Quote:
    One particular result called my attention: 9.999999999833 found in one Psion Diamond Mako/Revo Plus. The other three models, namely Psion Series 3a, Psion Series 5 (default application) and Psion Series 5 dCalc V. 2.02 (PocketIQ), are all in the range 8.99999999983xxx. any chance that the we have a typo in the Psion Diamond Mako/Revo Plus result?

    This is definitely a typo. All Psions running EPOC (now known as Symbian OS), i. e. the Series 5, 5mx, Revo/Mako, Series7, netBook, ...
    use IEEE arithmetic. The routines for the transcendental functions are built into the operating system, so any calculator software using the operating system services for math will give identical results on all these machines.

    #34

    Here are the results for n=10, 11 and 12:

    					    n=10		     n=11		     n=12

    sin(9) 0.1564344650 0.15643446504 0.156434465040
    cos(sin(9)) 0.9999962727 0.99999627274 0.999996272743
    tan(cos(sin(9))) 0.01745499985 0.017454999855 0.0174549998555
    atan(tan(cos(sin(9)))) 0.9999962724 0.99999627271 0.999996272744
    acos(atan(tan(cos(sin(9))))) 0.1564416604 0.15643515514 0.156434441642
    asin(acos(atan(tan(cos(sin(9)))))) 9.000417403 9.0000400327 8.99999864267

    Quote:
    Look through the results at the Forensics web site, and see if any calculators meet the goal. Which ones are they?


    10-digit calculators:

    HP-19C, HP-27, HP-29C, HP-41C/CV, HP-67, HP-91, HP-97, HP-97S HP-10C, HP-11C, HP-15C, HP-31E, HP-32E, HP-33E, HP-34C and Canon Canola F-11

    12-digit calculators:

    HP-19BII, HP-20S, HP-22S, HP-27S, HP-28C/S, HP-32S, HP-32SII, HP-38G, HP-39G, HP-42S, HP-48S, HP-48G/G+/GX, HP-48GII, HP-49G, HP-49G+, HP-50G, HP-71B, HP-75C/D and HP-87

    Interestingly CASIO FX-501P gives 9.0000400303, which is close to what we should expect on the theoretical 11-digit calculator.

    Regards,

    Gerson


    Edited: 27 Nov 2006, 3:48 p.m.


    #35

    The result for 8 digits is
    8.9911614.

    For 9 digits, it seems to be, barring any mistakes on my part,
    8.99725190

    Several calculators have the result for 8 digits; none seem to have the result for 9 digits.


    #36

    Nice work Gerson and Crawl! I didn't get a chance to finish my table of results before you had it all done. I was using an HP 30S to figure this out since you need to have a calculator that has greater precision than the required result, what did you use, not a PC I hope! :)


    #37

    Quote:
    I was using an HP 30S to figure this out since you need to have a calculator that has greater precision than the required result, what did you use, not a PC I hope! :)

    Shame on me! :-(

    I was at work when I prepared the tables so I used the Windows calculator. I might as well have used Free42 as it was also available:

    0.15643446504
    9.99996272743e-1
    1.74549998555e-2
    9.99996272744e-1
    1.56434441642e-1
    8.99999864267

    Just a matter of copying and pasting to and from Notepad.

    Gerson.

    #38

    Quote:
    ...since you need to have a calculator that has greater precision than the required result

    I don't think there is a problem with the current example, but this is not sufficient in general :-)

    A contrived counter-example:

        1.49      true value
    1.5 rounded to 2 digits

    1. true value rounded to 1 digit
    2. 2 digit rounded value rounded to 1 digit

    It isn't too difficult to detect when this might be occuring.


    - Pauli

    #39

    My purpose here was to show that if a BCD calculator does the "best" than can be done with n-digit BCD arithmetic, there is a specific result that you should get for a particular calculation, and that result is not necessarily what "exact" arithmetic would get.

    If the calculator doesn't get that "best" result, then the calculator is committing "avoidable error".

    A very simple example: Type "6" into a calculator and press 1/x twice. Should you get exactly "6" back? No. You should get 5.99999999... to n digits. If you see exactly 6, then the calculator has some "hidden" digits (such as TI calculators, or any calculator with the display set to show fewer digits than are present in final results), OR the calculator is looking for special cases and perturbing the answer to "please" the user with the expected exact answer (as some Casios do, and the HP30 sometimes does), which is not a good thing.

    The results of my mini challenge for n = 13 to 16 are:

       n         Result
    13 8.999999860032
    14 8.9999999817692
    15 8.99999999759468
    16 9.000000000029361

    Gerson and Crawl have the correct results for the cases they worked out.

    For the n=16 result, only the HP100LX and HP200LX solver get the right result. The HP95LX solver is very close; probably made a very tiny error somewhere along the way.

    The result for n=14 ought to apply to a lot of the TI machines, but only the Canon F-710, TI-89, TI-92 and TI Voyager 200 get it right.

    Note that the result for n=13 is suspiciously close to what the HP33S gets. That makes me wonder if it is doing 13 digit arithmetic internally.

    There is a subtlety involved in this "Forensic" calculation. Most calculators compute trig functions in radians, and multiply the input argument (for SIN, COS, and TAN) by PI/180 for an argument in degrees (and postmultiply the result by 180/PI for ARCSIN, ARCCOS and ARCTAN). If the multiplication by PI/180 is done with more than n digits internally (in a calculator that displays a maximum of n digits), the result may be different than if the multiplication by PI/180 is done with only n digits.

    For example, the HP41, HP15, etc., are known to do internal arithmetic with 13 digits although final results are only 10 digits. That is, a 13 digit approximation to PI/180 is used to convert the argument from degrees to radians.

    If only 10 digits were used for internal calculations, then the Forensic calculations would go like this:

    sin(9)                                  0.1564344651
    cos(sin(9)) 0.9999962727
    tan(cos(sin(9))) 0.01745499986
    atan(tan(cos(sin(9)))) 0.9999962729
    acos(atan(tan(cos(sin(9))))) 0.1564311679
    asin(acos(atan(tan(cos(sin(9)))))) 8.999808733

    In the list of Forensic results, only the Casio CFX watches appear to be doing this (maybe!) for 9 digit results, with (probably) 10 digit internal calculations.

    For some values of n, a different result is obtained if the PI/180 factor is used internally to only the same number of decimal places that the calculator uses to display its results. Those values are:

       n         Result
    8 8.9789683
    10 8.999808733
    13 9.000000712181
    14 9.0000000548100
    15 8.99999999759462

    There are a number of results that match the n = 8 case, although off by 1 count in the LSD.

    It's interesting to do the Forensic calculation one step at a time on a TI machine such as the TI-86 and see where it makes its mistake. It's usually only one count in the LSD, but that's all it takes to get a substantially different final result.

    So, the result to be desired when this "Forensic" test is done by any calculator is not a result of exactly 9. If this happens, the calculator is not doing correct BCD floating point math (I'm not including systems that can do exact arithmetic in this discussion). No n-digit calculator, 8 <= n <= 16, gets exactly 9 if it's doing "correct" BCD floating point arithmetic. The result to be desired is the result from this mini challenge for an n-digit calculator (within an LSD or so).

    Isn't it interesting that the 10 and 12 digit HP machines since about the time of the HP29C get the correct result? Excepting those not designed in Corvallis, I think.

    (I posted something a while back about the HP30 trickery. Do the Forensic calculation starting with 9 and you get exactly 9 back. Do it starting with 9.1 and you can see the errors in the final result.)


    #40

    Quote:
    If the calculator doesn't get that "best" result, then the calculator is committing "avoidable error".

    I had already noticed good forensic results did not always imply the best algorithms had been used. A slight error in the next to the last result can lead to a better final result, for instance. While I don't agree with tricks to get some "user-pleasing" answers, I appreciate the use of internal digits to avoid cumulative errors. Early calculators could not afford more than a couple of them, though.

    In my program for the HP-12C Platinum, which uses polynomial approximations, I had to decide between two sets of coefficients. They gave these results:

                                                set 1             set 2

    sin(9) 0.156434465041 0.156434465041
    cos(sin(9)) 0.999996272744 0.999996272742
    tan(cos(sin(9))) 0.0174549998556 0.0174549998555
    atan(tan(cos(sin(9)))) 0.999996272743 0.999996272738
    acos(atan(tan(cos(sin(9))))) 0.156434461500 0.156434566400
    asin(acos(atan(tan(cos(sin(9)))))) 8.99999979459 9.00000587986

    I chose the second set because the constants needed less digits. Since it's a ten-digit calculator with two extra internal digits, I thought this wouldn't affect the displayed results.

    Interestingly those results are close to those of Elektronika C3-15. Early version gives 8.99999979443 and later version gives 9.00000588129. Isn't this an indication they were using the same method?

    Gerson.

    Edited: 28 Nov 2006, 6:50 p.m.


    #41

    Hi Gerson,

    This result looks strange to me:

    acos(0.999996272743)=0.156434461500

    To 12 sig digits the answer has 2627 instead of the terminal 1500 above. The 1500 is therefore 1127 too light - for the second set of coefficients the 6400 is 1153 too light. On the face of it it looks like your acos(x) for x near 1 truncates? I guess using acos(x)=90-asin(x) might explain it - yup the difference between 2 numbers up near 90 will certainly create the "00" for the 11th and 12th d.p. of the diff.

    Cheers,
    Tony


    Edited: 28 Nov 2006, 8:19 p.m.


    #42

    You're right!

    Quote:
    I guess using acos(x)=90-asin(x) might explain it.

    But only because I used asin(x) = atan(x/(sqrt(1 - x2)).

    I've been able to avoid this problem in cos(x) by using cos(x) = sin(pi/2 - x) instead of cos(x)=sqrt(1 - sin2(x)). I ought to have used a better expression for asin(x). I remember for a while I became aware of this problem, but then I simply forgot about it.

    Thanks for pointing this out. I will have to change the program or at least the title: Fast and Accurate... except for acos(x) when x > 0.99999

    :-)

    Cheers,

    Gerson.


    Edited: 28 Nov 2006, 9:06 p.m.


    #43

    Hi Gerson,

    Quote:
    But only because I used asin(x) = atan(x/(sqrt(1 - x^2)).

    That won't get around the "big" 90 which pushes away a couple of digits. This might help for x near 1:

    acos(x) = atan((sqrt(1 - x^2)/x)

    "Tricky Trigonometrics" indeed :-)

    Tony


    #44

    Thanks, Tony!

    Too bad this might make an already large program even larger. I think my older version didn't have this bug, but the range for acos(x) was limited to positive x. I'll check this later.

    Regards,

    Gerson.

    Edited: 30 Nov 2006, 8:47 a.m.

    #45

    That's the formula that I used in my BASIC astronomy programs years ago where the only inverse function that was available was ATN.

    tm


    #46

    My first computer had only SIN, COS and ATAN. I remember the manual presented a list with BASIC expressions for all derived functions.

    Gerson.

    ---------------------

    Also TAN was available. Found an emulator:

    http://codigolivre.org.br/frs/?group_id=1367&release_id=1598

    5  HOME
    10 K = ATN (1) / 45
    20 F = ATN ( TAN ( COS ( SIN (9 * K) * K) * K)) / K
    30 X = F: GOSUB 200:F = X
    40 X = F: GOSUB 100:F = X
    50 PRINT F
    60 END
    100 X = ATN (X / SQR (1 - X * X)) / K
    110 RETURN
    200 X = ATN ( SQR (1 - X * X) / X) / K
    210 RETURN

    9.00015733


    Edited: 30 Nov 2006, 7:35 p.m.


    Forum Jump: