The Museum of HP Calculators

HP Articles Forum

[Return to the Index ]
[ Previous | Next ]


Fast and Accurate Trigonometric Functions on the HP-12C Platinum

Posted by Gerson W. Barbosa on 6 Aug 2006, 3:10 p.m.

The program

001- 2 023- CHS 045- g x<=y 067- / 089- g GTO 266 002- STO 6 024- g x<=y 046- g GTO 017 068- + 090- 2 003- Rv 025- g GTO 031 047- Rv 069- RCL 8 091- STO 6 004- RCL .0 026- 2 048- x<>y 070- * 092- Rv 005- / 027- * 049- STO 8 071- ENTER 093- RCL .0 006- x<>y 028- x<>y 050- g x^2 072- g x^2 094- / 007- STO 4 029- - 051- ENTER 073- 4 095- CHS 008- x<>y 030- ENTER 052- ENTER 074- * 096- 9 009- ENTER 031- Rv 053- ENTER 075- CHS 097- 0 010- g x^2 032- ENTER 054- RCL .2 076- 3 098- + 011- g SQRT 033- g x^2 055- * 077- + 099- g GTO 006 012- 3 034- g SQRT 056- RCL .3 078- * 100- RCL .0 013- 6 035- ENTER 057- - 079- RCL 6 101- / 014- 0 036- ENTER 058- * 080- 2 102- STO 4 015- g x<=y 037- 9 059- RCL .4 081- g x<=y 103- x<>y 016- g GTO 120 038- 0 060- + 082- g GTO 087 104- STO 7 017- Rv 039- - 061- * 083- Rv 105- 1 018- Rv 040- g x=0 062- RCL .5 084- g x=0 106- STO 6 019- 9 041- g GTO 047 063- - 085- g GTO 114 107- RCL 4 020- 0 042- Rv 064- * 086- g GTO 109 108- g GTO 009 021- g x<=y 043- 9 065- RCL .1 087- Rv 109- CLx 022- g GTO 026 044- 0 066- 3 088- x<>y 110- STO 6

111- x<>y 152- 0 193- g x<=y 234- RCL .7 275- Rv 112- RCL 4 153- * 194- g GTO 197 235- + 276- g GTO 000 113- g GTO 095 154- RCL .0 195- 1/x 236- * 114- x<>y 155- * 196- 9 237- RCL .8 115- RCL 4 156- g GTO 266 197- 0 238- - 116- x<>y 157- x<>y 198- RCL .1 239- * 117- / 158- STO 4 199- * 240- RCL .9 277- 1.745329252 118- RCL 7 159- x<>y 200- STO 6 241- + 288- EEX 119- g GTO 267 160- 2 201- x<>y 242- * 289- 2 120- x<>y 161- STO 7 202- 2 243- 3 290- CHS 121- EEX 162- Rv 203- ENTER 244- 1/x 291- STO .1 122- 1 163- 1 204- 3 245- - 292- 2.1 123- 2 164- x<>y 205- g SQRT 246- * 295- EEX 124- g x<=y 165- + 206- STO 9 247- 1 296- 26 125- g GTO 146 166- g x=0 207- - 248- + 298- CHS 126- Rv 167- g GTO 269 208- x<>y 249- RCL 9 299- STO .2 127- x<>y 168- g LSTx 209- g x<=y 250- * 300- 4.47569 128- / 169- 1 210- g GTO 222 251- RCL 8 307- EEX 129- g LSTx 170- x<>y 211- ENTER 252- + 308- 20 130- x<>y 171- - 212- ENTER 253- RCL 6 310- CHS 131- g INTG 172- g LSTx 213- RCL 9 254- g x=0 311- STO .3 132- * 173- 1 214- * 255- x<>y 312- 5.55391606 133- g x<=y 174- + 215- 1 256- x<>y 322- EEX 134- CHS 175- / 216- - 257- - 323- 14 135- + 176- g SQRT 217- x<>y 258- RCL 5 325- CHS 136- g GTO 019 177- g GTO 184 218- RCL 9 259- * 326- STO .4 137- x<>y 178- x<>y 219- + 260- RCL 7 327- 3.281837614 138- STO 4 179- STO 4 220- / 261- * 338- EEX 139- x<>y 180- x<>y 221- 3 262- RCL .1 339- 8 140- ENTER 181- 1 222- 0 263- / 340- CHS 141- g x^2 182- STO 7 223- RCL .1 264- RCL .0 341- STO .5 142- 1 183- x<>y 224- * 265- * 342- .0784 143- - 184- ENTER 225- STO 8 266- RCL 4 347- CHS 144- g x=0 185- g x^2 226- x<>y 267- x<>y 348- STO .6 145- g GTO 150 186- g SQRT 227- STO 9 268- g GTO 000 349- .110351 146- CHS 187- g x=0 228- g x^2 269- RCL 7 356- STO .7 147- g SQRT 188- g GTO 266 229- ENTER 270- g GTO 151 357- .142838 148- / 189- / 230- ENTER 271- 1 364- STO .8 149- g GTO 181 190- STO 5 231- ENTER 272- g GTO 274 365- .199999839 150- x<>y 191- 1 232- RCL .6 273- RCL .1 375- STO .9 151- 9 192- g LSTx 233- * 274- STO .0 376- g GTO 271

Lines 277 through 376 are optional. This is an initialization routine that loads the constants used by the program and sets the angular to degrees. This has to be done before running the program for the first time or whenever the registers in the range R.1 through R.9 are accidentally cleared or changed (un unlikely possibility since at least three key presses would be required: STO . 0..9; f CLEAR REG). In case these lines are not keyed in the constants should be entered manually:

R.1: 1.745329252E-02 (pi/180)

R.2: 2.1E-26 R.3: 4.47569E-20 R.4: 5.55391606E-14 R.5: 3.281837614E-08 (sine polynomial coefficients)

R.6: -0.0784 R.7: 0.110351 R.8: 0.142838 R.9: 0.199999839 (arctangent polynomial coefficients)

The next arctangent polynomial coefficient, 1/3, is built into the program (lines 243 and 244). The registers 4 through .0 are used by the program but they can be overwritten anytime. The registers 0 through 3 and the financial registers are always available to the user.

Usage

SIN: R/S (*) COS: g GTO 090 R/S TAN: g GTO 100 R/S ASIN: g GTO 137 R/S ACOS: g GTO 157 R/S ATAN: g GTO 178 R/S

DEG: g GTO 271 R/S RAD: g GTO 273 R/S

g GTO 277 R/S: If available, this loads all constants needed by the program.

(*) When just coming from an error condition due to an invalid argument, use g GTO 001 R/S instead, as the calculator will not return to program line 000.

Another way to set the angular mode to DEGREES is just storing 1 in register R.0: 1 STO .0. Likewise, .9 1/x STO .0 sets the angular mode to GRADS. As a convenience, register .1 (pi/180) may be used to convert to and from RADIANS.

Input ranges

SIN, COS, TAN: |x| < 1E12 (degrees)

ASIN, ACOS: |x| <= 1

ATAN: |x| <= 9.99999999E49

Notes:

1) Although in DEGREES mode the answers of SIN, COS and ATAN functions are accurate at least to 11 significant figures all through the valid ranges, in RADIANS mode a loss of accuracy occurs for larger absolute values of x.

2) The program will return 9.999999E99 for input arguments with absolute values equal or greater than 1E50. Arguments out of the allowed ranges will cause the program to return Error 0. Arguments for which the tangent function are not defined, such as 90 degrees, will also cause an Error 0 message to appear.

3) Pay attention to what appears to be a bug in the HP-12C Platinum: negating a number just entered by pressing the CHS (change sign) key does not work the way it used to on the classic HP-12C. For instance:

Keystrokes Display

45 ENTER 45.00000000 CHS -45.00000000 The number displayed by the calculator is negative, 5 5. + 50.00000000 yet it has been interpreted as a positive number.

15 ENTER 15.00000000 we've just entered 15, CHS -15.00000000 but then we rememember we want to compute sin(-15) and press CHS R/S 0.258819045 then R/S. Wrong result, the calculator shows sin(15) instead!

So, either we enter the correct number the first time or we remember to press ENTER after pressing CHS:

15 CHS R/S -0.258819045 sine(-15) 15 ENTER 15.00000000 CHS ENTER -15.00000000 R/S -0.258819045 sine(-15), as expected.

Running Times

SIN, COS: 2.5 s (|x| <= 180) TAN: 4.9 s (|x| <= 180) ASIN, ACOS, ATAN: 2.6 s

Sample Calculations

First set the display to FIX 9 and DEG modes if not already set:

f 9; g GTO 271 R/S

1) asin(acos(atan(tan(cos(sin(9)))))):

Keystrokes Display 9 9. R/S 0.156434465 g GTO 090 R/S 0.999996273 g GTO 100 R/S 0.017455000 g GTO 178 R/S 0.999996273 g GTO 157 R/S 0.156434568 g GTO 137 R/S 9.000005947

2) sin2(30) + cos2(30):

Notice that the stack register X is always preserved, so there's no need to store intermediate results in simple chain calculations (one-level only), as in this example.

Keystrokes Display

30 R/S 0.500000000 g x2 0.250000000 30 g GTO 090 R/S 0.866025404 g x2 0.750000000 + 1.000000000

3) sin(-1234567.8) + acos(-0.77)

Keystrokes Display

1234567.8 CHS -1,234,567.8 R/S -0.790155012 .77 CHS g GTO 157 R/S 140.3538889 + 139.5637338

4) tan(-6.2 rad)/asin(0.15 rad)

Set the angular mode to RADIANS:

g GTO 273 R/S

Keystrokes Display

6.2 CHS g GTO 100 R/S 0.083377715 0.15 g GTO 137 R/S 0.150568273 / 0.553753545

g GTO 271 R/S 0.553753545 ; back to DEG mode

Accuracy comparison with other HP calculators

Sin(x):

x (deg) HP-12CP HP-50G (2006) HP-15C (1982) HP-35 (1972) ---------------------------------------------------------------------------------- 0.000000 0.00000000000E+00 0.00000000000E+00 0.000000000E+00 0.000000000E+00 0.000010 1.74532925200E-07 1.74532925199E-07 1.745329252E-07 1.745000000E-07 0.000110 1.91986217720E-06 1.91986217719E-06 1.919862177E-06 1.919800000E-06 0.022000 3.83972426005E-04 3.83972426004E-04 3.839724260E-04 3.839723910E-04 3.330000 5.80867495981E-02 5.80867495977E-02 5.808674960E-02 5.808674960E-02 14.44000 2.49366025116E-01 2.49366025115E-01 2.493660251E-01 2.493660250E-01 25.55000 4.31298587031E-01 4.31298587031E-01 4.312985870E-01 4.312985871E-01 36.66000 5.97065256389E-01 5.97065256389E-01 5.970652564E-01 5.970652561E-01 47.77000 7.40452782677E-01 7.40452782677E-01 7.404527827E-01 7.404527828E-01 58.88000 8.56086728293E-01 8.56086728292E-01 8.560867283E-01 8.560867285E-01 69.99000 9.39632912698E-01 9.39632912698E-01 9.396329127E-01 9.396329127E-01 81.11000 9.87986852782E-01 9.87986852778E-01 9.879868528E-01 9.879868527E-01 88.88000 9.99808950037E-01 9.99808950038E-01 9.998089500E-01 9.998089499E-01 89.99900 9.99999999848E-01 9.99999999848E-01 9.999999998E-01 9.999999998E-01 89.99990 9.99999999998E-01 9.99999999998E-01 1.000000000E+00 1.000000000E+00 90.00000 1.00000000000E+00 1.00000000000E+00 1.000000000E+00 1.000000000E+00

Tan(x):

x (deg) HP-12CP HP-50G (2006) HP-15C (1982) HP-35 (1972) ---------------------------------------------------------------------------------- 0.000000 0.00000000000E+00 0.00000000000E+00 0.000000000E+00 0.000000000E+00 0.000010 1.74532925200E-07 1.74532925199E-07 1.745329252E-07 1.745000000E-07 0.000110 1.91986217721E-06 1.91986217720E-06 1.919862177E-06 1.919800000E-06 0.022000 3.83972454310E-04 3.83972454309E-04 3.839724543E-04 3.839724542E-04 3.330000 5.81849926709E-02 5.81849926706E-02 5.818499267E-02 5.818499260E-02 14.44000 2.57500649118E-01 2.57500649118E-01 2.575006491E-01 2.575006490E-01 25.55000 4.78047179798E-01 4.78047179799E-01 4.780471798E-01 4.780471798E-01 36.66000 7.44291588299E-01 7.44291588300E-01 7.442915883E-01 7.442915880E-01 47.77000 1.10168657755E+00 1.10168657755E+00 1.101686578E+00 1.101686578E+00 58.88000 1.65641139080E+00 1.65641139081E+00 1.656411391E+00 1.656411391E+00 69.99000 2.74598611677E+00 2.74598611678E+00 2.745986117E+00 2.745986119E+00 81.11000 6.39316645081E+00 6.39316645080E+00 6.393166451E+00 6.393166426E+00 88.88000 5.11504299317E+01 5.11504299320E+01 5.115042993E+01 5.115042860E+01 89.99900 5.72957795071E+04 5.72957795073E+04 5.729577951E+05 5.729655162E+04 89.99990 5.72957795128E+05 5.72957795130E+05 5.729577951E+05 5.730193057E+05 90.00000 Error 0 TAN Error 9.999999999E+99 9.999999999E+99

ArcTan(x):

x HP-12CP HP-50G (2006) HP-15C (1982) HP-35 (1972) ---------------------------------------------------------------------------------- 0.00000 0.00000000000E+00 0.00000000000E+00 0.000000000E+00 0.000000000E+00 0.00011 6.30253572098E-03 6.30253572102E-03 6.302535721E-03 6.302535688E-03 0.15500 8.81073298598E+00 8.81073298598E+00 8.810732986E+00 8.810732984E+00 0.26795 1.50000431708E+01 1.50000431708E+01 1.500004317E+01 1.500004317E+01 0.41421 2.24998257819E+01 2.24998257819E+01 2.249982578E+01 2.249982579E+01 0.57735 2.99999884325E+01 2.99999884324E+01 2.999998843E+01 2.999998843E+01 0.77700 3.78472067673E+01 3.78472067673E+01 3.784720677E+01 3.784720676E+01 0.88800 4.16050764577E+01 4.16050764579E+01 4.160507646E+01 4.160507646E+01 1.00000 4.50000000000E+01 4.50000000000E+01 4.500000000E+01 4.500000000E+01 1.22200 5.07054870168E+01 5.07054870169E+01 5.070548702E+01 5.070548702E+01 1.48880 5.61114572284E+01 5.61114572285E+01 5.611145723E+01 5.611145722E+01 2.11100 6.46526573517E+01 6.46526573518E+01 6.465265735E+01 6.465265735E+01 4.88800 7.84378235856E+01 7.84378235855E+01 7.843782359E+01 7.843782360E+01 7.55500 8.24600068273E+01 8.24600068272E+01 8.246000683E+01 8.246000679E+01 99.9990 8.94270555733E+01 8.94270555733E+01 8.942705557E+01 8.942705555E+01 3333333 8.99999828113E+01 8.99999828113E+01 8.999998281E+01 8.999998281E+01 1.1E+12 8.99999999999E+01 8.99999999999E+01 9.000000000E+01 9.000000000E+01

Of course, the HP-12C Platinum cannot display all those twelve significant digits simultaneously. They were obtained in most cases by multiplying the number by a power of 10 and taking the fractional part of the result. As we can notice in these examples, all answers are accurate to at least 11 digits, so the results displayed by the HP-12C will most of the time match those displayed by all scientific HP 10-series calculators (HP-10C, HP-11C and HP-15C). On some occasions there might be a difference of one unit in the last significant displayed digit. This happens because these calculators round the answers to 10 digits and on rare occasions there may be a difference on one unit in the last digit of the mantissa, as we can see in the following examples (DEG mode assumed):

1) HP-12CP HP-15C ----------------------------------------------------------------------- Keystrokes Display Keystrokes Display

32.888 32.888 32.888 32.888 R/S 0.542998588 SIN 0.542998589 f CLEAR PREFIX 5429985885 f CLEAR PREFIX 5429985885

sin(32.888) to 16 digits: 0.5429985884655923 rounding to 12 digits: 0.542998588466 HP-12CP answer: 0.542998588466 displayed 9-digit answer: 0.542998588 ; correct on the HP-12CP, because of the leading '4' 10-digit mantissa: 5429985885 ; correct on the HP-12CP, because of the leading '6'

rounding to 10 digits: 0.5429985885 displayed 9-digit answer: 0.542998589 ; correct on the HP-15C, because of the leading '5'

In this case, the answer of the HP-15C is right, according to its rounding method.

2) HP-12CP HP-15C ----------------------------------------------------------------------- Keystrokes Display Keystrokes Display

34.444 34.444 34.444 34.444 R/S 0.565600478 SIN 0.565600479 f CLEAR PREFIX 5656004784 f CLEAR PREFIX 5656004785

sin(34.444) to 16 digits: 0.5656004784499306 rounding to 12 digits: 0.565600478450 HP-12CP answer: 0.565600478449 displayed 9-digit answer: 0.565600478 ; correct on the HP-12CP, because of the leading '4' 10-digit mantissa: 5656004784 ; correct on the HP-12CP, because of the leading '4'

rounding to 10 digits: 0.5656004784 displayed 9-digit answer: 0.565600479 ; incorrect on the HP-15C, because of the leading '4' 10-digit mantissa: 5656004785 ; incorrect on the HP-15C, because of the leading '4'

In this case, according to its rounding method, the HP-15C should return 0.5656004784. Unless, of course, a slight error in the thirteenth significant digit has turned those '499' into '500' or greater.

Since the HP-12CP does not round the answers to 10 digits it doesn't make much sense to implement this feature in the program. Anyway, this is possible by means of a slight modification:

1) Change the line 268 to

268-	g GTO 377

2) Key in these lines:

	
377-	ENTER
378-	g x^2
379-	g SQRT
380-	g x=0
381-	g GTO 000
382-	g LN
383-	1
384-	0
385-	g LN
386-	/
387-	g INTG
388-	1
389-	-
390-	1
391-	0
392-	x<>y
393-	y^x
394-	/
395-	g LSTx
396-	x<>y
397-	f RND
398-	x<>y
399-	*

This will round the answers to the number of the decimal places of the display plus one. However results with absolute value less than 1E-50 will be rounded to zero.

After this change, let's do the first example again, actually a standard test created by Mike Sebastian to check trigonometric algorithms in calculators, called Calculator Forensics. The angular mode should be set do DEG and the display to FIX 9 mode.

1) asin(acos(atan(tan(cos(sin(9)))))):

Keystrokes Display

9 9. R/S 0.156434465 g GTO 090 R/S 0.999996273 g GTO 100 R/S 0.017455000 g GTO 178 R/S 0.999996272 g GTO 157 R/S 0.156441660 g GTO 137 R/S 9.000417403

The result exactly matches that of the HP-15C and other HP 10-digit calculators. Because of the difference in the last digit that occasionally may happen, as we have seen, this rounding won't assure the answers will always be the same on both calculators.

Don't forget to revert to the original program by editing the line 268:

g GTO 267
f P/R
g GTO 000
f P/R

The method

In order to achieve results that accurate under relatively short times, we used minimax polynomial approximations and range reduction. This technique yields maximum accuracy throughout a desired range by approximating a function by means of a polynomial of as low degree as possible. In this case, the sine function is approximated by a 9th-degree polynomial in the range [0..pi/6] and the arctangent function is approximated by an 11th-degree polynomial in the range [0..2-sqrt(3)]. These are not the full polynomials with all coefficients of all available powers we would obtain by default when calling minimax in most numerical math packages. Instead, by means of an initial work, we have obtained polynomials with odd powers only and with lowest-order coefficients equal to one, more akin to the classical infinite series for these functions.

Programs that approximate trigonometric functions using the plain Taylor series are generally more compact. Nevertheless, to obtain the same accuracy throughout these ranges using Taylor series we would need to evaluate up to 20 percent more terms in the sine polynomial and up to 50 percent more terms in the arctangent polynomial, not to mention that at the end of each iteration a test should be made to verify the accuracy. Furthermore, computing the coefficients on the fly would make the program even slower. Of course, in the older version of the HP-12C calculator, with only 99 programming steps and very few general use storage registers, there were few options left, if any. However, the HP-12C Platinum, with four time as much programming memory and 20 general use registers always available, is quite suitable for this solution. Since there's plenty of memory, the program is not so optimized as it could be, so there's still room for improvement. Anyway, the results obtained so far are better than what we might have expected. In fact, the program turns the HP-12C Platinum into the most accurate HP 10C-series calculator, regarding trigonometric functions.

The following expressions have been used:

(i) sin(3x) = sin(x)*(3 - 4 sin2(x))

(ii) sin(x) = x(a1 + x2(a2 + x2(a3 + x2(a4 + a5x2)))), x in radians;

These are the original coefficients, optimized for the range [0..pi/6]:

a1: 1.00000000000E+00 a2: -1.66666666666E-01 a3: 8.33333320429E-03 a4: -1.98410347969E-04 a5: 2.74201854577E-06

They are multiplied by successive odd powers of pi/540 for radians to degrees conversion and to avoid an initial division by three:

a1: 5.81776417333E-03 a2: -3.28183761371E-08 a3: 5.55391605878E-14 a4: -4.47566022234E-20 a5: 2.09351175453E-26

This yields a maximum absolute error of 8E-14. As we don't need this level of accuracy, the coefficients used in the program have been slightly modified to use as few digits as possible while keeping the maximum absolute error under 1.2E-12.

(iii) atan(x)=(pi/6 + atan((|x|*srqt(3) - 1)/(|x| + sqrt(3)))) * x/|x|, (2 - sqrt(3)) < |x| <= 1;

(iv) arctan(x) = (pi/2 - arctan(1/|x|) * x/|x|, |x| > 1;

(v) atan(x) = x(a1 + x2(a2 + x2(a3 + x2(a4 + x2(a5 + a6x2))))), x in radians;

These are the original coefficients, optimized for the range [0..2-sqrt(3)]:

a1: 1.00000000000E+00 a2: -3.33333333089E-01 a3: 1.99999829275E-01 a4: -1.42837953536E-01 a5: 1.10350743154E-01 a6: -7.83982197548E-02

The maximum absolute error with the original coefficients is less than 5E-12. Again, in order to minimize the number of digits, the following set of coefficients has been chosen. Although this would in theory give us a maximum absolute error of 7.8E-12, the maximum observed error is 5E-12. It is possible the internal floating point errors have compensated the maximum error for better. Also, a small part of the error is due to the use of 1.745329252E-02 as the angle conversion factor instead of 1.74532925199E-02. Entering the exact value to 12 digits has not produced any noticeable improvement in the overall accuracy though.

a1: 1.00000000000E+00 a2: -3.33333333333E-01 a3: 1.99999839000E-01 a4: -1.42838000000E-01 a5: 1.10351000000E-01 a6: -7.84000000000E-02

(vi) cos(x) = sin(pi/2 - x)

(vii) tan(x) = sin(x)/cos(x)

(viii) arcsin(x) = arctan(x/(sqrt(1 - x2)) |x| <= 1;

(ix) arccos(x) = 2 * arctan(sqrt((1 - x)/(1 + x))) 0 < x <= 1;

References

1) For some methods to approximate trigonometric functions, refer to the following paper, by Robin Green of Sony Computer Entertainment America:

Faster Math Functions

http://www.research.scea.com/research/pdfs/RGREENfastermath_GDC02.pdf

2) HP-71B Minimax Polynomial Fit, by Valentin Albillo:

http://membres.lycos.fr/albillo/calc/pdf/DatafileVA016.pdf

Edited: 30 Dec 2006, 9:15 p.m.

Password:

[ Return to the Message Index ]

Go back to the main exhibit hall