跳到主要内容

斐波那契数列

· 2 分钟阅读
Josh Cena

这是 C 社算法团队的第一次活动。我们介绍了斐波那契数列的第nn项求解方法,主要运用了矩阵快速幂算法。

题目#

斐波那契数列:
Fn={0,n=01,n=1Fn2+Fn1,n>1F_{n}= \begin{cases} 0,&n=0\\ 1,&n=1\\ F_{n-2}+F_{n-1},&n>1 \end{cases}

给定nn,求Fn mod 109+7F_{n}\text{ mod }10^9+7

数据规模内存限制运行时间
0n10190\le n\le 10^{19}64MB1.0s

题解#

101910^{19}显然灭掉了所有用循环解决的想法。有没有比简单的O(n)\mathcal{O}(n)更好一点的方法?用矩阵快速幂,可以达到O(logn)\mathcal{O}(\log n)。观察到:

(Fn+1Fn+2)=(Fn+1Fn+Fn+1)=(0111)(FnFn+1)\begin{pmatrix}F_{n+1}\\F_{n+2}\end{pmatrix}=\begin{pmatrix}F_{n+1}\\F_{n}+F_{n+1}\end{pmatrix}=\begin{pmatrix}0&1\\1&1\end{pmatrix}\begin{pmatrix}F_{n}\\F_{n+1}\end{pmatrix}

这一步对于所有递推数列都是适用的,因此在有经验之后应该非常容易得到。一般地,对于Fn+2=aFn+bFn+1F_{n+2}=aF_{n}+bF_{n+1},有

(Fn+1Fn+2)=(Fn+1aFn+bFn+1)=(01ab)(FnFn+1)\begin{pmatrix}F_{n+1}\\F_{n+2}\end{pmatrix}=\begin{pmatrix}F_{n+1}\\aF_{n}+bF_{n+1}\end{pmatrix}=\begin{pmatrix}0&1\\a&b\end{pmatrix}\begin{pmatrix}F_{n}\\F_{n+1}\end{pmatrix}

从递推式中有

(Fn+mFn+m+1)=(0111)m(FnFn+1)\begin{pmatrix}F_{n+m}\\F_{n+m+1}\end{pmatrix}=\begin{pmatrix}0&1\\1&1\end{pmatrix}^m\begin{pmatrix}F_{n}\\F_{n+1}\end{pmatrix}

n=0n=0,得到

(FmFm+1)=(0111)m(F0F1)\begin{pmatrix}F_{m}\\F_{m+1}\end{pmatrix}=\begin{pmatrix}0&1\\1&1\end{pmatrix}^m\begin{pmatrix}F_0\\F_1\end{pmatrix}

因此把问题转化成了如何求矩阵mm次方的问题。如果设m=20a0+21a1+22a2+m=2^0a_0+2^1a_1+2^2a_2+\dots(也就是把mm用二进制表示),那么有

(0111)m=((0111)1)a0×((0111)2)a1×((0111)4)a2\begin{pmatrix}0&1\\1&1\end{pmatrix}^m=\left(\begin{pmatrix}0&1\\1&1\end{pmatrix}^{1}\right)^{a_0}\times \left(\begin{pmatrix}0&1\\1&1\end{pmatrix}^{2}\right)^{a_1}\times \left(\begin{pmatrix}0&1\\1&1\end{pmatrix}^{4}\right)^{a_2}\dots

而这些矩阵的2k2^k次方,完全可以预处理。当mm的数量级为101910^{19}时,k<log21019<64k<\log_2 10^{19}<64,最多只需要存储 63 个矩阵。并且

(0111)2k=(0111)2k1×(0111)2k1\begin{pmatrix}0&1\\1&1\end{pmatrix}^{2^k}=\begin{pmatrix}0&1\\1&1\end{pmatrix}^{2^{k-1}}\times \begin{pmatrix}0&1\\1&1\end{pmatrix}^{2^{k-1}}

这些乘方可以在O(logm)\mathcal{O}(\log m)时间内得到。这便是快速幂的思想:计算所有的2k2^k次方,然后把其中需要的那些组合起来即可。

程序#

下面是 C++ 代码,其中最繁琐的部分是实现矩阵乘法:

#include <iostream>#include <cmath>
using namespace std;
struct mat {    unsigned long long a[4];    mat operator *(mat o) {        mat t;        t.a[0] = (this->a[0] * o.a[0] + this->a[1] * o.a[2]) % 1000000007;        t.a[1] = (this->a[0] * o.a[1] + this->a[1] * o.a[3]) % 1000000007;        t.a[2] = (this->a[2] * o.a[0] + this->a[3] * o.a[2]) % 1000000007;        t.a[3] = (this->a[2] * o.a[1] + this->a[3] * o.a[3]) % 1000000007;        return t;    }};
// 预处理的矩阵 2^k 次幂mat mat_pow[64];
int fib(unsigned long long k) {    // 临时矩阵,每次在此上面乘以 mat_pow 中的某项    mat tmp;    tmp.a[0] = 1;    tmp.a[1] = 0;    tmp.a[2] = 0;    tmp.a[3] = 1;    for (int i = 0; i < 64; i++) {        // 如果 a_i 为 1        if (k & (1ull << i)) {            tmp = tmp * mat_pow[i];        }    }    return tmp.a[1];}
int main() {    mat_pow[0].a[0] = 0;    mat_pow[0].a[1] = 1;    mat_pow[0].a[2] = 1;    mat_pow[0].a[3] = 1;    for (int i = 1; i < 64; i++) {        mat_pow[i] = mat_pow[i-1] * mat_pow[i-1];    }    unsigned long long n;    cin >> n;    cout << fib(n) << endl;    return 0;}

补充一下矩阵的乘法公式:

(a0a1a2a3)×(b0b1b2b3)=(a0b0+a1b2a0b1+a1b3a2b0+a3b2a2b1+a3b3)\begin{pmatrix}a_0&a_1\\a_2&a_3\end{pmatrix}\times\begin{pmatrix}b_0&b_1\\b_2&b_3\end{pmatrix}=\begin{pmatrix}a_0b_0+a_1b_2&a_0b_1+a_1b_3\\a_2b_0+a_3b_2&a_2b_1+a_3b_3\end{pmatrix}

(a0a1a2a3)×(b0b1)=(a0b0+a1b1a2b0+a3b1)\begin{pmatrix}a_0&a_1\\a_2&a_3\end{pmatrix}\times\begin{pmatrix}b_0\\b_1\end{pmatrix}=\begin{pmatrix}a_0b_0+a_1b_1\\a_2b_0+a_3b_1\end{pmatrix}