Faster Computation of Fibonacci Numbers

Summer 2024 · CS691-XVIII · Applications of Linear Algebra

Course Project loosely based on Miniature 1 and 2 from the book “Thirty-three Miniatures: Mathematical and Algorithmic Applications of Linear Algebra” by Jiřì Matoušek
By Sourabh Warrier & Cherian George

Consider the following sequence, Sn=a1Snk+a2Snk+1+...+akSn1aiRS_n = a_1S_{n-k} + a_2S_{n-k+1} +... + a_kS_{n-1} \ni a_i \in \mathbb{R}

Any term is a linear combination of the previous kk terms and first kk terms S1S_1 to SkS_k are given by the base cases b1b_1 to bkb_k of the recursion. Given a set of kk terms in such a sequence, the process of obtaining subsequent terms can be represented as linear transformation of vectors in Rk\mathbb{R}^k. Consider the following vectors,

v0=[Si+k1Si+k2...Si]v_0 = \begin{bmatrix} S_{i+k-1} \\ S_{i+k-2}\\.\\.\\.\\ S_{i} \end{bmatrix} and v1=[Si+kSi+k1...Si+1]v1=Av0v_1 = \begin{bmatrix} S_{i+k} \\ S_{i+k-1}\\.\\.\\.\\ S_{i+1} \end{bmatrix} \ni v_1 = Av_0, where AA is a linear transformation. The entries of AA would depend on the coefficients a1a_1 to aka_k specific to the sequence. In this case A=[akak1ak2...a11...000...0...1000...010]A = \begin{bmatrix} a_k & a_{k-1} & a_{k-2} & ... & a_1 \\ 1 & ...& 0 & 0 & 0 \\.\\.\\.\\ 0 & ...& 1 & 0 & 0 \\ 0 & ...& 0 & 1 & 0 \end{bmatrix}

Eigenbasis of AA

The eigenvalues of AA are the roots of the polynomial P(λ)=AλI=0P(\lambda) = |A - \lambda{I}| = 0. Suppose by some means, eigenvalues λ1\lambda_1 to λk\lambda_k and their corresponding eigenvectors ϕ1\phi_1 to ϕk\phi_k could be obtained. The transformation AA can be written as PDP1PDP^{-1}, where P=[ϕ1ϕ2...ϕk]P = \begin{bmatrix} \phi_1 & \phi_2 & ... & \phi_k \end{bmatrix} and D=[λ10...00λ2...0...00...λk]D = \begin{bmatrix} \lambda_1 & 0 & ... & 0 \\ 0 & \lambda_2 & ... & 0 \\ .\\.\\.\\ 0 & 0 & ... & \lambda_k \end{bmatrix}. This gives us vn=Anv0=PDnP1v0v_n=A^nv_0 = PD^nP^{-1}v_0.

Examples

{Sn=Sn2+Sn1,{S1=S2=1}S_n = S_{n-2} + S_{n-1}, \hspace{10 mm}\{S_1 = S_2 = 1\}}

1, 1, 2, 3, 5, 8, 13, 21, 34, 55, … The first example we’ll look at is the Fibonacci sequence (OEIS \href{https://oeis.org/A000045}{A000045}). From the definition, we have a1=a2=1a_1 = a_2 = 1 and S1=S2=1S_1 = S_2= 1. Therefore, A=[1110]A = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}, v0=[11]v_0 = \begin{bmatrix} 1 \\ 1 \end{bmatrix} and P(λ)=AλI=0    λ2λ1=0P(\lambda) = |A - \lambda{I} = 0| \implies \lambda^2 - \lambda -1 = 0, giving us λ1=1+52\lambda_1 = \frac{1+\sqrt{5}}{2} and λ2=152\lambda_2 = \frac{1-\sqrt{5}}{2}. With the augmented matrix [1λ101λ0]\left[\hspace{-5pt}\begin{array}{cc|c} 1-\lambda & 1 & 0 \\ 1 & -\lambda & 0 \end{array}\hspace{-5pt}\right], written in REF as [1λ100λ2λ11λ0]    ϕ=[1λ1]\left[\hspace{-5pt}\begin{array}{cc|c} 1-\lambda & 1 & 0 \\ 0 & \frac{\lambda^2 - \lambda-1}{1-\lambda} & 0 \end{array}\hspace{-5pt}\right] \implies \phi = \begin{bmatrix} 1 \\ \lambda-1 \end{bmatrix}. Plugging in λ1\lambda_1 and λ2\lambda_2 we find that ϕ1=[1λ11]\phi_1 = \begin{bmatrix} 1 \\ \lambda_1 -1 \end{bmatrix} and ϕ2=[1λ21]\phi_2 = \begin{bmatrix} 1 \\ \lambda_2 -1 \end{bmatrix}. Here, P=[11λ11λ21]P = \begin{bmatrix} 1 & 1 \\ \lambda_1 -1 & \lambda_2 -1 \end{bmatrix}, D=[λ100λ2]D = \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix} and P1=1λ2λ1[λ2111λ11]P^{-1} = \frac{1}{\lambda_2-\lambda_1}\begin{bmatrix} \lambda_2 -1 & -1 \\ 1-\lambda_1 & 1 \end{bmatrix}. Since vn=Anv0=PDnP1v0v_n=A^nv_0 = PD^nP^{-1}v_0, we look at the first coordinate of PDn2P1v0PD^{n-2}P^{-1}v_0 to obtain the nthn^{th} term of the sequence. This gives us [SnSn1]=1λ2λ1[11λ11λ21][λ1n200λ2n2][λ2111λ11][11]=1λ2λ1[λ1n2(λ22)+λ2n2(2λ1)λ1n2(λ11)(λ22)+λ2n2(λ21)(2λ1)]\begin{bmatrix} S_n \\ S_{n-1} \end{bmatrix}=\frac{1}{\lambda_2-\lambda_1}\begin{bmatrix} 1 & 1 \\ \lambda_1 -1 & \lambda_2 -1 \end{bmatrix} \begin{bmatrix} \lambda_1^{n-2} & 0 \\ 0 & \lambda_2^{n-2} \end{bmatrix}\begin{bmatrix} \lambda_2 -1 & -1 \\ 1-\lambda_1 & 1 \end{bmatrix}\begin{bmatrix} 1 \\ 1 \end{bmatrix} = \frac{1}{\lambda_2-\lambda_1}\begin{bmatrix} \lambda_1^{n-2}(\lambda_2-2) + \lambda_2^{n-2}(2-\lambda_1)\\ \lambda_1^{n-2}(\lambda_1-1)(\lambda_2-2) + \lambda_2^{n-2}(\lambda_2-1)(2-\lambda_1) \end{bmatrix}. Extracting SnS_n, from where Sn=1λ2λ1(λ1n2(λ22)+λ2n2(2λ1))=1λ1λ2(λ1n2(λ1+1)λ2n2(λ2+1)){λ2=1λ1}S_n = \frac{1}{\lambda_2-\lambda_1}(\lambda_1^{n-2}(\lambda_2-2) + \lambda_2^{n-2}(2-\lambda_1)) = \frac{1}{\lambda_1-\lambda_2}(\lambda_1^{n-2}(\lambda_1+1)-\lambda_2^{n-2}(\lambda_2+1)) \{\because \lambda2 = 1-\lambda_1\}. Since from the characteristic equation, λ+1=λ2\lambda +1 = \lambda^2, this simplifies to 1λ1λ2(λ1nλ2n)\frac{1}{\lambda_1-\lambda_2}(\lambda_1^n-\lambda_2^n). Plugging in the values of λ1\lambda_1 and λ2\lambda_2 produces Sn=15(1+52)n15(152)nS_n=\frac{1}{\sqrt{5}}\left(\frac{1+\sqrt{5}}{2}\right)^n - \frac{1}{\sqrt{5}}\left(\frac{1-\sqrt{5}}{2}\right)^n

{Sn=Sn2+2Sn1,{S1=S2=1}S_n = S_{n-2} + 2S_{n-1}, \hspace{10 mm}\{S_1 = S_2 = 1\}}

1,1,3,7,17,41,99,239,577,1393,...1, 1, 3, 7, 17, 41, 99, 239, 577, 1393, ...

This is sequence \href{https://oeis.org/A001333}{A001333} in the OEIS. Each term is the sum of twice the previous term and once the term before that. The base cases are identical to that of the Fibonacci sequence. The transformation for this sequence is A=[2110]A = \begin{bmatrix} 2 & 1 \\ 1 & 0 \end{bmatrix} and the characteristic polynomial λ22λ1=0\lambda^2-2\lambda -1 = 0. The eigenvalues and eigenvectors of AA are λ1=2+82\lambda_1 = \frac{2 + \sqrt{8}}{2}, λ2=282\lambda_2 = \frac{2 - \sqrt{8}}{2}, ϕ1=[1λ12]\phi_1 = \begin{bmatrix} 1 \\ \lambda_1-2 \end{bmatrix} and ϕ2=[1λ22]\phi_2 = \begin{bmatrix} 1 \\ \lambda_2-2 \end{bmatrix}.

P=[11λ12λ22]P = \begin{bmatrix} 1 & 1 \\ \lambda_1 -2 & \lambda_2 -2 \end{bmatrix}, D=[λ100λ2]D = \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix} and P1=1λ2λ1[λ2212λ11]P^{-1} = \frac{1}{\lambda_2-\lambda_1}\begin{bmatrix} \lambda_2 -2 & -1 \\ 2-\lambda_1 & 1 \end{bmatrix} and [SnSn1]=1λ2λ1[11λ12λ22][λ1n200λ2n2][λ2212λ11][11]    Sn=1λ2λ1(λ1n2(λ23)+λ2n2(3λ1))=1λ1λ2(λ1n2(3λ2)λ2n2(3λ1))\begin{bmatrix} S_n \\ S_{n-1} \end{bmatrix}=\frac{1}{\lambda_2-\lambda_1}\begin{bmatrix} 1 & 1 \\ \lambda_1 -2 & \lambda_2 -2 \end{bmatrix} \begin{bmatrix} \lambda_1^{n-2} & 0 \\ 0 & \lambda_2^{n-2} \end{bmatrix}\begin{bmatrix} \lambda_2 -2 & -1 \\ 2-\lambda_1 & 1 \end{bmatrix}\begin{bmatrix} 1 \\ 1 \end{bmatrix} \implies S_n = \frac{1}{\lambda_2-\lambda_1}(\lambda_1^{n-2}(\lambda_2-3)+\lambda_2^{n-2}(3-\lambda_1)) = \frac{1}{\lambda_1-\lambda_2}(\lambda_1^{n-2}(3-\lambda_2)-\lambda_2^{n-2}(3-\lambda_1)), which simplifies to the formula Sn=12+8(1+2)n+128(12)nS_n = \frac{1}{2+\sqrt{8}}\left(1+\sqrt{2}\right)^{n} + \frac{1}{2-\sqrt{8}}\left(1-\sqrt{2}\right)^{n}

{General solution for k = 2}

Sn=a1Sn2+a2Sn1,{S1=b1,S2=b2}S_n = a_1S_{n-2} + a_2S_{n-1}, \hspace{10 mm}\{S_1 = b_1,S_2 = b_2\}

In the general case when the nthn^{th} term is a linear combination of the previous two terms with coefficients a1a_1 and a2a_2 and base cases b1b_1 and b2b_2, we represent the corresponding transformation by A=[a2a110]A = \begin{bmatrix} a_2 & a_1 \\ 1 & 0 \end{bmatrix}. The eigenvalues of AA are obtained by solving λ2a2λa1=0\lambda^2-a_2\lambda - a_1 = 0, from where λ1=a2+a22+4a12\lambda_1 = \frac{a_2 + \sqrt{a_2^2+4a_1}}{2} and λ2=a2a22+4a12\lambda_2 = \frac{a_2 - \sqrt{a_2^2+4a_1}}{2}. Solving [a2λa101λ0][a2λa100λ2a2λa1a2λ0]    ϕ1=[1λ1a2a1]\left[\hspace{-5pt}\begin{array}{cc|c} a_2-\lambda & a_1 & 0 \\ 1 & -\lambda & 0 \end{array}\hspace{-5pt}\right] \rightsquigarrow \left[\hspace{-5pt}\begin{array}{cc|c} a_2-\lambda & a_1 & 0 \\ 0 & \frac{\lambda^2-a_2\lambda-a_1}{a_2-\lambda} & 0 \end{array}\hspace{-5pt}\right] \implies \phi_1 = \begin{bmatrix} 1 \\ \frac{\lambda_1-a_2}{a_1} \end{bmatrix} and ϕ2=[1λ2a2a1]\phi_2 = \begin{bmatrix} 1 \\ \frac{\lambda_2-a_2}{a_1} \end{bmatrix}. P=[11λ1a2a1λ2a2a1]P = \begin{bmatrix} 1 & 1 \\ \frac{\lambda_1-a_2}{a_1} & \frac{\lambda_2-a_2}{a_1} \end{bmatrix}, D=[λ100λ2]D = \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix}, and P1=a1λ2λ1[λ2a2a11a2λ1a11]P^{-1} = \frac{a_1}{\lambda_2-\lambda_1}\begin{bmatrix} \frac{\lambda_2-a_2}{a_1} & -1 \\ \frac{a_2-\lambda_1}{a_1} & 1 \end{bmatrix}. We obtain vn2v_{n-2} from the expression

[SnSn1]=a1λ2λ1[11λ1a2a1λ2a2a1][λ1n200λ2n2][λ2a2a11a2λ1a11][b2b1]\begin{bmatrix} S_n \\ S_{n-1} \end{bmatrix}=\frac{a_1}{\lambda_2-\lambda_1}\begin{bmatrix} 1 & 1 \\ \frac{\lambda_1-a_2}{a_1} & \frac{\lambda_2-a_2}{a_1} \end{bmatrix} \begin{bmatrix} \lambda_1^{n-2} & 0 \\ 0 & \lambda_2^{n-2} \end{bmatrix}\begin{bmatrix} \frac{\lambda_2-a_2}{a_1} & -1 \\ \frac{a_2-\lambda_1}{a_1} & 1 \end{bmatrix}\begin{bmatrix} b_2 \\ b_1 \end{bmatrix}, from where Sn=a1λ2λ1(λ1n2(b2(λ2a2)a1b1)+λ2n2(b2(a2λ1)a1+b1))S_n = \frac{a_1}{\lambda_2-\lambda_1}\left(\lambda_1^{n-2}\left(\frac{b_2(\lambda_2-a_2)}{a_1} - b_1\right) + \lambda_2^{n-2}\left(\frac{b_2(a_2-\lambda_1)}{a_1} + b_1\right)\right)

=1λ1λ2(λ1n2(a1b1+b2λ1)λ2n2(a1b1+b2λ2))= \frac{1}{\lambda_1-\lambda_2}\left(\lambda_1^{n-2}\left(a_1b_1+b_2\lambda_1\right) - \lambda_2^{n-2}\left(a_1b_1+b_2\lambda_2\right)\right). This produces the following formula TBC…

Written on June 22, 2024