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

一、何为斐波那契数列？

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

具体定义如下

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

二、解法一

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

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

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

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

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

整理得
k2−ak−b=0k^2-ak-b=0 k2−ak−b=0

假设该方程的根已求出，为 k1,k2k_1,k_2k1,k2，由韦达定理
a=k1+k2,b=−k1k2a=k_1+k_2,\,\,b=-k_1k_2 a=k1+k2,b=−k1k2

将 k=k1k=k_1k=k1 和上式代入 (∗)(*)(∗) 式，得
f(n)−k1f(n−1)=k2(f(n−1)−k1f(n−2))f(n)-k_1f(n-1)=k_2\big(f(n-1)-k_1f(n-2)\big) f(n)−k1f(n−1)=k2(f(n−1)−k1f(n−2))

则 an=k2an−1a_n=k_2 a_{n-1}an=k2an−1，所以 an=k2n−2a2a_n=k_2^{n-2}a_2an=k2n−2a2.
f(n)=k1f(n−1)+k2n−2a2=k12f(n−2)+k1k2n−3a2+k2n−2a2=k1n−1f(1)+a2(k1n−2k20+k1n−3k21+⋯+k11k2n−3+k2n−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} f(n)=k1f(n−1)+k2n−2a2=k12f(n−2)+k1k2n−3a2+k2n−2a2=k1n−1f(1)+a2(k1n−2k20+k1n−3k21+⋯+k11k2n−3+k2n−2)

若 k1=k2=kk_1=k_2=kk1=k2=k，则
f(n)=kn−1f(1)+a2(n−1)kn−2=kn(f(1)k−a2k2+a2k2n)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) f(n)=kn−1f(1)+a2(n−1)kn−2=kn(kf(1)−k2a2+k2a2n)

取
c1=f(1)k−a2k2,c2=a2k2c_1=\frac{f(1)}{k}-\frac{a_2}{k^2},\,\,c_2=\frac{a_2}{k^2} c1=kf(1)−k2a2,c2=k2a2

则
f(n)=(c1+nc2)knf(n)=(c_1+nc_2)k^n f(n)=(c1+nc2)kn
若 k1≠k2k_1\neq k_2k1=k2，则
f(n)=k1n−1f(1)+a2k1n−2(1−(k2/k1)n−1)1−k2/k1=k1n−1f(1)+a2k1n−1−k2n−1k1−k2=1k1(f(1)+a2k1−k2)k1n+1k2(−a2k1−k2)k2n\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} f(n)=k1n−1f(1)+a21−k2/k1k1n−2(1−(k2/k1)n−1)=k1n−1f(1)+a2k1−k2k1n−1−k2n−1=k11(f(1)+k1−k2a2)k1n+k21(−k1−k2a2)k2n

取

c1=1k1(f(1)+a2k1−k2),c2=1k2(−a2k1−k2)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}) c1=k11(f(1)+k1−k2a2),c2=k21(−k1−k2a2)

则
f(n)=c1k1n+c2k2nf(n)=c_1k_1^n+c_2k_2^n f(n)=c1k1n+c2k2n

所以
f(n)={c1k1n+c2k2n,k1≠k2(c1+nc2)kn,k1=k2=kf(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} f(n)={c1k1n+c2k2n,(c1+nc2)kn,k1=k2k1=k2=k

其中 c1,c2c_1,c_2c1,c2 的值通过将初始值 f(1),f(2)f(1),f(2)f(1),f(2) 代入求得.

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

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

解得
k1=1+52,k2=1−52k_1=\frac{1+\sqrt{5}}{2},\,\,k_2=\frac{1-\sqrt{5}}{2} k1=21+5,k2=21−5

k1≠k2k_1\neq k_2k1=k2，所以
f(n)=c1k1n+c2k2nf(n)=c_1{k_1}^n+c_2{k_2}^n f(n)=c1k1n+c2k2n

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

{f(1)=c1k1+c2k2=1f(2)=c1k12+c2k22=1\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} {f(1)f(2)=c1k1+c2k2=1=c1k12+c2k22=1

解得
c1=15,c2=−15c_1=\frac{1}{\sqrt{5}},\,\,c_2=-\frac{1}{\sqrt{5}} c1=51,c2=−51

所以
f(n)=15(1+52)n−15(1−52)nf(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 f(n)=51(21+5)n−51(21−5)n

三、解法二

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

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

表示成矩阵形式
(f(n)f(n−1))=(1110)(f(n−1)f(n−2))\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} (f(n)f(n−1))=(1110)(f(n−1)f(n−2))

为表达简便，用下标表示项数，并记 A=(1110)\displaystyle A=\begin{pmatrix}1 & 1\\1 & 0\end{pmatrix}A=(1110)，补充定义 f(0)=f0=0f(0)=f_0=0f(0)=f0=0，则
(fnfn−1)=A(fn−1fn−2)=A2(fn−2fn−3)=⋯=An−1(f1f0)=An−1(10)\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} (fnfn−1)=A(fn−1fn−2)=A2(fn−2fn−3)=⋯=An−1(f1f0)=An−1(10)

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

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

嘶，对了，用 特征分解：

任何一个 nnn 阶方阵 AAA，都可以表示为以下形式：

A=VDV−1A=VDV^{-1} A=VDV−1

其中 D=diag(λ1,λ2,⋯,λn)D=diag(\lambda_1,\lambda_2,\cdots,\lambda_n)D=diag(λ1,λ2,⋯,λn)，λi\lambda_iλi 为矩阵 AAA 的特征值，VVV 是由 AAA 的特征向量拼起来的矩阵. 则
Am=(VDV−1)m=VDV−1VDV−1⋯VDV−1=VD(V−1V)D(V−1V)⋯DV−1=VDmV−1\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} Am=(VDV−1)m=VDV−1VDV−1⋯VDV−1=VD(V−1V)D(V−1V)⋯DV−1=VDmV−1

由于 DDD 为对角阵，所以 Dm=diag(λ1m,λ2m,⋯,λnm)D^m=diag(\lambda_1^m,\lambda_2^m,\cdots,\lambda_n^m)Dm=diag(λ1m,λ2m,⋯,λnm).

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

得到
λ2−λ−1=0\lambda^2-\lambda-1=0 λ2−λ−1=0

细心的读者肯定会发现，这和前面的 k2−k−1=0k^2-k-1=0k2−k−1=0 实际上是同一个方程(其实这并不是巧合，具体原因会在文末分析)，所以解也是一样的.

λ1=1+52,λ2=1−52\lambda_1=\frac{1+\sqrt{5}}{2},\,\,\lambda_2=\frac{1-\sqrt{5}}{2} λ1=21+5,λ2=21−5

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

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

则
V−1=1λ1−λ2(1−λ2−1λ1)V^{-1}=\frac{1}{\lambda_1-\lambda_2} \begin{pmatrix}1 & -\lambda_2\\ -1 & \lambda_1\end{pmatrix} V−1=λ1−λ21(1−1−λ2λ1)

于是
(fnfn−1)=An−1(10)=VDn−1V−1(10)=1λ1−λ2(λ1λ211)(λ1n−100λ2n−1)(1−λ2−1λ1)(10)=1λ1−λ2(λ1nλ2nλ1n−1λ2n−1)(1−1)=1λ1−λ2(λ1n−λ2nλ1n−1−λ2n−1)\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} (fnfn−1)=An−1(10)=VDn−1V−1(10)=λ1−λ21(λ11λ21)(λ1n−100λ2n−1)(1−1−λ2λ1)(10)=λ1−λ21(λ1nλ1n−1λ2nλ2n−1)(1−1)=λ1−λ21(λ1n−λ2nλ1n−1−λ2n−1)

所以
f(n)=1λ1−λ2(λ1n−λ2n)=15((1+52)n−(1−52)n)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) f(n)=λ1−λ21(λ1n−λ2n)=51((21+5)n−(21−5)n)

四、合二为一

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

x2−x−1=0x^2-x-1=0 x2−x−1=0

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

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

表示成矩阵形式
(fnfn−1)=(ab10)(fn−1fn−2)=A(fn−1fn−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} (fnfn−1)=(a1b0)(fn−1fn−2)=A(fn−1fn−2)

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

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

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

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

五、实际实现

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

上文已提到，
(fnfn−1)=An−1(10),A=(1110)\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} (fnfn−1)=An−1(10),A=(1110)

根据矩阵乘法的定义，fnf_nfn 即为 An−1\small A^{n-1}An−1 首行首列的元素.

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

先考虑 xnx^nxn 的快速算法，将 nnn 利用二进制进行表示，
n=(dkdk−1⋯d1d0)2=dk2k+dk−12k−1+⋯+d12+d0n=(d_kd_{k-1}\cdots d_1d_0)_2=d_k2^k+d_{k-1}2^{k-1}+\cdots+d_12+d_0 n=(dkdk−1⋯d1d0)2=dk2k+dk−12k−1+⋯+d12+d0

di=0or1,0≤i≤k−1,dk=1\small d_i=0\,or\, 1,0\leq i\leq k-1,\,d_k=1di=0or1,0≤i≤k−1,dk=1，则 xnx^nxn 可以用连乘形式进行表示.
xn=xdk2k+dk−12k−1+⋯+d12+d0=xdk2kxdk−12k−1⋯xd12xd0x^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} xn=xdk2k+dk−12k−1+⋯+d12+d0=xdk2kxdk−12k−1⋯xd12xd0

考虑其中的某一项：

若 di=0d_i=0di=0，xdi2i=1x^{d_i2^i}=1xdi2i=1；
若 di=1d_i=1di=1，xdi2i=x2i=(x2i−1)(x2i−1)x^{d_i2^i}=x^{2^i}=(x^{2^{i-1}})(x^{2^{i-1}})xdi2i=x2i=(x2i−1)(x2i−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;
}

上述算法是不完整的，因为没有给出 did_idi 的求法，给出求法后，得到

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)\small O(\log n)O(logn). 类比得到 Am\small A^{m}Am 的快速算法.

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 > 0if mod(b,2)S = S*T;        endT = T*T;b = floor(b/2); end    f = S(1,1);
end