【洛谷P1939】矩阵加速（数列）

Description

定义序列a的通项公式为$a_i=\left\{\begin{aligned}1 && i \leq 3 \\a_{i-1}+a_{i-3} && i \geq 4\end{aligned}\right.$

求序列a的第n项对1e9+7取模的值

Solution

由于n的值很大，所以普通的递推无法解决，我们考虑使用矩阵乘法

假设$F(n)=\begin{bmatrix}a_n & a_{n-1} & a_{n-2}\end{bmatrix}$，我们考虑将$F(n-1)$乘上一个矩阵base之后变成$F(n)$

也就是说$\begin{bmatrix}a_{n-1} & a_{n-2} & a_{n-3}\end{bmatrix} \times base = \begin{bmatrix}a_n & a_{n-1} & a_{n-2}\end{bmatrix}$

显然，我们可以轻易计算出$base=\begin{bmatrix}1 & 1 & 0 \\0 & 0 & 1 \\ 1 & 0 & 0 \end{bmatrix}$

我们定义初始矩阵$a=\begin{bmatrix}1 & 1 & 1 \end{bmatrix}$，然后利用矩阵快速幂计算$a \times base^{n-3}$即可

时间复杂度为$O(3^3log_2n)$

Code

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 1e9 + 7;
int T, n;
struct Matrix {
    ll m[4][4];
    Matrix() {memset(m, 0, sizeof(m));}
    Matrix operator *(const Matrix &x) const {
        Matrix ret;
        for (register int i = 1; i <= 3; ++i)
            for (register int j = 1; j <= 3; ++j)
                for (register int k = 1; k <= 3; ++k)
                    ret.m[i][j] = (ret.m[i][j] + m[i][k] * x.m[k][j] % mod) % mod; 
        return ret;
    }
} a, base;
inline void init() {
    base.m[1][1] = 1; base.m[1][2] = 1; base.m[1][3] = 0;
    base.m[2][1] = 0; base.m[2][2] = 0; base.m[2][3] = 1;
    base.m[3][1] = 1; base.m[3][2] = 0; base.m[3][3] = 0;
    a.m[1][1] = 1; a.m[1][2] = 1; a.m[1][3] = 1;
}
void qpow(int p) {
    while (p) {
        if (p & 1) a = a * base;
        base = base * base;
        p >>= 1;
    }
}
int main() {
    scanf("%d", &T);
    while (T--) {
        init();
        scanf("%d", &n);
        if (n <= 3) {
            puts("1"); continue ;
        }
        qpow(n - 3);
        printf("%lld\n", a.m[1][1]);
    }
    return 0;
}