HP Forums

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.

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

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

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.

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.

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

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

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.

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.

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.

Gerson W. Barbosa

tony (nz)

Gerson W. Barbosa

tony (nz)

Gerson W. Barbosa

tony (nz)

Rodger Rosenbaum

Gerson W. Barbosa

Rodger Rosenbaum

Gerson W. Barbosa