递推公式斐波那契数列的几种求法

什么是递推公式

递推公式就是形如斐波那契数列那样，每一项都由前面几项运算求得

下面以斐波那契数列为例讲解递推公式的几种求解方法

斐波那契数列

斐波那契数列形式如下：

F(n)={1n=11n=2F(n−1)+F(n−2)n≥3F(n)=\left\{\begin{array}{rcl} 1~~~~~~~~~~~~~~~~~~~~~ n=1\\ 1~~~~~~~~~~~~~~~~~~~~~ n=2 \\ F(n-1)+F(n-2)~~~~n \geq 3 \end{array}\right.F(n)=⎩⎨⎧1 n=11 n=2F(n−1)+F(n−2) n≥3

斐波那契数列的前缀和公式为s(n)=f(n+2)−1s(n)=f(n+2)-1s(n)=f(n+2)−1

学习C语言后我们都知道斐波那契数列的递归求法，很好理解：

1.递归求法

时间复杂度O(2n)O(2^n)O(2n)

int Fibonacci(int n){return n > 2 ? Fibonacci(n-1) + Fibonacci(n-2) : 1;
}

根据上面的递归方法来看，由F(n)F(n)F(n)我们递归到F(n−1)F(n-1)F(n−1)和F(n−2)F(n-2)F(n−2)，接下来要求F(n−1)F(n-1)F(n−1)时，我们又递归到F(n−2)F(n-2)F(n−2)和F(n−3)F(n-3)F(n−3)；要求F(n−2)F(n-2)F(n−2)我们又需要递归到F(n−3)F(n-3)F(n−3)和F(n−4)F(n-4)F(n−4):

由上可知，F(n−2)F(n-2)F(n−2)、F(n−3)F(n-3)F(n−3)、F(n−4)F(n-4)F(n−4)…我们都求了两遍，这也是我们浪费时间的地方，稍微有点算法知识的都知道我们可以存起来每次求的值，办法有两种，一个是数组预处理，另一个是一次循环求

2.线性求法

数组预处理：

由于内存或者时间的限制此方法有一定的局限性

int F[maxn];
int num;    //num为我们需要求的最大的斐波那契数,num<N
void Fibonacci() {F[1] = F[2] = 1;for (int i = 3; i <= num; i++) {F[i] = F[i - 1] + F[i - 2];}
}

循环求:

通过中间变量保存F(n−1)F(n-1)F(n−1)和F(n−2)F(n-2)F(n−2)，由于每次ansansans都要更新，所以我们必须将ansansans初始值设为F(2)F(2)F(2)，然后令x=F(0)=0x=F(0)=0x=F(0)=0，y=F(1)y=F(1)y=F(1)，每次先更新xxx再更新ansansans再更新yyy，这样即可得出答案，见下面表格分析：

n	\	3	4	5	…
x	0	1	2	3	…
ans	1	2	3	5	…
y	1	1	2	3	…

代码：

int Fibonacci(int n) {if (n == 1 || n == 2) return 1;int x = 0, y = 1, ans = 1;for (int i = 3; i <= n; i++) {x = ans;ans = ans + y;y = x;}return ans;
}

3. 公式求法

斐波那契数列的通项公式为:

F(n)=15[(1+52)n−(1−52)n]F(n)= \frac{1}{\sqrt{5}}[(\frac{1+\sqrt{5}}{2})^n-(\frac{1-\sqrt{5}}{2})^n]F(n)=51[(21+5)n−(21−5)n]

对于斐波那契数列这样的二阶线性递推数列，可采用特征方程法：

特征方程

对于递推公式为xn=axn−1+bxn−2x_n=ax_{n-1}+bx_{n-2}xn=axn−1+bxn−2且a,b,x1,x2a,b,x_1,x_2a,b,x1,x2已知的递推数列，尝试待定系数法，验证{xn−kxn−1}\{x_n-kx_{n-1}\}{xn−kxn−1}是否为等比数列。设存在r,sr,sr,s使得xn−rxn−1=s(xn−1−rxn−2)x_n-rx_{n-1}=s(x_{n-1}-rx_{n-2})xn−rxn−1=s(xn−1−rxn−2)

所以xn=(s+r)xn−1−srxn−2x_n=(s+r)x_{n-1}-srx_{n-2}xn=(s+r)xn−1−srxn−2，那么可以得到a=s+r,b=−sra = s+r,b=-sra=s+r,b=−sr，根据韦达定理（x1+x2=−ba,x1x2=cax_1+x_2 = -\frac{b}{a}, x_1x_2 = \frac{c}{a}x1+x2=−ab,x1x2=ac，得到s,rs,rs,r满足方程：y2−ay−b=0y^2-ay-b=0y2−ay−b=0。因此{xn−rxn−1},{xn−sxn−1}\{x_n-rx_{n-1}\}, \{x_n-sx_{n-1}\}{xn−rxn−1},{xn−sxn−1}均是等比数列，分别列出等比数列的通项公式然后作差即可得到xnx_nxn。

如下是结论：

若方程有两相异根p、qp、qp、q，则：

xn=Apn+Bqn,A=x2−qx1p(p−q),B=px1−x2q(p−q)x_n = Ap^n+Bq^n,A=\frac{x_2-qx_1}{p(p-q)},B=\frac{px_1-x_2}{q(p-q)}xn=Apn+Bqn,A=p(p−q)x2−qx1,B=q(p−q)px1−x2
若方程有两等根ppp，则：

xn=(A+Bn)pn,A=2px1−x2p2,B=x2−px1p2x_n = (A+Bn)p^n,A=\frac{2px_1-x_2}{p^2},B=\frac{x_2-px_1}{p^2}xn=(A+Bn)pn,A=p22px1−x2,B=p2x2−px1

对于斐波那契数列的特征方程r2−r−1=0r^2-r-1=0r2−r−1=0，解得两不等根p=1+52,q=1−52p=\frac{1+\sqrt{5}}{2},q=\frac{1-\sqrt{5}}{2}p=21+5,q=21−5，可以直接代入上式求解A,BA,BA,B，也可以根据x1,x2x_1,x_2x1,x2已知直接利用结论，联立如下两个方程x1=Ap+Bq,x2=Ap2+Bq2x_1=Ap+Bq,x_2=Ap^2+Bq^2x1=Ap+Bq,x2=Ap2+Bq2求解A,BA,BA,B

最后我们得出其通项公式：F(n)=15[(1+52)n−(1−52)n]F(n)= \frac{1}{\sqrt{5}}[(\frac{1+\sqrt{5}}{2})^n-(\frac{1-\sqrt{5}}{2})^n]F(n)=51[(21+5)n−(21−5)n]

由于nnn增大时涉及大量浮点运算，会导致精度损失，因此无法通过浮点数计算斐波那契数列

int Fibonacci(int n){return (pow((1 + sqrt(5)) / 2, n) - pow((1 - sqrt(5)) / 2, n)) / sqrt(5);
}

但是通项公式还是可用的，只是要结合很多数论知识在模意义下使用

在ICPCICPCICPC中求大范围的斐波那契数一定都是需要取模的，而且取模的数一般都是大质数。假设为1e9+91e9+91e9+9，那么我们能求出x2≡5(mod1e9+9)x^2 \equiv 5(mod~~1e9+9)x2≡5(mod 1e9+9)是有解的，解出的两个解为616991993616991993616991993和383008016383008016383008016，也就是在这个取模的意义下5\sqrt{5}5和这两个解是相等的，那么我们就能直接通过公式计算F(n)F(n)F(n)了，设二次剩余方程的解存在且为xxx，那么有：

F(n)=inv(x)[((1+x)∗inv(2)%p)n−((1−x+p)%p∗inv(2))n]F(n)=inv(x)[((1+x)*inv(2)\%p)^n-((1-x+p)\%p*inv(2))^n]F(n)=inv(x)[((1+x)∗inv(2)%p)n−((1−x+p)%p∗inv(2))n]

4. 矩阵快速幂

对于一般的递推公式: ①F(n)=a∗F(n−1)+b∗F(n−2)F(n) = a*F(n-1) + b*F(n-2)F(n)=a∗F(n−1)+b∗F(n−2)

我们可以得到 ② F(n−1)=F(n−1)+0∗F(n−2)F(n-1) = F(n-1) + 0*F(n-2)F(n−1)=F(n−1)+0∗F(n−2)

根据线性代数知识，联立①和②，我们求得：

[F(n)F(n−1)]=[ab10][F(n−1)F(n−2)]\left[ \begin{array}{ccc} F(n) \\ F(n-1)\end{array} \right]=\left[ \begin{array}{ccc} a & b \\ 1 & 0 \end{array} \right]\left[ \begin{array}{ccc} F(n-1) \\ F(n-2) \end{array} \right][F(n)F(n−1)]=[a1b0][F(n−1)F(n−2)]

其中的常数矩阵称为递推式的特征矩阵

upd 2020.10.26：实际上上述的推导特征矩阵的过程不是通用的

假设Xn−1=[F(n−1)F(n−2)]X_{n-1}=\left[ \begin{array}{ccc} F(n-1) \\ F(n-2)\end{array} \right]Xn−1=[F(n−1)F(n−2)]，显然有Xn=[F(n)F(n−1)]X_n=\left[ \begin{array}{ccc} F(n) \\ F(n-1)\end{array} \right]Xn=[F(n)F(n−1)]，现在已知的是我们可以通过构造一个二阶方阵mmm使得Xn=m∗Xn−1X_n=m*X_{n-1}Xn=m∗Xn−1，首先假设方阵的所有项都是未知数，即

依次将等号右边的第二个式子代入，可以得到：m=[pqrs]m = \left[ \begin{array}{ccc} p & q \\ r & s \end{array} \right]m=[prqs]，那么有[F(n)F(n−1)]=[pqrs][F(n−1)F(n−2)]\left[ \begin{array}{ccc} F(n) \\ F(n-1)\end{array} \right]=\left[ \begin{array}{ccc} p & q \\ r & s \end{array} \right]\left[ \begin{array}{ccc} F(n-1) \\ F(n-2)\end{array} \right][F(n)F(n−1)]=[prqs][F(n−1)F(n−2)]，根据矩阵乘法的规律，我们不难一一求出p,q,r,sp,q,r,sp,q,r,s，最后即可得到特征方阵。通过递推规律，我们可以得到XnX_nXn和前两项的规律：

[F(n)F(n−1)]=[ab10]n−2[F(2)F(1)]\left[ \begin{array}{ccc} F(n) \\ F(n-1)\end{array} \right]=\left[ \begin{array}{ccc} a & b \\ 1 & 0 \end{array} \right]^{n-2}\left[ \begin{array}{ccc} F(2) \\ F(1) \end{array} \right][F(n)F(n−1)]=[a1b0]n−2[F(2)F(1)]

这样我们可以通过构造矩阵，然后使用矩阵快速幂来求解斐波那契数列的第nnn项了：

typedef long long ll;
const int Mod = 1e9 + 7;struct Matrix {ll matrix[10][10];int n, m;Matrix() {}Matrix(int x, int y) {n = x, m = y;memset(matrix, 0, sizeof(matrix));}
};Matrix mul(Matrix a, Matrix b) {Matrix ans(a.n, b.m);for (int i = 1; i <= ans.n; i++)for (int j = 1; j <= ans.m; j++)for (int k = 1; k <= a.m; k++) {   //前行乘后列ans.matrix[i][j] += (a.matrix[i][k] * b.matrix[k][j]) % Mod;ans.matrix[i][j] %= Mod;}return ans;
}Matrix qkp(Matrix mx, ll n) {Matrix ans(mx.n, mx.m);for (int i = 1; i <= mx.n; i++) ans.matrix[i][i] = 1;while (n) {if (n & 1) ans = mul(ans, mx);mx = mul(mx, mx);n >>= 1;}return ans;
}ll solve(int x) {if (x == 1 || x == 2) return 1;Matrix a(2, 2);  //特征矩阵及初始化a.matrix[1][1] = 1, a.matrix[1][2] = 1;a.matrix[2][1] = 1, a.matrix[2][2] = 0;Matrix b(2, 1);b.matrix[1][1] = 1, b.matrix[2][1] = 1;Matrix ans = mul(qkp(a, x - 2), b);return ans.matrix[1][1];
}

构造方阵解法：

对于[F(n)F(n−1)]=[ab10]n−2[F(2)F(1)]\left[ \begin{array}{ccc} F(n) \\ F(n-1)\end{array} \right]=\left[ \begin{array}{ccc} a & b \\ 1 & 0 \end{array} \right]^{n-2}\left[ \begin{array}{ccc} F(2) \\ F(1) \end{array} \right][F(n)F(n−1)]=[a1b0]n−2[F(2)F(1)]，我门可以继续递推，然后变成下面的式子：

[F(n)F(n−1)]=[ab10]n[F(1)F(−1)]\left[ \begin{array}{ccc} F(n) \\ F(n-1)\end{array} \right]=\left[ \begin{array}{ccc} a & b \\ 1 & 0 \end{array} \right]^{n}\left[ \begin{array}{ccc} F(1) \\ F(-1) \end{array} \right][F(n)F(n−1)]=[a1b0]n[F(1)F(−1)]

同理可得：[F(n+1)F(n)]=[ab10]n[F(1)F(0)]\left[ \begin{array}{ccc} F(n+1) \\ F(n)\end{array} \right]=\left[ \begin{array}{ccc} a & b \\ 1 & 0 \end{array} \right]^{n}\left[ \begin{array}{ccc} F(1) \\ F(0) \end{array} \right][F(n+1)F(n)]=[a1b0]n[F(1)F(0)]

联立上述两个矩阵，可以得到最终的方阵运算通式：

[F(n+1)F(n)F(n)F(n−1)]=[ab10]n[F(1)F(0)F(0)F(−1)]\left[ \begin{array}{ccc} F(n+1) & F(n) \\ F(n) & F(n-1)\end{array} \right]=\left[ \begin{array}{ccc} a & b \\ 1 & 0 \end{array} \right]^{n}\left[ \begin{array}{ccc} F(1) & F(0) \\ F(0) &F(-1) \end{array} \right][F(n+1)F(n)F(n)F(n−1)]=[a1b0]n[F(1)F(0)F(0)F(−1)]

由于有的递推式初始项没有给出F(0)F(0)F(0)和F(−1)F(-1)F(−1)，但我们由F(2)=F(1)+F(0)F(2)=F(1)+F(0)F(2)=F(1)+F(0)和F(1)=F(0)+F(−1)F(1)=F(0)+F(-1)F(1)=F(0)+F(−1)可以求出。下面的代码以斐波那契数列的递推矩阵为例：

typedef long long ll;
const int Mod = 1e9 + 7;struct Matrix {ll matrix[10][10];
};int n;  //矩阵的阶数Matrix mul(Matrix a, Matrix b) {Matrix res;memset(res.matrix, 0, sizeof res.matrix);for (int i = 1; i <= n; i++)for (int j = 1; j <= n; j++)for (int k = 1; k <= n; k++) {   //对于新矩阵的(i,j)位置，即前行乘后列res.matrix[i][j] += (a.matrix[i][k] * b.matrix[k][j]) % Mod;res.matrix[i][j] %= Mod;}return res;
}Matrix quick_pow(Matrix mx, ll x) {Matrix ans;   //快速幂我们初始化为1，那么这里就是矩阵的单位“1”，即单位矩阵memset(ans.matrix, 0, sizeof ans.matrix);for (int i = 1; i <= n; i++) ans.matrix[i][i] = 1;while (x) {if (x & 1) ans = mul(ans, mx);mx = mul(mx, mx);x >>= 1;}return ans;
}ll solve(int m = 2, int x) {n = m;Matrix mx;  //特征矩阵及初始化mx.matrix[1][1] = 1, mx.matrix[1][2] = 1;mx.matrix[2][1] = 1, mx.matrix[2][2] = 0;ll ans = quick_pow(mx, x);return mx.matrix[2][1];
}

5. 快速倍增法

斐波那契数列有如下性质：

F2n=Fn(2Fn+1−Fn)F_{2n}=F_n(2F_{n+1}-F_n)F2n=Fn(2Fn+1−Fn)
F2n+1=Fn+12+Fn2F_{2n+1}=F_{n+1}^2+F_n^2F2n+1=Fn+12+Fn2

那么可以通过递归增的方法求解二元组{F(n),F(n+1)}\{F(n),F(n+1)\}{F(n),F(n+1)}：

typedef long long ll;
typedef pair<ll, ll> pii;
const int Mod = 1e9 + 7;pii fib(ll n) {if (n == 0) return {0, 1};auto p = fib(n >> 1);ll c = (2 * p.second % Mod - p.first + Mod) % Mod * p.first % Mod;ll d = (p.first * p.first % Mod + p.second * p.second % Mod) % Mod;if (n & 1) return {d, c + d};else return {c, d};
}