(50g) Nth Fibonacci Number

[Eddie W. Shore]

nth Fibonaaci Number. n is on Level 1 of the stack.

Code:

<< → N
<< 1 5 √ + 2 / N
^ 1 5 √ - 2 / N ^ 
- 5 √ / EVAL 0 RND >> >>

[Gerald H]

Another way to do it for integer input:

::
CK1&Dispatch
# FF
::
::
DUP
Z2_
Z<
case
::
BINT4
NDUPN
DROP
;
FPTR2 ^Z>ZH
Z-2_
Z0Z1_
4PICK
FPTR2 ^ZBits
SWAPDROP#1-_
ZERO_DO
FPTR2 ^ZSQ_
SWAP
FPTR2 ^ZSQ_
2DUP
FPTR2 ^RSUBext
DUP
FPTR2 ^RADDext
3PICK
FPTR2 ^RADDext
3UNROLL
FPTR2 ^RADDext
SWAPROT
FPTR2 ^RADDext
3PICK
ISTOP-INDEX
#1-
FPTR2 ^ZBit?
SWAPDROP
ITE
::
DUPUNROT
FPTR2 ^QAdd
Z-2_
;
Z2_
3UNROLL
LOOP
;
4UNROLL3DROP
;
;

[Joe Horn]

Here's my favorite exact Nth Fibonacci Number program for the 50g:

Code:

<<
  [[ 1 1 ]
   [ 1 0 ]] SWAP ^ 2. GET
>>

BYTES: 45.0 #CD5Ah

Exact mode, of course. Short & sweet. Pretty fast for a URPL program: 512th Fibonacci number (107 digits) in 0.75 seconds. 1024th Fibonacci number (214 digits) in 1.63 seconds.

[Offroad]

(02-22-2015 09:06 PM)Joe Horn Wrote:  Here's my favorite exact Nth Fibonacci Number program for the 50g:

Code:

<<
  [[ 1 1 ]
   [ 1 0 ]] SWAP ^ 2. GET
>>

BYTES: 45.0 #CD5Ah

Exact mode, of course. Short & sweet. Pretty fast for a URPL program: 512th Fibonacci number (107 digits) in 0.75 seconds. 1024th Fibonacci number (214 digits) in 1.63 seconds.

Where's the LIKE button?

[rprosperi]

(02-22-2015 09:06 PM)Joe Horn Wrote:  Here's my favorite exact Nth Fibonacci Number program for the 50g:

Code:

<<
  [[ 1 1 ]
   [ 1 0 ]] SWAP ^ 2. GET
>>

BYTES: 45.0 #CD5Ah

Exact mode, of course. Short & sweet. Pretty fast for a URPL program: 512th Fibonacci number (107 digits) in 0.75 seconds. 1024th Fibonacci number (214 digits) in 1.63 seconds.

Wow, completely unexpected. So, can you explain that? I follow the code, trivial as it is, so I guess I don't know the math behind the matrix power operation.

[Han]

(02-26-2015 01:42 PM)rprosperi Wrote:  

(02-22-2015 09:06 PM)Joe Horn Wrote:  Here's my favorite exact Nth Fibonacci Number program for the 50g:

Code:

<<
  [[ 1 1 ]
   [ 1 0 ]] SWAP ^ 2. GET
>>

BYTES: 45.0 #CD5Ah

Exact mode, of course. Short & sweet. Pretty fast for a URPL program: 512th Fibonacci number (107 digits) in 0.75 seconds. 1024th Fibonacci number (214 digits) in 1.63 seconds.

Wow, completely unexpected. So, can you explain that? I follow the code, trivial as it is, so I guess I don't know the math behind the matrix power operation.

Very nice! The way it works is that matrix multiplication is essentially a dot product of row i (of the left matrix) with column j (of the right matrix) to form the i,j entry in the resulting matrix.
\[ \left[ \begin{array}{cc} a & b \\
c & d \end{array} \right] \cdot \left[ \begin{array}{cc} r & s \\
t & u \end{array} \right] = \left[ \begin{array}{cc} ar+bt & as+bu \\
cr+dt & cs+du \end{array} \right] \]

So now let
\[ A = \left[ \begin{array}{cc} F_2 & F_1 \\
F_1 & F_0 \end{array} \right] \]
where \( F_0 = 0\) and \(F_1 = 1 \). Notice \( F_2 = F_1 + F_0 = 1 \).
\[ A^2 = A \cdot \left[ \begin{array}{cc} 1 & 1 \\
1 & 0 \end{array} \right] = \left[ \begin{array}{cc} F_2 & F_1 \\
F_1 & F_0 \end{array} \right] \cdot \left[ \begin{array}{cc} 1 & 1 \\
1 & 0 \end{array} \right] =
\left[ \begin{array}{cc} F_1+F_2 & F_2+0 \\
F_1+F_0 & F_1+0 \end{array} \right]
= \left[ \begin{array}{cc} F_3 & F_2 \\
F_2 & F_1 \end{array} \right]\]
\[ A^3 = A^2 \cdot A = \left[ \begin{array}{cc} F_3 & F_2 \\
F_2 & F_1 \end{array} \right] \cdot \left[ \begin{array}{cc} 1 & 1 \\
1 & 0 \end{array} \right] =
\left[ \begin{array}{cc} F_4 & F_3 \\
F_3 & F_2 \end{array} \right] \]

(It also works if you prefer to do \( A^n = A \cdot A^{n-1} \) since \( A \) and its powers are symmetric.)

Very neat indeed.

[rprosperi]

(02-26-2015 07:39 PM)Han Wrote:  Very neat indeed.

Indeed! Brilliant!

And very nicely illustrated and explained Han, thanks for taking the time to walk thru that. Me thinks I'm not the only one who didn't see this. And like so many clever solutions, the math itself is indeed simple, but the insights to approach it this way is remarkable.

Joe - Do you recall who crafted such an insightful solution? Am I honoring you?

[Joe Horn]

(02-26-2015 08:23 PM)rprosperi Wrote:  Joe - Do you recall who crafted such an insightful solution?

No, I don't know who found it first, but I enjoyed serendipitously rediscovering it while playing with numbers on HP calculators, and immediately placed it into my personal mathematical petting zoo.

[Thomas Klemm]

This is the diagonalization of this matrix:



When \(M\) is squared the inner product of \(S^{-1} \cdot S\) cancel:
\(M^2=(S \cdot J \cdot S^{-1})^2=S \cdot J \cdot S^{-1} \cdot S \cdot J \cdot S^{-1}=S \cdot J \cdot J \cdot S^{-1}=S \cdot J^2 \cdot S^{-1}\)
This can be generalized for all powers so that we end up with:
\(M^n=S \cdot J^n \cdot S^{-1}\)
But the power of a diagonalized matrix consists just of the power of its elements on the diagonal (i.e. the eigenvalues).
This brings us back to Eddie's initial solution. Okay, maybe after a little exercise in algebra.

Cheers
Thomas

[Thomas Klemm]

(02-22-2015 09:06 PM)Joe Horn Wrote:  

Code:

<<
  [[ 1 1 ]
   [ 1 0 ]] SWAP ^ 2. GET
>>

Gerald H posted recently a program to calculate the Matrix Integer Powers which can be used to calculate this with the HP 42S.

Cheers
Thomas

[Han]

Another way to produce the constants obtained in Eddie's solution is to use formal series and generating functions to solve recurrence equations. Let \( a_n \) denote the \( n \)th term, and let
\[ f(x) = \sum_{n=0}^\infty a_n x^n \]
\[ a_{n+2} = a_{n+1} + a_{n} \quad \Longrightarrow \quad
a_{n+2}\cdot x^{n+2} = x \cdot a_{n+1}\cdot x^{n+1} + x^2 \cdot a_n \cdot x^n \]
Summing each side with \( n = 0, 1, 2, \dotsm \) produces
\[ \sum_{n=0}^\infty a_{n+2} x^{n+2} = x \sum_{n=0}^\infty a_{n+1}x^{n+1}
+ x^2 \sum_{n=0}^\infty a_n x^n \quad
\Longrightarrow \quad f(x)-a_1 x - a_0 = x \cdot [ f(x)-a_0 ] + x^2f(x) \]
Solve for \( f(x) \) to obtain
\[ f(x) = \frac{a_1 x -a_0x + a_0}{1-x-x^2} = \frac{1}{1-x-x^2} \]
In the last equality, we (for simplicity's sake) assume that \( a_0 = a_1 = 1 \). Using a partial fraction decomposition with \( \alpha = -\frac{1}{2}(1-\sqrt{5}) \) and \( \beta = -\frac{1}{2}(1+\sqrt{5}) \) gives
\[ f(x) = \frac{-1}{(x-\alpha)(x-\beta)} = \frac{A}{x-\alpha} + \frac{B}{x-\beta} \]
or
\[ -1 = A\cdot (x-\beta) + B\cdot (x-\alpha) \]
Letting \( x= \alpha \) and then \( x=\beta \) produces \(A = -\frac{1}{\alpha-\beta} = -\frac{1}{\sqrt{5}} \) and \( B = \frac{1}{ \sqrt{5}} \) respectively. Now rewrite \( f(x) \) as a sum of power series:
\[ f(x) = \frac{1}{\sqrt{5}} \left( -\frac{1}{x-\alpha} + \frac{1}{x-\beta} \right)
= \frac{1}{\sqrt{5}}\left( \frac{1}{\alpha-x} - \frac{1}{\beta-x}\right)
= \frac{1}{\sqrt{5}} \left( \frac{1}{\alpha\cdot (1-\frac{x}{\alpha})} -
\frac{1}{\beta \cdot (1-\frac{x}{\beta})} \right)
= \frac{1}{\sqrt{5}} \left( \frac{1}{\alpha} \sum_{n=0}^\infty \left(\frac{x}{\alpha}\right)^n - \frac{1}{\beta} \sum_{n=0}^\infty \left( \frac{x}{\beta}\right)^n \right) \]
by making use of the identity
\[ \frac{1}{1-t} = \sum_{n=0}^\infty t^n \]
So
\[ f(x) = \sum_{n=0}^\infty \underbrace{\frac{1}{\sqrt{5}}\left( \frac{1}{\alpha^{n+1}}-\frac{1}{\beta^{n+1}}\right)}_{\text{ this is } a_n} x^n \]
The exponent \( n+1 \) (as opposed to \( n \) ) is due to us setting \( a_0 = 1 \) and not \( a_0 = 0 \).

[Thomas Klemm]

(02-26-2015 01:42 PM)rprosperi Wrote:  So, can you explain that?

For \(x_{n+2}=x_n+x_{n+1}\) you can define a 2nd sequence \(y_n=x_{n+1}\).
Thus the relation becomes:

\(y_{n+1}=x_{n+2}=x_n+y_n\).

Now you combine both sequences into a vector:

\(z_n=\begin{bmatrix}
x_n \\
y_n
\end{bmatrix}\)

This allows to merge both equations into a single one:

\(z_{n+1}=\begin{bmatrix}
x_{n+1} \\
y_{n+1}
\end{bmatrix}=\begin{bmatrix}
y_n \\
x_n+y_n
\end{bmatrix}=\begin{bmatrix}
0 & 1 \\
1 & 1
\end{bmatrix}\cdot\begin{bmatrix}
x_n \\
y_n
\end{bmatrix}=M\cdot z_n\)

But this is just the recursive definition of a geometric sequence leading to:

\(z_n=M^n\cdot z_0\) with \(z_0=\begin{bmatrix}
1 \\
0
\end{bmatrix}\)

HTH
Thomas

[Thomas Klemm]

A poor man's approach to solve the recurrence would be to use the ansatz

\(x_n=\alpha\cdot\phi^n\)

Insert this into the recursive definition \(x_{n+2}=x_{n+1}+x_n\):

\(\alpha\cdot\phi^{n+2}=\alpha\cdot\phi^{n+1}+\alpha\cdot\phi^n\)

Divide both sides by \(\alpha\cdot\phi^n\):

\(\phi^2=\phi+1\)

This is the quadratic equation for the golden ratio which has two solutions: \(\phi_1,\phi_2\)

Now use the linear combination of both:

\(x_n=\alpha_1\cdot\phi_1^n+\alpha_2\cdot\phi_2^n\)

The initial values \(x_0=0\) and \(x_1=1\) can be used to determine the unknown factors \(\alpha_1\) and \(\alpha_2\).

Definitely not as cool as using generating functions and partial fraction decomposition and stuff.

Cheers
Thomas

PS: I feel a little bit like cheating as there's no explanation for why this ansatz works.

[rprosperi]

(02-27-2015 03:29 AM)Han Wrote:  Another way to produce the constants obtained in Eddie's solution is to use formal series and generating functions to solve recurrence equations....
[snip]

So, you are a Math Teacher/Professor !

I have long suspected, but this removes all doubt.

My head hurts just following this, never mind composing it, but it seemed the least I could do. Thanks for the derivation. Reading Eddie's original did indeed lead me to wonder where the heck those constants were from, but honestly not thinking I would ever actually know.

Thanks Han.

[rprosperi]

(02-27-2015 05:10 AM)Thomas Klemm Wrote:  

(02-26-2015 01:42 PM)rprosperi Wrote:  So, can you explain that?

For \(x_{n+2}=x_n+x_{n+1}\) you can define a 2nd sequence \(y_n=x_{n+1}\).
Thus the relation becomes:

\(y_{n+1}=x_{n+2}=x_n+y_n\).

Now you combine both sequences into a vector:

\(z_n=\begin{bmatrix}
x_n \\
y_n
\end{bmatrix}\)

This allows to merge both equations into a single one:

\(z_{n+1}=\begin{bmatrix}
x_{n+1} \\
y_{n+1}
\end{bmatrix}=\begin{bmatrix}
y_n \\
x_n+y_n
\end{bmatrix}=\begin{bmatrix}
0 & 1 \\
1 & 1
\end{bmatrix}\cdot\begin{bmatrix}
x_n \\
y_n
\end{bmatrix}=M\cdot z_n\)

But this is just the recursive definition of a geometric sequence leading to:

\(z_n=M^n\cdot z_0\) with \(z_0=\begin{bmatrix}
1 \\
0
\end{bmatrix}\)

HTH
Thomas

Another one! There are mathematicians everywhere I look here. And thank god for that, as us mere mortals need folks like you to explain stuff like this so we can go about implementing the algorithms on our machines, without having to derive the equations ourselves.

More seriously, thanks Thomas for these nice derivations. Your patient step-by-step approach is easy to follow and instructional. I certainly make no claim about having any of the inspiration behind the derivation, which is of course where the magic lies, but its still interesting to follow along.

[rprosperi]

(02-27-2015 05:58 AM)Thomas Klemm Wrote:  A poor man's approach to solve the recurrence would be to use the ansatz
[snip]
PS: I feel a little bit like cheating as there's no explanation for why this ansatz works.

So it appears an ansatz is a postulated theorem for some behavior which works, but without knowing why?

[Han]

(02-27-2015 01:48 PM)rprosperi Wrote:  

(02-27-2015 05:58 AM)Thomas Klemm Wrote:  A poor man's approach to solve the recurrence would be to use the ansatz
[snip]
PS: I feel a little bit like cheating as there's no explanation for why this ansatz works.

So it appears an ansatz is a postulated theorem for some behavior which works, but without knowing why?

My suspicion is that it is tied to the fact that the Fibonacci sequence is a linear recurrence. All linear recurrences will have rational generating functions (i.e. the \( f(x) \) in my earlier post will always be of the form \( \frac{p(x)}{q(x)} \) for polynomials \( p(x) \) and \( q(x) \)). In the case of lower order recurrences, one can often obtain a partial fraction decomposition of the rational function which then leads to solutions of the aforementioned forms.

[Gerald H]

(02-27-2015 01:48 PM)rprosperi Wrote:  

(02-27-2015 05:58 AM)Thomas Klemm Wrote:  A poor man's approach to solve the recurrence would be to use the ansatz
[snip]
PS: I feel a little bit like cheating as there's no explanation for why this ansatz works.

So it appears an ansatz is a postulated theorem for some behavior which works, but without knowing why?

A regular translation of Ansatz is English "approach" to a problem. It could be a first approximation.

[Thomas Klemm]

(02-27-2015 01:48 PM)rprosperi Wrote:  So it appears an ansatz is a postulated theorem for some behavior which works, but without knowing why?

It's easy to verify that the sequence \(x_n=\alpha\cdot\phi^n\) satisfies the equation \(x_{n+2}=x_{n+1}+x_n\): that's exactly how we calculate \(\phi\).

However just from looking at the recurrence it's not obvious to use a geometric sequence.
In this case the ansatz is a hint that allows you to solve this problem without too much linear algebra. In other cases, it might be a trick that has been proved successfully elsewhere.

My discomfort was related to the fact that I probably couldn't motivate the ansatz enough.
However you seem to estimate the post all the same.

Cheers
Thomas