hp41programs

KFibonacci

K-Fibonacci Numbers for the HP-41


Overview
 

 1°)  Fibonacci Numbers ( k = 2 )
 2°)  Tribonacci Numbers ( k = 3 )
 3°)  Tetranacci Numbers ( k = 4 )
 4°)  General Case
 5°)  Multiprecision

   a)  Fibonacci Numbers
   b)  Tribonacci Numbers
   c)  General Case
 

-The programs listed in paragraph 1°) are also listed in "Difference Equations for the HP41"
 
 

1°)  Fibonacci Numbers ( k = 2 )
 

-This sequence is defined by

  u0 = 0  ,   u1 = 1  &   un+1 = un + un-1
 

Data Registers: /
Flags: /
Subroutines: /
 
 

 01  LBL "FIB"
 02  ENTER^
 03  SIGN
 04  ENTER^
 05  CHS
 06  X<>Y
 07  ST+ Z
 08  LBL 01
 09  ST+ Y
 10  X<>Y
 11  DSE Z
 12  GTO 01
 13  END

 
( 25 bytes / SIZE 000 )
 
 

      STACK        INPUTS      OUTPUTS
           Y             /           un-1
           X             n            un
           L             /            n

  where n is a non-negative integer

Examples:

  49  XEQ "FIB"  >>>>   u49  = 7778742049           X<>Y   u48  = 4807526976
 480      R/S        >>>>   u480 = 9.216845715 E99   X<>Y   u479 = 5.696323921 E99

Notes:

-This routine uses a dummy  u -1 = 1
-The difference equation  un+1 = un + un-1  which may be re-written  un+1 - un - un-1 = 0   has a characteristic equation

    r2 - r - 1 = 0   which has 2 real roots  r1 = ( 1 + sqrt(5) ) / 2  and  r2 = ( 1 - sqrt(5) ) / 2

-So, the solution may be written   un = c1 r1n + c2 r2n       u0 = 0  &  u1 = 1 lead to   c1 = 1/sqrt(5)  &  c2 = -1/sqrt(5)
-And we could write the program
 

Data Registers:  R00 = n , R01 = sqrt(5)
Flags: /
Subroutines: /
 
 

 01  LBL "FIB"
 02  STO 00
 03  5
 04  SQRT
 05  STO 01
 06  1
 07  +
 08  2
 09  /
 10  X<>Y
 11  Y^X
 12  1
 13  RCL 01
 14  -
 15  2
 16  /
 17  RCL 00
 18  Y^X
 19  -
 20  RCL 01
 21  /
 22  END

 
( 30 bytes / SIZE 002 )
 
 

      STACK        INPUTS      OUTPUTS
           X             n            un

  where n is a non-negative integer

Examples:

 49  XEQ "FIB"  >>>>   u49 = 7778742108
 480      R/S        >>>>   OUT OF RANGE

Remarks:

-Though the second version is faster, there are also larger roundoff-errors.
-Moreover, an OUT OF RANGE occurs before un overflows !
 

2°)  Tribonacci Numbers ( k = 3 )
 

-Tribonacci numbers are defined by     u0 = 0  ,   u1 = 0  ,   u2 = 1    &   un+1 = un + un-1 + un-2
 

Data Registers: /
Flags: /
Subroutines: /
 
 

 01 LBL "TRIB"
 02 ENTER
 03 SIGN
 04 ENTER
 05 CHS
 06 X<>Y
 07 0
 08 LBL 01
 09 ST+ Z
 10 X<>Y
 11 ST+ Z
 12 X<> Z
 13 DSE T
 14 GTO 01
 15 END

 
( 29 bytes / SIZE 000 )
 
 

      STACK        INPUTS      OUTPUTS
           Z             /           un-2
           Y             /           un-1
           X             n            un

  where n is a non-negative integer

Examples:

    7  XEQ "TRIB"  >>>>   u7 =  13
   40       R/S          >>>>  u40 =  7046319384                     ---Execution time =  8s---
  380      R/S          >>>> u380 ~  6.741757257 E99            ---Execution time = 74s---
 

Notes:

-The difference equation  un+1 = un + un-1 + un-2    which may be re-written  un+1 - un - un-1 - un-2  = 0   has a characteristic equation

    r3 - r2 - r - 1 = 0   which has 2 real roots  r1 ~ 1.839286755  ,  r2 ~  -0.4196433776 + 0.6062907292 i  ,  r3 ~  -0.4196433776 - 0.6062907292 i

-So, the solution may be written   un = c1 r1n + c2 r2n + c3 r3n

    with   c1 ~ 0.1828035330  ,   c2 ~ -0.01940176648 + 0.3405465308 i   ,   c3 ~ -0.01940176648 - 0.3405465308 i

-This leads to the following program:
 

Data Registers: /
Flags: /
Subroutines: /
 
 

 01 LBL "TRIB"
 02 DEG
 03 STO Y
 04 124.6889974
 05 *
 06 105.0239517
 07 +
 08 COS
 09 .7051984758
 10 *
 11 .7373527058
 12 RCL Z
 13 Y^X
 14 *
 15 1.839286755
 16 R^
 17 Y^X
 18 .182803533
 19 *
 20 +
 21 END

 
( 91 bytes / SIZE 000 )
 
 

      STACK        INPUT      OUTPUT
           X             n            un

  where n is a non-negative integer

Examples:

    7  XEQ "TRIB"  >>>>   u7 ~  12.99999999
   40       R/S          >>>>  u40 ~   7046319352
  377      R/S          >>>> u377 ~  1.083489631 E99
 

Notes:

-Roundoff-errors are larger than with the 1st version
-But the routine is much faster for large numbers
-If n > 377 , we get OUT OF RANGE

-A simpler formula may be used ( cf reference [1] ):

   un = RND0 [ 3.b ( a + c + 1 )n-1 / ( b2 - 2.b + 4 ) ]

  where  a = ( 19 + 3 sqrt(33) )1/3  ,  c = ( 19 - 3 sqrt(33) )1/3  ,  b = ( 589 + 102 sqrt(33) )1/3
 

Data Registers: /
Flags: /
Subroutines: /
 
 

 01 LBL "TRIB"
 02 1.839286755
 03 X<>Y
 04 1
 05 -
 06 Y^X
 07 .336228117
 08 *
 09 FIX 0
 10 RND
 11 FIX 9
 12 END

 
( 42 bytes / SIZE 000 )
 
 

      STACK        INPUT      OUTPUT
           X             n            un

  where n is a non-negative integer

Examples:

    7  XEQ "TRIB"  >>>>   u7 =    13
   40       R/S          >>>>  u40 ~   7046319354
  378      R/S          >>>> u378 ~  1.992848127 E99
 

Note:

-If n > 378 , we get an OUT OF RANGE message
 

3°)  Tetranacci Numbers ( k = 4 )
 

-Similarly, tetranacci numbers are defined by  u0 = 0  ,   u1 = 0  ,   u2 = 0  ,   u3 = 1     &   un+1 = un + un-1 + un-2 + un-3
 

Data Register:  R00: temp
Flags: /
Subroutines: /
 
 

 01 LBL "TETRN"
 02 STO 00
 03 CLST
 04 SIGN
 05 ENTER
 06 CHS
 07 X<>Y
 08 R^
 09 LBL 01
 10 R^
 11 R^
 12 ST+ Y
 13 R^
 14 ST+ Z
 15 R^
 16 ST+ T
 17 R^
 18 DSE 00
 19 GTO 01
 20 END

 
( 35 bytes / SIZE 001 )
 
 

      STACK        INPUTS      OUTPUTS
           T             /           un-3
           Z             /           un-2
           Y             /           un-1
           X             n            un

  where n is a non-negative integer

Examples:

    7  XEQ "TETRN"  >>>>  u12 =  208
   38       R/S              >>>>  u38 =   5350220959                     ---Execution time = 10s---
  354      R/S              >>>> u354 ~  6.180313110 E99             ---Execution time = 94s---
 

Note:

-There is an OUT OF RANGE if n > 354
 

4°)  General Case
 

-We could write several programs for k = 5 , 6 , 7 , ....  but it's simpler to create a general program for all k > 1
 

Data Registers:      R00 = n   R01 = u(k)n-k+1  R02 = u(k)n-k+2  .........  Rkk = u(k)n
Flags: /
Subroutines: /
 
 

 01 LBL "KFIB" 
 02 STO 00
 03 X<>Y
 04 STO M
 05 .1
 06 %
 07 ISG X
 08 CLRGX
 09 SIGN		  10 STO IND M
 11 RCL 00
 12 +
 13 RCL M
 14 -
 15 STO N
 16 X>0?
 17 GTO 01
 18 LASTX		  19 +
 20 RCL IND X
 21 GTO 10
 22 LBL 01
 23 RCL M
 24 STO O
 25 RCL IND X 
 26 ENTER
 27 DSE O		  28 LBL 02
 29 X<> IND O
 30 ST+ Y
 31 DSE O
 32 GTO 02
 33 X<>Y
 34 STO IND M
 35 DSE N		  36 GTO 01
 37 LBL 10
 38 RCL 00       
 39 SIGN
 40 X<> M
 41 X<>Y
 42 CLA
 43 END

 
     ( 74 bytes / SIZE kkk+1 )
 
 

      STACK        INPUTS      OUTPUTS
           Y           k > 1            k
           X             n           u(k)n
           L             /            n

  where n is a non-negative integer

Examples:

     4  ENTER^  12   XEQ "KFIB"  >>>>    u(4)12 =  208                                   ---Execution time =  7s---
          CLX        38         R/S          >>>>    u(4)38 =   5350220959                     ---Execution time = 27s---
          CLX       354        R/S          >>>>   u(4)354 ~   6.180313114 E99            ---Execution time = 4m25s---
 

Notes:

-This program is of course slower than the corresponding  "FIB"  "TRIB"  "TETRN"  ...
-But it gives the same results
-The alpha "register" is cleared at the end.

-A simpler and faster method may be used ( cf reference [2] ):

   u(k)n = RND0 [ ( r-1 ) rn-k+1 / { 2 + (k+1).(r-2) }    where  r is the positive root of the characteristic equation

-The following variant employs Newton's method to find the real root of the characteristic equation.
-After finding this root, simply press R/S for another u(k)n with the same k-value.
 

Data Registers:      R00 to R03: temp  ( R01 = root )
Flags: /
Subroutines: /
 
 
 

 01 LBL "KFIB"
 02 STO 00
 03 X<>Y
 04 STO 03
 05 2
 06 STO 01
 07 LBL 00
 08 RCL 01
 09 RCL 01
 10 RCL 03
 11 Y^X
 12 STO 02
 13 LASTX		  14 1
 15 +
 16 *
 17 RCL 01      
 18 RCL 03
 19 DSE X
 20 Y^X
 21 RCL 03
 22 ST+ X
 23 *
 24 -
 25 X<> 02
 26 ST* Y		  27 ST+ X
 28 -
 29 1
 30 +
 31 RCL 02      
 32 /
 33 ST- 01
 34 ABS
 35  E-7
 36 X<Y?
 37 GTO 00
 38 RCL 00
 39 LBL 10		  40 STO 00
 41 RCL 01
 42 X<>Y
 43 RCL 03      
 44 DSE X
 45 -
 46 Y^X
 47 RCL 01
 48 1
 49 -
 50 ST* Y
 51 LASTX
 52 -		  53 RCL 03       
 54 1
 55 +
 56 *
 57 2
 58 +
 59 /
 60 FIX 0
 61 RND
 62 FIX 9
 63 RTN
 64 GTO 10
 65 END

 
    ( 88 bytes / SIZE 004 )
 
 

    STACK     INPUT1   OUTPUT1     INPUT2   OUTPUT2
        Y       k > 1          /           /          /
        X         n        u(k)n          m        u(k)m

  where n , m are non-negative integers

Examples:

    4   ENTER^   12    XEQ "KFIB"    >>>>   u(4)12 =  208
             38                        R/S            >>>>   u(4)38 ~  5350221004
            353                       R/S            >>>>  u(4)353 ~  3.206285314 E99

Notes:

-The method is quite good, but since the HP41 only works with 10 digits,
  roundoff-errors may be much larger than with the 1st version.
-It's also much faster, especially for large n-values

-Instead of solving the characteristic equation  f(r) =  rk - rk-1 - ..... - r - 1 = 0  ,  the HP41 solves  ( r - 1 ) f(r) = 0 = rk+1 - 2.rk + 1 = 0
-The root is always close to 2, so the 1st estimation is 2 ( line 05 )
 

5°)  Multiprecision
 

     a)  Fibonacci Numbers
 
 
 
 

Data Registers:      R00 = n , n-1 , .... , 0

                                   R01-R02-............-Rmm = un
                                   Rm+1-................-R2m = un-1
Flags: /
Subroutines: /
 
 
 

 01 LBL "FIBM"
 02 STO 00
 03  E9
 04 STO O
 05 CLX
 06 STO 02
 07 SIGN
 08 ST- 00
 09 STO 01
 10 STO N
 11 LBL 01
 12 RCL N
 13 STO M
 14 ST+ X
 15 ENTER		  16 CLX
 17 LBL 02
 18 RCL IND Y
 19 +
 20 RCL IND M
 21 STO IND Z
 22 +
 23 RCL X
 24 RCL O
 25 MOD
 26 STO IND M
 27 -
 28 RCL O
 29 /
 30 DSE Y		  31 DSE M
 32 GTO 02
 33 X=0?
 34 GTO 00      
 35 X<>Y
 36  E3
 37 /
 38 RCL N
 39 +
 40 3
 41 +
 42  E3
 43 /
 44 RCL N
 45 1		  46 +
 47 STO N
 48 +
 49 REGMOVE
 50 LASTX
 51 .1
 52 %
 53 +
 54 -
 55 1
 56 +
 57 REGMOVE
 58 RCL N
 59 LASTX
 60 +		  61 0
 62 STO IND Y
 63 R^
 64 STO 01
 65 LBL 00
 66 DSE 00
 67 GTO 01
 68 RCL N
 69 .1
 70 %
 71 1
 72 +
 73 CLA
 74 END

 
    ( 115 bytes / SIZE 2m+1 )
 
 

      STACK        INPUT      OUTPUT
           X          n > 1        1.mmm

  And  un  is stored in registers  R01 R02 .... Rmm

Examples:

    49  XEQ "FIBM"   >>>>   1.002                               ---Execution time = 38s---

    R01 = 7
    R02 = 778742049

-So,    u49 = 7778742049

-Likewise, u48 = 4807526976   is in registers R03-R04

   100   R/S   >>>>   1.003                               ---Execution time = 1m46s---

   R01 = 354
   R02 = 224848179
   R03 = 261915075

-So,    u100 = 354224848179261915075

-Likewise, u99 = 218922995834555169026  is in registers R04-R05-R06

   1000  R/S  >>>>  1.024                                 ---Execution time = 10s---                        with V41 in turbo mode

    R01 = 43                   R07 = 371780402    R13 = 295922593    R19 = 906533187
    R02 = 466557686     R08 = 481729089    R14 = 080322634    R20 = 938298969
    R03 = 937456435     R09 = 536555417    R15 = 775209689    R21 = 649928516
    R04 = 688527675     R10 = 949051890    R16 = 623239873    R22 = 003704476
    R05 = 040625802     R11 = 403879840    R17 = 322471161    R23 = 137795166
    R06 = 564660517     R12 = 079255169    R18 = 642996440    R24 = 849228875

-So,  u1000 = 43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169
                  295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875

and  u999  is in R25 to R48
 

Notes:

-Add zeros on the left if there are less than 9 digits in R02 ... Rmm
-For example, if  R02 = 123 , read  R02 = 000000123
 

     b)  Tribonacci Numbers
 
 

Data Registers:      R00 = n , n-1 , .... , 0

                                   R01 R02 ............ Rmm = un
                                   Rm+1 ................. R2m = un-1
                                   R2m+1 ............... R3m = un-2
Flags: /
Subroutines: /
 
 

 01 LBL "TRIBM"
 02 STO 00
 03 CLX
 04 STO 02
 05 STO 03
 06 SIGN
 07 STO 01
 08 STO N
 09 2
 10 ST- 00
 11 LBL 01
 12 RCL N
 13 ENTER
 14 STO M
 15 ST+ X
 16 STO O
 17 +		  18 ENTER
 19 CLX
 20 LBL 02
 21 RCL IND Y 
 22 +
 23 RCL IND O
 24 STO IND Z
 25 +
 26 RCL IND M
 27 STO IND O
 28 +
 29 RCL X
 30  E9
 31 MOD
 32 STO IND M
 33 -
 34  E9		  35 /
 36 DSE Y
 37 DSE O
 38 DSE M
 39 GTO 02
 40 X=0?
 41 GTO 00       
 42 X<>Y
 43 .1
 44 %
 45 +
 46 1.004
 47 +
 48 RCL N
 49  E6
 50 /
 51 +		  52 REGMOVE
 53 RCL N
 54 .1
 55 %
 56 +
 57  E-3
 58 +
 59 -
 60 REGMOVE 
 61 LASTX
 62 -
 63 REGMOVE
 64 X<>Y
 65 STO 01
 66 RCL N
 67 1
 68 +		  69 STO N
 70 ENTER
 71 ST+ X
 72 1
 73 ST+ Y
 74 ST+ Z
 75 CLX
 76 STO IND Y 
 77 STO IND Z
 78 LBL 00
 79 DSE 00
 80 GTO 01
 81 RCL N
 82 .1
 83 %
 84 ISG X
 85 CLA
 86 END

 
    ( 142 bytes / SIZE 3m+1 )
 
 

      STACK        INPUT      OUTPUT
           X          n > 2        1.mmm

  And  un  is stored in registers  R01 R02 .... Rmm

Examples:

•  40  XEQ "TRIBM"   >>>>  1.002                                 ---Execution time = 44s---

    R01 = 7
    R02 = 46319384      whence   u40 = 7046319384

    R03 = 3
    R04 = 831006429    whence   u39 = 3831006429

    R05 = 2
    R06 = 82876103     whence   u38 = 2082876103

•  100   R/S   >>>>   1.003                                 ---Execution time = 2m58s---

    R01 = 53324762
    R02 = 928098149
    R03 = 64722658     whence   u100 = 53324762928098149064722658

>>>  u99 & u98 are in R04-R05-R06 & R07-R08-R09  respectively.

    u99 = 28992087708416717612934417  &  u98 = 15762679542071167858843489
 

Note:

-Add zeros on the left if there are less than 9 digits in R02 ...
 

     c)  General Case
 
 

Data Registers:      R00 = n-k , n-k-1 , .... , 0

                                   R01 R02 ............ Rmm = un
                                   Rm+1 ................. R2m = un-1
                                   R2m+1 ............... R3m = un-2
Flags: /
Subroutines: /

-Line 105 is a three-byte GTO 01
 
 

 01 LBL "KFIBM"
 02 STO 00
 03 STO a
 04 X<>Y
 05 ST- 00
 06 STO M
 07 .1
 08 %
 09 ISG X
 10 CLRGX
 11 SIGN
 12 STO 01
 13 STO N
 14 ISG 00
 15 LBL 01
 16 RCL M
 17 RCL N
 18 ST* Y
 19 STO O
 20  E5
 21 /
 22 +
 23 STO Q		  24 CLX
 25 LBL 02
 26 RCL Q
 27 STO P
 28 X<>Y
 29 RCL IND Y  
 30 +
 31 DSE Y
 32 LBL 03
 33 RCL IND Y
 34 STO IND P
 35 +
 36 DSE P
 37 DSE Y
 38 GTO 03
 39 RCL X
 40 PI
 41 X^2
 42 INT
 43 10^X
 44 MOD
 45 STO IND P
 46 ST- Y		  47 X<> L
 48 /
 49 PI
 50 SIGN
 51 ST- Q
 52 X<>Y
 53 DSE O
 54 GTO 02        
 55 X=0?
 56 GTO 00
 57 RCL N
 58 RCL M
 59 STO O
 60 DSE X
 61 *
 62 1
 63 +
 64 .1
 65 %
 66 +
 67 RCL M
 68  E3
 69 /		  70 +
 71 RCL N
 72  E6
 73 /
 74 +
 75 LBL 04
 76 REGMOVE  
 77 RCL N
 78 .1
 79 %
 80 +
 81 -
 82  E-3
 83 -
 84 DSE O
 85 GTO 04
 86 ISG N
 87 CLX
 88 CLX
 89 RCL N
 90  E5
 91 /
 92 RCL M		  93 RCL N
 94 *
 95  E3
 96 /
 97 +
 98 1
 99 +
100 CLRGX
101 X<>Y
102 STO 01        
103 LBL 00
104 DSE 00
105 GTO 01
106 RCL M
107 RCL N
108  E3
109 /
110 ISG X
111 RCL a
112 SIGN
113 RDN
114 CLA
115 END

 
    ( 187 bytes / SIZE k.m+1 )
 
 

      STACK        INPUTS      OUTPUTS
           Y      1<  k < 11             k
           X        n > k-1        1.mmm
           L             /             n

  And  u(k)n  is stored in registers  R01 R02 .... Rmm

Example1:

     3     ENTER^
   100   XEQ "KFIBM"   >>>>   1.003                                 ---Execution time = 4m37s---

   R01 = 53324762
   R02 = 928098149
   R03 = 64722658     whence   u(3)100 = 53324762928098149064722658

>>>  u(3)99 & u(3)98 are stored in R04-R05-R06 & R07-R08-R09  respectively:

    u(3)99 = 28992087708416717612934417  &  u(3)98 = 15762679542071167858843489

Example2:

    4     ENTER^
  100       R/S       >>>>   1.004                                 ---Execution time = 5m33s---

  R01 = 2
  R02 = 505471397
  R03 = 838180985
  R04 = 96739296   whence   u(4)100 = 2505471397838180985096739296

>>>  u(4)99 , u(4)98 & u(4)97  are stored in R05-R06-R07-R08 , R09-R10-R11-R12 & R13-R14-R15-R16  respectively:

   u(4)99 = 1299813666022576563302979257
   u(4)98 =   674330414562637002982620399
   u(4)97 =   349835918709535855778683552
 

Notes:

-Add zeros on the left if there are less than 9 digits in R02 ...

-Synthetic registers P & Q are used, so don't stop this routine ( no risk of crash, however ).
-Synthetic register a is also employed, so "KFIBM" must not be called as more than a first level subroutine.
-Register "a" is only used to save n in L-register. Otherwise, simply delete lines 111-112-113 and line 03.

-There is a risk of roundoff-error if k > 10
 

References:

[1]  https://en.wikipedia.org/w/index.php?title=Generalizations_of_Fibonacci_numbers&oldid=933134595
[2]  Gregory P. B. Dresden & Zhaohui Du - "A Simplified Binet Formula for k-Generalized Fibonacci Numbers"