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 XY  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"