[模板][题解][Luogu1939]矩阵乘法加速递推（详解）

题目传送门

题目大意：计算数列a的第n项，其中：

\[a[1] = a[2] = a[3] = 1\]

\[a[i] = a[i-3] + a[i - 1]\]

\[(n ≤ 2 \times 10^9)\]

一般的递推是O(n)的，显然时间和空间都不能承受。

由于每一步递推都是相同的。这句话包含了2个层面：首先，递推式是相同的；其次，递推的条件也要是相同的。综合来说，每一步的递推都是相同的。这是应用矩阵加速递推的充分条件。

那么怎么进行矩阵加速呢？我们首先观察，第\(i\)项和哪些项有关？ 与\(i-3\)和\(i-1\)优化，也就是和前3项有关。为了能够**仅仅利用矩阵就能推出\(a[i]\)， 我们需要记录前3项，于是，我们构造一个3*3的矩阵：
\[ A=\begin{bmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{bmatrix} \]
有同学会问：这个矩阵是怎么构造出来的呢？

我们首先要构造出类似于DP的状态来存下所有计算过程中可能会用到的信息，对于这道题，我们需要记录：(假设我们要从\(a[i]\)推到\(a[i+1]\))
\[ B=\begin{bmatrix} a[i] \\ a[i-1] \\ a[i-2] \\ \end{bmatrix} \]
这个矩阵要推到：
\[ C=\begin{bmatrix} a[i+1] \\ a[i] \\ a[i-1] \\ \end{bmatrix} \]
也就是说，我们需要构造一个矩阵\(A\)，使得\(A*B = C\)，根据线性代数的相关定义，A一定是一个\(3*3\)的矩阵，没错吧。

那好，我们已经得到：
\[ \begin{bmatrix} ? & ? & ? \\ ? & ? & ? \\ ? & ? & ? \\ \end{bmatrix}*\begin{bmatrix} a[i] \\ a[i-1] \\ a[i-2] \\ \end{bmatrix} =\begin{bmatrix} a[i+1] \\ a[i] \\ a[i-1] \\ \end{bmatrix} \\A*B=C \]
我们只需要根据递推式，把矩阵\(A\)中的数填满就可以了，比如说：

由于$a[i+1] = a[i-2] +a[i] \(，而根据矩阵，\)a[i+1] = (0,0) * a[i] + (0,1) * a[i-1] + (0,2) * a[i-2]$，所以，矩阵的第一行是\(1,0,1\)，以此类推，就可以吧矩阵填满了。

然后，我们可以得到：
\[ \begin{bmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{bmatrix}*\begin{bmatrix} a[i] \\ a[i-1] \\ a[i-2] \\ \end{bmatrix} =\begin{bmatrix} a[i+1] \\ a[i] \\ a[i-1] \\ \end{bmatrix} \\A*B=C \]
可是，有了这个，怎么从\(a[1]\)推到\(a[n]\)呢？

我们知道：\(a[1] = a[2] = a[3] = 1\),如果把它们代入矩阵\(B\)（就是中间的那个矩阵），我们会得到：
\[ \begin{bmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{bmatrix}*\begin{bmatrix} a[3]=1 \\ a[2]=1 \\ a[1]=1 \\ \end{bmatrix} =\begin{bmatrix} a[4]=1*a[3] + 1 * a[1] = 2 \\ a[3] = 1 \\ a[2] = 1 \\ \end{bmatrix} \\A * B = C \]
一开始我们只知道\(a[1], a[2], a[3]\)，但是两个矩阵相乘后，我们得到了一个新的值\(a[4] = 2\)，很开心有木有。如果我们用矩阵\(A\)去乘矩阵\(C\)，会得到一个新的矩阵,暂且叫\(C'\)，你会得到有一个新的值\(a[5]\)，我们有点兴奋起来，这有点想多米诺骨牌，推到第一个，会一直向前倒，知道最后一个。我相信你脑子一定有了这样一个式子：
\[ …… \begin{bmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{bmatrix}* \begin{bmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{bmatrix}* \begin{bmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{bmatrix}*\begin{bmatrix} a[3]=1 \\ a[2]=1 \\ a[1]=1 \\ \end{bmatrix} =\begin{bmatrix} a[n] \\ a[n-1] \\ a[n-2] \\ \end{bmatrix} \]
矩阵乘法有结合律（但没有交换律，不要问我为什么），所以左边一堆相同的矩阵（不妨叫系数矩阵）可以用一个括号括起来，就像这样：
\[ \left(…… \begin{bmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{bmatrix}* \begin{bmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{bmatrix}* \begin{bmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{bmatrix}\right)*\begin{bmatrix} a[3]=1 \\ a[2]=1 \\ a[1]=1 \\ \end{bmatrix} =\begin{bmatrix} a[n] \\ a[n-1] \\ a[n-2] \\ \end{bmatrix} \]
想到了什么？
\[ \begin{bmatrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{bmatrix}^k*\begin{bmatrix} a[3]=1 \\ a[2]=1 \\ a[1]=1 \\ \end{bmatrix} =\begin{bmatrix} a[n] \\ a[n-1] \\ a[n-2] \\ \end{bmatrix} \]
我们可以得到\(k = n - 3\)（想想为什么？），由于n很大，能不能快速求矩阵k次方呢？

在mod p意义下？矩阵乘法满足结合律？

快速幂！

为什么可以用快速幂这个黑科技呢？

普通的快速幂是用来求\(b^k mod\ p\)的，其原理是把\(k\)二进制拆分成\(k=2^{a_1}+2^{a_2}+ ……+2^{a_k}\),从而得到\(b^k mod \ p = b^{2^{a_1}} * b^{2^{a_2}} * ……*b^{2^{a_k}} mod \ p = ((b^{2^{a_1}} mod \ p) *(b^{2^{a_2}} mod \ p) * …… * (b^{2^{a_k}} mod \ p))\ mod\ p\) ，只要满足乘法结合律，就可以进行快速幂。

矩阵快速幂通常用来加速递推。比如说斐波那契数列的第n项mod p的值也可以用矩阵快速幂来求。但是并不是所有的递推都可以用矩阵快速幂，只有那些转移具有周期性的递推才能使用。

代码模块

1、矩阵的定义（结构体）

struct Square{
    int mat[3][3];
    void clear(){
        memset(mat, 0, sizeof(mat));
    }
} Base, Result;

2、方阵的乘法

void Times(Square &A, Square B){
    Square C; C.clear();
    for (int i = 0; i <= 2; ++i)
        for (int j = 0; j <= 2; ++j)
            for (int k = 0; k <= 2; ++k)
                (C.mat[i][j] += (LL)A.mat[i][k] * B.mat[k][j] % p) %= p;
    A = C;
}

3、矩阵快速幂

void SquareQpow(Square Base, int k){
    if (k == 1){
        Result = Base;
        return;
    }
    Result.clear();
    SquareQpow(Base, k / 2);
    Times(Result, Result);
    if (k % 2 == 1) Times(Result, Base);    
}

4、矩阵初始化

void init(){
    Base.clear();
    Base.mat[0][0] = 1; Base.mat[0][2] = 1;
    Base.mat[1][0] = 1; Base.mat[2][1] = 1;
}

易错点

1. 乘法时需进行强制类型转换：(C.mat[i][j] += (LL)A.mat[i][k] * B.mat[k][j] % p) %= p;
2. C++数组从0开始的哦QAQ
3. 计算答案时注意要加3项，不要只加2项

完整代码

#include<iostream>
#include<cstdio>
#include<cstring>
#define LL long long
using namespace std;
const int p = 1e9 + 7;
struct Square{
    int mat[3][3];
    void clear(){
        memset(mat, 0, sizeof(mat));
    }
} Base, Result;

void init(){
    Base.clear();
    Base.mat[0][0] = 1; Base.mat[0][2] = 1;
    Base.mat[1][0] = 1; Base.mat[2][1] = 1;
}

void Times(Square &A, Square B){
    Square C; C.clear();
    for (int i = 0; i <= 2; ++i)
        for (int j = 0; j <= 2; ++j)
            for (int k = 0; k <= 2; ++k)
                (C.mat[i][j] += (LL)A.mat[i][k] * B.mat[k][j] % p) %= p;
    A = C;
}

void SquareQpow(Square Base, int k){
    if (k == 1){
        Result = Base;
        return;
    }
    Result.clear();
    SquareQpow(Base, k / 2);
    Times(Result, Result);
    if (k % 2 == 1) Times(Result, Base);    
}

int main(){
    int T; scanf("%d", &T);
    while (T--){
        int n;
        scanf("%d", &n);
        if (n <= 3) printf("1\n");
        else{
            init();
            SquareQpow(Base, n-3);
            printf("%d\n", ((Result.mat[0][0] + Result.mat[0][1]) % p + Result.mat[0][2]) % p);
        }
    }
    return 0;
}

转载于:https://www.cnblogs.com/YJZoier/p/9442086.html