HP-33S trigonometric inaccuracy -- new findings



#19

All --

Recently, Gerson Barbosa made me aware of a trigonometric-function bug in the HP-33S in this archived post:

http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv016.cgi?read=103065#103065

His example was the tangent of 89.999... degrees, but I've tried to get to the heart of the matter in that thread and subsequently.

It was opined in the thread that tan(x) might be computed using sin(x)/cos(x). Although not unreasonable, I doubt that it is correct, because timing tests using the following example program on the HP-32SII and HP-33S showed that computation of tangent is faster than that of sine and cosine, for values between 0 (noninclusive) and at least 85 degrees:

LBL L
30 (or other value)
TAN (or SIN or COS)
DSE L
GTO L
RTN

(Store an integer value in L -- say, 100 -- then XEQ L.)

Perhaps the tried-and-true CORDIC routines are still being used, in which tangent or cotangent is computed first, and then sine and cosine are calculated therefrom:

sin(x) = (+/-)tan(x) / sqrt[1 + tan2(x)]

cos(x) = (+/-)cot(x) / sqrt[1 + cot2(x)]

The above equations are printed in an article from the June 1977 Hewlett-Packard Journal, "Personal Calculator Algorithms II: Trigionometric Functions". (A scan of this article is in "77JUNAL2.PDF" from the MoHPC CD/DVD set.)


Here are some calculations to 12 significant digits using an angle of 0.0001 radians:

104 * tan(0.0001)

HP-32SII 1.00000000333
HP-33S 1.00000000330
actual 1.00000000333

105 * sin(0.0001)

HP-32SII 9.99999998333
HP-33S 9.99999998300
actual 9.99999998333

So, we can see that one significant digit is lost for the scaled tangent in the HP-33S, and two digits are lost for the scaled sine. The additional lost digit could be a result of the higher multiplier needed to reveal the last digit, or might be related to calculating sine from tangent.

My "actual" values came from the Windows XP Calculator, but one can easily confirm the results using Taylor series with terms less than the 15th significant digit truncated:

sin(x) = x - x3/3! + x5/5! - ...

cos(x) = 1 - x2/2! + x4/4! - ...

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

For x = 0.0001, the first two terms will suffice for 15 significant digits:

sin(0.0001) ~= 10-4 - 10-12/6
~= 0.0000999999998333
~= 9.99999998333 x 10-5

cos(0.0001) ~= 1 - 10-8/2
~= 0.9999999950000000

sin(0.0001)
tan(0.0001) = -----------
cos(0.0001)

10-4 - 10-12/6
~= -------------
1 - 10-8/2

Multiplying through by (1 + 10-8/2),

10-4 - 10-12/6 + 10-12/2 - 10-20/12
~= --------------------------------
1 - 10-16/4

Truncating the smallest terms and combining,

~= 10-4 + 10-12/3

~= 0.000100000000333

~= 1.00000000333 x 10-4

And there you have it -- the HP-32SII gave correct answers to 12 significant digits for sine and tangent with small arguments, but the HP-33S does not. Chalk up another for the "HP-33S bug list".

-- KS

Edited: 10 Dec 2006, 1:27 a.m.


#20

Quote:
Perhaps the tried-and-true CORDIC routines are still being used, in which tangent or cotangent is computed first, and then sine and cosine are calculated therefrom

Hi Karl,

I thought modern calculators didn't use the CORDIC algorithm anymore. By what I know the CORDIC algorithms were used in early calculators because their processors lacked multiplying instructions. They had however left and right shifts instructions (that is, simple multiplication and dive by 2 instructions) which were required by the method. But it's possible HP have sticked to CORDIC because the algorithms were tested and stable. I really don't know which algorithms are used on various calculators.

Quote:
the HP-32SII gave correct answers to 12 significant digits for sine and tangent with small arguments, but the HP-33S does not. Chalk up another for the "HP-33S bug list".

Even my program gives correct answers to 12 significant digits to your examples on the 12C Platinum, although my commitment is to 11 digits. By the way, now that the major bugs have been fixed, what is the updated "HP-33S bug list"?

Regards,

Gerson.


Edited: 10 Dec 2006, 1:57 p.m.


#21

Hi, Gerson --

Quote:
I thought modern calculators didn't use the CORDIC algorithm anymore. By what I know the CORDIC algorithms were used in early calculators because their processors lacked multiplying instructions.

I've been trying to learn more about CORDIC, having read the original 1959 article and a 1971 article provided courtesy of Jacques Laporte, as well as the 1977 Hewlett-Packard Journal article is cited in my last post. I still don't really understand the procedure.

I actually hope that you're right, in that more-direct methods of calculation would be employed today. Any microprocessor for a handheld calcuator that didn't have built-in commands for multiplication and division seems unfathomable.

It's clear that tan(x) is not generated by sin(x)/cos(x), based on the timing tests. Obviously, a more-efficient method is being used, but I haven't really researched what it is.

A 1984 article in Hewlett-Packard Journal ("84JUL71.PDF") about the HP-71B contains this blurb about its microprocessor:

Quote:

"The instruction set is based on the style of architecture used in the HP-41 computer and is oriented toward the arithmetic nature of handheld calculators with the capability of both hexadecimal and binary-coded decimal (BCD) arithmetic."

No specifics about methods, however...

Quote:
By the way, now that the major bugs have been fixed, what is the updated "HP-33S bug list"?

I was referring more to a comprehensive list of bugs on the originally-released version, many of which -- including mine -- are still out there. I know of no computational bugs other than the three that were fixed and this trigonometric one, but if other users keep finding 'em, I'll analyze 'em!

Regards,

-- KS


#22

Hello Karl,

Quote:
I've been trying to learn more about CORDIC, having read the original 1959 article and a 1971 article provided courtesy of Jacques Laporte, as well as the 1977 Hewlett-Packard Journal article is cited in my last post. I still don't really understand the procedure.

Neither do I. I don't have much intereste in the CORDIC because it is adequate only at the hardware level. I should read more about it though. Once I found a TI program that illustrated the method. Here is a version that will run on the HP-32SII (and on the HP-33S):

Cordic Algorithm - tan(x)

A01 LBL A
AO2 STO T
A03 1
A04 STO X
A05 STO Q
A06 STO i
A07 0
A08 STO Y
B01 LBL B
B02 RCL T
B03 7.E-13
B04 x>y?
B05 GTO E
C01 LBL C
C02 RCL(i)
C03 RCL T
C04 x>=y
C05 GTO D
C06 1
C07 STO+ i
C08 10
C09 STO÷ Q
C10 GTO C
D01 LBL D
D02 RCL(i)
D03 STO- T
D04 RCL X
D05 STO W
D06 RCL Q
D07 RCLx Y
D08 STO- X
D09 RCL Q
D10 RCLx W
D11 STO+ Y
D12 GTO B
E01 LBL E
E02 RCL Y
E03 RCL X
E04 ÷
E05 STO Z
E06 RTN

Constants Table

F01 LBL F
F02 1.015
F03 STO i
G01 LBL G
G02 RCL i
G03 IP
G04 1
G05 -
G06 +/-
G07 10^x
G08 ATAN
G09 STO(i)
G10 ISG i
611 GTO G
G12 RTN

HP-32SII:

LBL CK Size
A C117 012.0
B 9B62 015.5
C 914E 015.0
D 9C61 018.0
E 1A1B 009.0
F 4351 012.5
G 333A 018.0

HP-33S:

LBL CK Size
A FDA5 48
B 788F 27
C C143 54
D 6D9D 36
E 2D80 18
F D079 21
G F6A3 48

The program doesn't work for atan(0.0001) unless line B03 is changed to 5.E-15. In this case your first calculation will be 1.00000000330 on the HP-32SII and 1.00000000340 on the HP-33S.

Of course, on a real system the constants are hardcoded in ROM and the multiplications and divisions are by 2 instead of 10. Here they are generated by a program (XEQ F). If the calculator is set to DEGREES when F is executed, the arguments will have to be in degrees befor executing A.

Quote:
but if other users keep finding 'em, I'll analyze 'em!

Keep up the good work, is all I have to say. Let's hope HP read 'em! :-)

Regards,

Gerson.

Edited: 10 Dec 2006, 8:27 p.m.


#23

Gerson --

Thanks for the program and compliment. I surmise that someone forwarded the previous discussions of HP-33S bugs in the Forum to KinHPo, because they did indeed get addressed.

-- KS

#24

Quote:
I actually hope that you're right, in that more-direct methods of calculation would be employed today

Your phrasing almost makes it sound like you think there's something fundamentally wrong with CORDIC, but in point of fact there's no compelling reason NOT to use CORDIC in a calculator. It is a very efficient means of computing various elementary functions. Use of other methods, such as power series, rational approximations, etc. are not particularly better; any of these methods can be used to calculate results to any desired precision and accuracy. Depending on the hardware available, one method might be faster than another. As with any kind of engineering, there are tradeoffs involved.

The thing that must be understood by the engineers is that regardless of which method is chosen, an error analysis must be used to determine the number of terms and the precision to which the calculations must be performed.

Many of the known bugs/inaccuracies in calculator elementary functions are due to not doing such an error analysis, or perhaps deliberately ignoring it to increase performance. (Using fewer terms or lower precision will get worse results faster.)

Quote:
Any microprocessor for a handheld calcuator that didn't have built-in commands for multiplication and division seems unfathomable.

It is not unfathomable when you understand that they are optimizing for low manufacturing cost, and processors with hardware multiply and divide cost more than those without. Hardware divide is just about useless for a BCD calculator anyhow, and hardware muliply is only of marginal utility.


#25

Eric --

That was a good, thoughtful response. I admit that I am not very knowledgeable about computer hardware or the specific low-level computational methods -- it's not my area.

As long as a method such as CORDIC returns accurate answers quickly, and is well-suited to the microprocessor and format of data, then it is perfectly acceptable, and may even be the best available.

It just seems odd to me that methods developed a long time ago for far-less-capable hardware, and which I find arcane and almost-inscrutable with my -- well, respectable --background in mathematics, may still be the de facto standard.

Regards,

-- KS

Edited: 14 Dec 2006, 11:53 p.m. after one or more responses were posted


#26

The reason CORDIC seems inscrutable is that many places where the CORDIC algorithms are explained actually do a fairly poor job of it. I'll try to summarize here, but the Wikipedia article is better.

CORDIC operates in two modes, called rotation and vectoring. Rotation mode is used to compute the sine and cosine of an angle, and vectoring mode is used to compute an arctangent or do a rectangular to polar conversion.

Suppose we want to use rotation mode to compute sin(theta) and cos( theta) in binary arithmetic. (The calculator does this in base 10 rather than straight binary, but the fundamentals of the algorithm are the same.) We start with a unit vector on the x axis, (1, 0). We'll call the components of the vector x0 and y0, and the angle, which is zero, theta0.

We'll perform a series of rotations of the vector by decreasing angles, until thetan is approximately equal to theta. At that point the vector components xn and yn are approximations of cos(theta) and sin(theta). Obviously the closer thetan approaches theta, the better the results will approximate the true values of the sine and cosine of theta.

The usual way of performing 2D vector rotation involves computing the sine and cosine of the rotation angle, and doing a matrix multiply. In this case, obviously we can't do that, since the sine and cosine functions are what we're trying to implement. Instead, we do rotations only by a predetermined series of angles rthetan, chosen such that if we factor a cos(thetan) out of the rotation matrix, the elements on the main diagonal of the matrix are 1, and the other elements are of the form +/- 2-n. Note that binary multiplication by an integral power of two, such as -n, is simply shifting. Thus the matrix multiply can be performed by only shift and add, and we multiply the results by 1/cos(rthetan).

Since the rthetan is a constant for each n, we can defer this multiplication until all n iterations have completed, then multiply by the reciprocal of the product of the cosines of the rtheta values.

In vectoring mode, we start with a vector at an unknown angle, and theta0 of zero. We apply the same series of rotations, choosing the rotation directions to make the vector approach the unit vector on the x axis, (1, 0). As we do the rotations, we keep track of thetan, such that at the conclusion thetan is an approximation of theta, the value we were seeking.

(error corrected on 15 DEC 2006)

Edited: 15 Dec 2006, 2:00 p.m. after one or more responses were posted


#27

Eric --

Thanks for taking the time to prepare the quick summary of CORDIC. There are indeed many detailed documents available through the Wikipedia links. I just might dabble in these...

Regards,

-- KS

#28

Hello, why shouldn't CORDIC be used in modern calculators ?

Historically there were two methods used : polynomial approximations, and CORDIC. Both scale to higher precisions, but IMHO the latter is much faster. Recall that originally Volder worked on airplane radar, and CORDIC made it just simply feasible, beating polynomial approximations which exist since the 17th century at least.

Also, CORDIC allows you to get trigonometrics, hyperbolics, inverse, square root, and others, in basically one same loop construct. I wouldn't be surprised if only one code were producing all scientific functions on tight code implementations like the HP35.

Could anyone in the know tell us what other methods exist ? I kno NO efficient polynomial approximation.

Also, on BCD calculators all shifting is by powers of ten, very easy and precise.


#29

Hello,

Quote:
Historically there were two methods used : polynomial approximations, and CORDIC. Both scale to higher precisions, but IMHO the latter is much faster.

From wikipedia:

"CORDIC is generally faster than other approaches when a hardware multiplier is unavailable (e.g. in a microcontroller), or when the number of gates required to implement one needs to be minimized (e.g. in an FPGA). On the other hand, when a hardware multiplier is available (e.g. in a DSP microprocessor), table-lookup methods and power series are generally faster than CORDIC.

So CORDIC is the fastest solution even for the HP-33S (Sunplus microcontroller SPLB31A), a modern calculator, and you are right!

Quote:
Could anyone in the know tell us what other methods exist ? I kno NO efficient polynomial approximation.

I am not in the know but I know other methods do exist, suitable for certain applications, like some described here:

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

Regards,

Gerson.

---------------


This is problably known to you:

http://www.jacques-laporte.org/LeSecretDesAlgorithmes.htm

Very interesting indeed, this and the other articles therein!


Edited: 11 Dec 2006, 7:56 p.m.

#30

Here are additional (or more-severe) manifestions of previously-identified HP-33S bugs.


TRIGONOMETRIC BUG:

109 * cos(89.9999999 degrees) 

HP-32SII 1.74532925199
HP-33S 1.74532000000
actual 1.74532925199

With six significant digits lost by the HP-33S! Every calculator of mine with a 10-digit display (except the HP-35 from 1971) will return the correct value of 1.745329252 -- even the shoddy "Le World" and Sharp EL-520R.


COMBINATION/PERMUTATION BUG:

This bug, in which results of the combination function that exceeded 10500 were displayed instead as values within the displayable range, also extended to the permutation function.

Perm(500,250)

HP-32SII Overflow
HP-33S 5.92335407928 E 499
actual 3.77417592222 E 641

Can anyone with a 2005 HP-33S confirm that this bug was fixed?

Thanks,

-- KS



Edited: 11 Dec 2006, 2:51 a.m.


#31

I have an HP33S CNA6080xxxx.

500 ENTER 250 leftshift nPr makes the calc work a bit (about a second or two), gives the OVERFLOW notice, and leaves 1e500 in the x register.

I gather this is the desired behaviour.

My unit gives the cos(89.9999999) bug you mention.

This is really disappointing! In so many ways I really like this calculator....

Les

#32

It is even more interesting:

Model     mantissa

89.9->
32sii: 174532836590
33s: 174532836590

89.99->
32sii: 174532924313
33s: 174532924306

89.999->
32sii: 174532925191
33s: 174532925091

89.9999->
32sii 174532925199
33s 174532925000

89.999 99 (five nines)->
32sii: 174532925199
33s: 174532920000

89.999 999 (6 nines)->
32sii: 174532925199
33s: 174532900000

89.999 999 9 (seven nines)->
32sii: 174532925199
33s: 174532000000

89.999 999 99 (8 nines)->
32sii: 174532925199
33s: 174532925199

89.999 999 999 (9 nines)->
32sii: 174532925199
33s: 174532925199

89.999 999 999 9 (ten nines)->
32sii: 147532925199
33s: 174532925199


Note that it returns the correct result at extremely high significant figures for the angle approuaching the limit of 90.

#33

Quote:
Can anyone with a 2005 HP-33S confirm that this bug was fixed?

Mine from 2005 (CNA 51500163) still has all the bugs! Keyboard and display are ok but the ROM seems to be the original one.

Marcus


#34

My brand new hp33s has all the same bugs! This is very disappointing !
Are there patches/downloads available?
Are any of these bugs in the more expensive hp graphing calculators also?


#35

Steve --

What is the serial # of your HP-33S? (First 6 characters will do.)

BTW, there are no patches/downloads available for fixing the "buggy" ROM, because there is no way to upload them into the calculator.

I seriously doubt that the same bugs are present on the HP-50G, HP-49G+, and HP-48GII, because they use different microprocessors, for which different code was developed.

The HP-33S' errors were introduced during the "porting" from the predecessor HP-32SII, which used the legacy Saturn mircoprocessor, originally developed in 1983.

-- KS

Edited: 16 Dec 2006, 1:49 p.m.

#36

My CNA532 displays OVERFLOW then 1 E500 for the above permutation.


Possibly Related Threads...
Thread Author Replies Views Last Post
  Trigonometric Functions (HP-17BII) Gerson W. Barbosa 2 414 04-09-2013, 11:57 PM
Last Post: Gerson W. Barbosa
  Degrees conversion for trigonometric function Beth Snipes 2 275 08-30-2011, 03:32 PM
Last Post: Beth Snipes
  CORDIC Trigonometric Algorithm on an HP17BII Charles 4 365 03-14-2008, 12:14 AM
Last Post: Eric Smith
  The trigonometric bug is spreading !? Lyuka 33 1,707 09-07-2007, 04:57 AM
Last Post: DaveJ
  HP 35s trigonometric functions alternative Lyuka 5 390 08-30-2007, 10:01 PM
Last Post: Gerson W. Barbosa
  Fast and Accurate Trigonometric Functions on the HP-17BII Gerson W. Barbosa 2 274 01-14-2007, 02:37 AM
Last Post: Gerson W. Barbosa
  Even Faster and More Accurate Trigonometric Functions on the HP-12C Platinum (Long) Gerson W. Barbosa 8 478 08-02-2006, 10:50 PM
Last Post: Gerson W. Barbosa
  HP-12CP Trigonometric functions program still not tested. Anyone? Gerson W. Barbosa 8 552 02-05-2006, 08:40 AM
Last Post: Gerson W. Barbosa
  Integration Times "Old" 33s vs "New" 33s John Smitherman 21 1,326 12-14-2005, 12:04 AM
Last Post: Karl Schneider
  New Trigonometric Functions Program for the New 12C Platinum Gerson W. Barbosa 2 275 11-11-2005, 11:21 AM
Last Post: Gerson W. Barbosa

Forum Jump: