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:

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 )

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

-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

( 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
Flags: /
Subroutine:   a program that returns

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

-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(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
Flags: /
Subroutine:  a program that returns

T-register = Im(bn)
Z-register = Re(bn)
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 83  X

( 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:

[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