斐波那契数列通项的两种求法

一、何为斐波那契数列？

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

具体定义如下

$$\{\begin{array}{l}f(n)=f(n-1)+f(n-2),n\ge 3\\ f(1)=f(2)=1\end{array}$$

二、解法一

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

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

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

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

将左端看作是等比数列 { a n } \{a_n\} {an​} 中的 $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_1{k}_2$$

将 $k=k_1$ 和上式代入 $(\ast )$ 式，得
$$f(n)-k_{1}f(n-1)=k_2(f(n-1)-k_{1}f(n-2))$$

则 $a_n=k_2{a}_{n-1}$，所以 $a_n=k_2^{n-2}{a}_2$.
$$\begin{array}{rl}f(n)& =k_{1}f(n-1)+k_2^{n-2}{a}_2\\ & =k_1^{2}f(n-2)+k_1{k}_2^{n-3}{a}_2+k_2^{n-2}{a}_2\\ & =k_1^{n-1}f(1)+a_2(k_1^{n-2}{k_2}^0+k_1^{n-3}{k}_2^1+\cdots +k_1^1{k}_2^{n-3}+k_2^{n-2})\end{array}$$

若 $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+n{c}_2)k^n$$
若 $k_1\ne k_2$，则
$$\begin{array}{rl}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}(f(1)+\frac{a_2}{k_1-k_2})k_1^n+\frac{1}{k_2}(-\frac{a_2}{k_1-k_2})k_2^n\end{array}$$

取

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

则
$$f(n)=c_1{k}_1^n+c_2{k}_2^n$$

所以
$$f(n)=\{\begin{array}{ll}{c}_1{k}_1^n+c_2{k}_2^n,& k_1\ne k_2\\ (c_1+n{c}_2)k^n,& k_1=k_2=k\end{array}$$

其中 $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\ne k_2$，所以
$$f(n)=c_1{k_1}^n+c_2{k_2}^n$$

将 $f(1),f(2)$ 的值代入，得

$$\{\begin{array}{l}\begin{array}{rl}f(1)& =c_1{k}_1+c_2{k}_2=1\\ f(2)& =c_1{k_1}^2+c_2{k_2}^2=1\end{array}\end{array}$$

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

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

三、解法二

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

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

表示成矩阵形式
$$\left(\begin{array}{c}f(n)\\ f(n-1)\end{array}\right)=\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)\left(\begin{array}{c}f(n-1)\\ f(n-2)\end{array}\right)$$

为表达简便，用下标表示项数，并记 $A=\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)$，补充定义 $f(0)=f_0=0$，则
$$\left(\begin{array}{c}{f}_n\\ f_{n-1}\end{array}\right)=A\left(\begin{array}{c}{f}_{n-1}\\ f_{n-2}\end{array}\right)=A^2\left(\begin{array}{c}{f}_{n-2}\\ f_{n-3}\end{array}\right)=\cdots =A^{n-1}\left(\begin{array}{c}{f}_1\\ f_0\end{array}\right)=A^{n-1}\left(\begin{array}{c}1\\ 0\end{array}\right)$$

根据矩阵乘法的定义， $f_n$ 即为 $A^{n-1}$ 首行首列的元素.

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

嘶，对了，用 特征分解：

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

$$A=VD{V}^{-1}$$

其中 $D=diag({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n)$，${\lambda }_i$ 为矩阵 $A$ 的特征值，$V$ 是由 $A$ 的特征向量拼起来的矩阵. 则
$$\begin{array}{rl}{A}^m=(VD{V}^{-1}{)}^m& =VD{V}^{-1}VD{V}^{-1}\cdots VD{V}^{-1}\\ & =VD(V^{-1}V)D(V^{-1}V)\cdots D{V}^{-1}\\ & =V{D}^m{V}^{-1}\end{array}$$

由于 $D$ 为对角阵，所以 $D^m=diag({\lambda }_1^m,{\lambda }_2^m,\cdots ,{\lambda }_n^m)$.

首先求矩阵 $A$ 的特征值，根据定义，有
$$\left|A-\lambda I\right|=\left|\begin{array}{cc}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=\left(\begin{array}{cc}{\lambda }_1& {\lambda }_2\\ 1& 1\end{array}\right)$$

则
$$V^{-1}=\frac{1}{{\lambda }_1-{\lambda }_2}\left(\begin{array}{cc}1& -{\lambda }_2\\ -1& {\lambda }_1\end{array}\right)$$

于是
$$\begin{array}{rl}\left(\begin{array}{c}{f}_n\\ f_{n-1}\end{array}\right)& =A^{n-1}\left(\begin{array}{c}1\\ 0\end{array}\right)\\ & =V{D}^{n-1}{V}^{-1}\left(\begin{array}{c}1\\ 0\end{array}\right)\\ & =\frac{1}{{\lambda }_1-{\lambda }_2}\left(\begin{array}{cc}{\lambda }_1& {\lambda }_2\\ 1& 1\end{array}\right)\left(\begin{array}{cc}{\lambda }_1^{n-1}& 0\\ 0& {\lambda }_2^{n-1}\end{array}\right)\left(\begin{array}{cc}1& -{\lambda }_2\\ -1& {\lambda }_1\end{array}\right)\left(\begin{array}{c}1\\ 0\end{array}\right)\\ & =\frac{1}{{\lambda }_1-{\lambda }_2}\left(\begin{array}{cc}{\lambda }_1^n& {\lambda }_2^n\\ {\lambda }_1^{n-1}& {\lambda }_2^{n-1}\end{array}\right)\left(\begin{array}{c}1\\ -1\end{array}\right)\\ & =\frac{1}{{\lambda }_1-{\lambda }_2}\left(\begin{array}{c}{\lambda }_1^n-{\lambda }_2^n\\ {\lambda }_1^{n-1}-{\lambda }_2^{n-1}\end{array}\right)\end{array}$$

所以
$$f(n)=\frac{1}{{\lambda }_1-{\lambda }_2}({\lambda }_1^n-{\lambda }_2^n)=\frac{1}{\sqrt{5}}\left(\right(\frac{1+\sqrt{5}}{2}{)}^n-\left(\frac{1-\sqrt{5}}{2}{)}^n\right)$$

四、合二为一

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

$$x^2-x-1=0$$

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

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

表示成矩阵形式
$$\left(\begin{array}{c}{f}_n\\ f_{n-1}\end{array}\right)=\left(\begin{array}{cc}a& b\\ 1& 0\end{array}\right)\left(\begin{array}{c}{f}_{n-1}\\ f_{n-2}\end{array}\right)=A\left(\begin{array}{c}{f}_{n-1}\\ f_{n-2}\end{array}\right)$$

解法一中方程的由来我们已经讲得很清楚了，两边同时减去 $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}{cc}a-k& b\\ 1& -k\end{array}\right|=0$$

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

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

五、实际实现

  虽然已经求出完美的通项公式，理论上将可以求得数列的任意项. 但实际上，由于通项公式中含有无理数，计算机进行计算时，是用的有限位浮点数来表示无理数，所以当指数 $n$ 较大时，会产生较大的运算误差，这是我们不想要的.

  上文已提到，
$$\left(\begin{array}{c}{f}_n\\ f_{n-1}\end{array}\right)=A^{n-1}\left(\begin{array}{c}1\\ 0\end{array}\right),A=\left(\begin{array}{cc}1& 1\\ 1& 0\end{array}\right)$$

根据矩阵乘法的定义，$f_n$ 即为 $A^{n-1}$ 首行首列的元素.

  由于是整数的乘法和加法运算，所以这种方法不存在运算误差，但是需要计算 $A^{n-1}$，下面介绍一种快速算法，能够将时间复杂度控制在 $O(\log n)$ 量级.

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

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

考虑其中的某一项：

  若 $d_i=0$，$x^{d_i{2}^i}=1$；
  若 $d_i=1$，$x^{d_i{2}^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;
}

时间复杂度为 $O(\log n)$. 类比得到 $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