原题链接
本题需要用到矩阵加速
先介绍一下矩阵加速!
矩阵加速
1.从矩阵快速幂说起
矩阵乘法的定义是:设C = A * B,n为两个矩阵长度相等一边的长度,则$C_{ij} = \sum^{n}{k=1} A{ik} * B_{kj}。由此易得出矩阵乘法满足结合律,即(A * B) * C=A * (B * C)。那么矩阵的幂运算也可以使用快速幂,只是乘换为了矩阵乘。
此外,在矩阵乘法中,有一种特殊的矩阵,它的值为:\begin{bmatrix}
1 & 0 & … & 0 \
0 & 1 & … & 0 \
… & … & … & …\
0 & 0 & … & 1
\end{bmatrix}$
即除了对角线上的数字为1外其余皆为0的方阵,任何矩阵乘该矩阵结果均为原矩阵。我们称这个矩阵为单位矩阵,它的作用类似于整数乘法中的1.
下面贴代码:
struct Matrix {
int m[110][110];
int size;
Matrix() {
size = 0;
memset(m, 0, sizeof(m));
}
Matrix(int mode, int s) {
if (mode == 0) {
size = s;
memset(m, 0, sizeof(m));
}
if (mode == 1) {
size = s;
memset(m, 0, sizeof(m));
for (int i = 1; i <= size; ++i) {
m[i][i] = 1;
}
} //构造单位矩阵
}
Matrix operator*(Matrix& rhs) {
Matrix t;
int n = this->size;
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= n; ++j) {
for (int k = 1; k <= n; ++k) {
t.m[i][j] += this->m[i][k] * rhs.m[k][j];
}
}
}
t.size = this->size;
return t;
} //矩阵乘
void read(int s) {
size = s;
for (int i = 1; i <= size; ++i) {
for (int j = 1; j <= size; ++j) {
scanf("%lld", &this->m[i][j]);
}
}
}
void print() {
for (int i = 1; i <= size; ++i) {
for (int j = 1; j <= size; ++j) {
printf("%lld ", this->m[i][j]);
}
printf("\n");
}
}
}; //输出矩阵
Matrix fast_power_matrix(Matrix x, int n) {
Matrix ans(1, x.size);
while (n) {
if (n & 1) {
ans = ans * x;
}
x = x * x;
n >>= 1;
}
return ans;
} //矩阵快速幂
(直接复制粘贴的送你见祖宗哦=ω)
矩阵加速
讲完了矩阵快速幂,那么什么是矩阵加速呢?
对于a_{i}=a_{i-x}+a_{i-y}+a_{i-z}+…这样的递推式,我们可以构造一个矩阵\begin{bmatrix}
a_{i-x}\
a_{i-y}\
a_{i-z}\
…
\end{bmatrix}
联系上面所说到的知识,显然这里可以构造一个矩阵使上面的矩阵和构造的矩阵相乘得出以下矩阵:\begin{bmatrix}
a_{i}\
a_{i-x}\
a_{i-y}\
…
\end{bmatrix}
由此,我们可以用递推式开始的已知或易求得的几项组成的矩阵不断乘上一个构造的矩阵得到答案。又矩阵乘法的结合律可得,我们可以对构造的矩阵进行求幂运算,再乘上初始矩阵直接得到答案。
现在,以Luogu P1939为例,讲解矩阵如何进行构造。
本题的递推式为a_{i}=a_{i-1}+a_{i-3},a_{1}=a_{2}=a_{3}=1,则得初始矩阵\begin{bmatrix}
a_{i-3}\
a_{i-2}\
a_{i-1}
\end{bmatrix}
现在我们要求的一个矩阵使初始矩阵乘上它为:\begin{bmatrix}
a_{i-2}\
a_{i-1}\
a_{i-3}+a_{i-1}
\end{bmatrix}
易得该矩阵为:\begin{bmatrix}
0&1&0\
0&0&1\
1&0&1
\end{bmatrix}
解释:
设\begin{bmatrix}
a_{i-3}\
a_{i-2}\
a_{i-1}
\end{bmatrix}与\begin{bmatrix}
0&1&0\
0&0&1\
1&0&1
\end{bmatrix}相乘得矩阵C,则
C_{11}=a_{i-3} * 0+a_{i-2} * 1+a_{i-1} * 0=a_{i-2}
C_{12}=a_{i-3} * 0+a_{i-2} * 0+a_{i-1} * 1=a_{i-1}
C_{13}=a_{i-3} * 1+a_{i-2} * 0+a_{i-1}*1=a_{i-3}+a_{i-1}=a_{i}
所以得到的矩阵C为\begin{bmatrix}
a_{i-2}\
a_{i-1}\
a_{i}
\end{bmatrix}。
接下来即可对结果进行计算,如求a_{i}:\begin{bmatrix}
1\
1\
1
\end{bmatrix} * \begin{bmatrix}
0&1&0\
0&0&1\
1&0&1
\end{bmatrix}^{i-3} = \begin{bmatrix}
a_{i-2}\
a_{i-1}\
a_{i}
\end{bmatrix}
附AC代码:
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
const int MOD = 1000000007;
struct Matrix {
long long m[110][110];
int size;
Matrix() {
size = 0;
memset(m, 0, sizeof(m));
}
Matrix(int mode, int s) {
if (mode == 0) {
size = s;
memset(m, 0, sizeof(m));
}
if (mode == 1) {
size = s;
memset(m, 0, sizeof(m));
for (int i = 1; i <= size; ++i) {
m[i][i] = 1;
}
}
}
Matrix operator*(Matrix& rhs) {
Matrix t;
int n = this->size;
for (int i = 1; i <= n; ++i) {
for (int j = 1; j <= n; ++j) {
for (int k = 1; k <= n; ++k) {
t.m[i][j] += (this->m[i][k] % MOD * rhs.m[k][j] % MOD) % MOD;
t.m[i][j] %= MOD;
}
}
}
t.size = this->size;
return t;
}
void read(int s) {
size = s;
for (int i = 1; i <= size; ++i) {
for (int j = 1; j <= size; ++j) {
scanf("%lld", &this->m[i][j]);
}
}
}
void print() {
for (int i = 1; i <= size; ++i) {
for (int j = 1; j <= size; ++j) {
printf("%lld ", this->m[i][j]);
}
printf("\n");
}
}
};
Matrix fast_power_matrix(Matrix x, long long n) {
Matrix ans(1, x.size);
while (n) {
if (n & 1) {
ans = ans * x;
}
x = x * x;
n >>= 1;
}
return ans;
}
long long T, n;
int t[2][30] = { { 0 } }, res[2][30];
Matrix mt;
int main() {
scanf("%lld", &T);
mt.size = 3;
for (int i = 1; i <= 3; i++) {
t[1][i] = 1;
}
for (int i = 1; i <= T; i++) {
memset(res, 0, sizeof(res));
memset(mt.m, 0, sizeof(mt.m));
scanf("%lld", &n);
for (int i = 1; i <= 3; i++) {
for (int j = 1; j <= 3; j++) {
if (j == i + 1 || (i == 3 && j != 2)) {
mt.m[i][j] = 1;
}
}
}
if (n > 3) {
mt = fast_power_matrix(mt, n - 3);
for (int i = 1; i <= 1; ++i) {
for (int j = 1; j <= 3; ++j) {
for (int k = 1; k <= 3; ++k) {
res[i][j] += ((t[i][k] % MOD) * (mt.m[j][k] % MOD)) % MOD;
res[i][j] %= MOD;
}
}
}
printf("%d\n", res[1][3]);
} else {
printf("1\n");
}
}
return 0;
}