ContFrac

Continued Fractions for the HP-41

Overview

1°) Real Arguments

a) Wallis method
b) Modified Lenz's algorithm

2°) Complex Arguments

a) Wallis method
b) Modified Lenz's algorithm

-The 4 following programs compute a function f(x) defined by a continued fraction:

                            a₁       a₂       a₃                               a_n
      f(x) = b₀ + -----   -----   -----   .....................   ----- .......................     where a_n & b_n are functions of x and n
                          b₁ +    b₂ +    b₃ +                            b_n +

a) The Wallis method calculates A_n & B_n via the formulae A_n = b_n A_n-1 + a_n A_n-2 , B_n = b_n B_n-1 + a_n B_n-2 with A_-1 = 1 , A₀ = b₀ , B_-1 = 0 , B₀ = 1
f(x) is the limit A_n/B_n as n tends to infinity and these programs stop when 2 consecutive fractions are equal.

b) In themodified Lenz's algorithm,
we set f₀ = b₀ if b₀ # 0 or f₀ = tiny if b₀ = 0 ( tiny = E-30 in the following programs )
C₀ = f₀ and D₀ = 0

then, for n = 1 , 2 , ....... C_n = b_n + a_n/C_n-1 ( C_n is replaced by tiny if C_n = 0 )
D_n = 1/ [ b_n + a_n D_n-1 ] ( D_n = 1/tiny if the denominator = 0 )

and f_n = C_n D_n f_n-1 until | C_n D_n - 1 | <= a "small" number ( E-9 in the following programs ) f(x) = the last f_n

1°) Real Arguments

a) Wallis Method

Data Registers: • R00 = Subroutine name ( Register R00 is to be initialized before executing "CFR" )

R01 = x R03 = A_n R05 = A_n-1 R07 = the successive A_n/B_n and f(x) at the end
R02 = n R04 = B_n R06 = B_n-1

Flags: / Y-register = b_n
Subroutine: this program must produce X-register = a_n assuming x is in R01 and n is in R02 ( n > 0 )

-To avoid an OUT OF RANGE   add the following instuctions after line 14
RCL 03   ABS   RCL 04   ABS   X<Y?   X<>Y   X#0?   LOG   INT   10^X   ST/ 03   ST/ 04   ST/ 05   ST/ 06
-To avoid a possible infinite loop, you can replace line 39 by   ST- Y   X=0?   SIGN   /   ABS   E-9   X<Y?

01 LBL "CFR"
02 STO 01
03 X<>Y
04 STO 03
05 STO 07
06 CLX
07 STO 02
08 STO 06
09 SIGN

10 STO 04
11 STO 05
12 LBL 01
13 ISG 02
14 CLX
15 XEQ IND 00
16 X=0?
17 GTO 02
18 ST* 06

19 RCL 05
20 *
21 RCL Y
22 RCL 03
23 STO 05
24 *
25 +
26 STO 03
27 RCL 06

28 RCL 04
29 STO 06
30 R^
31 *
32 +
33 STO 04
34 X=0?
35 GTO 01
36 /

37 ENTER^
38 X<> 07
39 X#Y?
40 GTO 01
41 LBL 02
42 RCL 07
43 END

( 59 bytes / SIZE 008 )

STACK	INPUTS	OUTPUTS
Y	b₀	/
X	x	f(x)

Example1: tanh x may be computed by a continued fraction:

                   x      x²     x²     x²
   tanh x = ---- ---- ---- ---- .....................
                  1+    3+   5+    7+

-Load the short routine ( let's call it "T" , any global LBL , maximum 6 characters )

LBL "T"        GTO 01      -               X^2
RCL 02        RCL 01      RCL 02     RTN
1                   RTN           +
X#Y?           LBL 01       RCL 01

T ASTO 00 and, for instance, to get tanh 1

0 ENTER^ ( here b₀ = 0 )
1 XEQ "CFR" >>>> tanh 1 = 0.761594156 ( in 9 seconds )

Example2: If f is defined by: b₀ = 1 ; a_n = 2.x + n , b_n = x² + n² ( n > 0 ) ; evaluate f(2) & f(Pi)

LBL "T"        X^2           RCL 02
RCL 01         +               +
X^2              RCL 01     RTN
RCL 02        ST+ X

T ASTO 00

1 ENTER^
2 XEQ "CFR" >>>> f(2) = 1.876576535 ( 8 seconds )

1 ENTER^
PI R/S >>>> f(Pi) = 1.636265662 ( 9 seconds )

b) Modified Lenz's Algorithm

Data Registers: • R00 = Subroutine name ( Register R00 is to be initialized before executing "LCFR" )

R01 = x R03 = f_n R05 = D_n R07 = tiny = E-30
R02 = n R04 = C_n , b_n R06 = a_n

Flags: / Y-register = b_n
Subroutine: this program must produce X-register = a_n assuming x is in R01 and n is in R02 ( n > 0 )

-Lines 38 thru 41 may be replaced by X#Y? but with a risk of infinite loop.

01 LBL "LCFR"
02 STO 01
03 E-30
04 STO 07
05 0
06 STO 02
07 STO 05
08 R^
09 X=0?

10 RCL 07
11 STO 03
12 STO 04
13 LBL 01
14 ISG 02
15 CLX
16 XEQ IND 00
17 STO 06
18 RCL 04

19 /
20 X<>Y
21 STO 04
22 +
23 X=0?
24 RCL 07
25 X<> 04
26 RCL 05
27 RCL 06

28 *
29 +
30 X=0?
31 RCL 07
32 1/X
33 STO 05
34 RCL 04
35 *
36 ST* 03

37 1
38 -
39 ABS
40 E-9
41 X<Y?
42 GTO 01
43 RCL 03
44 END

( 63 bytes / SIZE 008 )

STACK	INPUTS	OUTPUTS
Y	b₀	/
X	x	f(x)

Example: b₀ = 1 ; a_n = 2.x + n , b_n = x² + n² ( n > 0 ) ; evaluate f(2)

LBL "T"        X^2           RCL 02
RCL 01         +               +
X^2              RCL 01     RTN
RCL 02        ST+ X

T ASTO 00

1 ENTER^
2 XEQ "LCFR" >>>> f(2) = 1.876576535 ( 10 seconds )

-"CFR" and "LCFR" are almost equivalent.

2°) Complex Arguments

a) Wallis Method

-Here we have to compute f(z) = Z = X + i.Y with z = x + i.y

Data Registers: • R00 = Subroutine name ( Register R00 is to be initialized before executing "CFRZ" )

                                         R01 = x       R04 R05 = A_n    R08 R09 = A_n-1       R12-R13: temp           R14 = X
                                         R02 = y       R06 R07 = B_n    R10 R11 = B_n-1                                           R15 = Y
                                         R03 = n
Flags: /
Subroutine:   a program that returns

     T-register = Im(b_n)
     Z-register = Re(b_n)
     Y-register = Im(a_n)              assuming x is in R01 , y is in R02 and n is in R03   ( n > 0 )
     X-register = Re(a_n)

-To avoid an OUT OF RANGE   add the following instuctions after line 24
RCL 04   RCL 06   X<Y?   X<>Y   X#0?   LOG   INT   10^X   ST/ 04   ST/ 06   ST/ 08   ST/ 10
-To avoid a possible infinite loop, replace line 95 by   E-8     X<Y?      ( or another "small" number )

01 LBL "CFRZ"
02 STO 01
03 RDN
04 STO 02
05 X<> Z
06 STO 15
07 X<>Y
08 STO 14
09 R-P
10 STO 04
11 X<>Y
12 STO 05
13 CLX
14 STO 03
15 STO 07
16 STO 09
17 STO 10
18 STO 11
19 SIGN
20 STO 06

21 STO 08
22 LBL 01
23 ISG 03
24 CLX
25 XEQ IND 00
26 R-P
27 X=0?
28 GTO 02
29 ST* 08
30 ST* 10
31 RDN
32 ST+ 09
33 ST+ 11
34 RDN
35 R-P
36 STO 12
37 X<>Y
38 STO 13
39 RCL 05
40 +

41 X<>Y
42 RCL 04
43 *
44 P-R
45 RCL 09
46 RCL 08
47 P-R
48 X<>Y
49 ST+ T
50 RDN
51 +
52 R-P
53 X<> 04
54 STO 08
55 X<>Y
56 X<> 05
57 STO 09
58 RCL 07
59 RCL 13
60 +

61 RCL 06
62 RCL 12
63 *
64 P-R
65 RCL 11
66 RCL 10
67 P-R
68 X<>Y
69 ST+ T
70 RDN
71 +
72 R-P
73 X<> 06
74 STO 10
75 X<>Y
76 X<> 07
77 STO 11
78 RCL 05
79 RCL 07
80 -

81 RCL 04
82 RCL 06
83 X=0?
84 GTO 01
85 /
86 P-R
87 ENTER^
88 X<> 14
89 -
90 X<>Y
91 ENTER^
92 X<> 15
93 -
94 R-P
95 X#0?
96 GTO 01
97 LBL 02
98 RCL 15
99 RCL 14
100 END

( 127 bytes / SIZE 016 )

STACK	INPUTS	OUTPUTS
T	Im(b₀)	/
Z	Re(b₀)	0
Y	y	Y
X	x	X

with f(z) = f(x+i.y) = X + i.Y

Example1: Let's compute tanh ( 1+2.i )

LBL "T"       X#Y?          RCL 01       RCL 03       R-P           X^2
CLST          GTO 01       RTN           +                 X<>Y        P-R
RCL 03       CLX            LBL 01       RCL 02       ST+ X       RTN
1                  RCL 02       -                 RCL 01       X<>Y

T ASTO 00

       0   ENTER^    ENTER^      ( b₀ = 0 = 0 + 0.i )
       2   ENTER^
       1   XEQ "CFRZ"   >>>>    1.166736258   X<>Y -0.243458200   whence   tanh ( 1+2.i ) = 1.166736258 - 0.243458200 i    ( in 96 seconds )

Example2: If f(z) is defined by: b₀ = 0.2 + 0.3 i ; a_n = 2.z + n , b_n = z² + n² ( n > 0 ) ; evaluate f(1+2.i)

LBL "T"
RCL 02
RCL 01
XEQ "Z^2" ( cf "Complex Functions for the HP-41" )
RCL 03
X^2
+
RCL 01
ST+ X
RCL 03
+
RCL 02
ST+ X
X<>Y
RTN

T ASTO 00

0.3 ENTER^
0.2 ENTER^
2 ENTER^
1 XEQ "CFRZ" >>>> 1.084596569 X<>Y -0.749791581 whence f(1+2.i) = 1.084596569 - 0.749791581 i ( in 74 seconds )

b) Modified Lenz's Algorithm

f(z) = Z = X + i.Y with z = x + i.y

Data Registers: • R00 = Subroutine name ( Register R00 is to be initialized before executing "LCFRZ" )

                                         R01 = x       R04 R05 = f_n     R08 R09 = D_n       R12-R13 = b_n
                                         R02 = y       R06 R07 = C_n    R10 R11 = a_n       R14 = tiny = E-30
                                         R03 = n
Flags: /
Subroutine: a program that returns

     T-register = Im(b_n)
     Z-register = Re(b_n)
     Y-register = Im(a_n)              assuming x is in R01 , y is in R02 and n is in R03   ( n > 0 )
     X-register = Re(a_n)

01 LBL "LCFRZ"
02 STO 01
03 RDN
04 STO 02
05 CLX
06 E-30
07 STO 14
08 RDN
09 R-P
10 X=0?
11 R^
12 STO 04
13 STO 06
14 X<>Y
15 STO 05
16 STO 07
17 CLX
18 STO 03

19 STO 08
20 STO 09
21 LBL 01
22 ISG 03
23 CLX
24 XEQ IND 00
25 R-P
26 STO 10
27 RDN
28 STO 11
29 RDN
30 STO 12
31 RDN
32 STO 13
33 X<> Z
34 RCL 07
35 -
36 X<>Y

37 RCL 06
38 /
39 P-R
40 X<> Z
41 +
42 X<>Y
43 RCL 12
44 +
45 R-P
46 X=0?
47 RCL 14
48 STO 06
49 X<>Y
50 STO 07
51 RCL 09
52 RCL 11
53 +
54 RCL 08

55 RCL 10
56 *
57 P-R
58 RCL 13
59 ST+ Z
60 CLX
61 RCL 12
62 +
63 R-P
64 X=0?
65 RCL 14
66 1/X
67 STO 08
68 RCL 06
69 *
70 ST* 04
71 X<>Y
72 CHS

73 STO 09
74 RCL 07
75 +
76 ST+ 05
77 X<>Y
78 P-R
79 1
80 -
81 R-P
82 E-9
83 X<Y?
84 GTO 01
85 RCL 05
86 RCL 04
87 P-R
88 END

( 111 bytes / SIZE 015 )

STACK	INPUTS	OUTPUTS
T	Im(b₀)	/
Z	Re(b₀)	/
Y	y	Y
X	x	X

with f(z) = f(x+i.y) = X + i.Y

Example: f(z) is defined by: b₀ = 0.2 + 0.3 i ; a_n = 2.z + n , b_n = z² + n² ( n > 0 ) ; evaluate f(1+2.i)

LBL "T"
RCL 02
RCL 01
XEQ "Z^2" ( cf "Complex Functions for the HP-41" )
RCL 03
X^2
+
RCL 01
ST+ X
RCL 03
+
RCL 02
ST+ X
X<>Y
RTN

T ASTO 00

0.3 ENTER^
0.2 ENTER^
2 ENTER^
1 XEQ "LCFRZ" >>>> 1.084596568 X<>Y -0.749791576 whence f(1+2.i) = 1.084596568 - 0.749791576 i ( in 60 seconds )

-"LCFRZ" is shorter than "CFRZ"
-Moreover, it seems faster ( there are less R-P & P-R instructions ).
-The modified Lenz's algorithm also avoids the risk of "OUT OF RANGE"

References:

[1] Abramowitz & Stegun - "Handbook of Mathematical Functions" - Dover Publications - ISBN 0-486-61272-4
[2] Press , Vetterling , Teukolsky , Flannery - "Numerical Recipes in C++" - Cambridge University Press - ISBN 0-521-75033-4

hp41programs

Continued Fractions for the HP-41