The Museum of HP Calculators

HP Forum Archive 16

[ Return to Index | Top of Index ]

Fast and Accurate Trigs (HP-12C Platinum) - Revised version.
Message #1 Posted by Gerson W. Barbosa on 3 Dec 2006, 12:46 a.m.

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.

      
Re: Fast and Accurate Trigs (HP-12C Platinum) - Revised version.
Message #2 Posted by tony (nz) on 3 Dec 2006, 3:39 p.m.,
in response to message #1 by Gerson W. Barbosa

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

            
Re: Fast and Accurate Trigs (HP-12C Platinum) - Revised version.
Message #3 Posted by Gerson W. Barbosa on 3 Dec 2006, 8:07 p.m.,
in response to message #2 by tony (nz)

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

                  
Re: Fast and Accurate Trigs (HP-12C Platinum) - Revised version.
Message #4 Posted by tony(nz) on 4 Dec 2006, 5:34 p.m.,
in response to message #3 by Gerson W. Barbosa

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.

                        
Re: Fast and Accurate Trigs (HP-12C Platinum) - Revised version.
Message #5 Posted by Gerson W. Barbosa on 4 Dec 2006, 9:35 p.m.,
in response to message #4 by tony(nz)

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.

                              
Re: Fast and Accurate Trigs (HP-12C Platinum) - Revised version.
Message #6 Posted by tony (nz) on 5 Dec 2006, 1:23 a.m.,
in response to message #5 by Gerson W. Barbosa

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

                              
Re: Fast and Accurate Trigs (HP-12C Platinum) - Revised version.
Message #7 Posted by Rodger Rosenbaum on 5 Dec 2006, 6:20 a.m.,
in response to message #5 by Gerson W. Barbosa

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).

                                    
Re: Fast and Accurate Trigs (HP-12C Platinum) - Revised version.
Message #8 Posted by Gerson W. Barbosa on 5 Dec 2006, 10:00 a.m.,
in response to message #7 by Rodger Rosenbaum

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.

                                          
Re: Fast and Accurate Trigs (HP-12C Platinum) - Revised version.
Message #9 Posted by Rodger Rosenbaum on 5 Dec 2006, 5:14 p.m.,
in response to message #8 by Gerson W. Barbosa

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.

                                                
Re: Fast and Accurate Trigs (HP-12C Platinum) - Revised version.
Message #10 Posted by Gerson W. Barbosa on 5 Dec 2006, 6:10 p.m.,
in response to message #9 by Rodger Rosenbaum

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.


[ Return to Index | Top of Index ]

Go back to the main exhibit hall