斐波那契 python 斐波那契数列通项

目录：

• 一、何为斐波那契数列？
• 二、解法一
• 三、解法二
• 四、合二为一
• 五、实际实现

一、何为斐波那契数列？

$1,1,2,3,5,8,13,\cdots$

具体定义如下

$\begin{cases} f(n)=f(n-1)+f(n-2),\quad n \geq 3\\ f(1)=f(2)=1 \end{cases}$

二、解法一

记得高中时，数学老师曾让我尝试一下这个问题. 奈何苦思良久没有思路，前段时间算法课上重遇，解法自己从脑子中蹦出来了，也算是一种缘分吧！

对于下面这种形式的递推公式
$f(n)=af(n-1)+bf(n-2) \quad(b\neq0)$

有一种通用的解法：看着跟等比数列有点儿像，所以想办法 构造出一个等比数列.

两边同时减去 $kf(n-1)$，原式变为
$f(n)-kf(n-1)=(a-k)f(n-1)+bf(n-2) \tag{*}$

将左端看作是等比数列 $\{a_n\}$ 中的 $a_n$，右端看作是 $(a-k)\,a_{n-1}$，则有对应系数成比例
$\frac{1}{a-k}=\frac{-k}{b}$

整理得
$k^2-ak-b=0$

假设该方程的根已求出，为 $k_1,k_2$，由韦达定理
$a=k_1+k_2,\,\,b=-k_1k_2$

将 $k=k_1$ 和上式代入 $(*)$ 式，得
$f(n)-k_1f(n-1)=k_2\big(f(n-1)-k_1f(n-2)\big)$

则 $a_n=k_2 a_{n-1}$，所以 $a_n=k_2^{n-2}a_2$.
$\begin{aligned} f(n)&=k_1f(n-1)+k_2^{n-2}a_2\\&= k_1^2f(n-2)+k_1k_2^{n-3}a_2+k_2^{n-2}a_2\\&= k_1^{n-1}f(1)+a_2\big(k_1^{n-2}{k_2}^{0}+k_1^{n-3}k_2^{1}+\cdots+k_1^1k_2^{n-3}+k_2^{n-2}\big) \end{aligned}$

若 $k_1=k_2=k$，则
$f(n)=k^{n-1}f(1)+a_2(n-1)k^{n-2}=k^n(\frac{f(1)}{k}-\frac{a_2}{k^2}+\frac{a_2}{k^2}n)$

取
$c_1=\frac{f(1)}{k}-\frac{a_2}{k^2},\,\,c_2=\frac{a_2}{k^2}$

则
$f(n)=(c_1+nc_2)k^n$
若 $k_1\neq k_2$，则
$\begin{aligned} f(n)&=k_1^{n-1}f(1)+a_2\frac{k_1^{n-2}(1-(k_2/k_1)^{n-1})}{1-k_2/k_1}\\&= k_1^{n-1}f(1)+a_2\frac{k_1^{n-1}-k_2^{n-1}}{k_1-k_2}\\&= \frac{1}{k_1}\big(f(1)+\frac{a_2}{k_1-k_2}\big)k_1^{n}+\frac{1}{k_2}(-\frac{a_2}{k_1-k_2})k_2^{n} \end{aligned}$

取

$c_1=\frac{1}{k_1}\big(f(1)+\frac{a_2}{k_1-k_2}\big),\,\,c_2=\frac{1}{k_2}(-\frac{a_2}{k_1-k_2})$

则
$f(n)=c_1k_1^n+c_2k_2^n$

所以
$f(n)= \begin{cases} c_1k_1^n+c_2k_2^n, & k_1\neq k_2 \\ (c_1+nc_2)k^n, & k_1=k_2=k \end{cases}$

其中 $c_1,c_2$ 的值通过将初始值 $f(1),f(2)$

用此方法来求解斐波那契数列的通项公式，数列的递推公式如下：
$f(n)=f(n-1)+f(n-2)$

首先，求解 $k$，对应方程为
$k^2-k-1=0$

解得
$k_1=\frac{1+\sqrt{5}}{2},\,\,k_2=\frac{1-\sqrt{5}}{2}$

$k_1\neq k_2$，所以
$f(n)=c_1{k_1}^n+c_2{k_2}^n$

将 $f(1),f(2)$

$\begin{cases} \begin{aligned} f(1)&=c_1k_1+c_2k_2=1 \\ f(2)&=c_1{k_1}^2+c_2{k_2}^2=1 \end{aligned} \end{cases}$

解得
$c_1=\frac{1}{\sqrt{5}},\,\,c_2=-\frac{1}{\sqrt{5}}$

所以
$f(n)=\frac{1}{\sqrt{5}}\big(\frac{1+\sqrt{5}}{2}\big)^n-\frac{1}{\sqrt{5}}\big(\frac{1-\sqrt{5}}{2}\big)^n$

三、解法二

这个解法也是算法课上提到的，不过我把后续的特征分解部分添了上去，这还要得益于高等代数打下的基础！

$\begin{aligned} f(n)&=f(n-1)+f(n-2)\\ f(n-1)&=f(n-1) \end{aligned}$

表示成矩阵形式
$\begin{pmatrix}f(n)\\f(n-1)\end{pmatrix}=\begin{pmatrix}1 & 1\\1 & 0\end{pmatrix}\begin{pmatrix}f(n-1)\\f(n-2)\end{pmatrix}$

为表达简便，用下标表示项数，并记 $\displaystyle A=\begin{pmatrix}1 & 1\\1 & 0\end{pmatrix}$，补充定义 $f(0)=f_0=0$，则
$\begin{pmatrix}f_n\\f_{n-1}\end{pmatrix}= A\begin{pmatrix}f_{n-1}\\f_{n-2}\end{pmatrix}= A^2\begin{pmatrix}f_{n-2}\\f_{n-3}\end{pmatrix}= \cdots = A^{n-1}\begin{pmatrix}f_{1}\\f_{0}\end{pmatrix}= A^{n-1}\begin{pmatrix}1\\0\end{pmatrix}$

根据矩阵乘法的定义， $f_n$ 即为 $A^{n-1}$

问题转化为求 $A^{n-1}$，emm，看起来好像不是很简单.

嘶，对了，用 特征分解：

$n$ 阶方阵 $A$，都可以表示为以下形式：

$A=VDV^{-1}$

其中 $D=diag(\lambda_1,\lambda_2,\cdots,\lambda_n)$，$\lambda_i$ 为矩阵 $A$ 的特征值，$V$ 是由 $A$ 的特征向量拼起来的矩阵. 则
$\begin{aligned} A^m=(VDV^{-1})^m&=VDV^{-1}VDV^{-1}\cdots VDV^{-1}\\&= VD(V^{-1}V)D(V^{-1}V)\cdots DV^{-1}\\&=VD^mV^{-1} \end{aligned}$

由于 $D$ 为对角阵，所以 $D^m=diag(\lambda_1^m,\lambda_2^m,\cdots,\lambda_n^m)$.

首先求矩阵 $A$ 的特征值，根据定义，有
$\left|A-\lambda I\right|= \left|\begin{array}{cccc} 1-\lambda & 1 \\ 1 & -\lambda \end{array}\right|=0$

得到
$\lambda^2-\lambda-1=0$

细心的读者肯定会发现，这和前面的 $k^2-k-1=0$

$\lambda_1=\frac{1+\sqrt{5}}{2},\,\,\lambda_2=\frac{1-\sqrt{5}}{2}$

根据特征向量的定义
$(A-\lambda I)x=0$

求得一个解 $x=(\lambda,1)^T$. 所以 $\lambda_1$ 的一个特征向量为 $(\lambda_1,1)^T$，$\lambda_2$ 的一个特征向量为 $(\lambda_2,1)^T$，取
$V=\begin{pmatrix}\lambda_1 & \lambda_2 \\1 & 1\end{pmatrix}$

则
$V^{-1}=\frac{1}{\lambda_1-\lambda_2} \begin{pmatrix}1 & -\lambda_2\\ -1 & \lambda_1\end{pmatrix}$

于是
$\begin{aligned} \begin{pmatrix}f_n\\f_{n-1}\end{pmatrix}&= A^{n-1}\begin{pmatrix}1\\0\end{pmatrix}\\&= VD^{n-1}V^{-1}\begin{pmatrix}1\\0\end{pmatrix}\\&= \frac{1}{\lambda_1-\lambda_2} \begin{pmatrix}\lambda_1 & \lambda_2 \\1 & 1\end{pmatrix} \begin{pmatrix}\lambda_1^{n-1} & 0\\0 & \lambda_2^{n-1}\end{pmatrix} \begin{pmatrix}1 & -\lambda_2\\ -1 & \lambda_1\end{pmatrix} \begin{pmatrix}1\\0\end{pmatrix}\\&= \frac{1}{\lambda_1-\lambda_2} \begin{pmatrix}\lambda_1^{n} & \lambda_2^{n}\\\lambda_1^{n-1} & \lambda_2^{n-1}\end{pmatrix} \begin{pmatrix}1\\-1\end{pmatrix}\\&= \frac{1}{\lambda_1-\lambda_2} \begin{pmatrix}\lambda_1^{n}-\lambda_2^{n}\\\lambda_1^{n-1}-\lambda_2^{n-1}\end{pmatrix} \end{aligned}$

所以
$f(n)=\frac{1}{\lambda_1-\lambda_2}(\lambda_1^{n}-\lambda_2^{n})=\frac{1}{\sqrt{5}}\Big(\big(\frac{1+\sqrt{5}}{2}\big)^n-\big(\frac{1-\sqrt{5}}{2}\big)^n\Big)$

四、合二为一

下面来解释一下为什么特征根 $\lambda$ 满足的方程与解法一中 $k$

$x^2-x-1=0$

考虑一般形式的递推公式：

$f(n)=af(n-1)+bf(n-2)$

表示成矩阵形式
$\begin{pmatrix}f_n\\f_{n-1}\end{pmatrix}= \begin{pmatrix}a & b\\1 & 0\\\end{pmatrix} \begin{pmatrix}f_{n-1}\\f_{n-2}\end{pmatrix}= A\begin{pmatrix}f_{n-1}\\f_{n-2}\end{pmatrix}$

解法一中方程的由来我们已经讲得很清楚了，两边同时减去 $kf(n-1)$，则原式变为
$f(n)-kf(n-1)=(a-k)f(n-1)+bf(n-2)$

为使 $f(n)-kf(n-1)$ 与 $f(n-1)-kf(n-2)$ 为等比数列，需要求上式两边对应系数成比例，而这也可以用行列式表示，两行成比例等价于行列式为零，即
$\left| \begin{array}{cccc}a-k & b \\1 & -k \end{array} \right|=0$

这正是矩阵 $A$ 的特征值 $\lambda$ 需满足的方程 $|A-\lambda I|=0$.

所以，此文介绍的两种方法本质上是一样的.

五、实际实现

$n$

  上文已提到，
$\begin{pmatrix}f_n\\f_{n-1}\end{pmatrix}= A^{n-1}\begin{pmatrix}1\\0\end{pmatrix},A=\begin{pmatrix}1&1\\1&0\end{pmatrix}$

根据矩阵乘法的定义，$f_n$ 即为 $\small A^{n-1}$

$\small A^{n-1}$，下面介绍一种快速算法，能够将时间复杂度控制在 $\small O(\log n)$

先考虑 $x^n$ 的快速算法，将 $n$ 利用二进制进行表示，
$n=(d_kd_{k-1}\cdots d_1d_0)_2=d_k2^k+d_{k-1}2^{k-1}+\cdots+d_12+d_0$

$\small d_i=0\,or\, 1,0\leq i\leq k-1,\,d_k=1$，则 $x^n$ 可以用连乘形式进行表示.
$x^n=x^{d_k2^k+d_{k-1}2^{k-1}+\cdots+d_12+d_0}=x^{d_k2^k}x^{d_{k-1}2^{k-1}}\cdots x^{d_12}x^{d_0}$

考虑其中的某一项：

$d_i=0$，$x^{d_i2^i}=1$；
  若 $d_i=1$，$x^{d_i2^i}=x^{2^i}=(x^{2^{i-1}})(x^{2^{i-1}})$.

所以具体程序可以这样来写

number power(x,n)
{
	s = 1; // 存放结果
	t = x; // 记录x^{2^i}
	for(i = 0; i <= k; i++)
	{
		if(d_i != 0) // d_i表示二进制系数
		{
			s = s*t;
		}
		t = t*t;
	}
	return s;
}

上述算法是不完整的，因为没有给出 $d_i$

number power(x,n)
{
	s = 1; // 存放结果
	a = n; // 记录幂次
	t = x; // 记录x^{2^i}
	while(a > 0)
	{
		if(a mod 2) // d_i i=0,1,2,...,k
		{
			s = s*t;
		}
		t = t*t;
		a = a/2;    // 除2取整
	}
	return s;
}

时间复杂度为 $\small O(\log n)$. 类比得到 $\small A^{m}$

matrix power(A,m)
{
	S = I; // 存放结果，I为单位矩阵
	a = m; // 记录幂次
	T = A; // 记录A^{2^i}
	while(a > 0)
	{
		if(a mod 2) // d_i i=0,1,2,...,k
		{
			S = S*T;
		}
		T = T*T;
		a = a/2;    // 除2取整
	}
	return S;
}

所以求斐波那契数列通项的程序(MATLAB版)如下：

function f = fibs(n)
% 函数说明：求斐波那契数列第 n 项的快速算法
    S = eye(2);    % 单位矩阵 
    T = [1,1;1,0]; % 递推矩阵
    b = n-1;
    % 迭代计算
    while b > 0
        if mod(b,2)
            S = S*T;        
        end
        T = T*T;
        b = floor(b/2); 
    end    
    f = S(1,1);
end