 The Museum of HP Calculators

# Continued Fractions for the HP-41

This program is Copyright © 2006-2007 by Jean-Marc Baillard and is used here by permission.

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

1°) Real Arguments

a) Wallis method   ( improved )
b) Modified Lenz's algorithm    ( new )

2°) Complex Arguments

a) Wallis method  ( improved )
b) Modified Lenz's algorithm    ( new )

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

a1       a2       a3                               an
f(x) =  b0 +  -----   -----   -----   .....................   -----  .......................     where  an &  bn are functions of  x and n
b1 +    b2 +    b3 +                            bn +

a) The Wallis method calculates  An  &  Bn  via the formulae   An = bn An-1 + an An-2  ,   Bn = bn Bn-1 + an Bn-2   with  A-1 = 1 ,  A0 = b0 ,  B-1 = 0 ,  B0 = 1
f(x) is the limit  An/Bn as n tends to infinity and these programs stop when 2 consecutive fractions are equal.

b) In the modified Lenz's algorithm,
we set   f0 = b0  if  b0 # 0  or  f0 = tiny  if  b0 = 0    ( tiny = E-30 in the following programs )
C0 = f0 and D0 = 0

then, for  n = 1 , 2 , .......      Cn = bn + an/Cn-1              (  Cn is replaced by tiny if Cn = 0 )
Dn = 1/ [ bn + an Dn-1 ]      (  Dn = 1/tiny  if the denominator = 0 )

and  fn = Cn Dn fn-1   until  | Cn Dn - 1 |  <=  a "small" number  (  E-9 in the following programs )    f(x) = the last  fn

1°) Real Arguments

a) Wallis Method

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

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

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

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                                        To avoid an  OUT OF RANGE   add  the following instuctions after line 14
13  ISG 02
14  CLX                                            RCL 03   ABS   RCL 04   ABS   X<Y?   X<>Y   X#0?   LOG   INT   10^X   ST/ 03   ST/ 04   ST/ 05   ST/ 06
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?                                         to avoid an infinite loop,  replace line 39  by   ST- Y   X=0?   SIGN   /   ABS   E-9   X<Y?
40  GTO 01
41  LBL 02
42  RCL 07
43  END

( 59 bytes / SIZE 008 )

 STACK INPUTS OUTPUTS Y b0 / X x f(x)

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

x      x2     x2     x2
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 b0 = 0 )
1  XEQ "CFR"                     >>>>   tanh 1 =  0.761594156  ( in 9 seconds )

Example2:    If  f is defined by:    b0 = 1  ;   an = 2.x + n  ,  bn =  x2 + n2  ( 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 = fn            R05 = Dn        R07 = tiny = E-30
R02 = n       R04 = Cn , bn    R06 = an

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

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  -                       lines 38 thru 41 may be replaced by  X#Y?  but with a risk of infinite loop
39  ABS
40   E-9                 or another "small" number like 2 E-9
41  X<Y?
42  GTO 01
43  RCL 03
44  END

( 63 bytes / SIZE 008 )

 STACK INPUTS OUTPUTS Y b0 / X x f(x)

Example:       b0 = 1  ;   an = 2.x + n  ,  bn =  x2 + n2  ( 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 = An    R08 R09 = An-1       R12-R13: temp           R14 = X
R02 = y       R06 R07 = Bn    R10 R11 = Bn-1                                           R15 = Y
R03 = n
T-register = Im(bn)
Flags: /                                                     Z-register = Re(bn)
Subroutine:  this program must produce   Y-register = Im(an)   assuming  x is in R01 ,  y is in R02  and n is in R03   ( n > 0 )
X-register = Re(an)

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                                      To avoid an  OUT OF RANGE   add  the following instuctions after line 24
23  ISG 03
24  CLX                                          RCL 04   RCL 06   X<Y?   X<>Y   X#0?   LOG   INT   10^X   ST/ 04   ST/ 06   ST/ 08   ST/ 10
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?                         to avoid an infinite loop,  replace line 95 by   E-8     X<Y?      ( or another "small" number )
96  GTO 01
97  LBL 02
98  RCL 15
99  RCL 14
100  END

( 127 bytes / SIZE 016 )

 STACK INPUTS OUTPUTS T Im(b0) / Z Re(b0) 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^      (  b0 = 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:    b0 = 0.2 + 0.3 i  ;   an = 2.z + n  ,  bn =  z2 + n2  ( 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 = fn     R08 R09 = Dn       R12-R13 = bn
R02 = y       R06 R07 = Cn    R10 R11 = an       R14 = tiny = E-30
R03 = n
T-register = Im(bn)
Flags: /                                                     Z-register = Re(bn)
Subroutine:  this program must produce   Y-register = Im(an)   assuming  x is in R01 ,  y is in R02  and n is in R03   ( n > 0 )
X-register = Re(an)

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                 or another "small" number like 2 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(b0) / Z Re(b0) / Y y Y X x X

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

Example:           f(z)  is defined by:    b0 = 0.2 + 0.3 i  ;   an = 2.z + n  ,  bn =  z2 + n2  ( 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:

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