Posts: 2,761
Threads: 100
Joined: Jul 2005
Some days ago, while I was concerned about the HP-33S 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.999999999E-01) returning 8.102847000E-04 rather than 8.102846845E-04. 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 C3-15, 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 E-13 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 12-digit 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 HP-32SII HP-25 HP-12CP HP-15C
--------------------------------------------------------------------------------------------------------
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 HP-32SII and HP-15C. The HP-25 has errors beginning in the fifth significant digits. In this range the HP-12CP program is closer to the HP-15C than the HP-25. Out of this critical range, the program has no problems. However the HP-25, 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 HP-32SII HP-25 HP-12CP HP-15C
--------------------------------------------------------------------------------------------------------
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 step-by-step 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 12-digit 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 step-by-step 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 step-by-step example:
acos(x) = atan(sqrt(1 - x^2)/x)
Try this formula for acos:
acos(x) = 2 * atan(SQRT(1-x)/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 built-in 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(1-x)/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.999999999E-01 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 HP-15C 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 HP-12C 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
1-x = 0.000000050000
1+x = 1.999999950000
(1-x)/(1+x) = 0.000000025000000625
sqrt((1-x)/(1+x)) = 0.000158113884985
atan(sqrt((1-x)/(1+x))) = 0.00905925821651
acos(x)=2*atan(sqrt((1-x)/(1+x))) = 0.0181185164330
acos(x)=2*atan(sqrt(1-x)/sqrt(1+x))
x = 0.9999999500
1-x = 0.000000050000
1+x = 1.99999995000
sqrt(1-x) = 0.000223606797750
sqrt(1+x) = 1.41421354470
sqrt(1-x)/sqrt(1+x) = 0.000158113884984
atan(sqrt(1-x)/sqrt(1+x)) = 0.00905925821645
2*atan(sqrt(1-x)/sqrt(1+x)) = 0.0181185164329
actual: 0.0181185164331
-------------------------------------------------------
acos(x) = 2 * atan(sqrt((1 - x)/(1 + x)))
x = 0.999999900000
1-x = 0.000000100000
1+x = 1.99999999000
(1-x)/(1+x) = 0.0000000500000025000
sqrt((1-x)/(1+x)) = 0.000223606803340
atan(sqrt((1-x)/(1+x))) = 0.0128117258820
acos(x)=2*atan(sqrt((1-x)/(1+x))) = 0.0256234517764
acos(x)=2*atan(sqrt(1-x)/sqrt(1+x))
x = 0.999999900000
1-x = 0.000000100000
1+x = 1.99999999000
sqrt(1-x) = 0.000316227766017
sqrt(1+x) = 1.41421352702
sqrt(1-x)/sqrt(1+x) = 0.000223606803340
atan(sqrt(1-x)/sqrt(1+x)) = 0.0128117258882
2*atan(sqrt(1-x)/sqrt(1+x)) = 0.0256234517764
actual: 0.0256234517765
----------------------------------------------------------
Edited: 5 Dec 2006, 6:12 p.m.
|