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

刺骨的言语ヽ痛彻心扉 2023-06-05 08:52 12阅读 0赞

### 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 ###

![ContractedBlock.gif][] ![ExpandedBlockStart.gif][]

1 #include <bits/stdc++.h>
     2 using namespace std;
     3 typedef long long ll;
     4 const int mod = 1e9 + 7;
     5 int T, n;
     6 struct Matrix {
     7     ll m[4][4];
     8     Matrix() {memset(m, 0, sizeof(m));}
     9     Matrix operator *(const Matrix &x) const {
    10         Matrix ret;
    11         for (register int i = 1; i <= 3; ++i)
    12             for (register int j = 1; j <= 3; ++j)
    13                 for (register int k = 1; k <= 3; ++k)
    14                     ret.m[i][j] = (ret.m[i][j] + m[i][k] * x.m[k][j] % mod) % mod; 
    15         return ret;
    16     }
    17 } a, base;
    18 inline void init() {
    19     base.m[1][1] = 1; base.m[1][2] = 1; base.m[1][3] = 0;
    20     base.m[2][1] = 0; base.m[2][2] = 0; base.m[2][3] = 1;
    21     base.m[3][1] = 1; base.m[3][2] = 0; base.m[3][3] = 0;
    22     a.m[1][1] = 1; a.m[1][2] = 1; a.m[1][3] = 1;
    23 }
    24 void qpow(int p) {
    25     while (p) {
    26         if (p & 1) a = a * base;
    27         base = base * base;
    28         p >>= 1;
    29     }
    30 }
    31 int main() {
    32     scanf("%d", &T);
    33     while (T--) {
    34         init();
    35         scanf("%d", &n);
    36         if (n <= 3) {
    37             puts("1"); continue ;
    38         }
    39         qpow(n - 3);
    40         printf("%lld\n", a.m[1][1]);
    41     }
    42     return 0;
    43 }

AC Code

转载于:https://www.cnblogs.com/shl-blog/p/11335910.html

[ContractedBlock.gif]: https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif
[ExpandedBlockStart.gif]: /images/20230601/200d5e26122741ca93d5c4a7f4062d06.png