Arctangent function program [HP-42S, Free42]

[Gerson W. Barbosa]

As I was evaluating ATAN(1/300) on Free42, I noticed the answer was approximately

1/(300 + 1/(900 + 1/(375 + 3/2800)))
=
0.00333332098773662486119811248858...

This is equivalent to

1/(300 + 1/(900 + 4/(1500 + 9/2100)))
=
0.00333332098773662486119811248858...

which suggests that

atan(1/x) = 1/(x + 1^2/(3x + 2^2/(5x + 3^2/(7x + ... ))))

I thought of writing a program for the HP-12C but I tried it on the HP-42S first, because it’s easier. I might write a 12C version later, but I don’t think it would be short and fast enough.
On the HP-42S replace 43 with 14, for arguments in the range ]0..1].

Code:


00 { 47-Byte Prgm }
01▸LBL "ATN"
02 1/X
03 ENTER
04 STO+ ST X
05 0
06 43
07 RCL× ST Z
08 STO+ ST T
09 X<> ST L 
10▸LBL 00
11 X<> ST T
12 +
13 LASTX
14 RCL- ST Z
15 X<> ST T
16 X↑2
17 X<>Y
18 STO÷ ST Y
19 X<> ST L
20 DSE ST X
21 GTO 00
22 X<> ST T
23 +
24 1/X
25 END

Example:

1 XEQ “ATN” 4 × ->

3.141592653589793238462643383279503

[Joe Horn]

What does the "PM" mean in line 09 of the program?

[Gerson W. Barbosa]

(09-14-2024 12:59 PM)Joe Horn Wrote:  What does the "PM" mean in line 09 of the program?

I wonder how that crept in. Anyway, gone for good now. Thank you!

Gerson.

[Thomas Klemm]

(09-14-2024 10:42 AM)Gerson W. Barbosa Wrote:  which suggests that

atan(1/x) = 1/(x + 1^2/(3x + 2^2/(5x + 3^2/(7x + ... ))))

From Inverse Tangent:

Quote:The inverse tangent has continued fraction representations



(Lambert 1770; Lagrange 1776; Wall 1948, p. 343; Olds 1963, p. 138)

[Gerson W. Barbosa]

(09-14-2024 01:43 PM)Thomas Klemm Wrote:  

(09-14-2024 10:42 AM)Gerson W. Barbosa Wrote:  which suggests that

atan(1/x) = 1/(x + 1^2/(3x + 2^2/(5x + 3^2/(7x + ... ))))

From Inverse Tangent:

Quote:The inverse tangent has continued fraction representations



(Lambert 1770; Lagrange 1776; Wall 1948, p. 343; Olds 1963, p. 138)

Yes, Thomas. That’s one of the many continued fraction representations for the arctangent(x) function l found when googling for them yesterday. What I have used is a continued fraction for arctangent(1/x). It appeared to me somewhat simpler to program, but I may be wrong. I would have to write programs for a few of those and see how they compare.

Cheers,

Gerson.

[Thomas Klemm]

Both formulas are essentially the same: just substitute \(x\) with \(\frac{1}{x}\).
My impression was that you came up with this representation by yourself.
So I just wanted to confirm that indeed it is correct.

[Albert Chan]

(09-14-2024 10:42 AM)Gerson W. Barbosa Wrote:  atan(1/x) = 1/(x + 1^2/(3x + 2^2/(5x + 3^2/(7x + ... ))))

Formula work well if user input = 1/x, is small (within ±1)

(05-31-2021 09:51 PM)Albert Chan Wrote:  \(\displaystyle\arctan(x) = 2\arctan\left( {x \over \sqrt{1+x^2}+1} \right)\)

RHS atan argument is always within ±1      plot(x/(1+sqrt(1+x^2)), x = -10 .. 10)

With atan "half-angle" formula, we can extend "ATN" code for all ranges.

Code:

00 { 59-Byte Prgm }
01▸LBL "ATN"
02 ENTER
03 X↑2
04 1
05 +
06 SQRT
07 1
08 +
09 X<>Y
10 ÷
11 ENTER
12 STO+ ST X
13 0
14 43
15 RCL× ST Z
16 STO+ ST T
17 X<> ST L
18▸LBL 00
19 X<> ST T
20 +
21 LASTX
22 RCL- ST Z
23 X<> ST T
24 X↑2
25 X<>Y
26 STO÷ ST Y
27 X<> ST L
28 DSE ST X
29 GTO 00
30 X<> ST T
31 +
32 1/X
33 STO+ ST X
34 END

1E99 XEQ "ATN" 2 ×      → 3.141592653589793238462643383279502

[Gerson W. Barbosa]

(09-14-2024 03:00 PM)Thomas Klemm Wrote:  Both formulas are essentially the same: just substitute \(x\) with \(\frac{1}{x}\).
My impression was that you came up with this representation by yourself.
So I just wanted to confirm that indeed it is correct.

You are right, Thomas. Thank you! A question still remains: would the program take up less steps if implemented according to the original formula? The initial 1/X instruction wouldn’t be necessary anymore, of course, but would that make the program shorter? A more challenging problem – at least for me – would be making the program HP-41 compatible. Recall arithmetic is quite handy!

[Thomas Klemm]

(09-14-2024 04:41 PM)Gerson W. Barbosa Wrote:  A question still remains: would the program take up less steps if implemented according to the original formula?

This is what I came up with:

Code:

00 { 39-Byte Prgm }
01▸LBL "ATN"
02 43
03 X<>Y
04 0
05▸LBL 00
06 2
07 RCL× ST T
08 +
09 1
10 +
11 RCL ST Z
12 RCL× ST Z
13 X↑2
14 X<>Y
15 ÷
16 DSE ST Z
17 GTO 00
18 1
19 +
20 ÷
21 END

[Thomas Klemm]

(09-14-2024 04:41 PM)Gerson W. Barbosa Wrote:  A more challenging problem – at least for me – would be making the program HP-41 compatible. Recall arithmetic is quite handy!

Without recall arithmetic we'd probably have to use storage arithmetic.
This leads to a lot of stack movements.
Therefore I suggest to use a register for the loop counter.
This gives us enough wiggle room to use more or less the same program:

Code:

01▸LBL "ATN"
02 14
03 STO 00
04 CLX
05▸LBL 00
06 2
07 RCL 00
08 ×
09 +
10 1
11 +
12 X<>Y
13 RCL 00
14 ×
15 X↑2
16 X<>Y
17 ÷
18 DSE 00
19 GTO 00
20 1
21 +
22 ÷
23 END

Instead of register 00 we could use register M.

[Gerson W. Barbosa]

(09-14-2024 05:37 PM)Thomas Klemm Wrote:  

Code:

01▸LBL "ATN"
02 14
03 STO 00
04 CLX
05▸LBL 00
06 2
07 RCL 00
08 ×
09 +[code]
10 1
11 +
12 X<>Y
13 RCL 00
14 ×
15 X↑2
16 X<>Y
17 ÷
18 DSE 00
19 GTO 00
20 1
21 +
22 ÷
23 END

Very nice!

Here’s my port to the HP-12C:

Code:


   001 {        1 } 1
   002 {        2 } 2
   003 {    44  0 } STO 0
   004 {    43 24 } g FRAC
   005 {        2 } 2
   006 {    45  0 } RCL 0
   007 {       20 } *
   008 {       40 } +
   009 {        1 } 1
   010 {       40 } +
   011 {       34 } X<=>Y
   012 {    45  0 } RCL 0
   013 {       20 } *
   014 {       36 } ENTER
   015 {       20 } *
   016 {       34 } X<=>Y
   017 {       10 } /
   018 {    45  0 } RCL 0
   019 {        1 } 1
   020 {       30 } -
   021 {    44  0 } STO 0
   022 {    43 35 } g x=0
   023 { 43 33 26 } g GTO 26
   024 {       33 } Rv
   025 { 43 33 05 } g GTO 05
   026 {       40 } +
   027 {        1 } 1
   028 {       40 } +
   029 {       10 } /
   030 { 43 33 00 } g GTO 00

1 g GTO 00 R/S 4 × -> 3.141592653
( ~ 17 seconds, almost instantly on the 12C+ )

(Edit to include a few missing prefixes in the listing)

[Thomas Klemm]

Here's the program for the HP-12C:

Code:

01 {       01 } 1
02 {       05 } 5
03 {    44 00 } STO O
04 {       35 } CLx
05 {       36 } ENTER
06 {       01 } 1
07 { 44 30 00 } STO - 0
08 {       40 } +
09 {    45 00 } RCL 0
10 {    43 35 } x=0
11 { 43 33 23 } GTO 23
12 {       02 } 2
13 {       20 } *
14 {       40 } +
15 {       34 } X<=>Y
16 {    45 00 } RCL 0
17 {       20 } *
18 {       36 } ENTER
19 {       20 } *
20 {       34 } X<=>Y
21 {       10 } /
22 { 43 33 06 } GTO 06
23 {       33 } Rv
24 {       10 } /

[Thomas Klemm]

And here for the HP-25C:

Code:

01: 31       : ENTER
02: 01       : 1
03: 05       : 5
04: 23 00    : STO 0
05: 34       : CLx
06: 31       : ENTER
07: 01       : 1
08: 23 41 00 : STO - 0
09: 51       : +
10: 24 00    : RCL 0
11: 15 71    : g x=0
12: 13 23    : GTO 23
13: 02       : 2
14: 61       : *
15: 51       : +
16: 21       : x<->y
17: 24 00    : RCL 0
18: 61       : *
19: 15 02    : g x^2
20: 21       : x<->y
21: 71       : /
22: 13 07    : GTO 07
23: 22       : Rv
24: 71       : /

[Gerson W. Barbosa]

(09-14-2024 10:37 PM)Thomas Klemm Wrote:  Here's the program for the HP-12C:

Code:

01 {       01 } 1
02 {       05 } 5
03 {    44 00 } STO O
04 {       35 } CLx
05 {       36 } ENTER
06 {       01 } 1
07 { 44 30 00 } STO - 0
08 {       40 } +
09 {    45 00 } RCL 0
10 {    43 35 } x=0
11 { 43 33 23 } GTO 23
12 {       02 } 2
13 {       20 } *
14 {       40 } +
15 {       34 } X<=>Y
16 {    45 00 } RCL 0
17 {       20 } *
18 {       36 } ENTER
19 {       20 } *
20 {       34 } X<=>Y
21 {       10 } /
22 { 43 33 06 } GTO 06
23 {       33 } Rv
24 {       10 } /

Significantly shorter, but again about 17 seconds on my older 12C. Do we really need 15 iterations? 13 appear to be just fine. Timing down to ~ 14.5 seconds.

[Thomas Klemm]

(09-14-2024 11:03 PM)Gerson W. Barbosa Wrote:  Do we really need 15 iterations? 13 appear to be just fine. Timing down to ~ 14.5 seconds.

Excellent!

[Albert Chan]

From wikipedia Gauss's continued fraction



Wikipedia solve for \(_2 F _1(1,b;c;z) \) CF with a=0, c=c-1, and simplify. Here is another way, with b=1



\( \displaystyle \begin{align}
_2 F _1 (a,1;c;z)
&= 1 + \frac{a\;z}{c} \; _2F _1(a+1,1;c+1;z)
\\ &= \frac{1}{1-a\;z \left( \frac{_2F _1(a+1,1;c+1;z)}
{c\; _2F _1(a,1;c;z)}\right)} \\ \\ &=
\frac{1}{1\;-}\;
\frac{a\,z}{c \;- }\;
\frac{(c-a)\,1z}{(c+1) \;- }\;
\frac{(c+0)(a+1)\,z}{(c+2) \;- }\;
\frac{(c-a+1)\,2z}{(c+3) \;- }\;
\frac{(c+1)(a+2)\,z}{(c+4) \;- }\;
\frac{(c-a+2)\,3z}{(c+5) \;- }\; ...
\end{align}\)


\(\displaystyle \begin{align}
\arctan(z) &=
z\;_2\!F_1 \left(\frac{1}{2}, 1; \frac{3}{2}; -z^2 \right)

\\ &=
\frac{z}{1\;+}\;
\frac{\frac{1}{2}\,z^2}{\frac{3}{2} \;+ }\;
\frac{1 \, z^2}{\frac{5}{2} \;+ }\;
\frac{\frac{9}{4} \, z^2}{\frac{7}{2}\;+ }\;
\frac{4 \, z^2}{\frac{9}{2}\;+ }\;
\frac{\frac{25}{4}\, z^2}{\frac{11}{2}\;+ }\;
\frac{9\, z^2}{\frac{13}{2} \;+ }\; ...

\\ \\ &=
\frac{z}{1\;+}\;
\frac{z^2}{3 \;+ }\;
\frac{2^2 \,z^2}{5 \;+ }\;
\frac{3^2 \, z^2}{7\;+ }\;
\frac{4^2 \, z^2}{9\;+ }\;
\frac{5^2 \, z^2}{11\;+ }\;
\frac{6^2 \, z^2}{13\;+ }\; ...
\end{align}\)

This may also be of interest:

Continued Fraction expansion of tan(x)
Continued Fraction of tan(n*x)

[Gerson W. Barbosa]

(09-14-2024 10:37 PM)Thomas Klemm Wrote:  Here's the program for the HP-12C:

Code:

01 {       01 } 1
02 {       05 } 5
03 {    44 00 } STO O
04 {       35 } CLx
05 {       36 } ENTER
06 {       01 } 1
07 { 44 30 00 } STO - 0
08 {       40 } +
09 {    45 00 } RCL 0
10 {    43 35 } x=0
11 { 43 33 23 } GTO 23
12 {       02 } 2
13 {       20 } *
14 {       40 } +
15 {       34 } X<=>Y
16 {    45 00 } RCL 0
17 {       20 } *
18 {       36 } ENTER
19 {       20 } *
20 {       34 } X<=>Y
21 {       10 } /
22 { 43 33 06 } GTO 06
23 {       33 } Rv
24 {       10 } /

Since we’re at it, here’s tangent using basically the same nice scheme of yours:

Code:

01 {       36 } ENTER
02 {       20 } *
03 {        8 } 8
04 {     44 0 } STO 0
05 {       35 } CLx
06 {       36 } ENTER
07 {        1 } 1
08 { 44 30  0 } STO - 0
09 {       40 } +
10 {    45  0 } RCL 0
11 {    43 35 } x=0
12 { 43 33 19 } GTO 19
13 {        2 } 2
14 {       20 } *
15 {       40 } +
16 {       10 } /
17 {       16 } CHS
18 { 43 33 07 } GTO 7
19 {       33 } Rv
20 {       22 } 1/x
21 {       34 } x<=>y
22 {    43 21 } SQRTx
23 {       20 } *
24 { 43 33 00 } GTO 00

3.141592654 ENTER 4 / g GTO 00 R/S -> 1.000000000

1.5707 g GTO 00 R/S -> 10381.32(037)

( ~ 7 seconds on the good old 12C )

———-

P.S.:

Or, if we want all of them on the stack:

Code:

...
24 {       36 } ENTER
25 {       36 } ENTER
26 {       36 } ENTER
27 {       20 } *
28 {        1 } 1
29 {       40 } +
30 {    43 21 } SQRTx
31 {       22 } 1/x
32 {       20 } *
33 {    43 36 } LASTx
34 {       34 } x<=>y
35 { 43 33 00 } GTO 00

g GTO 00 R/S ->

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

———-

Code:


36 { xx xx xx } XXXX 
37 { xx xx xx } XXXX
...
...  (*)

g GTO 36 R/S ->

Z: ATAN(x)
Y: ACOS(x)
X: ASIN(x)


(*) These are left as an exercise to the interested reader.

[Thomas Klemm]

Instead of using



we could use



This makes the program a bit shorter at the cost of an additional register:

Code:

01 {    44 01 } STO 1
02 {       09 } 9
03 {    44 00 } STO O
04 {       35 } CLx
05 {       36 } ENTER
06 {    45 00 } RCL 0
07 {    43 35 } x=0
08 { 43 33 20 } GTO 20
09 {       02 } 2
10 {       20 } *
11 {       01 } 1
12 { 44 30 00 } STO - 0
13 {       30 } -
14 {    45 01 } RCL 1
15 {       10 } /
16 {       34 } X<=>Y
17 {       30 } -
18 {       22 } 1/x
19 { 43 33 06 } GTO 06
20 {       33 } Rv

[Gerson W. Barbosa]

(09-15-2024 09:18 PM)Thomas Klemm Wrote:  This makes the program a bit shorter at the cost of an additional register:

Code:

01 {    44 01 } STO 1
02 {       09 } 9
03 {    44 00 } STO O
04 {       35 } CLx
05 {       36 } ENTER
06 {    45 00 } RCL 0
07 {    43 35 } x=0
08 { 43 33 20 } GTO 20
09 {       02 } 2
10 {       20 } *
11 {       01 } 1
12 { 44 30 00 } STO - 0
13 {       30 } -
14 {    45 01 } RCL 1
15 {       10 } /
16 {       34 } X<=>Y
17 {       30 } -
18 {       22 } 1/x
19 { 43 33 06 } GTO 06
20 {       33 } Rv

No problem! Stack-only is not important for the HP-12C. Size is. Thanks for the improvement!
Except for the borrowed code, the following is not size-optimized. But that’s just a test.

HP-12C:

Code:


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

Initialization:

3.141592654 ENTER 180 / STO 5

0 (DEG) or 1 (RAD) g GTO 47 R/S

Usage:

x R/S ->

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

First and second quadrants only. 90 degrees case not handled.

——

HP-15C listing the JRPN 15C Simulator:

Code:


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

[Thomas Klemm]

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