Posts: 2,761
Threads: 100
Joined: Jul 2005
Some days ago, while I was concerned about the HP33S trig bug I was warned of a bug in my own program (thanks, Tony!). In the next to last line below the program returned 0.156434566400, that is, a truncation in the last two significant digits. As Tony observed, this was because arccos(x) was computed as 90  arcsin(x). The fixed version uses the formula he has suggested, incidentally the same I had used used in an earlier version. However, the program is now 21 steps longer.
The bug was not so critical, the worst case being arccos(9.999999999E01) returning 8.102847000E04 rather than 8.102846845E04. I hope this was the only and the last one :)
sin(9) 0.156434465041
cos(sin(9)) 0.999996272742
tan(cos(sin(9))) 0.0174549998555
atan(tan(cos(sin(9)))) 0.999996272738
acos(atan(tan(cos(sin(9))))) 0.156434566425
asin(acos(atan(tan(cos(sin(9)))))) 9.00000588135
This result is pretty close to that of Elektronika C315, later version: 9.00000588129.
Gerson.
Edited: 3 Dec 2006, 12:48 a.m.
Posts: 93
Threads: 2
Joined: Jul 2005
Hi Gerson,
Congrats of your Magnum Opus!
My apologies for the extra lines :(
I was just idly wondering if you subtracted E13 from the pi/180 (which makes it very marginally more accurate) ... how would the forensic 9 degree test go? Maybe you'll get the 8.99999864267 like the 19BII (for example).
Cheers,
Tony
Posts: 2,761
Threads: 100
Joined: Jul 2005
Hello Tony,
Thanks for the compliment, but Magnum Opus is a bit exaggerated :)
For pi/180 I have used 0.01745329252 because it's close to the actual 12digit constant 0.0174532925199. I tried 0.0174532925198 as you suggested but I noticed no improvement (I obtained 9.00000588157).
Although the last two digits are not truncated anymore, there is a lost of accuracy beginning in the seventh or the eighth significant digit for arguments equal or greater than 0.999999. This appears to be due to floating point errors and doesn't have anything to do with the choice of contants. I have to repeat the same tests in the table below with the older version. In case there isn't any noticeable improvement, perhaps I'll get back to the older version because it's slightly shorter. Anyway, this was a nice exercise for a Saturday afternoon and evening :)
x HP32SII HP25 HP12CP HP15C

0.9999999999 0.000810284684548 0.0008102769138 0.000810284684566 0.0008102846845
0.999999999 0.00256234515652 0.002562324556 0.00256234515715 0.002562345157
0.99999999 0.00810284685217 0.008102826434 0.00810284687236 0.008102846852
0.99999995 0.0181185164331 0.01811852976 0.0181185166595 0.01811851643
0.9999999 0.0256234517765 0.02562347473 0.0256234524170 0.02562345178
0.999999 0.0810284752065 0.08102849351 0.0810284752060 0.08102847521
0.999998 0.114591578125 0.1145916163 0.114591578125 0.1145915781
0.999997 0.140345459308 0.1403455535 0.140345459308 0.1403454593
0.999995 0.181185239070 0.1811854695 0.181185239070 0.1811852391
0.999992 0.229183270841 0.2291830034 0.229183270841 0.2291832708
In this table the results are exact to 12 and 10 digits on the HP32SII and HP15C. The HP25 has errors beginning in the fifth significant digits. In this range the HP12CP program is closer to the HP15C than the HP25. Out of this critical range, the program has no problems. However the HP25, a commercial product in the mid 70's, presents an error in the tenth significant digit, when computing acos(0.9). Considering I was mostly interested in a fast running program, I am still pleased with the program. It's faster than the russian calculator in the link, if this can be used as a comparison. Of course I am also pleased with the overall accuracy, thanks to the extra digits in the Platinum.
Cheers,
Gerson.

_{ Edited to include this table:}
After debugging (by using acos(x) = 2*atan(sqrt((1  x)/(1+x))), as suggested by Tony):
x HP32SII HP25 HP12CP HP15C

0.9999999999 0.000810284684548 0.0008102769138 0.000810284684543 0.0008102846845
0.999999999 0.00256234515652 0.002562324556 0.00256234515651 0.002562345157
0.99999999 0.00810284685217 0.008102826434 0.00810284685213 0.008102846852
0.99999995 0.0181185164331 0.01811852976 0.0181185164330 0.01811851643
0.9999999 0.0256234517765 0.02562347473 0.0256234517764 0.02562345178
0.999999 0.0810284752065 0.08102849351 0.0810284752060 0.08102847521
0.999998 0.114591578125 0.1145916163 0.114591578125 0.1145915781
0.999997 0.140345459308 0.1403455535 0.140345459308 0.1403454593
0.999995 0.181185239070 0.1811854695 0.181185239069 0.1811852391
0.999992 0.229183270841 0.2291830034 0.229183270840 0.2291832708
_{The older version may be definitely discarded.}
Edited: 13 Dec 2006, 3:52 p.m. after one or more responses were posted
Posts: 93
Threads: 2
Joined: Jul 2005
Hi Gerson, I checked out the acos(.9999999) on the train this morning. This may be a coincidence but your HP12cp result is *exactly* what I get from a two term Taylor:
(x  (x^3)/3)*180/PI = 0.0256234524170
where x=sqrt(1.9999999^2)/.9999999
To get 17765 for the last five places we need to prolong the Taylor.
Hope this helps in diagnosing the problem.
Cheers,
Tony
Edited: 4 Dec 2006, 9:34 p.m.
Posts: 2,761
Threads: 100
Joined: Jul 2005
Hello Tony,
This is really a weird coincidence since I haven't used Taylor's series. I have used a polynomial approximation instead.
The problem is due to floating point errors, as I suspected. Just take a look at this stepbystep example:
acos(x) = atan(sqrt(1  x^2)/x)
12CP ACTUAL
x = 0.999999950000000
x^2 = 0.999999900000000 0.999999900000002500000000000
1  x^2 = 0.000000100000000 0.000000099999997500000000000
sqrt(1  x^2) = 0.000316227766017 0,000316227762063990833284121
sqrt(1  x^2)/x = 0.000316227781828 0.000316227777875379727053107
atan(sqrt(1  x^2)/x) = 0.0181185166595 0.018118516433109151619077930

atan(0.000316227781828) = 0.018118516659577588629936690
That is, the program gives the best answer given the 12CP is a 12digit calculator. It seems HP calculators use double precision in this critical range. Perhaps that's why they have pi to 20 or 24 four places built in.
This problem doesn't occur with asin(x) because the error is not critical for arguments in that range.
Perhaps I should go back to the older version, since it's shorter. After all, since the errors begin in the seventh digit there is no problem in truncating the last two digits, since they are wrong anyway.
acos(0.999996272743)
Old version: 0.156434461500
Revised version: 0.156434461498
actual: 0.156434462627
acos(0.999996272738)
Old version: 0.156434566400
Revised version: 0.156434566425
actual: 0.156434567553
Perhaps I should use those extra steps to implement the HP rounding philosophy as Rodger Rosenbaum has suggested (see my next reply) :)
Regards,
Gerson.
Posts: 93
Threads: 2
Joined: Jul 2005
Hi Gerson,
Quote:
This is really a weird coincidence since I haven't used Taylor's series. I have used a polynomial approximation instead.
The problem is due to floating point errors, as I suspected. Just take a look at this stepbystep example:
acos(x) = atan(sqrt(1  x^2)/x)
OK, try
acos(x) = 2*atan(sqrt((1  x)/(1+x)))
Cheers,
Tony
Posts: 305
Threads: 17
Joined: Jun 2007
Quote: The problem is due to floating point errors, as I suspected. Just take a look at this stepbystep example:
acos(x) = atan(sqrt(1  x^2)/x)
Try this formula for acos:
acos(x) = 2 * atan(SQRT(1x)/SQRT(1+x)), for x>0, and with appropriate changes for x<0.
I think that when x is near 1, this formula will retain more accurate digits (if the atan function is accurate).
Quote: It seems HP calculators use double precision in this critical range. Perhaps that's why they have pi to 20 or 24 four places built in.
The Saturn based machines have a builtin 31 digit pi for range reduction of trig functions, and do 31 digit arithmetic operations during the range reduction.
That's why they can get accurate results from trig functions in radian mode whose arguments are as large as 10^19. I know of no other hand held calculators that can do that (most of them just error out for arguments that large).
Posts: 2,761
Threads: 100
Joined: Jul 2005
Tony:
Quote:
try
acos(x) = 2*atan(sqrt((1  x)/(1+x)))
Rodger:
Quote:
Try this formula for acos:
acos(x) = 2 * atan(SQRT(1x)/SQRT(1+x)), for x>0, and with appropriate changes for x<0.
I think that when x is near 1, this formula will retain more accurate digits (if the atan function is accurate).
Thanks Rodger and Tony for your suggestions. I made some tests and the first expression appears to grant enough accuracy. If it doesn't I may use the second. The math functions are very fast on the Platinum, even LN and e^x. So, another SQRT won't delay the answer. And yes, the atan function is accurate to 11 digits all through the valid range (x <= 9.99999999E49), (no counter example found so far). The program gives the same answers of the 15C, except for acos(x) (when x >= 0.999999). If the acos function can be fixed, the forensic result will match exactly the one on the 15C.
Regards,
Gerson.
Posts: 305
Threads: 17
Joined: Jun 2007
Gerson,
Those two formulas are mathematically equivalent, of course. But you will see a difference in the errors associated with the values they return. Sometimes one formula will return an exactly correct value, and other times the other one will.
For example, for an input argument of .999999999966, the formula with only one square root gives the correct result. For an input argument of .999999999964, the formula with two square roots gives the correct result. These calculations were made on the HP50.
Both formulas should give an error of only 1 or 2 units in the LSD for arguments near 1.
Unless you can find a way to get a couple of extra digits (guard digits) in your computations, you won't be able to get perfectly exact (to 10 or 12 digits) results.
Posts: 2,761
Threads: 100
Joined: Jul 2005
Hello Rodger,
Since the program runs on the 12C Platinum, only 10 significant digits can normally be entered by the user. So arguments up to 9.999999999E01 will be entered. And considering two guard digits are available and the atan(x) function is accurate to 11 significant digits I am confident it will be possible to replicate the HP15C answers, according to the tests I've done so far. And yes, I am convinced it's better to have always 10 correct digits than sometimes 11 and sometimes 12, so I decided to use the HP rounding philosophy, which looks wonderful :)
The tests below have been performed on the HP12C Platinum. As you can see, the atan(x) function is accurate enough to display at least 10 correct digits.
Thanks for the suggestion and the past "What should we get" threads.
Regards,
Gerson.

acos(x) = 2 * atan(sqrt((1  x)/(1 + x)))
x = 0.9999999500
1x = 0.000000050000
1+x = 1.999999950000
(1x)/(1+x) = 0.000000025000000625
sqrt((1x)/(1+x)) = 0.000158113884985
atan(sqrt((1x)/(1+x))) = 0.00905925821651
acos(x)=2*atan(sqrt((1x)/(1+x))) = 0.0181185164330
acos(x)=2*atan(sqrt(1x)/sqrt(1+x))
x = 0.9999999500
1x = 0.000000050000
1+x = 1.99999995000
sqrt(1x) = 0.000223606797750
sqrt(1+x) = 1.41421354470
sqrt(1x)/sqrt(1+x) = 0.000158113884984
atan(sqrt(1x)/sqrt(1+x)) = 0.00905925821645
2*atan(sqrt(1x)/sqrt(1+x)) = 0.0181185164329
actual: 0.0181185164331

acos(x) = 2 * atan(sqrt((1  x)/(1 + x)))
x = 0.999999900000
1x = 0.000000100000
1+x = 1.99999999000
(1x)/(1+x) = 0.0000000500000025000
sqrt((1x)/(1+x)) = 0.000223606803340
atan(sqrt((1x)/(1+x))) = 0.0128117258820
acos(x)=2*atan(sqrt((1x)/(1+x))) = 0.0256234517764
acos(x)=2*atan(sqrt(1x)/sqrt(1+x))
x = 0.999999900000
1x = 0.000000100000
1+x = 1.99999999000
sqrt(1x) = 0.000316227766017
sqrt(1+x) = 1.41421352702
sqrt(1x)/sqrt(1+x) = 0.000223606803340
atan(sqrt(1x)/sqrt(1+x)) = 0.0128117258882
2*atan(sqrt(1x)/sqrt(1+x)) = 0.0256234517764
actual: 0.0256234517765

Edited: 5 Dec 2006, 6:12 p.m.
