▼
Posts: 305
Threads: 17
Joined: Jun 2007
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 ndigit 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 socalled 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 ndigit 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?
▼
Posts: 4,027
Threads: 172
Joined: Aug 2005
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 'outofrange' 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 PR55NC. 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 32bit 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)
▼
Posts: 1,792
Threads: 62
Joined: Jan 2005
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 32bit 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/cgisys/cgiwrap/hpmuseum/archv015.cgi?read=85973#85973
http://www.hpmuseum.org/cgisys/cgiwrap/hpmuseum/archv015.cgi?read=88282#88282
http://www.hpmuseum.org/cgisys/cgiwrap/hpmuseum/archv015.cgi?read=83934#83934
The "forensics" expression, while straightforward, always seemed a bit arbitrary to me.
The HP30S, with only two exponent digits, apparently uses an 80bit extendeddouble floatingpoint representation. It allows extra precision  up to 24 digits  for the mantissa. The Saturnprocessor 64bit binarycoded 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 HP30S returns the exact answer of 9 for the forensics test is that the HP30S 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 Saturnprocessor calculator (or their KinHPo descendants), take "sin 3.14159265358" in radians mode. You will get the numericallycorrect result of 9.79323846264 x 10^{12}  which notcoincidentally are the next 12 significant digits of pi, given that the input was not, in fact, exactly pi.
Now, try the same on the HP30S. Start with "sin 3.1415926535". The displayed result is 8.979323846 x 10^{11}  correct to its 10digit 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 wellsuited for calculators, which usually display a result immediately after each calculation. In programs, a "test for zero" can't be bollixed by a tiny floatingpoint argument that wasn't quite represented as zero. In HP's 64bit 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 HP30S. 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 HP30S will keep giving more incorrect digits if the procedure is continued.
HP30S  REPRESENTATION VS. COMPUTATION
The HP30S' representational accuracy of 24 significant digits does not extend to the same computational accuracy.
The digitrevealing procedure after computing e^(1) on the HP30S 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 HP30S 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 floatingpoint representation instead of BCD on a calculator (and particularly the slow, nonprogrammable, nongraphing HP30S), I'd like to see it.
 KS
Edited: 2 Dec 2006, 9:40 p.m. after one or more responses were posted
▼
Posts: 305
Threads: 17
Joined: Jun 2007
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/cgisys/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 HP30S, with only two exponent digits, uses a 64bit floatingpoint 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 HP30S 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 HP30S. 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 digitrevealing procedure after computing e^(1) on the HP30S 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 reconversion of a properly rounded 55bit 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.
▼
Posts: 1,792
Threads: 62
Joined: Jan 2005
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 loworder digits of pi, but we have the same procedure in mind.
I'm also not sure about internal rounding in the HP30S, or how many bits, but I think we both agree that the HP30S' calculations based on its 80bit internal representation are not as robust as they ought to be. Frankly, I'm not sure what purpose was served by abandoning the 64bit BCD and all the Saturnprocessor mathematical algorithms developed for it. The HP30S is no computer  it's a slow, nongraphing, nonprogrammable calculator. Whatever benefits an 80bit floatingpoint representation would offer were just not applicable here.
 KS
▼
Posts: 2,448
Threads: 90
Joined: Jul 2005
Do you think JeanYves would in fact know the answer to the question "why?"
Posts: 305
Threads: 17
Joined: Jun 2007
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."
▼
Posts: 801
Threads: 36
Joined: Nov 2005
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!
Posts: 4,027
Threads: 172
Joined: Aug 2005
Hi, Karl;
thanks for pointing that out... 32bit 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 64bit precision, that is even wrong when compared to the actual 80bit.
Thanks again.
Very good post this one of yours, as always.
Best regards.
Luiz (Brazil)
Posts: 3,283
Threads: 104
Joined: Jul 2005
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.
Posts: 2,761
Threads: 100
Joined: Jul 2005
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?
10digit calculators:
HP19C, HP27, HP29C, HP41C/CV, HP67, HP91, HP97, HP97S HP10C, HP11C, HP15C, HP31E, HP32E, HP33E, HP34C and Canon Canola F11
12digit calculators:
HP19BII, HP20S, HP22S, HP27S, HP28C/S, HP32S, HP32SII, HP38G, HP39G, HP42S, HP48S, HP48G/G+/GX, HP48GII, HP49G, HP49G+, HP50G, HP71B, HP75C/D and HP87
Interestingly CASIO FX501P gives 9.0000400303, which is close to what we should expect on the theoretical 11digit calculator.
Regards,
Gerson
Edited: 27 Nov 2006, 3:48 p.m.
▼
Posts: 306
Threads: 3
Joined: Sep 2009
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.
▼
Posts: 1,477
Threads: 71
Joined: Jan 2005
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! :)
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.99996272743e1
1.74549998555e2
9.99996272744e1
1.56434441642e1
8.99999864267
Just a matter of copying and pasting to and from Notepad.
Gerson.
Posts: 3,229
Threads: 42
Joined: Jul 2006
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 counterexample:
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
Posts: 305
Threads: 17
Joined: Jun 2007
My purpose here was to show that if a BCD calculator does the "best" than can be done with ndigit 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 F710, TI89, TI92 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 TI86 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 ndigit 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 ndigit 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.)
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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 "userpleasing" 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 HP12C 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 tendigit calculator with two extra internal digits, I thought this wouldn't affect the displayed results.
Interestingly those results are close to those of Elektronika C315. 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.
▼
Posts: 93
Threads: 2
Joined: Jul 2005
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)=90asin(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.
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
You're right!
Quote:
I guess using acos(x)=90asin(x) might explain it.
But only because I used asin(x) = atan(x/(sqrt(1  x^{2})).
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  sin^{2}(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.
▼
Posts: 93
Threads: 2
Joined: Jul 2005
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
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 406
Threads: 47
Joined: Jul 2005
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
▼
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
