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