(15C) Fibonacci Numbers

[Eddie W. Shore]

This program calculates the nth Fibonacci number. Store n in memory register 0 before running the program.

Formula: ( (1 + √5)^n – (1 - √5)^n ) / (2^n * √5)

Caution: Due to internal calculator calculation, you may not get an integer answer. You might have to round the answer manually. From Joe Horn, with thanks and appreciation for letting me post his comment:

"Hi, Eddie! I just keyed up your HP-15C program for Fibonacci numbers, and noticed that it gets wrong answers due to round-off errors much of the time. For example, with an input of 5, it should output 5, but it outputs 4.9999999998. Even rounding to the nearest integer doesn't fix the problem for inputs of 40 through 49. The only way to get the right answers for all inputs is with a loop (the usual method)." - Joe Horn

I would recommend running this program when n< 40.

Code:


Step    Key    Code
001    LBL D    42, 21, 14
002    1    1
003    ENTER    36
004    5    5
005    √     11
006    +    40
007    RCL 0    45, 0
008    Y^X    14
009    1    1
010    ENTER    36
011    5    5
012    √     11
013    -    30
014    RCL 0    45, 0
015    Y^X    14
016    -    30
017    2    2
018    RCL 0    45, 0
019    Y^X    14
020    ÷    10
021    5    5
022    √     11
023    ÷     10
024    RTN    43, 32

Example: R0 = n = 16, Output: 987

[Eddie W. Shore]

Loop Method – Joe Horn

Joe Horn provided a more accurate program which uses loops. It is slower, however should speed should not be a problem if a HP 15C Limited Edition or emulator is used. Full credit and thanks goes to Joe Horn for providing this program.

Code:


Step    Key    Code
001    LBL A    42, 21, 11
002    STO 0    44, 0
003    1    1
004    ENTER    36
005    0    0
006    LBL 1    42, 21, 1
007    +    40
008    LST X    43, 36
009    X<>Y    34
010    DSE 0    42, 5, 0
011    GTO 1    22, 1
012    RTN    43, 32

A nice part is that you don’t have to pre-store n in memory 0. This method is accurate for n ≤ 49.

Comparing Results
Here are some results of some selected n:

Code:


n    Approximation    Loop Method
6    8    8
15    610    610
16    987    987
22    17710.9999    17711
29    514228.9979    514229
36    14930351.92    14930352
40    102334154.4    102334155
44    701408728.7    701408733
49    7778741992    7778742049

Bug?

I manually calculated the approximation formula on my HP 15C LE (Limited Edition) and HP Prime for n = 44 and n = 49.

n = 44
HP15C LE: 701408728.7
HP Prime: 701408733.002

n = 49
HP 15C LE: 7778741992
HP Prime: 7778742049.02

I think we may have found a bug on the HP 15C LE.

[Thomas Klemm]

Computation by rounding

We can use rounding to the nearest integer:

\(
F_{n}=\left\lfloor {\frac {\varphi ^{n}}{\sqrt {5}}}\right\rceil ,\ n\geq 0.
\)

However we experience errors when \(n \gt 40\).

Splitting the Index

But since the biggest value for \(n\) is only \(49\) we can use the following formulas to reduce the index:

\(
\begin{aligned}
F_{2n-1}&={F_{n}}^{2}+{F_{n-1}}^{2}\\
F_{2n}&=(2F_{n-1}+F_{n})F_{n}\\
\end{aligned}
\)

Now even in the worst case of \(F_{49}\) we can use rounding to calculate \(F_{24}\) and \(F_{25}\).

Examples

\(
\begin{aligned}
F_{44}&=&(2F_{21}+F_{22})F_{22} &=& (2 \cdot 10,946 + 17,711) \cdot 17,711 &=& 701,408,733 \\
F_{49}&=&{F_{25}}^{2}+{F_{24}}^{2} &=& 75,025^2 + 46,386^2 &=& 7,778,742,049 \\
\end{aligned}
\)

Program for the HP-15C

Code:

   001 {    42 21 11 } f LBL A
   002 {       43 20 } g x=0
   003 {       43 32 } g RTN
   004 {           1 } 1
   005 {          30 } −
   006 {           2 } 2
   007 {          10 } ÷
   008 {          36 } ENTER
   009 {       43 44 } g INT
   010 {    43 30  5 } g TEST x=y
   011 {       22  0 } GTO 0
   012 {       32  1 } GSB 1
   013 {          34 } x↔y
   014 {           2 } 2
   015 {          20 } ×
   016 {          34 } x↔y
   017 {          40 } +
   018 {       43 36 } g LSTΧ
   019 {          20 } ×
   020 {       43 32 } g RTN
   021 {    42 21  0 } f LBL 0
   022 {       32  1 } GSB 1
   023 {       43 11 } g x²
   024 {          34 } x↔y
   025 {       43 11 } g x²
   026 {          40 } +
   027 {       43 32 } g RTN
   028 {    42 21  1 } f LBL 1
   029 {           5 } 5
   030 {          11 } √x̅
   031 {          36 } ENTER
   032 {          36 } ENTER
   033 {           1 } 1
   034 {          40 } +
   035 {           2 } 2
   036 {          10 } ÷
   037 {          36 } ENTER
   038 {       43 33 } g R⬆
   039 {          14 } yˣ
   040 {       43 33 } g R⬆
   041 {          10 } ÷
   042 {          20 } ×
   043 {       43 36 } g LSTΧ
   044 {       32  2 } GSB 2
   045 {          34 } x↔y
   046 {    42 21  2 } f LBL 2
   047 {          48 } .
   048 {           5 } 5
   049 {          40 } +
   050 {       43 44 } g INT
   051 {       43 32 } g RTN

Timing

For \(F_{49}\) it takes about 5 seconds while Joe's program takes about 30 seconds.
This was measured using the Nonpareil emulator on a MacBook.
But I assume that the timing will be similar on the real calculator.

References

• Fibonacci sequence
  

[Thomas Klemm]

Another idea is to calculate:

\(
\left[
\begin{matrix}
0 & 1\\
1 & 1
\end{matrix}
\right]^n
\)

Here's a Python program that uses exponentiation by squaring:

Code:

from sympy import Matrix

def fib(n):
    A = Matrix([[1, 0], [0, 1]])
    B = Matrix([[0, 1], [1, 1]])
    while n:
        if n % 2:
            A *= B
        B *= B
        n //= 2
    return A

For \(F_{49}\) we can run fib(49) to get:

\(
\left[
\begin{matrix}
4807526976 & 7778742049\\
7778742049 & 12586269025
\end{matrix}
\right]
\)

Apparently we have to pick \(7778742049\).

This program for the HP-15C uses matrix operations:

Code:

   001 {    42 21 11 } f LBL A
   002 {           2 } 2
   003 {          36 } ENTER
   004 {    42 23 11 } f DIM A
   005 {    42 23 12 } f DIM B
   006 {    42 26 13 } f RESULT C
   007 {    42 16  1 } f MATRIX 1
   008 {       43 35 } g CLx
   009 {    44 12  u } STO B
   010 {           1 } 1
   011 {    44 12  u } STO B
   012 {    44 12  u } STO B
   013 {    44 12  u } STO B
   014 {       43 32 } g RTN
   015 {    44 11  u } STO A
   016 {          34 } x↔y
   017 {    44 11  u } STO A
   018 {    44 11  u } STO A
   019 {          34 } x↔y
   020 {    44 11  u } STO A
   021 {       43 32 } g RTN
   022 {       43 33 } g R⬆
   023 {    42 21  0 } f LBL 0
   024 {       43 20 } g x=0
   025 {       22  2 } GTO 2
   026 {           2 } 2
   027 {          10 } ÷
   028 {          36 } ENTER
   029 {       43 44 } g INT
   030 {    43 30  5 } g TEST x=y
   031 {       22  1 } GTO 1
   032 {    45 16 11 } RCL MATRIX A
   033 {    45 16 12 } RCL MATRIX B
   034 {          20 } ×
   035 {    44 16 11 } STO MATRIX A
   036 {          33 } R⬇
   037 {    42 21  1 } f LBL 1
   038 {    45 16 12 } RCL MATRIX B
   039 {          36 } ENTER
   040 {          20 } ×
   041 {    44 16 12 } STO MATRIX B
   042 {          33 } R⬇
   043 {       22  0 } GTO 0
   044 {    42 21  2 } f LBL 2
   045 {    45 11  u } RCL A
   046 {    45 11  u } RCL A
   047 {       43 32 } g RTN

Example

It takes about 21 seconds to calculate \(F_{49}\).

49 A
7,778,742,049.

[Joe Horn]

(08-06-2023 11:06 AM)Thomas Klemm Wrote:  ,,, For \(F_{49}\) it takes about 5 seconds while Joe's program takes about 30 seconds.

Brilliant solution! On the 15C CE, your program returns \(F_{49}\) immediately.

[Albert Chan]

(03-23-2017 01:32 PM)Eddie W. Shore Wrote:  n = 44
HP15C LE: 701408728.7
HP Prime: 701408733.002

n = 49
HP 15C LE: 7778741992
HP Prime: 7778742049.02

I think we may have found a bug on the HP 15C LE.

√5 ≈ 2.23606797749979

Round to 10 digits, we have huge ulp error, almost 1/2 ULP
(1+√5) relative error = 0.49979E-9 / 3.23607 = 1.54E-10

(1 + ε)^n ≈ 1 + n*ε

F(44) corrections ≈ 44 * 1.54E-10 * 701408728.7 ≈ 4.8
F(44) ≈ 701408728.7 + 4.8 = 701408733.5 (error = -0.5)

F(49) corrections ≈ 49 * 1.54E-10 * 7778741992 ≈ 59
F(49) ≈ 7778741992 + 59 = 7778742051 (error = -2)

---

√20 is better, ulp error = 2 * 0.49979 - 1 = -0.00042 ULP

F(n) ≈ (2+√20)^n / (2^(2n-1)*√20)

On my HP12C:

F(44) = 701,408,733.1 (error -0.1)
F(49) = 7,778,742,050 (error -1)

[Thomas Klemm]

That's nice!
We can use your formula to calculate \(F_{47}\) exactly by rounding.
This allows us to calculate:

\(
F_{93}=F_{47}^2+F_{46}^2
\)

---

Edit: Bummer, the HP-16C has no \(y^x\) or similar function.

[Werner]

(08-07-2023 10:09 PM)Albert Chan Wrote:  F(n) ≈ (2+√20)^n / (2^(2n-1)*√20)

Hello Albert!
This overflows already for n=123.
perhaps better

F(n) ≈ 2*(0.2+√0.2)^n / (0.4^n*√20)

which can also very easily be turned into a program, eg 41-style
(I don't have a 15C.. should arrive any minute now ;-)
can handle up till n=248, correct for n=0..44

Update: snipped older version, see two posts down for new

Cheers, Werner

[Albert Chan]

(08-08-2023 11:05 AM)Werner Wrote:  perhaps better

F(n) ≈ 2*(0.2+√0.2)^n / (0.4^n*√20)

We might as well combine the constant, with *same* relative error as √20

2 / √20 = √20 / 10

F(n) ≈ (0.2+√0.2)^n / 0.4^n * √0.2

Scale √0.2 last, within n = 0 .. 49, rounded, only F(48) was off by 1

[Werner]

Yes, I realized that too. The program now becomes (with 15C instructions only, I think - UPS is late):
(Update: now correct up to n=47. N=48 is off by 1, and N=49 is correct, too)

29 Bytes:

Code:

01 LBL"FIB"
02 .2
03 ENTER
04 SQRT
05 ENTER
06 RDN
07 +
08 X<>Y
09 Y^X
10 .4
11 LASTX
12 Y^X
13 /
14 *
15 ENTER
16 FRC
17 +
18 INT
19 END

Cheers, Werner