Trig algorithms on HP RPN?

[MinkLib]

Can someone tell me if the trig algorithms for HP (and other RPN calculators) are all done in radians (behind the scenes). The reason I ask is I've been trying out the Mike Sebastian Forensic test.
It is of course meant to find the underlying chip used, regardless of brand name. But it gives a (very) rough estimation of accuracy, and the results I get using 0.15707 radians are startlingly better than with 9 degrees.

HP-25 degrees 9.0040 error 0.05%
HP-25 radians 8.9999 error 0.00003%

Any Soviet Elektronika programmable RPN (eg B3-34) degrees 9.0881 error 1%
Any Soviet Elektronika programmable RPN (eg B3-34) radians 9.0000 error 0.0003%

Sinclair Scientific radians 8.9381 error 0.7% (no degree mode)

Thanks.

[Idnarn]

(11-30-2024 01:39 PM)MinkLib Wrote:  Can someone tell me if the trig algorithms for HP (and other RPN calculators) are all done in radians (behind the scenes).

While the answer cannot be generalized as it may depend on the particular calculator, if you want an example, see the paper "Algorithms and Accuracy in the HP-35" by David S. Cochran (example, here: https://www.keesvandersanden.nl/calculat..._HP-35.pdf) which describes trigonometry function calculation on the HP-35. I read somewhere (was it an interview with William Kahan?) that the HP 15C code derived from it.

You can also browse Jacques Laporte's excellent reverse engineering of the HP-35 here: https://archived.hpcalc.org/laporte/index.html

https://archived.hpcalc.org/laporte/Trigonometry.htm

[naddy]

See this thread where I asked the same seven weeks ago:
https://www.hpmuseum.org/forum/thread-22486.html

[Thomas Klemm]

(11-30-2024 01:39 PM)MinkLib Wrote:  Can someone tell me if the trig algorithms for HP (and other RPN calculators) are all done in radians (behind the scenes).

Although the HP-35 only uses degrees for trigonometric functions, it still seems to use the \(\tan^{-1}\) values ​​in radians:

Code:

0-0716 0730  load constant 7        reg['C'] = 0x00700000000000
0-0717 1030  load constant 8        reg['C'] = 0x00780000000000
0-0720 0530  load constant 5        reg['C'] = 0x00785000000000
0-0721 0330  load constant 3        reg['C'] = 0x00785300000000
0-0722 1130  load constant 9        reg['C'] = 0x00785390000000
0-0723 1030  load constant 8        reg['C'] = 0x00785398000000
0-0724 0130  load constant 1        reg['C'] = 0x00785398100000
0-0725 0630  load constant 6        reg['C'] = 0x00785398160000
0-0726 0330  load constant 3        reg['C'] = 0x00785398163000
0-0727 0530  load constant 5        reg['C'] = 0x00785398163500

And then a bit later:

Code:

0-0617 1030  load constant 8        reg['C'] = 0x00099668666666
0-0620 0630  load constant 6        reg['C'] = 0x00099668666666
0-0621 0530  load constant 5        reg['C'] = 0x00099668656666
0-0622 0230  load constant 2        reg['C'] = 0x00099668652666
0-0623 0430  load constant 4        reg['C'] = 0x00099668652466
0-0624 1130  load constant 9        reg['C'] = 0x00099668652496

This is output from running x11-calc-35 with the -t option:

Code:

  -t,                      trace

Compare this to (using RAD):

1 ATAN

0.78539816340

0.1 ATAN

0.09966865249

From my understanding of CORDIC we might as well use angles in degrees.
Honestly, I don't know why that wasn't used instead.

[Gerson W. Barbosa]

(11-30-2024 03:49 PM)Thomas Klemm Wrote:  Honestly, I don't know why that wasn't used instead.

Probably to take advantage of the fact that arctan(x) ~ ln(1 + x) for small values of x, which allows part of the arctangent table to be used in the computation of the natural logarithm.

For example, arctan(0.00001) ~ ln (1.00001). That’s the fifth element in the arctangent table starting at register 1, generated by the HP-15C program B below, when in RAD mode. Program A computes tan(x) using the CORDIC algorithm.

Code:


   001 { 42 21 11 } f LBL A
   002 {    44 .5 } STO 15
   003 {        1 } 1
   004 {    44 .7 } STO 17
   005 {    44 .4 } STO 14
   006 {    44 25 } STO I
   007 {        0 } 0
   008 {    44 .8 } STO 18
   009 { 42 21  0 } f LBL 0
   010 {    45 .5 } RCL 15
   011 {        7 } 7
   012 {       26 } EEX
   013 {       16 } CHS
   014 {        1 } 1
   015 {        1 } 1
   016 { 43 30  7 } g TEST 7
   017 {    22  3 } GTO 3
   018 { 42 21  1 } f LBL 1
   019 {    45 24 } RCL (i)
   020 {    45 .5 } RCL 15
   021 { 43 30  9 } g TEST 9
   022 {    22  2 } GTO 2
   023 {        1 } 1
   024 { 44 40 25 } STO + I
   025 {        1 } 1
   026 {        0 } 0
   027 { 44 10 .4 } STO / 14
   028 {    22  1 } GTO 1
   029 { 42 21  2 } f LBL 2
   030 {    45 24 } RCL (i)
   031 { 44 30 .5 } STO - 15
   032 {    45 .7 } RCL 17
   033 {    44 .6 } STO 16
   034 {    45 .4 } RCL 14
   035 { 45 20 .8 } RCL * 18
   036 { 44 30 .7 } STO - 17
   037 {    45 .4 } RCL 14
   038 { 45 20 .6 } RCL * 16
   039 { 44 40 .8 } STO + 18
   040 {    22  0 } GTO 0
   041 { 42 21  3 } f LBL 3
   042 {    45 .8 } RCL 18
   043 { 45 10 .7 } RCL / 17
   044 {    43 32 } g RTN
   045 { 42 21 12 } f LBL B
   046 {        1 } 1
   047 {       48 } .
   048 {        0 } 0
   049 {        1 } 1
   050 {        3 } 3
   051 {    44 25 } STO I
   052 { 42 21  4 } f LBL 4
   053 {    45 25 } RCL I
   054 {    43 44 } g INT
   055 {        1 } 1
   056 {       30 } -
   057 {       16 } CHS
   058 {       13 } 10^x
   059 {    43 25 } g TAN-1
   060 {    44 24 } STO (i)
   061 { 42  6 25 } f ISG I
   062 {    22  4 } GTO 4
   063 {    43 32 } g RTN

I cannot find a reference to the original TI program this one is based upon.

[Albert Chan]

(11-30-2024 03:49 PM)Thomas Klemm Wrote:  

Code:

0-0727 0530  load constant 5        reg['C'] = 0x00785398163500

Should last digit for pi/4 a 4?

lua> pi/4
0.7853981633974483

(11-30-2024 03:49 PM)Thomas Klemm Wrote:  Honestly, I don't know why that wasn't used instead.

Just a guess, if angles are in radian, we can just use tan code to get atan

tan sum formula: tan(x+y) = (tan(x) + tan(y)) / (1 - tan(x)*tan(y))
Flip for atan:          atan(x) = atan( y + (x-y) )

atan(x) = atan(y) + atan(z),     where z = (x-y)/(1+x*y)

if y = tan(x), (x-y) ≈ x^3/3 --> z shrink fast, cubic convergence!

But we can do better!

y can be set to anything. As long as y is close to x, z will shrink fast.
Below, y = tan(x2 = x/(1+x*x/3)), to get convergence rate of O(x^5)

Code:

function myatan(x, s)
    if abs(x)>1 then return myatan(-1/x, copysign(pi/2,x)) end
    s = s or 0
    print(s .. " + atan( " .. x .. " )")
    if 3+x*x==3 then return s + x end   -- atan(x) ≈ x
    local x2 = x / (1+x*x/3)
    local y = tan(x2)           -- y is close to x
    local z = (x-y) / (1+x*y)   -- atan(x) = atan(y) + atan(z)
    return myatan(z, s + x2)    -- atan(y) = x2
end

lua> myatan(0.1)
0 + atan( 0.1 )
0.09966777408637874 + atan( 8.784047832967151e-07 )
0.09966865249116204 + atan( -1.0587911840670585e-22 )
0.09966865249116204

lua> myatan(1)
0 + atan( 1 )
0.75 + atan( 0.03541295579818369 )
0.7853981584542247 + atan( 4.943223638603625e-09 )
0.7853981633974483

[Albert Chan]

(12-01-2024 05:17 PM)Gerson W. Barbosa Wrote:  For example, arctan(0.00001) ~ ln (1.00001). That’s the fifth element in the arctangent table starting at register 1, generated by the HP-15C program B below, when in RAD mode. Program A computes tan(x) using the CORDIC algorithm.

Can you explain in words (better, example) how tan and log1p is used to get atan?

[Thomas Klemm]

Your program illustrates my point very nicely: Run program B in DEG mode to use program A to calculate \(\tan(x)\) in degrees.

But then we end up with constants like these:

45.00000000
5.710593137
0.572938698
0.057295760
0.005729578
0.000572958
…

And ROM was a scarce resource.

---

(12-01-2024 05:17 PM)Gerson W. Barbosa Wrote:  I cannot find a reference to the original TI program this one is based upon.

Initially I considered the HP-15C for this article: Exploring the CORDIC algorithm with the WP-34S
But I'm not sure if that would have worked as well because you can't store complex numbers in registers.
And the WP-34C provides INC and decimal shift which turned out to be useful.

But you can still follow the example on the HP-15C manually:

1 ENTER 0.1 I
ENTER
ENTER
3 ENTER 4 I

3.000000000

×

2.600000000

×

2.170000000

×

1.714000000

×

1.236300000

×

0.741460000

×

0.234257000

×

-0.280360600

[Thomas Klemm]

(12-01-2024 07:27 PM)Albert Chan Wrote:  Just a guess, if angles are in radian, we can just use tan code to get atan.

This is already the case.
From A Unified Algorithm for Elementary Functions:

J. S. WALTHER Wrote:This paper describes a single unified algorithm for the calculation of elementary functions including multiplication, division, sin, cos, tan, arctan, sinh, cosh, tanh, arctanh, ln, exp and square-root.

[Gerson W. Barbosa]

(12-01-2024 07:30 PM)Albert Chan Wrote:  

(12-01-2024 05:17 PM)Gerson W. Barbosa Wrote:  For example, arctan(0.00001) ~ ln (1.00001). That’s the fifth element in the arctangent table starting at register 1, generated by the HP-15C program B below, when in RAD mode. Program A computes tan(x) using the CORDIC algorithm.

Can you explain in words (better, example) how tan and log1p is used to get atan?

Actually, atan and log1p lookup tables are used to get tan and log (and exp) functions (please see https://archived.hpcalc.org/laporte/TheS...rithms.htm). Sorry if my text was not clear enough.

.x arctan(1/10^x) ln(1+1/10^x)
00 0.785398163397 0.693147180560
01 0.099668652491 0.095310179804
02 0.009999666687 0.009950330853
03 0.000999999667 0.000999500333
04 0.000100000000 0.000099995000
05 0.000010000000 0.000009999950
06 0.000001000000 0.000000999999
07 0.000000100000 0.000000100000
08 0.000000010000 0.000000010000
09 0.000000001000 0.000000001000
10 0.000000000100 0.000000000100
11 0.000000000010 0.000000000010
12 0.000000000001 0.000000000001

As Thomas Klemm has correctly pointed out, tan(x) would accept arguments in degrees if the atan table elements were in degrees rather than radians. It appears they chose radians in order do save memory, as thus the second half of the elements would be common to both tables.