矩阵加速(一些题目)
要做矩阵加速首先要明白矩阵快速幂
明白了矩阵快速幂后就可以通过矩阵的乘法来加速递推
例1 斐波那契数列
题目大意
求第 n n n项斐波那契数模 1000000007 1000000007 1000000007的值。 n n n在 l o n g    l o n g long\,\,long longlong范围内
输入
一行， n n n
输出
f ( n ) % 1000000007 f(n) {\%} 1000000007 f(n)%1000000007的值
样例输入
10
样例输出
55
分析
首先，直接递推肯定超时，可以观察如下两个式子
f [ n ] = f [ n − 1 ] ∗ 1 + f [ n − 2 ] ∗ 1 f[n] = f[n-1] *1 + f[n-2] *1 f[n]=f[n−1]∗1+f[n−2]∗1
f [ n − 1 ] = f [ n − 1 ] ∗ 1 + f [ n − 2 ] ∗ 0 f[n-1] = f[n-1]*1 +f[n-2]*0 f[n−1]=f[n−1]∗1+f[n−2]∗0
把这两个线性方程表示成矩阵，提取系数，可得到如下表达式
[
1
1
1
0
]
×
[
f
(
n
−
1
)
f
(
n
−
2
)
]
=
[
f
(
n
)
f
(
n
−
1
)
]
\left[ \begin{matrix} 1 & 1 \\ 1 & 0 \end{matrix} \right] \times \left[ \begin{matrix} f(n-1) \\ f(n-2) \end{matrix} \right] =\left[ \begin{matrix} f(n) \\ f(n-1) \end{matrix} \right]
[1110]×[f(n−1)f(n−2)]=[f(n)f(n−1)]
令
A = [ 1 1 1 0 ] A =\left[ \begin{matrix} 1 & 1\\ 1 & 0 \end{matrix} \right] A=[1110]
可以得到
A n × [ f ( 1 ) f ( 0 ) ] = [ f ( n + 1 ) f ( n ) ] A^n \times \left[ \begin{matrix} f(1) \\ f(0) \end{matrix} \right] =\left[ \begin{matrix} f(n+1) \\ f(n) \end{matrix} \right] An×[f(1)f(0)]=[f(n+1)f(n)]
其中 A n A^n An可以通过矩阵快速幂得到，这样就可以通过矩阵 A A A来加速递推运算
代码
#include <iostream>
#include <vector>
#include <cstdio>
using namespace std;
typedef long long ll;
typedef vector<ll> vec;
typedef vector<vec> mat;
const ll mod = 1000000007;
inline mat mult (mat A, mat B) {
int row = A.size(), col = B[0].size(), l = B.size();
mat C(row, vec(col));
for(int i = 0; i < row; i++)
for(int j = 0; j < col; j++)
for(int k = 0; k < l; k++)
C[i][j] = C[i][j] % mod + A[i][k] * B[k][j] % mod;
return C;
}
mat Power(mat &A, ll m) { //求快速幂
ll n = A.size();
mat res(n, vec(n));
for(int i = 0; i < n; i++)
res[i][i] = 1;
while(m) {
if(m & 1) res = mult(res, A);
m >>= 1;
A = mult(A, A);
}
return res;
}
int main () {
ll n;
mat a(2, vec(2));
cin >> n;
a[0][0] = 1, a[0][1] = 1; \\通过矩阵a来进行加速
a[1][0] = 1, a[1][1] = 0;
mat t = Power(a, n);
cout << t[1][0] % mod << "\n";
return 0;
}
例二 数列
题目大意
a [ 1 ] = a [ 2 ] = a [ 3 ] = 1 a[1]=a[2]=a[3]=1 a[1]=a[2]=a[3]=1， a [ x ] = a [ x − 3 ] + a [ x − 1 ] ( x > 3 ) a[x]=a[x-3]+a[x-1] (x>3) a[x]=a[x−3]+a[x−1](x>3)
求 a a a数列的第 n n n项对 1000000007 （ 1 0 9 + 7 ） 1000000007（10^9+7） 1000000007（109+7）取余的值。
输入
第一行一个整数 T T T，表示询问个数。
以下 T T T行，每行一个正整数 n n n。
输出
每行输出一个非负整数表示答案。
输入样例
3
6
8
10
输出样例
4
9
19
分析
同样的用矩阵加速递推，关键在于找到用哪个矩阵来进行加速，由于这个递推公式跨度有 3 3 3个长度，所以可以写 3 3 3个递推式如下
f [ n ] = 1 × f [ n − 1 ] + 0 × f [ n − 2 ] + 1 × f [ n − 3 ] f[n] = 1 \times f[n-1] + 0 \times f[n-2] + 1 \times f[n-3] f[n]=1×f[n−1]+0×f[n−2]+1×f[n−3]
f [ n − 1 ] = 1 × f [ n − 1 ] + 0 × f [ n − 2 ] + 0 × f [ n − 3 ] f[n-1] = 1 \times f[n-1] + 0 \times f[n-2] +0 \times f[n-3] f[n−1]=1×f[n−1]+0×f[n−2]+0×f[n−3]
f [ n − 2 ] = 0 × f [ n − 1 ] + 1 × f [ n − 2 ] + 0 × f [ n − 3 ] f[n-2] = 0 \times f[n-1] + 1 \times f[n-2] +0 \times f[n-3] f[n−2]=0×f[n−1]+1×f[n−2]+0×f[n−3]
同样提取相应的系数，转换为矩阵乘法
[
1
0
1
1
0
0
0
1
0
]
×
[
f
(
n
−
1
)
f
(
n
−
2
)
f
(
n
−
3
)
]
=
[
f
(
n
)
f
(
n
−
1
)
f
(
n
−
2
)
]
\left[ \begin{matrix} 1 & 0 & 1\\ 1 & 0 & 0\\ 0 & 1 & 0 \end{matrix} \right] \times \left[ \begin{matrix} f(n-1) \\ f(n-2) \\ f(n-3) \end{matrix} \right] =\left[ \begin{matrix} f(n)\\ f(n-1)\\ f(n-2) \end{matrix} \right]
⎣⎡110001100⎦⎤×⎣⎡f(n−1)f(n−2)f(n−3)⎦⎤=⎣⎡f(n)f(n−1)f(n−2)⎦⎤
令
A
=
[
1
0
1
1
0
0
0
1
0
]
A = \left[ \begin{matrix} 1 & 0 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \end{matrix} \right]
A=⎣⎡110001100⎦⎤
因此
A
n
−
3
×
[
f
(
3
)
f
(
2
)
f
(
1
)
]
=
[
f
(
n
)
f
(
n
−
1
)
f
(
n
−
2
)
]
                
[
f
(
3
)
f
(
2
)
f
(
1
)
]
=
[
1
1
1
]
A^{n-3} \times \left[ \begin{matrix} f(3)\\ f(2)\\ f(1) \end{matrix} \right] =\left[ \begin{matrix} f(n)\\ f(n-1)\\ f(n-2) \end{matrix} \right] \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \left[ \begin{matrix} f(3)\\ f(2)\\ f(1) \end{matrix} \right] =\left[ \begin{matrix} 1\\ 1\\ 1 \end{matrix} \right]
An−3×⎣⎡f(3)f(2)f(1)⎦⎤=⎣⎡f(n)f(n−1)f(n−2)⎦⎤⎣⎡f(3)f(2)f(1)⎦⎤=⎣⎡111⎦⎤
代码
#include <iostream>
#include <cctype>
#include <vector>
#include <cstdio>
#define ll long long
#define in(a) a = read();
using namespace std;
inline int read () {
int f = 1, x = 0; char c = getchar();
while(!isdigit(c)) {if(c == '-') f = -1; c = getchar();}
while(isdigit(c)) {x = (x << 3) + (x << 1) + c - 48; c = getchar();}
return f*x;
}
typedef vector<ll> vec;
typedef vector<vec> mat;
const ll mod = 1e9+7;
inline mat mul(mat A, mat B) {
int row = A.size(), col = B[0].size(), l = B.size();
mat C(row, vec(col));
for(int i = 0; i < row; i++)
for(int j = 0; j < col; j++)
for(int k = 0; k < l; k++)
C[i][j] = C[i][j] % mod + A[i][k] * B[k][j] % mod;
return C;
}
inline mat powMod (mat A, int m) {
int n = A.size();
mat res(n, vec(n));
for(int i = 0; i < n; i++)
res[i][i] = 1;
while(m) {
if(m & 1) res = mul(res, A);
m >>= 1;
A = mul(A, A);
}
return res;
}
int main () {
int T1, n;
mat F(3, vec(3)), e(3, vec(1)); \\F矩阵用来加速，e表示初始的f(3)，f(2)，f(1)
F[0][0] = 1, F[0][1] = 0, F[0][2] = 1;
F[1][0] = 1, F[1][1] = 0, F[1][2] = 0;
F[2][0] = 0, F[2][1] = 1, F[2][2] = 0;
e[0][0] = 1, e[1][0] = 1, e[2][0] = 1;
in(T1);
while (T1 --) {
in(n);
n -= 3;
if(n <= 0) printf("1\n");
else {
mat t = powMod(F, n);
t = mul(t, e);
printf("%lld\n", t[0][0] % mod);
}
}
return 0;
}
例3 Blocks
题目大意
N N N个方块排成一列，有红、黄、绿、蓝四种颜色， 对方块染色，求当红色方块数和绿色方块数同为偶数个数的时候有多少种染色方案。结果对 10007 10007 10007取余数
输入
第一行一个 T T T，表示测试数据数量 ( 1 ≤ T ≤ 100 ) (1 \leq T \leq100) (1≤T≤100)
接下来 T T T行，每行一个 n n n，表示方块数量。 ( 1 ≤ N ≤ 1 0 9 ) (1 \leq N \leq 10^9) (1≤N≤109)
输出
对于每一个 n n n，输出一个结果，每个结果换一行
输入样例
2
1
2
输出样例
2
6
分析
首先要找出递推关系式，然后才能用矩阵加速
-
红绿都是偶数的方案数为 a i a_i ai
-
红绿其中有一个偶数另一个是奇数的方案数为 b i b_i bi
-
红绿都是奇数的方案数为 c i c_i ci
考虑 a i a_i ai，那么 a i + 1 = 2 × a i + b i a_{i+1} = 2 \times a_i + b_i ai+1=2×ai+bi
(递推 a i + 1 a_{i+1} ai+1块的时候，为保证红绿都为偶数。上一种红绿都是偶数的时候，新来一块可以选择黄蓝染色，同时还有上一种红绿只有一个为奇数，那么新来的一块只能染成这个奇数的颜色)
考虑 b i b_i bi，那么 b i + 1 = 2 × a i + 2 × b i + 2 × c i b_{i+1} = 2 \times a_i + 2 \times b_i + 2 \times c_i bi+1=2×ai+2×bi+2×ci
（递推 b i + 1 b_{i+1} bi+1块的时候，为了保证红绿只有一奇。上一种红绿都为偶数，新来一块可以选择红或者绿， 上一种红绿有一个为奇数，新来的一块可以染黄或者蓝色，如果上一种红绿都为奇数，新来的一块染红或者绿）
考虑 c i c_i ci，那么 c i + 1 = 2 × c i + b i c_{i+1} = 2 \times c_i+b_i ci+1=2×ci+bi
(递推 c i + 1 c_{i+1} ci+1的时候，为了保证红绿都为奇。上一种红绿都为偶数不满足，上一种红绿有一个为奇数的时候，新来的一块对红绿中偶数颜色的进行染色即可，上一种红绿都为奇数的时候，新来的一块可以染黄或者蓝色)
因此得到
a i + 1 = 2 × a i + 1 × b i + 0 × c i a_{i+1} = 2 \times a_i + 1 \times b_i + 0 \times c_i ai+1=2×ai+1×bi+0×ci
b i + 1 = 2 × a i + 2 × b i + 2 × c i b_{i+1} = 2 \times a_i + 2 \times b_i + 2 \times c_i bi+1=2×ai+2×bi+2×ci
c i + 1 = 0 × a i + 1 × b i + 2 × c i c_{i+1} = 0 \times a_i+1 \times b_i + 2 \times c_i ci+1=0×ai+1×bi+2×ci
提取相应系数，转换为矩阵乘法
[
2
1
0
2
2
2
0
1
2
]
×
[
a
(
n
−
1
)
b
(
n
−
1
)
c
(
n
−
1
)
]
=
[
a
(
n
)
b
(
n
)
c
(
n
)
]
\left[ \begin{matrix} 2 & 1 & 0\\ 2 & 2 & 2\\ 0 & 1 & 2 \end{matrix} \right] \times \left[ \begin{matrix} a(n-1)\\ b(n-1)\\ c(n-1) \end{matrix} \right]= \left[ \begin{matrix} a(n)\\ b(n)\\ c(n) \end{matrix} \right]
⎣⎡220121022⎦⎤×⎣⎡a(n−1)b(n−1)c(n−1)⎦⎤=⎣⎡a(n)b(n)c(n)⎦⎤
令
A
=
[
2
1
0
2
2
2
0
1
2
]
A = \left[ \begin{matrix} 2 & 1 & 0\\ 2 & 2 & 2\\ 0 & 1 & 2 \end{matrix} \right]
A=⎣⎡220121022⎦⎤
则
A
n
×
[
a
(
0
)
b
(
0
)
c
(
0
)
]
=
[
a
(
n
)
b
(
n
)
c
(
n
)
]
                 
[
a
(
0
)
b
(
0
)
c
(
0
)
]
=
[
1
0
0
]
A^n \times \left[ \begin{matrix} a(0)\\ b(0)\\ c(0) \end{matrix} \right]= \left[ \begin{matrix} a(n)\\ b(n)\\ c(n) \end{matrix} \right]\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \left[ \begin{matrix} a(0)\\ b(0)\\ c(0) \end{matrix} \right]= \left[ \begin{matrix} 1\\ 0\\ 0 \end{matrix} \right]
An×⎣⎡a(0)b(0)c(0)⎦⎤=⎣⎡a(n)b(n)c(n)⎦⎤⎣⎡a(0)b(0)c(0)⎦⎤=⎣⎡100⎦⎤
代码
#include <iostream>
#include <cctype>
#include <vector>
#include <cstdio>
#define ll long long
#define in(a) a = read();
using namespace std;
inline int read () {
int f = 1, x = 0; char c = getchar();
while(!isdigit(c)) {if(c == '-') f = -1; c = getchar();}
while(isdigit(c)) {x = (x << 3) + (x << 1) + c - 48; c = getchar();}
return f*x;
}
typedef vector<ll> vec;
typedef vector<vec> mat;
const ll mod = 10007;
inline mat mul(mat A, mat B) {
int row = A.size(), col = B[0].size(), l = B.size();
mat C(row, vec(col));
for(int i = 0; i < row; i++)
for(int j = 0; j < col; j++)
for(int k = 0; k < l; k++)
C[i][j] = C[i][j] % mod + A[i][k] * B[k][j] % mod;
return C;
}
inline mat powMod (mat A, int m) {
int n = A.size();
mat res(n, vec(n));
for(int i = 0; i < n; i++)
res[i][i] = 1;
while(m) {
if(m & 1) res = mul(res, A);
m >>= 1;
A = mul(A, A);
}
return res;
}
int main () {
int T1, n;
mat F(3, vec(3)), e(3, vec(1));
F[0][0] = 2, F[0][1] = 1, F[0][2] = 0;
F[1][0] = 2, F[1][1] = 2, F[1][2] = 2;
F[2][0] = 0, F[2][1] = 1, F[2][2] = 2;
e[0][0] = 1, e[1][0] = 0, e[2][0] = 0;
in(T1);
while (T1 --) {
in(n);
mat t = powMod(F, n);
t = mul(t, e);
printf("%lld\n", t[0][0] % mod);
}
return 0;
}
例4 Fibonacci Sum
题目大意
给出两个非负整数 N N N, M M M，计算 ∑ i = N M F i \sum_{i=N}^{M}F_i ∑i=NMFi， F i F_i Fi表示斐波那契数列第 i i i项， F 0 = 0 F_0 = 0 F0=0， F 1 = 1 F_1 = 1 F1=1
输入
第一行为一个整数 T T T,代表数据组数，接着是两个正整数 N N N, M M M
输入
T T T行，为每组数据的表达式结果，结果对 1000000007 1000000007 1000000007取模。
输入样例
3
0 3
3 5
10 19
输出样例
4
10
10857
分析
∑ i = N M F i = S ( M ) − S ( N − 1 ) \sum_{i=N}^{M}F_i = S(M)-S(N-1) ∑i=NMFi=S(M)−S(N−1)
同样要找到加速矩阵，观察下面3个表达式
s n + 2 = 1 × s n + 1 + 1 × f ( n ) + 1 × f ( n + 1 ) s_{n+2}=1 \times s_{n+1} + 1 \times f(n)+1 \times f(n+1) sn+2=1×sn+1+1×f(n)+1×f(n+1)
f ( n + 1 ) = 0 × S n + 1 + 0 × f ( n ) + 1 × f ( n + 1 ) f(n+1) = 0 \times S_{n+1}+0 \times f(n)+1 \times f(n+1) f(n+1)=0×Sn+1+0×f(n)+1×f(n+1)
f ( n + 2 ) = 0 × S n + 1 + 1 × f ( n ) + 1 × f ( n + 1 ) f(n+2)=0 \times S_{n+1} +1 \times f(n)+1 \times f(n+1) f(n+2)=0×Sn+1+1×f(n)+1×f(n+1)
同样的提取系数
[
1
1
1
0
0
1
0
1
1
]
×
[
s
n
+
1
f
(
n
)
f
(
n
+
1
)
]
=
[
s
n
+
2
f
(
n
+
1
)
f
(
n
+
2
)
]
\left[ \begin{matrix} 1 & 1 & 1\\ 0 & 0 & 1\\ 0 & 1 & 1 \end{matrix} \right] \times \left[ \begin{matrix} s_{n+1}\\ f(n)\\ f(n+1) \end{matrix} \right] =\left[ \begin{matrix} s_{n+2}\\ f(n+1)\\ f(n+2) \end{matrix} \right]
⎣⎡100101111⎦⎤×⎣⎡sn+1f(n)f(n+1)⎦⎤=⎣⎡sn+2f(n+1)f(n+2)⎦⎤
因此，用来加速的矩阵 A A A
A = [ 1 1 1 0 0 1 0 1 1 ] A= \left[ \begin{matrix} 1 & 1 & 1\\ 0 & 0 & 1\\ 0 & 1 & 1 \end{matrix} \right] A=⎣⎡100101111⎦⎤
可得
A
n
−
2
×
[
S
2
f
(
1
)
f
(
2
)
]
=
[
S
n
f
(
n
−
1
)
f
(
n
)
]
A^{n-2} \times \left[ \begin{matrix} S_2\\ f(1)\\ f(2) \end{matrix} \right]= \left[ \begin{matrix} S_n\\ f(n-1)\\ f(n) \end{matrix} \right]
An−2×⎣⎡S2f(1)f(2)⎦⎤=⎣⎡Snf(n−1)f(n)⎦⎤
坑点：ans = s(m) - s(n-1)，当用快速幂求出s(m)和s(n-1)的时候，直接用s(m) - s(n-1)会出错！因为在取模的过程中可能导致s(m) < s(n-1)，因此ans = (s(m) + mod - s(n-1) % mod) % mod
代码(使用结构体存数组)
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#define ll long long
#define in(a) a = read()
#define wr(a) write(a), putchar(10)
using namespace std;
const ll mod = 1000000007;
inline ll read () {
ll x = 0, f = 1; char c = getchar();
while(!isdigit(c)) {if(c == '-') f = -1; c = getchar();}
while(isdigit(c)) {x = (x << 3) + (x << 1) + c - 48; c = getchar();}
return f*x;
}
inline void write(ll x) {
if(x < 0) {
putchar('-');
x = -x;
}
if(x >= 10) write(x / 10);
putchar(x%10+48);
}
struct mat{
int row, col;
ll a[105][105];
mat () { //默认构造函数，初始化为0
for(int i = 0; i < 105; i++)
for(int j = 0; j < 105; j++)
a[i][j] = 0;
}
void init () { //设置为单位矩阵
for(int i = 0; i < row; i++)
a[i][i] = 1;
}
mat operator * (const mat &T) { //重装*运算符，使得能够直接结构体相乘
mat res;
res.row = row, res.col = T.col;
for(int i = 0; i < res.row; i++)
for(int j = 0; j < res.col; j++)
for(int k = 0; k < col; k++)
res.a[i][j] = res.a[i][j] % mod + a[i][k] * T.a[k][j] % mod;
return res;
}
};
mat Quick_Pow (mat A, ll b) {
mat res;
res.row = A.row, res.col = A.row;
res.init();
while (b) {
if (b & 1) res = res * A;
b >>= 1;
A = A * A;
}
return res;
}
int main () {
mat f, ini, ans1, ans2;
f.a[0][0] = 1, f.a[0][1] = 1, f.a[0][2] = 1;
f.a[1][0] = 0, f.a[1][1] = 0, f.a[1][2] = 1;
f.a[2][0] = 0, f.a[2][1] = 1, f.a[2][2] = 1;
f.row = f.col = 3;
ini.row = 3, ini.col = 1;
ini.a[0][0] = 2, ini.a[1][0] = 1, ini.a[2][0] = 1;
ll T, n, m, t1, t2;
in(T);
while (T--) {
in(n), in(m);
if(m > 2) { //求前m项的和
ans1 = Quick_Pow(f, m-2) * ini;
t1 = ans1.a[0][0] % mod;
}else
t1 = m;
if(n > 3) { //求前n-1项的和
ans2 = Quick_Pow(f, n - 3) * ini;
t2 = ans2.a[0][0];
}else
t2 = n-1;
if(t2 < 0) t2 = 0;
wr((t1 + mod - t2 % mod)% mod);
}
return 0;
}
例5 Matrix Power Series
题目大意
给定 n × n n \times n n×n的矩阵 A A A和正整数 k k k和 m m m。假设 S = A 1 + A 2 + . . . + A k S = A^1 + A^2 + ...+A^k S=A1+A2+...+Ak，输出 S S S的各元素对 M M M取余的答案
其中 k ≤ 1 0 9 , n ≤ 30 , m ≤ 1 0 4 k \leq 10^9,n \leq 30, m \leq 10^4 k≤109,n≤30,m≤104
输入
第一行包含 3 3 3个整数， n 、 k 、 m n、k、m n、k、m，接下来 n n n行每行包含 n n n个不超过 32768 32768 32768的非负整数
输出
输出 S S S模 m m m后的矩阵 S S S
样例输入
2 2 4
0 1
1 1
样例输出
1 2
2 3
分析
观察两个式子
S n = A + A 2 + A 3 + . . . + A n S_n = A + A^2 + A^3 + ... + A^n Sn=A+A2+A3+...+An
S n + 1 = ( A + A 2 + A 3 + . . . + A n ) × A + A S_{n+1} = (A + A^2 + A^3 + ... + A^n) \times A + A Sn+1=(A+A2+A3+...+An)×A+A
可得
S n + 1 = S n × A + A S_{n+1} = S_n \times A +A Sn+1=Sn×A+A
变换成矩阵为
[
A
A
0
1
]
×
[
S
n
−
1
1
]
=
[
S
n
1
]
      
令
  
F
=
[
A
A
0
1
]
\left[ \begin{matrix} A & A\\ 0 & 1 \end{matrix} \right] \times \left[ \begin{matrix} S_{n-1}\\ 1 \end{matrix} \right]= \left[ \begin{matrix} S_n\\ 1 \end{matrix} \right]\,\,\,\,\,\, 令\,\,F= \left[ \begin{matrix} A & A\\ 0 & 1 \end{matrix} \right]
[A0A1]×[Sn−11]=[Sn1]令F=[A0A1]
得到
F
k
−
1
×
[
S
1
1
]
=
[
S
k
1
]
F^{k-1} \times \left[ \begin{matrix} S_1\\ 1 \end{matrix} \right] = \left[ \begin{matrix} S_k\\ 1 \end{matrix} \right]
Fk−1×[S11]=[Sk1]
所以
F
F
F就是其中的加速的矩阵，但是关键在于
A
A
A也是矩阵，这就变成了矩阵里面套矩阵，其实并没有关系，可以将矩阵
A
A
A在矩阵
F
F
F中展开，由于
A
A
A的长宽是
n
n
n，所以
F
F
F的宽长应该是
2
×
n
2 \times n
2×n，其中
1
1
1变换为长宽为
n
n
n的单位矩阵，
0
0
0则是长宽为
n
n
n的全
0
0
0的矩阵
最后的答案只需要取长宽为 n n n的右上角的部分即可
举个例子
A
=
[
0
1
1
1
]
−
>
               
F
=
[
0
1
0
1
1
1
1
1
0
0
1
0
0
0
0
1
]
               
[
S
1
1
]
=
[
A
1
1
]
=
[
0
1
1
1
1
0
0
1
]
A = \left[ \begin{matrix} 0 & 1\\ 1 & 1 \end{matrix} \right] -> \,\,\,\,\,\,\,\,\,\,\,\,\,\,\, F = \left[ \begin{matrix} 0 & 1 & 0 & 1\\ 1 & 1 & 1 & 1\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{matrix} \right] \,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \left[ \begin{matrix} S_1\\ 1 \end{matrix} \right]= \left[ \begin{matrix} A_1\\ 1 \end{matrix} \right]= \left[ \begin{matrix} 0 & 1\\ 1 & 1\\ 1 & 0\\ 0 & 1 \end{matrix} \right]
A=[0111]−>F=⎣⎢⎢⎡0100110001101101⎦⎥⎥⎤[S11]=[A11]=⎣⎢⎢⎡01101101⎦⎥⎥⎤
因此只需要做好相应的变换即可
#include <iostream>
#include <cstdio>
#include <cstring>
#include <cctype>
#define N 35
#define ll long long
#define rep(i, j) for (int i = 0; i < j; i++)
using namespace std;
ll n, k, mod;
struct mat {
int row, col;
ll a[N*2][N*2];
mat () { //默认全为0
for(int i = 0; i < N*2; i++)
for(int j = 0; j < N*2; j++)
a[i][j] = 0;
}
void init () { //初始化乘单位矩阵
for(int i = 0; i < row; i++)
a[i][i] = 1;
}
mat operator * (const mat &b) {
mat c;
c.row = row, c.col = b.col;
for(int i = 0; i < c.row; i++)
for(int j = 0; j < c.col; j++)
for(int k = 0; k < col; k++)
c.a[i][j] = c.a[i][j] % mod + a[i][k] * b.a[k][j] % mod;
return c;
}
};
mat Quick_Pow (mat A, ll b) { //快速幂
mat res;
res.row = res.col = A.row;
res.init();
while(b) {
if(b & 1) res = res * A;
b >>= 1;
A = A * A;
}
return res;
}
int main () {
scanf("%lld %lld %lld", &n, &k, &mod);
mat F, S, ans; // ans = S*F^(k-1)
F.row = F.col = n*2;
S.row = n*2, S.col = n;
for(int i = 0; i < n; i++) //左上部分为A
for(int j = 0; j < n; j++)
scanf("%lld", &F.a[i][j]);
for(int i = 0; i < n; i++) //右上部分为A
for(int j = n; j < 2*n; j++)
F.a[i][j] = F.a[i][j-n];
for(int i = n; i < 2*n; i++) //左下部分默认为0不用管，右下部分置为单位矩阵
F.a[i][i] = 1;
for(int i = 0; i < n; i++)
for(int j = 0; j < n; j++)
S.a[i][j] = F.a[i][j];
for(int i = n; i < 2*n; i++)
S.a[i][i-n] = 1;
ans = Quick_Pow(F, k-1);
ans = ans * S;
for(int i = 0; i < ans.row / 2; i++) { //输出上半部分
for(int j = 0; j < ans.col; j++)
printf("%lld ", ans.a[i][j] % mod);
printf("\n");
}
return 0;
}
小结
一般地，对于 m m m项递推式，如果递推式为
a n + m = ∑ i = 0 m − 1 b i a n + i a_{n+m} = \sum_{i=0}^{m-1}b_ia_{n+i} an+m=∑i=0m−1bian+i
则可把递推式写成如下矩阵
[
a
n
+
m
a
n
+
m
−
1
.
.
.
a
n
+
1
]
=
[
b
m
−
1
.
.
.
b
1
b
0
1
.
.
.
0
0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0
.
.
.
1
0
]
[
a
n
+
m
−
1
a
n
+
m
−
2
.
.
.
a
n
]
\left[ \begin{matrix} a_{n+m}\\ a_{n+m-1}\\ .\\.\\.\\ a_{n+1} \end{matrix} \right] = \left[ \begin{matrix} b_{m-1} & ... & b_1 & b_0\\ 1 & ... & 0 & 0\\ . & ... & . & .\\ . & ... & . & .\\ . & ... & . & .\\ 0 & ... & 1 & 0 \end{matrix} \right] \left[ \begin{matrix} a_{n+m-1}\\ a_{n+m-2}\\ .\\.\\.\\ a_{n} \end{matrix} \right]
⎣⎢⎢⎢⎢⎢⎢⎡an+man+m−1...an+1⎦⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎡bm−11...0..................b10...1b00...0⎦⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎡an+m−1an+m−2...an⎦⎥⎥⎥⎥⎥⎥⎤
时间复杂度为
O
(
m
3
l
o
g
  
n
)
O(m^3log\,\,n)
O(m3logn)，如果递推式有常数项
c
c
c则变成
[
a
n
+
m
a
n
+
m
−
1
.
.
.
a
n
+
1
1
]
=
[
b
m
−
1
.
.
.
b
1
b
0
c
1
.
.
.
0
0
0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
0
.
.
.
1
0
0
0
.
.
.
0
0
1
]
[
a
n
+
m
−
1
a
n
+
m
−
2
.
.
.
a
n
1
]
\left[ \begin{matrix} a_{n+m}\\ a_{n+m-1}\\ .\\.\\.\\ a_{n+1}\\ 1 \end{matrix} \right] = \left[ \begin{matrix} b_{m-1} & ... & b_1 & b_0 & c\\ 1 & ... & 0 & 0 & 0\\ . & ... & . & . & .\\ . & ... & . & . & .\\ . & ... & . & . & .\\ 0 & ... & 1 & 0 & 0\\ 0 & ... & 0 & 0 & 1 \end{matrix} \right] \left[ \begin{matrix} a_{n+m-1}\\ a_{n+m-2}\\ .\\.\\.\\ a_{n}\\ 1 \end{matrix} \right]
⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡an+man+m−1...an+11⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡bm−11...00.....................b10...10b00...00c0...01⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡an+m−1an+m−2...an1⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
