The Fibonacci sequence is a very famous sequence which is defined as below: $$ F_0 = 0, \quad F_1 = 1, \quad F_n = F_{n-1} + F_{n-2} \;\; \textrm{for} \; n \geq 2 $$ The first values of the sequence are $0,1,1,2,3,5,8,13,\ldots$

Obtaining $F_n$ in a non-recursive manner

Since the sequence is defined in a recursive manner, computing its $n$-th element without using recursion is not an immediately obvious task. It can, however, be done if we look at the elements of the Fibonacci sequence in pairs: $$ \begin{eqnarray} F_{n+1} &=& F_{n} + F_{n-1} \nonumber \\[5pt] F_{n} &=& F_{n} \label{post_c772250d88e35665d3a793882a7b70e5_fibonacci_pairs} \end{eqnarray} $$ The second equality is trivial, but it is useful as it helps us interpret equation \eqref{post_c772250d88e35665d3a793882a7b70e5_fibonacci_pairs} using matrices as shown below: $$ \left(\begin{matrix} F_{n+1} \\ F_n \end{matrix}\right) = \left(\begin{matrix} F_{n} + F_{n-1} \\ F_{n} \end{matrix}\right) = \left(\begin{matrix} 1 & 1 \\ 1 & 0\end{matrix}\right) \left(\begin{matrix} F_{n} \\ F_{n-1} \end{matrix}\right) \label{post_c772250d88e35665d3a793882a7b70e5_fibonacci_pairs_matrix} $$ Equation \eqref{post_c772250d88e35665d3a793882a7b70e5_fibonacci_pairs_matrix} is very interesting since it yields a recursive formula which associates the pairs $(F_{n+1},F_n)$ and $(F_n,F_{n-1})$. We can use this equation to relate $(F_{n+1},F_n)$ directly to $(F_1,F_0)$. Defining: $$ M = \left(\begin{matrix} 1 & 1 \\ 1 & 0\end{matrix}\right) $$ we have that: $$ \left(\begin{matrix} F_{n+1} \\ F_n \end{matrix}\right) = M\left(\begin{matrix} F_{n} \\ F_{n-1} \end{matrix}\right) = M^2\left(\begin{matrix} F_{n-1} \\ F_{n-2} \end{matrix}\right) = \ldots = M^n\left(\begin{matrix} F_{1} \\ F_{0} \end{matrix}\right) $$ Since $F_0 = 0$ and $F_1 = 1$, we obtain: $$ \boxed{ \left(\begin{matrix} F_{n+1} \\ F_n \end{matrix}\right) = M^n\left(\begin{matrix} 1 \\ 0 \end{matrix}\right) } \label{post_c772250d88e35665d3a793882a7b70e5_eq_Fnp1_Fn_power_of_M} $$ Computing $M^n$ can be done easily if we first diagonalize it ($M$ is symmetric and therefore diagonalizable). Let's first compute its eigenvalues: $$ \det(M - \lambda I) = 0 \Longrightarrow \det\left(\begin{matrix} 1 - \lambda & 1 \\ 1 & - \lambda\end{matrix}\right) = \lambda^2 - \lambda - 1 = 0 \label{post_c772250d88e35665d3a793882a7b70e5_eq_det_M_minus_lambdaI} $$ The solutions of equation \eqref{post_c772250d88e35665d3a793882a7b70e5_eq_det_M_minus_lambdaI} are the eigenvalues of $M$: $$ \displaystyle \lambda_{\pm} = \frac{1 \pm \sqrt{5}}{2} \label{post_c772250d88e35665d3a793882a7b70e5_lambda_pm} $$ The eigenvalue $\lambda_+$ is known as the golden ratio and is usually represented as the Greek letter $\phi$ (we will however use $\lambda_+$ in this discussion since it will be more convenient).

Let's now compute the eigenvectors of $M$ associated with these eigenvalues: $$ \left(\begin{matrix} 1 & 1 \\ 1 & 0\end{matrix}\right) \left(\begin{matrix}\alpha_{\pm} \\ \beta_{\pm} \end{matrix}\right) = \left(\begin{matrix}\alpha_{\pm} + \beta_{\pm}\\ \alpha_{\pm} \end{matrix}\right) = \lambda_{\pm}\left(\begin{matrix}\alpha_{\pm} \\ \beta_{\pm} \end{matrix}\right) \label{post_c772250d88e35665d3a793882a7b70e5_eigenvec_M_general} $$ Equation \eqref{post_c772250d88e35665d3a793882a7b70e5_eigenvec_M_general} implies that the eigenvectors of $M$ satisfy: $$ \alpha_{\pm} = \lambda_{\pm}\beta_{\pm} \label{post_c772250d88e35665d3a793882a7b70e5_alpha_beta_pm} $$ Since we can choose $\beta_+$ and $\beta_-$ arbitrarily, we will take $\beta_+ = \beta_- = 1$ to compute two eigenvectors ${\bf w}_+$ and ${\bf w}_-$ of $M$ associated with $\lambda_+$ and $\lambda_-$ respectively. With these choices of values for $\beta_{\pm}$ and from equation \eqref{post_c772250d88e35665d3a793882a7b70e5_alpha_beta_pm}, we get: $$ {\bf w}_+ = \left(\begin{matrix}\lambda_+ \\ 1 \end{matrix}\right) \quad\quad {\bf w}_- = \left(\begin{matrix}\lambda_- \\ 1 \end{matrix}\right) \label{post_c772250d88e35665d3a793882a7b70e5_eigenvectors_M} $$ NOTE: since $M$ is symmetric, its eigenvectors must be orthogonal. Indeed: $$ {\bf w}_+^T{\bf w}_- = \lambda_+\lambda_- + 1 = \displaystyle \left(\frac{1 + \sqrt{5}}{2}\right)\left(\frac{1 - \sqrt{5}}{2}\right) + 1 =\frac{1 - 5}{4} + 1 = 0 $$

Equation \eqref{post_c772250d88e35665d3a793882a7b70e5_eigenvectors_M} allows us to diagonalize $M$. Letting $P = ({\bf w}_+, {\bf w}_-)$ be the $2 \times 2$ matrix whose columns are ${\bf w}_+$ and ${\bf w}_-$: $$ P = \left(\begin{matrix} \lambda_+ & \lambda_- \\ 1 & 1 \end{matrix}\right) $$ we have that: $$ MP = (M{\bf w}_+, M{\bf w}_-) = (\lambda_+{\bf w}_+, \lambda_-{\bf w}_-) = ({\bf w}_+, {\bf w}_-)\left(\begin{matrix} \lambda_+ & 0 \\ 0 & \lambda_-\end{matrix}\right) = PD $$ and therefore: $$ M = PDP^{-1}, \quad D = \left(\begin{matrix} \lambda_+ & 0 \\ 0 & \lambda_-\end{matrix}\right), \quad P^{-1} = \frac{1}{\lambda_+ - \lambda_-}\left(\begin{matrix} 1 & -\lambda_- \\ -1 & \lambda_+\end{matrix}\right) \label{post_c772250d88e35665d3a793882a7b70e5_diagonalization_M} $$ where $P^{-1}$ was computed using the standard formula for the inverse of $2 \times 2$ matrices ($P$ is invertible since $\det(P) = \lambda_+ - \lambda_- \neq 0$).

Equation \eqref{post_c772250d88e35665d3a793882a7b70e5_diagonalization_M} is the diagonalization of $M$. This equation is very useful since it turns the computation of $M^n$ into a trivial task: $$ M^n = (PDP^{-1})^n = PD^nP^{-1} \label{%INDEX_eq_Mn} $$ where $$ D^n = \left(\begin{matrix} \lambda_+^n & 0 \\ 0 & \lambda_-^n\end{matrix}\right) \label{%INDEX_eq_Dn} $$ Using equations \eqref{%INDEX_eq_Mn} and \eqref{%INDEX_eq_Dn} on equation \eqref{post_c772250d88e35665d3a793882a7b70e5_eq_Fnp1_Fn_power_of_M}, we obtain: $$ \left(\begin{matrix} F_{n+1} \\ F_n \end{matrix}\right) = PD^nP^{-1}\left(\begin{matrix} 1 \\ 0 \end{matrix}\right) = \frac{1}{\lambda_+ - \lambda_-}\left(\begin{matrix} \lambda_+^{n+1} - \lambda_-^{n+1} \\ \lambda_+^{n} - \lambda_-^{n} \end{matrix}\right) \label{post_c772250d88e35665d3a793882a7b70e5_Fnp1_and_Fn_final} $$ Using the values of $\lambda_{\pm}$ given on equation \eqref{post_c772250d88e35665d3a793882a7b70e5_lambda_pm}, we obtain the final expression for $F_n$: $$ \boxed{ F_n = \displaystyle\frac{1}{\sqrt{5}}\left[\left(\frac{1 + \sqrt{5}}{2}\right)^n - \left(\frac{1 - \sqrt{5}}{2}\right)^n\right] } \label{post_c772250d88e35665d3a793882a7b70e5_final_Fn} $$

The golden ratio

Since $F_{n+1} = F_{n} + F_{n-1} \geq F_{n} + 1$ for all $n \geq 2$, the sequence $F_n$ diverges as $n \rightarrow \infty$. However, the sequence: $$ a_n = \displaystyle\frac{F_{n+1}}{F_n} \;\; \textrm{for} \; n \geq 1 $$ does not. As a matter of fact, it converges to the the golden ratio. Proving this fact is easy (and there are many ways to do it). From equation \eqref{post_c772250d88e35665d3a793882a7b70e5_Fnp1_and_Fn_final}, we have that: $$ \displaystyle\frac{F_{n+1}}{F_n} = \frac{\lambda_+^{n+1} - \lambda_-^{n+1}}{\lambda_+^{n} - \lambda_-^{n}} = \lambda_+\left(\frac{1 - (\lambda_-/\lambda_+)^{n+1}} {1 - (\lambda_-/\lambda_+)^{n}}\right) $$ Since $|\lambda_+| \approx 1.62 \gt |\lambda_-| \approx 0.62$, both $(\lambda_-/\lambda_+)^{n+1}$ and $(\lambda_-/\lambda_+)^{n}$ vanish as $n \rightarrow \infty$. This implies that: $$ \boxed{ \lim_{n \rightarrow \infty}\displaystyle\frac{F_{n+1}}{F_n} = \lambda_+ = \displaystyle\frac{1 + \sqrt{5}}{2} = 1.61803\ldots } $$

Computing $F_n$ for large values of $n$

Although equation \eqref{post_c772250d88e35665d3a793882a7b70e5_final_Fn} is mathematically correct, it is not adequate for computing $F_n$ using a computer. The problem is that the evaluation of $F_n$ will always have issues due to the accuracy errors which are inherent in floating point arithmetic.

One can, however, use arbitrary precision arithmetic to compute $F_n$ using equation \eqref{post_c772250d88e35665d3a793882a7b70e5_eq_Fnp1_Fn_power_of_M} and exponentiation by squaring. A Python script which does exactly that can be found here.

Notice, however, that the time complexity of the method just described is not $O(\log_2 n)$! In spite of the fact that the time complexity of exponentiation by squaring is logarithmic, the use of arbitrary precision makes the time complexity of the overall algorithm become $O(n)$. This is because $\textrm{digits}(F_n)$, the number of digits in $F_n$, is proportional to $n$ as $n$ becomes large and arithmetic with arbitrary precision involves operations which are $O(n)$ or worse if $n$ is the typical length of the numbers involved (consider adding two large numbers with $n$ digits: the time complexity of that operation is $O(n)$). To see why $\textrm{digits}(F_n)$ grows linearly with $n$, notice first that since $|\lambda_+| \gt |\lambda_-|$, $F_n \approx \lambda_+^n / \sqrt{5}$ for large values of $n$. Therefore, when $n$ becomes large, we have that: $$ \textrm{digits}(F_n) \approx \log_{10}(F_n) \approx \log_{10}(\lambda_+^n) - \log_{10}(\sqrt{5}) \approx n\log_{10}(\lambda_+) \propto n $$ To better illustrate the time complexity of computing $F_n$ using exponentiation by squaring and arbitrary precision arithmetic, I used the Python script mentioned above to compute $F_n$ for $n = 2^i$, $i = 20, 21, \ldots, 30$ (those are quite large values of $n$). Figure 1 shows that the time taken for computing $F_n$ grows linearly with $n$.