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:

                            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 the modified 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 )
 
 

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      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 )
 

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:       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
                                                                  T-register = Im(b_n)
Flags: /                                                     Z-register = Re(b_n)
Subroutine:  this program must produce   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 "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^      (  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
                                                                  T-register = Im(b_n)
Flags: /                                                     Z-register = Re(b_n)
Subroutine:  this program must produce   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                 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:    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
 
 
 

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall