Native mode of trig functions?

[naddy]

The talk of BCD vs. binary in calculators has made me think of another question:

What is the native mode for calculating trigonometric functions on, say, the HP-11C/15C?

Are degrees first converted to radians?
Or vice versa?
Or are there even separate implementations?

(This is just idle curiosity.)

[Johnh]

Its a good question. Id expect that radians would be native, since the usual algorithms for calculating trig functions from a series work directly with radians

[naddy]

I feel like I should be able to tell from the Appendix "Accuracy of Numerical Calculations" in the HP-15C Advanced Functions Handbook—but I can't.

Does the fact that the trig functions are more accurate for real arguments in degrees (and grads) mode than in radians mode imply that the modes use separate algorithms? Or is it due to internal d>r conversion using 13 digit accuracy?

[brouhaha]

TL;DR: the trigonometric and inverse trigonometric functions in the HP-15C, and other HP calculators up through the Pioneer series, work natively in radians.

The HP-15C microcode for trig and inverse trig is very close to that of the 11C, which was very close to that of the 41C, which was very close to that of the Spice series (3xE and 3xC series). If you trace it back further, there are some changes, as the HP engineers made significant improvements to the code in the six years between the HP 35 and the 30 series. I'll describe the algorithm specifically with regard to the HP-15C implementation, based on my own reverse-engineering work. Any mistakes in my description are obviously mine.

The separate entry points for sin, cos, and tan set distinct combinations of the s10 and s12 processor flags to track which function was requested. All three paths then merge together. A subroutine is called to read the angle mode from one of the internal status registers, and set up flags s4 and s5.

Argument range reduction is done first, in the current angle mode setting (based on s4 and s5). Grads are dealt with by converting them into degrees. Beyond this point, grads no longer have to be specifically handled. The path forks to get the appropriate constant for range reduction in degree or radian mode, then recombines to use restoring division by repeated subtraction, and appropriate shifting and exponent update. The remainder after the division is the reduced angle.

After range reduction, if the mode is degrees (or was in grads), then degrees are converted to radians.

Once the angle is in radians, a CORDIC rotation is performed to simultaneously computer the sin and cos, by rotating a unit vector (1, 0) through the angle. The resulting x and y are proportionate to the sin and cos, though they have been multiplied by a non-constant factor (because the decimal CORDIC rotation multiplies the vector length by a non-constant scale factor). However this non-constant factor is identical for the x and y components. This allows the tangent to be computed by a simple division.

Once the tangent has been computed, if sin or cos were requested (as recorded by s10 and s12), they are derived by trig identities from the tangent.

David S. Cochran, the microcode programmer of the HP-35, described the trig algorithm in "Algorithms and Accuracy in the HP-35", in Hewlett Packard Journal Vol. 23 No. 10, June 1972. He references the shift-and-add logarithm algorithm devised by Heny Briggs in 1624, and CORDIC method by Volder (1956, published 1959), and an equivalent formulation referred to a psuedo-division and pseudo-multiplication by Meggitt (1962). Cochran explains the algorithm in terms of Meggitt's pseudo-division and pseudo-multiplication, which I personally find to be somewhat more confusing terminology than that normally used to explain CORDIC. If you want to understand how the CORDIC algorithm works, I recommend the Wikipedia CORDIC article. Most descriptions of CORDIC focus on base 2 implementation, but the base 10 implementation is not much different.

Also see Jacques Laporte's explanation of how the HP-35 implemented trigonometric functions.

[brouhaha]

(10-10-2024 11:37 PM)naddy Wrote:  Does the fact that the trig functions are more accurate for real arguments in degrees (and grads) mode than in radians mode imply that the modes use separate algorithms? Or is it due to internal d>r conversion using 13 digit accuracy?

If you're talking about angles outside +/-90 degrees (+/- pi/2 radians), the answer is that the argument range reduction, performed before any conversion of degrees to radians, is inherently inaccurate because of the limited accuracy of the internal constant for pi/4 used in the range reduction for the radians case. Range reduction for degrees and grads are "exact", at least in so far as the constants used for the reduction can be represented exactly.

Within +/- pi/2 radians, the trig functions should be very slightly more accurate in radians than in degrees or grads.

[Thomas Klemm]

(10-11-2024 12:31 AM)brouhaha Wrote:  Once the angle is in radians, a CORDIC rotation is performed to simultaneously computer the sin and cos, by rotating a unit vector (1, 0) through the angle.

For those interested in doing this by themselves:
Exploring the CORDIC algorithm with the WP-34S

[Thomas Klemm]

(10-11-2024 12:38 AM)brouhaha Wrote:  Within +/- pi/2 radians, the trig functions should be very slightly more accurate in radians than in degrees or grads.

It is probably difficult to define what that means since the same number can't be entered exactly in both modes DEG and RAD.
Therefore you have to convert the entry one or the other anyway.

Example

\(
\begin{align}
\cos(36^{\circ})
&= \frac{1+\sqrt{5}}{4} \\
&\approx 0.80901699437\cdots \\
\end{align}
\)

On an HP-41C I get the following results:

DEG
36
COS

0.8090169944

RAD
36
D-R
COS

0.8090169944

5
SQRT
1
+
4
/

0.8090169943

Both results using COS are rounded correctly while the "exact" value using SQRT is not.

[Werner]

(10-11-2024 11:57 AM)Thomas Klemm Wrote:  5
SQRT
1
+
4
/

0.8090169943
[/font]
Both results using COS are rounded correctly while the "exact" value using SQRT is not.

No, because you have two rounding operations introducing errors.
If you calculate it as SQRT(5/16) + 0.25, you only have one, and get 0.8090169944 as well.

Cheers, Werner

[dm319]

(10-11-2024 12:31 AM)brouhaha Wrote:  I'll describe the algorithm specifically with regard to the HP-15C implementation, based on my own reverse-engineering work.

Very impressed you can figure that out from the ROM!

If I ever get that ROM out of the Rogue One I know who to ask.

PS Do you happen to know how the TVM functions are calculated on the 12c? I'm curious as to the formula and method used for the solve-for-i.

[AnnoyedOne]

(10-11-2024 01:56 PM)dm319 Wrote:  If I ever get that ROM out of the Rogue One I know who to ask.

Good luck with that! And if it is a NUT emulation and not just native 6502 code.

A1

[Albert Chan]

(10-11-2024 11:57 AM)Thomas Klemm Wrote:  

(10-11-2024 12:38 AM)brouhaha Wrote:  Within +/- pi/2 radians, the trig functions should be very slightly more accurate in radians than in degrees or grads.

It is probably difficult to define what that means since the same number can't be entered exactly in both modes DEG and RAD.

We define "more accurate" in the statistical sense.
If there is no angle reduction, DEG mode need a conversion back to RAD, thus suffered an extra rounding.

(10-11-2024 12:59 PM)Werner Wrote:  No, because you have two rounding operations introducing errors.
If you calculate it as SQRT(5/16) + 0.25, you only have one, and get 0.8090169944 as well.

It has less to do with two roundings, more to be with size of ulp

error(SQRT(5)) < 0.5 ULP of 1 = 5 ULP of 0.1

Even if all other operations are exact, error((SQRT(5)+1)/4) < 5/4 = 1.25 ULP of 0.1

SQRT(5) = 2.2360 67977 49978 96964 ...

For 10 decimal digits, SQRT(5) is almost a half-way case, above '<' turns to '≈'
It is this 1.25 ULP error that make phi result not accurate.

If we expand it out, SQRT(5/16) + 0.25, error from SQRT < 0.5 ULP of 0.1
Since sum still below 1, whole expression also < 0.5 ULP of 0.1, i.e. correctly rounded.

[Thomas Klemm]

(10-11-2024 04:30 PM)Albert Chan Wrote:  We define "more accurate" in the statistical sense.

Can you elaborate in that?

My point is that you can't enter \(\frac{\pi}{5}\) into a calculator.
So we always start with an approximation.
But in case of \(36^{\circ}\) we can start with the exact value.

Interestingly I get the following result on an HP-41C:

RAD
PI
5
/
COS

0.8090169943

Which is not the correctly rounded result of \(\cos\left(\frac{\pi}{5}\right)\).

[Albert Chan]

(10-11-2024 05:29 PM)Thomas Klemm Wrote:  Interestingly I get the following result on an HP-41C:

RAD
PI
5
/
COS

0.8090169943

Which is not the correctly rounded result of \(\cos\left(\frac{\pi}{5}\right)\).

Accuracy is only in statistical sense, assuming user input is correct.
What is user input of a function? Value on stack before we execute it.

Above 10-digits calculator example, value = 3.141592654 / 5 = 0.6283185308

COS(0.6283185308) = 0.8090169943 2672472713 ...

HP-41C COS result are correctly rounded, the best we can ask for.
We cannot input exact value of pi/5, but this does not count against accuracy of COS.

[J-F Garnier]

(10-11-2024 07:37 PM)Albert Chan Wrote:  COS(0.6283185308) = 0.8090169943 2672472713 ...
HP-41C COS result are correctly rounded, the best we can ask for.

The internal 13-digit value of the HP-41C COS result is actually 0.8090169943 276
then (here correctly) rounded to the user 10-digit result = 0.8090169943
However, since the 3 guard digits are far from being all correct, the rounding will be wrong by +/-1 in the last place from time to time, depending on the function argument.

J-F

[brouhaha]

(10-11-2024 01:56 PM)dm319 Wrote:  

(10-11-2024 12:31 AM)brouhaha Wrote:  I'll describe the algorithm specifically with regard to the HP-15C implementation, based on my own reverse-engineering work.

Very impressed you can figure that out from the ROM!

I've spent a LOT of time working through disassembly listings, originally just to turn the disassemblies of various calculator ROMs back into an approximation of source code, but later to actually develop an understanding of how various parts of it actually work.

At one time HP actually asked ME for the HP 12C source code, but I've actually never had the actual source code for any of the Voyager or Spice models, nor for most of the pre-Saturn models. Fortunately HP did not long after that manage to recover the 12C source code, though one of their engineers has said that the source code for the Voyager (and earlier) calculators is basically unusable. My interpretation is that they don't want to hack on Nut assembly code any more than absolutely necessary for simulation in the ARM-based Voyager hardware. I don't think that's an unreasonable position; only some of us Nut-jobs are willing to do it. :-)

Quote:If I ever get that ROM out of the Rogue One I know who to ask.

It's probably written in C compiled to 6502 (or 6502-like) machine code, and might be more difficult to extract useful knowledge from than Nut code, but I'd definitely be interested to look it over. I have reverse-engineered a lot of 6502 code over the years, but none that was compiled from a high-level language like C. The 6502 is rather ill-suited to C, so I expect the result to be rather convoluted, if the compiler is actually any good. If the compiler doesn't optimize well, it might contain a huge amount of boilerplate, which might actually make reverse-engineering easier.

Quote:PS Do you happen to know how the TVM functions are calculated on the 12c? I'm curious as to the formula and method used for the solve-for-i.

As are we all! Of course, that's the 12C's "secret sauce", so while HP has described aspects of the algorithm in general terms, they have not published any of the details.

While I have reverse-engineered some of the 12C code, I have not had time to really investigate the actual interest computation, which of course is the most complicated code it contains. I think it may be easier to do such reverse-engineering initially from the 37E. As far as I'm aware, the 12C algorithm is very close to the one the 37E, 38E, and 38C use, and the 37E only has 2.5 Kwords of instructions to decipher, vs. 6 Kwords in the 12C. Also, the first 1K of the 37E ROM, in the 1820-2122 Spice CPU chip, is shared with the 38E, and the next 0.5K is nearly identical between the two.

[dm319]

(10-12-2024 02:58 AM)brouhaha Wrote:  At one time HP actually asked ME for the HP 12C source code

That's hilarious!! Also worrying...

(10-12-2024 02:58 AM)brouhaha Wrote:  The 6502 is rather ill-suited to C, so I expect the result to be rather convoluted, if the compiler is actually any good. If the compiler doesn't optimize well, it might contain a huge amount of boilerplate, which might actually make reverse-engineering easier.

AnnoyedOne may have identified the IC - but looks like it might have some sort of decimal/nibble co-processor on it, which I suppose could make understanding it even harder.

(10-12-2024 02:58 AM)brouhaha Wrote:  As are we all! Of course, that's the 12C's "secret sauce", so while HP has described aspects of the algorithm in general terms, they have not published any of the details.

Seems so silly for everyone to have been re-inventing the wheel all this time.

(10-12-2024 02:58 AM)brouhaha Wrote:  While I have reverse-engineered some of the 12C code, I have not had time to really investigate the actual interest computation, which of course is the most complicated code it contains. I think it may be easier to do such reverse-engineering initially from the 37E.

I think they are similar:

Code:

|calculator                  |1b   |4    |5A   |6A   |7    |8    |11   |12   |
|:---------------------------|:----|:----|:----|:----|:----|:----|:----|:----|
|HP-12c Platinum             |11   |11.1 |ERR  |ERR  |21.4 |12.2 |11   |2.8  |
|HP-19BII                    |10.7 |11.1 |ERR  |ERR  |11.7 |12.2 |NA   |NA   |
|HP-27S                      |10.7 |11.1 |ERR  |ERR  |NA   |12.2 |10.9 |2.8  |
|HP-17B                      |10.7 |11.1 |ERR  |ERR  |11.9 |12.2 |10.9 |2.8  |
|HP-92                       |8.5  |Inf  |ERR  |ERR  |9.7  |11.8 |4.3  |1.2  |
|HP-12c                      |8.6  |9.1  |ERR  |ERR  |19.7 |11.5 |4.6  |1    |
|HP-37E                      |8.5  |9    |ERR  |ERR  |9.7  |11.8 |6    |1    |
|HP-38E                      |8.5  |9    |ERR  |ERR  |9.7  |11.8 |6    |1    |
|HP-12c Rogue One            |12.4 |ERR  |12.2 |13.2 |ERR  |12.2 |ERR  |ERR  |
|HP-41CX Financial Decisions |6.4  |Inf  |ERR  |ERR  |9.7  |7.6  |2    |0    |
|HP-27                       |9.5  |9    |ERR  |ERR  |ERR  |ERR  |4.3  |ERR  |
|HP-22                       |8.6  |9    |ERR  |ERR  |ERR  |ERR  |4.7  |ERR  |
|HP-70                       |6.6  |9    |ERR  |ERR  |ERR  |ERR  |-1   |ERR  |
|HP-80                       |5.1  |ERR  |ERR  |ERR  |ERR  |ERR  |-1   |ERR  |

The HP-12c and HP-12c platinum seem to do very well on problem 7:

Code:

| #  | Ref        | N            | I%YR      | PV       | PMT          | FV        | P/YR | Mode  |
|----|------------|--------------|-----------|----------|--------------|-----------|------|-------|
| 7  | A Chan     | 10           | ?         | -100     | 10           | 1e-10     | 12   | end   |

But strangely the 12c does better than what came before and what came after... But I'm sure the HP-37E would be a good start!

D