矩阵快速幂
矩阵快速幂(入门)
定义及模板
对于方阵存在求幂运算的概念。对于单个数值，我们通过将指数进行二进制拆分的思想将求幂运算的复杂度降至log(n)，这种思想同样可以应用到方阵的求幂运算中。事实上，方阵的快速幂与单个数值的快速幂的写法完全一致，只需要对于乘法、取模进行运算符重载即可。
下面给出矩阵快速幂的模板
#include <cstdio>
#include <cstdlib>
using namespace std;
#define N 5+1
// 最大矩阵规模,根据题目修改
typedef long long ll;
const ll mod = 1e9+7; // 取模
struct Matrix{
int x;
int y;
ll a[N][N];
Matrix(){};
Matrix(int m,int n){
x = m;
y = n;
for(int i = 1; i <= x; i++){
for(int j = 1; j <= y; j++){
a[i][j] = 0;
}
}
}
ll* operator[](int x){
return a[x];
}
void ones(){
for(int i = 1; i <= x; i++){
a[i][i] = 1;
}
}
Matrix operator *(Matrix b){
Matrix ans(x,b.y);
for(int i = 1; i <= x; i++){
for(int j = 1; j <= b.y; j++){
for(int k = 1; k <= y; k++){
ans[i][j] += a[i][k] * b[k][j];
}
ans[i][j] %= mod;
}
}
return ans;
}
Matrix operator %(ll md){
Matrix ans = *this;
for(int i = 1; i <= x; i++){
for(int j = 1; j <= y; j++){
ans[i][j] %= md;
}
}
return ans;
}
Matrix operator ^(ll p){
Matrix ans(x,x);
Matrix base = *this;
ans.ones();
while(p){
if(p&1){
ans = ans * base;
ans = ans % mod;
}
base = base * base;
base = base % mod;
p >>= 1;
}
return ans;
}
};
int main(){
return 0;
}
应用:矩阵加速
矩阵加速用于线性方程的加速递推
level 1:简单线性递推
- e.g1 : 求Fibonacci(n)%mod,mod=1e9+7
考虑到斐波那契数列的递推公式
f(n) = f(n-1) + f(n-2);
我们可以将数列的第n项、第n-1项、第n-2项组合成矩阵形式
// 数列第n、n-1、n-2的组合
[f(n) f(n-1) f(n-2)]
// 数列第n+1、n、n-1的组合
[f(n+1) f(n) f(n-1)]
这两个矩阵满足关系式子
\[\left[
\begin{matrix}
a_{n}&a_{n-1}&a_{n-2}
\end{matrix}
\right]
*
\left[
\begin{matrix}
1 & 1 & 0\\
1 & 0 & 1\\
0 & 0 & 0
\end{matrix}
\right]
=
\left[
\begin{matrix}
a_{n+1}&a_{n}&a_{n-1}
\end{matrix}
\right]
\]
进而得出
\[\left[
\begin{matrix}
2 & 1 & 1
\end{matrix}
\right]
*
\left[
\begin{matrix}
1 & 1 & 0\\
1 & 0 & 1\\
0 & 0 & 0
\end{matrix}
\right]^{n-2}
=
\left[
\begin{matrix}
a_{n}&a_{n-1}&a_{n-2}
\end{matrix}
\right]
\]
因此,要求fibonacci(n),通过矩阵加速可以将O(n)降至O(log(n))
- e.g2: 矩阵加速-数列
本题关系式为
\[\left[
\begin{matrix}
2&1&1&1
\end{matrix}
\right]
*
\left[
\begin{matrix}
1&0&0&0\\
0&1&1&0\\
1&0&0&1\\
0&1&0&0
\end{matrix}
\right]^{n-4}
=
\left[
\begin{matrix}
a_{n}&a_{n-1}&a_{n-2}&a_{n-3}
\end{matrix}
\right]
\]
#include <cstdio>
#include <cstdlib>
using namespace std;
#define N 5+1
// 最大矩阵规模,根据题目修改
typedef long long ll;
const ll mod = 1e9+7;
struct Matrix{
int x;
int y;
ll a[N][N];
Matrix(){};
Matrix(int m,int n){
x = m;
y = n;
for(int i = 1; i <= x; i++){
for(int j = 1; j <= y; j++){
a[i][j] = 0;
}
}
}
ll* operator[](int x){
return a[x];
}
void ones(){
for(int i = 1; i <= x; i++){
a[i][i] = 1;
}
}
Matrix operator *(Matrix b){
Matrix ans(x,b.y);
for(int i = 1; i <= x; i++){
for(int j = 1; j <= b.y; j++){
for(int k = 1; k <= y; k++){
ans[i][j] += a[i][k] * b[k][j];
}
ans[i][j] %= mod;
}
}
return ans;
}
Matrix operator %(ll md){
Matrix ans = *this;
for(int i = 1; i <= x; i++){
for(int j = 1; j <= y; j++){
ans[i][j] %= md;
}
}
return ans;
}
Matrix operator ^(ll p){
Matrix ans(x,x);
Matrix base = *this;
ans.ones();
while(p){
if(p&1){
ans = ans * base;
ans = ans % mod;
}
base = base * base;
base = base % mod;
p >>= 1;
}
return ans;
}
};
int main(){
int t;
scanf("%d",&t);
Matrix k(4,4);
k[1][1] = 1;
k[2][2] = 1;
k[2][3] = 1;
k[3][1] = 1;
k[3][4] = 1;
k[4][2] = 1;
Matrix s(1,4);
s[1][1] = 2;
s[1][2] = 1;
s[1][3] = 1;
s[1][4] = 1;
while(t--){
ll n;
scanf("%lld",&n);
if(n <= 3){
printf("1\n");
}else{
Matrix ans = s * (k^(n-4));
printf("%lld\n",ans[1][1]);
}
}
return 0;
}
level 2: 带有求和项的简单递推
由于矩阵加速题难点在于关系矩阵的寻找，而不在于代码实现，为了节省篇幅，之后的题解不再出现具体代码实现。
e.g1:HDU-3306 Another kind of Fibonacci
本题的关系矩阵略难推导，应该注意到
\[\begin{align}
& S_{n+1} = S_{n} + A_{n+1}^2\\
& A_{n+2}^2 = (xA_{n+1} + yA_{n})^2 = x^2A_{n+1}^2 + y^2A_n^2 + 2xyA_{n+1}A_n \\
& A_{n+2}A_{n+1} = (xA_{n+1}+yA_{n})A_{n+1} = xA_{n+1}^2 + yA_{n+1}A_n
\end{align}
\]
因此可以得到
\[\left[
\begin{matrix}
S_n & A_{n+1}^2 & A_n^2 & A_{n+1}A_n
\end{matrix}
\right]
*
\left[
\begin{matrix}
1 & 0 & 0 & 0\\
1 & x^2 & 1 & x\\
0 & y^2 & 0 & 0\\
0 & 2xy & 0 & y
\end{matrix}
\right]
=
\left[
\begin{matrix}
S_{n+1}^2 & A_{n+2}^2 & A_{n+1}^2 & A_{n+2}A_{n+1}
\end{matrix}
\right]
\]
本题小结
- 对于求和项有 \(S_{n+1}-S_n = ...\)
- 如果递推公式中出现了别的项,无法通过已有项的线性组合得出(如\(A_{n+1}A_n\)项无法通过\(A_{n+1}\)和\(A_n\)的线性组合得出)，那么考虑更改递推公式或在向量中添加新项
level 3: 打表寻找递推公式
e.g1:luogu5004 专心OI - 跳房子
本题不难想到二维DP的解法(虽然会爆)
伪代码如下,其中i表示在格子上停留的次数，j表示当前走在第j个格子上
procedure(DP):
int dp[N][N]
int n,m
read n,m
int maxstep = (n+1)/(m+1)
for i:= 1 to n
dp[1][i] = 1
for i:= 2 to n
for j:= 1 to n
int begin = 1
int end = j-1+m
for k:= begin to end
dp[i][j] += dp[i-1][k]
ans = 1
for i:= 1 to maxstep
for j:= 1 to n
ans += dp[i][j]
return ans
用这个思路来运行较小的n，m值，打表寻找规律:
- m=1
n = 0 m = 1 ans = 1
n = 1 m = 1 ans = 2
n = 2 m = 1 ans = 3
n = 3 m = 1 ans = 5
n = 4 m = 1 ans = 8
n = 5 m = 1 ans = 13
n = 6 m = 1 ans = 21
n = 7 m = 1 ans = 34
n = 8 m = 1 ans = 55
n = 9 m = 1 ans = 89
n = 10 m = 1 ans = 144
n = 11 m = 1 ans = 233
n = 12 m = 1 ans = 377
n = 13 m = 1 ans = 610
n = 14 m = 1 ans = 987
n = 15 m = 1 ans = 1597
n = 16 m = 1 ans = 2584
n = 17 m = 1 ans = 4181
n = 18 m = 1 ans = 6765
n = 19 m = 1 ans = 10946
n = 20 m = 1 ans = 17711
可以发现m=1时ans符合斐波那契数列的增长规律
\[m = 1 \to
\left\{
\begin{align}
& A_{0} = 1\\
& A_{1} = 2\\
& A_{2} = 3\\
& A_{n} = A_{n-1}+A_{n-2} (n \geqslant 3)
\end{align}
\right.
\]
- m = 2
n = 0 m = 2 ans = 1
n = 1 m = 2 ans = 2
n = 2 m = 2 ans = 3
n = 3 m = 2 ans = 4
n = 4 m = 2 ans = 6
n = 5 m = 2 ans = 9
n = 6 m = 2 ans = 13
n = 7 m = 2 ans = 19
n = 8 m = 2 ans = 28
n = 9 m = 2 ans = 41
n = 10 m = 2 ans = 60
n = 11 m = 2 ans = 88
n = 12 m = 2 ans = 129
n = 13 m = 2 ans = 189
n = 14 m = 2 ans = 277
n = 15 m = 2 ans = 406
n = 16 m = 2 ans = 595
n = 17 m = 2 ans = 872
n = 18 m = 2 ans = 1278
n = 19 m = 2 ans = 1873
n = 20 m = 2 ans = 2745
\[m = 2 \to
\left\{
\begin{align}
& A_{0} = 1\\
& A_{1} = 2\\
& A_{2} = 3\\
& A_{3} = 4\\
& A_{n} = A_{n-1}+A_{n-3} (n \geqslant 4)
\end{align}
\right.
\]
- m = 4
n = 0 m = 4 ans = 1
n = 1 m = 4 ans = 2
n = 2 m = 4 ans = 3
n = 3 m = 4 ans = 4
n = 4 m = 4 ans = 5
n = 5 m = 4 ans = 6
n = 6 m = 4 ans = 8
n = 7 m = 4 ans = 11
n = 8 m = 4 ans = 15
n = 9 m = 4 ans = 20
n = 10 m = 4 ans = 26
n = 11 m = 4 ans = 34
n = 12 m = 4 ans = 45
n = 13 m = 4 ans = 60
n = 14 m = 4 ans = 80
n = 15 m = 4 ans = 106
n = 16 m = 4 ans = 140
n = 17 m = 4 ans = 185
n = 18 m = 4 ans = 245
n = 19 m = 4 ans = 325
n = 20 m = 4 ans = 431
\[m = 4 \to
\left\{
\begin{align}
& A_{0} = 1\\
& A_{1} = 2\\
& A_{2} = 3\\
& A_{3} = 4\\
& A_{4} = 5\\
& A_{5} = 6\\
& A_{n} = A_{n-1}+A_{n-5} (n \geqslant 6)
\end{align}
\right.
\]
找到规律
\[\forall m \to
\left\{
\begin{align}
& A_n = n+1 \ \ \ \ \ \ \ (0 \leqslant n \leqslant m+1)\\
& A_{n} = A_{n-1}+A_{n-m-1} \ \ \ \ \ \ \ \ (n \geqslant m+2)
\end{align}
\right.
\]
因此有
本题关系式为
\[\left[
\begin{matrix}
A_n & A_{n-1} &\cdots & A_{n-m}
\end{matrix}
\right]*
\left[
\begin{matrix}
1 & 1 & 0 & 0 &\cdots & 0 & 0\\
0 & 0 & 1 & 0 &\cdots & 0 & 0\\
0 & 0 & 0 & 1 & \cdots & 0&0\\
0 & 0 & 0 & 0 & 1 & \cdots &0\\
\vdots & \vdots & \vdots & \vdots & \vdots & \ddots & 1& \\
1 & 0 & 0 & 0 & \cdots & 0 & 0\\
\end{matrix}
\right]
=
\left[
\begin{matrix}
A_{n+1} & A_{n} \cdots & A_{n-m+1}
\end{matrix}
\right]
\]
\[\R_{1\times (m+1)} * \R_{(m+1)\times(m+1)} = \R_{1 \times (m+1)}
\]
代码
#include <cstdio>
#include <cstdlib>
using namespace std;
#define N 18+1
// 最大矩阵规模,根据题目修改
typedef long long ll;
const ll mod = 1e9+7; // 取模
struct Matrix{
int x;
int y;
ll a[N][N];
Matrix(){};
Matrix(int m,int n){
x = m;
y = n;
for(int i = 1; i <= x; i++){
for(int j = 1; j <= y; j++){
a[i][j] = 0;
}
}
}
ll* operator[](int x){
return a[x];
}
void ones(){
for(int i = 1; i <= x; i++){
a[i][i] = 1;
}
}
Matrix operator *(Matrix b){
Matrix ans(x,b.y);
for(int i = 1; i <= x; i++){
for(int j = 1; j <= b.y; j++){
for(int k = 1; k <= y; k++){
ans[i][j] += a[i][k] * b[k][j];
}
ans[i][j] %= mod;
}
}
return ans;
}
Matrix operator %(ll md){
Matrix ans = *this;
for(int i = 1; i <= x; i++){
for(int j = 1; j <= y; j++){
ans[i][j] %= md;
}
}
return ans;
}
Matrix operator ^(ll p){
Matrix ans(x,x);
Matrix base = *this;
ans.ones();
while(p){
if(p&1){
ans = ans * base;
ans = ans % mod;
}
base = base * base;
base = base % mod;
p >>= 1;
}
return ans;
}
};
int dp[N][N];
int main(){
ll n,m;
scanf("%lld%lld",&n,&m);
Matrix a(1,m+1);
a[1][m+1] = 1;
for(int i = m; i >= 1; i--){
a[1][i] = a[1][i+1]+1;
}
Matrix b(m+1,m+1);
b[1][1] = 1;
b[m+1][1] = 1;
for(int i = 1; i <= m; i++){
b[i][i+1] = 1;
}
if(n < m){
printf("%lld\n",a[1][m+1-n]);
}else{
Matrix ans = a * (b^(n-m));
printf("%lld\n",ans[1][1]);
}
// system("pause");
return 0;
}
level4: 数论+DP+矩阵加速
e.g1: P3746 [六省联考2017]组合数问题
本题如果具备组合数递推的前置知识就不算难题，关于组合数的前置知识请看:浅谈组合数相关性质
根据组合数的递推性质
\[C_{m}^{n} = C_{m-1}^{n} + C_{m-1}^{n-1}
\]
它其实就是DP方程，表示前m个物品选择n个物品的方案数。根据最后一个物品选或不选分别对应出\(C_{m-1}^n\)和\(C_{m-1}^{n-1}\)
写成DP形式即为
\[dp[i][j] = dp[i-1][j] + dp[i-1][j-1]
\]
本题要求的为
\[\sum_{i=0}^{+\infin}C_{nk}^{nk+r}
\]
它所表示的是从nk个物品中选择m个物品，其中m%k = r，求对于所有的m的方案数的总和
我们设\(dp[i][j]\)表示从前i个物品中选择出 m % k = j 个物品的方案数
\[dp[i][0] = dp[i-1][0] + dp[i-1][k-1] \\
dp[i][j] = dp[i-1][j] + dp[i-1][j-1]
\]
因此递推关系式为：
\[\left[
\begin{matrix}
1 & 0 & 0 & 0 & \cdots & \cdots & 0
\end{matrix}
\right]
*
\left[
\begin{matrix}
1 & 1 & 0 & 0 & 0 & \cdots & 0 \\
0 & 1 & 1 & 0 & 0 & \cdots & 0 \\
0 & 0 & 1 & 1 & 0 & \cdots & 0 \\
0 & 0 & 0 & 1 & 1 & \cdots & 0 \\
\vdots & \vdots&\vdots& \vdots&\vdots &\cdots & 0 \\
\vdots & \vdots&\vdots& \vdots&\vdots &1 & 0 \\
\vdots & \vdots&\vdots& \vdots&\vdots &1 & 1 \\
1 & 0 & 0 & 0 & 0 & \cdots & 1
\end{matrix}
\right]^{nk}
=
\left[
\begin{matrix}
f_{nk,0} & f_{nk,1} & f_{nk,2} & \cdots & \cdots & f_{nk,k-1}
\end{matrix}
\right]
\]
最后，k=1需要特判，k=1时应该为
\[\left[
\begin{matrix}
1
\end{matrix}
\right]
*
\left[
\begin{matrix}
2
\end{matrix}
\right]^{n}
=
\left[
\begin{matrix}
f_{n,0}
\end{matrix}
\right]
\]
最终答案为
ans[1][r+1]; // ans: Matrix(1,k)
#include <cstdio>
#include <cstdlib>
using namespace std;
#define N 50+5
// 最大矩阵规模,根据题目修改
typedef long long ll;
ll mod; // 取模
struct Matrix{
int x;
int y;
ll a[N][N];
Matrix(){};
Matrix(int m,int n){
x = m;
y = n;
for(int i = 1; i <= x; i++){
for(int j = 1; j <= y; j++){
a[i][j] = 0;
}
}
}
ll* operator[](int x){
return a[x];
}
void ones(){
for(int i = 1; i <= x; i++){
a[i][i] = 1;
}
}
Matrix operator *(Matrix b){
Matrix ans(x,b.y);
for(int i = 1; i <= x; i++){
for(int j = 1; j <= b.y; j++){
for(int k = 1; k <= y; k++){
ans[i][j] += a[i][k] * b[k][j] % mod;
ans[i][j] %= mod;
}
}
}
return ans;
}
Matrix operator %(ll md){
Matrix ans = *this;
for(int i = 1; i <= x; i++){
for(int j = 1; j <= y; j++){
ans[i][j] %= md;
}
}
return ans;
}
Matrix operator ^(ll p){
Matrix ans(x,x);
Matrix base = *this;
base = base % mod;
ans.ones();
while(p){
if(p&1){
ans = ans * base;
ans = ans % mod;
}
base = base * base;
base = base % mod;
p >>= 1;
}
return ans;
}
};
int main(){
ll n,p,k,r;
scanf("%lld%lld%lld%lld",&n,&p,&k,&r);
mod = p;
Matrix a(1,k); // j from 0 to k-1
a[1][1] = 1; // dp[0][0] = 1
Matrix b(k,k);
for(int i = 1; i <= k; i++){
b[i][i] = 1;
b[i-1][i] = 1;
}
b[k][1] = 1;
if(k == 1){
b[1][1] = 2;
}
Matrix ans = a * (b^(n*k));
printf("%lld\n",ans[1][r+1]);
// system("pause");
return 0;
}
---- suffer now and live the rest of your life as a champion ----
