Arctangent function program [HP-42S, Free42]

[Gerson W. Barbosa]

(09-17-2024 01:56 AM)Thomas Klemm Wrote:  Minor typo in:

(09-16-2024 05:10 PM)Gerson W. Barbosa Wrote:  

Code:

   002 { 43 33 25 } GTO 2

Should be:

Code:

   002 { 43 33 25 } GTO 25

Thanks! At least the key codes are correct. Anyway, I’ve found out that more than 9 iterations are needed. Apparently no accuracy improvement beyond 11 or 12. Better accuracy should be expected on the 12C Platinum, but I don’t have 30 steps free in any of mine.

Typo correction and patch:

Code:

   002 { 43 33 25 } g GTO 25
   ...
   006 {        1 } 1
   007 {        1 } 1
   008 {    44  0 } STO 0
   009 {    43 24 } g FRAC

Edited to fix a typo

[Albert Chan]

(09-17-2024 04:23 AM)Gerson W. Barbosa Wrote:  Anyway, I’ve found out that [for tan(x)] more than 9 iterations are needed.
Apparently no accuracy improvement beyond 11 or 12.

It depends on domain. see https://demonstrations.wolfram.com/Conti...ntFunction
If tan argument limited to ±pi/4, 8 iterations already have 16+ digits accuracy.

Instead of getting CF result bottom-up, with fixed iterations, we can also do this top-down.
see https://www.ams.org/journals/mcom/1952-0...9650-3.pdf, Method III

Note: there is a typo in formula, it should be 1+p(i) = 1 / (1 + r(i)*(1+p(i-1)))
This is why I use variable p1 = 1+p instead.

Code:

function my_tan(x)
    local a, b = x, 1    -- CF = x/(1 - x^2/(3 - x^2(5 - x^2(7 - ...
    local t = a/b
    local s = t  
    local p1 = 1        -- = p+1
    for k=1,inf do      -- other CF terms
        a, b, b0 = -x^2, b+2, b
        p1 = b0 / (b0 + p1*a/b)
        t = t * (p1-1)
        s = s + t
        print(k, t, s)
        if s==s+t then return s end
    end
end

lua> my_tan(pi/4)
1      0.2032910765367568      0.9886892399342051
2      0.011098440980763355    0.9997876809149685
3      0.00021018750072653323  0.999997868415695
4      2.1181106602376684e-06  0.9999999865263552
5      1.341497025916972e-08   0.9999999999413255
6      5.84877225566155e-11    0.9999999999998133
7      1.8641515888296674e-13  0.9999999999999997
8      4.534798657058108e-16   1.0000000000000002
9      8.698065146274166e-19   1.0000000000000002
1.0000000000000002

[Albert Chan]

tan(x) continued fraction formula is good, but sin(x) taylor series is better. (*)

Sine x^(2k+1) / (2k+1)! term shrink very fast!
sin(x) argument to within ±pi/4, 7 iterations get 16+ digits accuracy.
Even if domain expanded to ±pi/2, only 3 more iterations are needed.

Here, I convert sin taylor series to continued fraction, just to check previous post CF code.



Euler's continued fraction formula proof is in code comment!
We could do taylor term directly, t = t * -a , bypass continued fraction calculations.

Code:

function my_sin(x)     
    local a, b = x, 1   -- sin(x) = x*(1 - x^2/6*(1 - x^2/20*(1 - x^2/42*(1 - ...
    local t = a/b
    local s = t  
    local p1 = 1        -- = p+1
    for k=1,inf do      -- other CF terms
        a = x*x/(2*k*(2*k+1))
        b, b0 = 1-a, b          -- p1 = b0
        p1 = b0/(b0 + p1*a/b)   -- p1 = b0/(b0 + b0*a/b) = b/(b+a) = b
        t = t * (p1-1)          -- p1-1 = b-1 = b-(a+b) = -a
        s = s + t
        print(k, t, s)
        if s==s+t then return s end
    end
end

lua> my_sin(pi/4)
1    -0.08074551218828074      0.7046526512091675
2    0.002490394570192714      0.7071430457793603
3    -3.657620418217711e-05    0.7071064695751781
4    3.1336168903781285e-07    0.7071067829368671
5    -1.7572476734434049e-09   0.7071067811796194
6    6.948453273886762e-12     0.7071067811865679
7    -2.0410263396642417e-14   0.7071067811865475
8    4.628704628834906e-17     0.7071067811865475
0.7071067811865475
lua> _ / sqrt(1 - _*_) -- = tan(pi/4)
0.9999999999999999


(*) Correction, tan vs sin convergence rate probably too close to call.
But tan(x) has bigger range, getting sin/cos from it should be more accurate.
Bonus, we can map all trig. functions with Tangent half-angle formula, without square root.

[Thomas Klemm]

(09-17-2024 04:32 PM)Albert Chan Wrote:  Instead of getting CF result bottom-up, with fixed iterations, we can also do this top-down.

Example

We can use the program in Convergents of a Continued Fraction to calculate:

\(
\tan(0.1) = [0; 10, -30, 50, -70, 90, -110, …]
\)

0 1 0 10 XEQ 00

0.1

-30 R/S

0.1003344481605351170568561872909699

50 R/S

0.1003346720214190093708165997322624

-70 R/S

0.1003346720854403773884482176487636

90 R/S

0.100334672085450544030807774009714

-110 R/S

0.100334672085450545058008197616473


0.1 TAN

0.1003346720854505450580800457811115

[Thomas Klemm]

This program for HP-42S calculates \(\sin(x)\) and \(\cos(x)\) based on \(\tan(\frac{x}{2})\):

Code:

00 { 63-Byte Prgm }
01 2
02 X<>Y
03 ÷
04 STO 01
05 STO 02
06 8
07 STO 00
08 1
09 STO 04
10 STO 05
11 CLX
12 STO 03
13 STO 06
14▸LBL 00
15 RCL 01
16 RCL 02
17 STO 01
18 2
19 ×
20 +
21 +/-
22 STO 02
23 RCL 04
24 RCL 03
25 STO 04
26 RCL 01
27 ×
28 +
29 STO 03
30 RCL 06
31 RCL 05
32 STO 06
33 RCL 01
34 ×
35 +
36 STO 05
37 DSE 00
38 GTO 00
39 ÷
40 ENTER
41 +
42 LASTX
43 ENTER
44 ×
45 1
46 X<>Y
47 -
48 1
49 LASTX
50 +
51 ÷
52 X<>Y
53 LASTX
54 ÷
55 X<>Y
56 END

Hitting the [÷] button calculates \(\tan(x)\).

Example

1 R/S

y: 0.8414709848078965065362209309116153
x: 0.54030230586813971758203414358187

÷

x: 1.557407724654902229769750316481308

[Gerson W. Barbosa]

(09-19-2024 04:25 AM)Thomas Klemm Wrote:  Example

1 R/S

y: 0.8414709848078965065362209309116153
x: 0.54030230586813971758203414358187

÷

x: 1.557407724654902229769750316481308

With the patched HP-12C program above I get

1 ENTER 47 R/S ; RAD mode
R/S ->

HP-12C:

Z: 1.557407724
Y: 0.5403023061
X: 0.8414709848

HP-12C Platinum:

Z: 1.557407725
Y: 0.5403023059
X: 0.8414709848

I had made some tests using the tangent half-angle formula in the second half of the first quadrant. That would need six instead of eleven iterations, but I don’t think it’s worth the additional steps since current 12C calculators are extremely fast.
The program loses accuracy near critical points like 90° and 180°, so it’s meant for everyday trigs only.
Anyway, when maximum accuracy is needed this program should be used instead (HP-12C Platinum only). Because of page width change, my formatting has been completely ruined. That is the same program in the old article forum, except for the one-step reduction to make the ATAN entry-point more consistent with the ones for ACOS and ASIN.

[Gerson W. Barbosa]

Updated programs

HP-12C:

Code:


   001 {    43 35 } g x==0
   002 { 43 33 24 } GTO 24
   003 {    45  6 } RCL 6
   004 {       10 } /
   005 {    44  1 } STO 1
   006 {        7 } 7
   007 {    44  0 } STO 0
   008 {    43 24 } g FRAC
   009 {    45  0 } RCL 0
   010 {    43 35 } g x==0
   011 { 43 33 23 } g GTO 23
   012 {        2 } 2
   013 {       20 } *
   014 {        1 } 1
   015 { 44 30  0 } STO - 0
   016 {       30 } -
   017 {    45  1 } RCL 1
   018 {       10 } /
   019 {       34 } X<=>Y
   020 {       30 } -
   021 {       22 } 1/x
   022 { 43 33 09 } g GTO 09
   023 {       33 } Rv
   024 {       22 } 1/x
   025 {    43 36 } g LSTx
   026 {       30 } -
   027 {        2 } 2
   028 {       34 } X<=>Y
   029 {       10 } /
   030 {       36 } ENTER
   031 {       36 } ENTER
   032 {       36 } ENTER
   033 {       20 } *
   034 {    43 21 } g sqrt(x)
   035 {    43 35 } g x==0
   036 {    43  3 } g n!
   037 {       10 } /
   038 {    43 35 } g x==0
   039 {    43  3 } g n!
   040 {       34 } X<=>Y
   041 {       36 } ENTER
   042 {       20 } *
   043 {        1 } 1
   040 {       40 } +
   045 {    43 21 } g sqrt(x)
   046 {       22 } 1/x
   047 {       20 } *
   048 {       20 } *
   049 {    43 36 } g LSTx
   050 {       34 } X<=>Y
   051 { 43 33 00 } g GTO 00
   052 {    45  5 } RCL 5
   053 {       34 } X<=>Y
   054 {    43 35 } g x==0
   055 {       40 } +
   056 {    44  6 } STO 6
   057 { 43 33 00 } g GTO 0

Initialization:

360 ENTER 3.141592654 / STO 5

0 (DEG) or 2 (RAD) g GTO 52 R/S

Usage:

x R/S ->

T: TAN(x)
Z: TAN(x)
Y: COS(x)
X: SIN(x)

First and second quadrants only.

——

HP-15C listing for the JRPN 15C Simulator:

Code:


   001 { 42 21 11 } f LBL A
   002 {    43 20 } g x==0
   003 {    22  2 } GTO 2
   004 {    45  6 } RCL 6
   005 {       10 } /
   006 {    44  1 } STO 1
   007 {        7 } 7
   008 {    44  0 } STO 0
   009 {    42 44 } f FRAC
   010 { 42 21  0 } f LBL 0
   011 {    45  0 } RCL 0
   012 {    43 20 } g x==0
   013 {    22  1 } GTO 1
   014 {        2 } 2
   015 {       20 } *
   016 {        1 } 1
   017 { 44 30  0 } STO - 0
   018 {       30 } -
   019 {    45  1 } RCL 1
   020 {       10 } /
   021 {       34 } X<=>Y
   022 {       30 } -
   023 {       15 } 1/x
   024 {    22  0 } GTO 0
   025 { 42 21  1 } f LBL 1
   026 {       33 } Rv
   027 { 42 21  2 } f LBL 2
   028 {       15 } 1/x
   029 {    43 36 } g LSTx
   030 {       30 } -
   031 {        2 } 2
   032 {       34 } X<=>Y
   033 {       10 } /
   034 {       36 } ENTER
   035 {       36 } ENTER
   036 {       36 } ENTER
   037 {       20 } *
   038 {       11 } sqrt(x)
   039 {    43 20 } g x==0
   040 {    42  0 } f x!
   041 {       10 } /
   042 {    43 20 } g x==0
   043 {    42  0 } f x!
   044 {       34 } X<=>Y
   045 {       36 } ENTER
   046 {       20 } *
   047 {        1 } 1
   048 {       40 } +
   049 {       11 } sqrt(x)
   050 {       15 } 1/x
   051 {       20 } *
   052 {       20 } *
   053 {    43 36 } g LSTx
   054 {       34 } X<=>Y
   055 {    43 32 } g RTN
   056 { 42 21  3 } f LBL 3
   057 {    45  5 } RCL 5
   058 {       34 } X<=>Y
   059 {    43 20 } g x==0
   060 {       40 } +
   061 {    44  6 } STO 6
   062 {    43 32 } g RTN

[Thomas Klemm]

For the sake of completeness I used only 3 variables \(p\), \(q\) and \(w\), similar to what I've done here.

Here is a Python program for a function that returns both \(\sin(x)\) and \(\cos(x)\):

Code:

def sc(x, n):
    a = b = 2 / x
    p = 1
    q = w = 0
    for _ in range(n):
        a, b = b, -(a + 2 * b)
        q = 1 / (a + q)
        p, w = w * q, (a * w + p) * q
    r = 1 + w**2
    return 2*w/r, (1 - w**2)/r

This is the corresponding program for the HP-42S:

Code:

00 { 64-Byte Prgm } ; x
01 2                ; 2     x
02 X<>Y             ; x     2
03 ÷                ; x/2
04 STO 00           ; a = x/2
05 STO 01           ; b = a
06 1                ; 1
07 STO 02           ; p = 1
08 CLX              ; 0
09 STO 03           ; q = 0
10 STO 04           ; w = 0
11 8                ; 8
12 STO 05           ; n = 8
13▸LBL 00           ; for loop
14 RCL 00           ; a
15 RCL 01           ; b     a
16 STO 00           ; a = b
17 2                ; 2     b       a
18 ×                ; 2b    a
19 +                ; a+2b
20 +/-              ; -(a+2b)
21 STO 01           ; b = -(a+2b)
22 RCL 02           ; p
23 RCL 04           ; w     p
24 STO 02           ; p = w
25 RCL 00           ; a     w       p
26 ×                ; a*w   p
27 +                ; a*w+p
28 STO 04           ; w = a*w+p
29 RCL 00           ; a
30 RCL 03           ; q      a
31 +                ; a+q
32 1/X              ; 1/(a+q)
33 STO 03           ; q = 1/(a+q)
34 STO× 02          ; p = w*q
35 STO× 04          ; w = (a*w+p)*q
36 DSE 05           ; n = n-1
37 GTO 00           ; repeat
38 RCL 04           ; t = tan(x/2)
39 ENTER            ; t     t
40 +                ; 2t
41 LASTX            ; t     2t
42 ENTER            ; t     t       2t
43 ×                ; t^2   2t
44 1                ; 1     t^2     2t
45 X<>Y             ; t^2   1       2t
46 -                ; 1-t^2 2t
47 1                ; 1     1-t^2   2t
48 LASTX            ; t^2   1       1-t^2   2t
49 +                ; 1+t^2 1-t^2   2t      2t
50 ÷                ; cos(x)=(1+t^2)/(1-t^2)  2t
51 X<>Y             ; 2t    cos(x)
52 LASTX            ; 1+t^2 2t      cos(x)
53 ÷                ; sin(x)=2t/(1+t^2)  cos(x)
54 X<>Y             ; cos(x) sin(x)
55 END              ; return

But it looks like we don't gain much with this approach.

[Albert Chan]

Slightly improved Thomas code (previous post)
Note: no changes done to tan(x/2) calculations

sin(x) * tan(x/2) = 2*sin(x/2)*cos(x/2) * sin(x/2)/cos(x/2) = 2*sin(x/2)^2 = versine(x) = 1 - cos(x)

Code:

00 { 63-Byte Prgm } ; x
01 2                ; 2     x
02 X<>Y             ; x     2
03 ÷                ; x/2
04 STO 00           ; a = x/2
05 STO 01           ; b = a
06 1                ; 1
07 STO 02           ; p = 1
08 CLX              ; 0
09 STO 03           ; q = 0
10 STO 04           ; w = 0
11 8                ; 8
12 STO 05           ; n = 8
13▸LBL 00           ; for loop
14 RCL 00           ; a
15 RCL 01           ; b     a
16 STO 00           ; a = b
17 2                ; 2     b       a
18 ×                ; 2b    a
19 +                ; a+2b
20 +/-              ; -(a+2b)
21 STO 01           ; b = -(a+2b)
22 RCL 02           ; p
23 RCL 04           ; w     p
24 STO 02           ; p = w
25 RCL 00           ; a     w       p
26 ×                ; a*w   p
27 +                ; a*w+p
28 STO 04           ; w = a*w+p
29 RCL 00           ; a
30 RCL 03           ; q      a
31 +                ; a+q
32 1/X              ; 1/(a+q)
33 STO 03           ; q = 1/(a+q)
34 STO× 02          ; p = w*q
35 STO× 04          ; w = (a*w+p)*q
36 DSE 05           ; n = n-1
37 GTO 00           ; repeat
38 RCL 04           ; t = tan(x/2)
39 ENTER

40 ENTER
41 +
42 LASTX
43 ENTER
44 ×             
45 1                ; 1    t^2  2t   t
46 +                
47 ÷                ; sin(x) t            t   t
48 STO× ST Y        ; sin(x) versine(x)   t   t
49 1                
50 RCL- ST Z        ; cos(x) sin(x) versine(x) tan(x/2)
51 END

PI 6 / R/S                  ; x = pi/6

T = 0.267949192431 ; tan(x/2)
Z = 0.133974596216 ; versine(x)
Y = 0.5                     ; sin(x)
X = 0.866025403784 ; cos(x)

Update: replaced "X↑2" by "Enter ×" (reason, see next post)

[Thomas Klemm]

(09-21-2024 05:31 PM)Albert Chan Wrote:  Slightly improved Thomas code

Nice.
But please note that the reason for the following lines is that operations like X↑2 or R↑ are missing on the HP-12C:

Code:

42 ENTER            ; t     t       2t
43 ×                ; t^2   2t

Therefore I avoided them to make a translation easier.
And with this of course goes all the arithmetic stack operations.
What about DSE, you may ask.
As often this is left as an exercise to the reader.