矩阵快速幂（矩阵套矩阵版）
「一本通 6.5 例 3」Fibonacci 前 n 项和
求 \(Fibonacci\) 数列前 \(n\) 项的和对 \(m\) 取模后的结果
输入 \(n\) 和 \(m\)
#include<bits/stdc++.h>
#define int long long
#define ll long long
#define I 1.0000000
#define fd(i,a,b) for(int i=a,_i=b;i<=_i;i=-~i)
#define bd(i,a,b) for(int i=a,_i=b;i>=_i;i=~-i)
#define Db(a,b,ans) cout<<a<<';'<<b<<','<<ans<<endl;>
using namespace std;
const int N=35/*+509*/,M=5;
int mod,m,n,k,l;
struct node
{
int m[N][N];
friend node operator+(node x,node y);
friend node operator*(node x,node y);
friend node operator%(node x,int mod);
friend node operator^(node x,int y);
node() {fd(i,1,N-1) fd(j,1,N-1) m[i][j]=0;}
node(bool flag) {if(!flag) return; fd(i,1,N-1) fd(j,1,N-1) m[i][j]=0;fd(i,1,N-1) m[i][i]=1;}
node operator+=(node x) {return (*this)=(*this)+x;}
node operator*=(node x) {return (*this)=(*this)*x;}
node operator%=(int mod) {return (*this)=(*this)%mod;}
node operator^=(int x) {return (*this)=(*this)^x;}
}A,O,E(1),B;
node operator+(node x,node y) {fd(i,1,l) fd(j,1,l) x.m[i][j]=(x.m[i][j]+y.m[i][j])%mod; return x;}
node operator*(node x,node y) {node res; fd(i,1,l) fd(j,1,l) fd(k,1,l) res.m[i][j]=(res.m[i][j]+x.m[i][k]*y.m[k][j]%mod+mod)%mod; return res;}
node operator%(node x,int mod) {fd(i,1,l) fd(j,1,l) x.m[i][j]=x.m[i][j]%mod; return x;}
node operator^(node x,int y) {node ans(1),base=x;for(;y;y>>=1) {if(y&1)ans=(ans*base)%mod;base=(base*base)%mod;} return ans%mod;}
struct node2
{
node m[M][M];
friend node2 operator*(node2 x,node2 y);
friend node2 operator^(node2 x,int y);
node2() {fd(i,1,M-1) fd(j,1,M-1) m[i][j]=O;}
node2(bool flag) {fd(i,1,M-1) fd(j,1,M-1) m[i][j]=O; if(!flag) return; fd(i,1,M-1) m[i][i]=E;}
node2 operator*=(node2 x) {return (*this)=(*this)*x;}
node2 operator^=(int x) {return (*this)=(*this)^x;}
}ans,base;
node2 operator*(node2 x,node2 y) {node2 res; fd(i,1,2) fd(j,1,2) fd(k,1,2) res.m[i][j]=(res.m[i][j]+x.m[i][k]*y.m[k][j])%mod; return res;}
node2 operator^(node2 x,int y) {node2 ans(1),base=x;for(;y;y>>=1) {if(y&1)ans*=base;base*=base;} return ans;}
namespace quicker
{
inline int read()
{
int num=0; bool flag=1; char c=getchar();
while(!isdigit(c)){if(c=='-') flag=0;c=getchar();}
while(isdigit(c)){num=(num<<1)+(num<<3)+(c-'0');c=getchar();}
return flag?num:-num;
}
inline void write(int x)
{
if(x<0) putchar('-'),x=-x;
if(x>9) write(x/10);
putchar((x%10)^48);
}
inline void tie0()
{
ios::sync_with_stdio(0);
cin.tie(0); cout.tie(0);
}
}
using namespace quicker;
signed main()
{
#ifdef FJ
freopen(".in","m",stdin);
freopen(".out","w",stdout);
quicker::tie0();
#endif
cin>>m>>mod;
l=3;
A.m[1][1]=1;
A.m[2][1]=1;
A.m[2][2]=1;
A.m[2][3]=1;
A.m[3][2]=1;
A^=(m-1);
// fd(i,1,n)
// {
// fd(j,1,n) cout<<A.m[i][j]<<' ';
// cout<<endl;
// }
B.m[1][1]=1;
B.m[1][2]=1;
B.m[1][3]=1;
B*=A;
cout<<B.m[1][1];
return 0;
}
Luogu P10502
Matrix Power Series
题面翻译
【题目描述】
给定一个 \(n×n\) 矩阵 \(A\) 和一个正整数 \(k\)，找出和 \(S=A+A^2 +A^3 +...+A^k\)。
【输入格式】
输入包含一个测试用例。输入的第一行包含三个正整数 \(n\)（\(n \le 30\)）、\(k\)（\(k \le 10^9\)）和 \(m\)（\(m < 10^4\)）。接下来的 \(n\) 行每行包含 \(n\) 个小于 32,768 的非负整数，按行主序给出 \(A\) 的元素。
【输出格式】
以与给定 \(A\) 相同的方式输出 \(S\) 的元素对 \(m\) 取模。
翻译来自于：ChatGPT
样例 #1
样例输入 #1
2 2 4
0 1
1 1
样例输出 #1
1 2
2 3
代码：
#include<bits/stdc++.h>
#define int long long
#define ll long long
#define I 1.0000000
#define fd(i,a,b) for(int i=a,_i=b;i<=_i;i=-~i)
#define bd(i,a,b) for(int i=a,_i=b;i>=_i;i=~-i)
#define Db(a,b,ans) cout<<a<<';'<<b<<','<<ans<<endl;>
using namespace std;
const int N=35/*+509*/,M=5;
int mod,n,k,l;
struct node
{
int m[N][N];
friend node operator+(node x,node y);
friend node operator*(node x,node y);
friend node operator%(node x,int mod);
friend node operator^(node x,int y);
node() {fd(i,1,N-1) fd(j,1,N-1) m[i][j]=0;}
node(bool flag) {if(!flag) return; fd(i,1,N-1) fd(j,1,N-1) m[i][j]=0;fd(i,1,N-1) m[i][i]=1;}
node operator+=(node x) {return (*this)=(*this)+x;}
node operator*=(node x) {return (*this)=(*this)*x;}
node operator%=(int mod) {return (*this)=(*this)%mod;}
node operator^=(int x) {return (*this)=(*this)^x;}
}A,O,E(1);
node operator+(node x,node y) {fd(i,1,l) fd(j,1,l) x.m[i][j]=(x.m[i][j]+y.m[i][j])%mod; return x;}
node operator*(node x,node y) {node res; fd(i,1,l) fd(j,1,l) fd(k,1,l) res.m[i][j]=(res.m[i][j]+x.m[i][k]*y.m[k][j]%mod+mod)%mod; return res;}
node operator%(node x,int mod) {fd(i,1,l) fd(j,1,l) x.m[i][j]=x.m[i][j]%mod; return x;}
node operator^(node x,int y) {node ans(1),base=x;for(;y;y>>=1) {if(y&1)ans=(ans*base)%mod;base=(base*base)%mod;} return ans%mod;}
struct node2
{
node m[M][M];
friend node2 operator*(node2 x,node2 y);
friend node2 operator^(node2 x,int y);
node2() {fd(i,1,M-1) fd(j,1,M-1) m[i][j]=O;}
node2(bool flag) {fd(i,1,M-1) fd(j,1,M-1) m[i][j]=O; if(!flag) return; fd(i,1,M-1) m[i][i]=E;}
node2 operator*=(node2 x) {return (*this)=(*this)*x;}
node2 operator^=(int x) {return (*this)=(*this)^x;}
}ans,base;
node2 operator*(node2 x,node2 y) {node2 res; fd(i,1,2) fd(j,1,2) fd(k,1,2) res.m[i][j]=(res.m[i][j]+x.m[i][k]*y.m[k][j])%mod; return res;}
node2 operator^(node2 x,int y) {node2 ans(1),base=x;for(;y;y>>=1) {if(y&1)ans*=base;base*=base;} return ans;}
namespace quicker
{
inline int read()
{
int num=0; bool flag=1; char c=getchar();
while(!isdigit(c)){if(c=='-') flag=0;c=getchar();}
while(isdigit(c)){num=(num<<1)+(num<<3)+(c-'0');c=getchar();}
return flag?num:-num;
}
inline void write(int x)
{
if(x<0) putchar('-'),x=-x;
if(x>9) write(x/10);
putchar((x%10)^48);
}
inline void tie0()
{
ios::sync_with_stdio(0);
cin.tie(0); cout.tie(0);
}
}
using namespace quicker;
signed main()
{
#ifdef FJ
freopen(".in","m",stdin);
freopen(".out","w",stdout);
quicker::tie0();
#endif
cin>>n>>k>>mod;
l=n+1;
fd(i,1,n)
fd(j,1,n)
cin>>A.m[i][j];
ans.m[1][2]=base.m[1][1]=E;
base.m[2][1]=base.m[2][2]=A;
ans*=base^k;
fd(i,1,n)
{
fd(j,1,n) cout<<ans.m[1][1].m[i][j]<<' ';
cout<<endl;
}
// A.m[1][1]=1;
// A.m[2][1]=1;
// A.m[1][2]=1;
// A^=n;
// cout<<A.m[1][2];
return 0;
}
upd on 2024.11.28
【模板】矩阵快速幂
可以直接用上面的板子，但是我重新板刷了一下，写了一份动态开空间的
#include<bits/stdc++.h>
#define int long long
#define ll long long
#define fd(i,a,b) for(int i=(a);i<=(b);i=-~i)
#define bd(i,a,b) for(int i=(a);i>=(b);i=~-i)
#define db(x) cout<<"DEBUG "<<#x<<" = "<<x<<endl;
#define endl '\n'
using namespace std;
//#define SIZE (1<<20)
//char In[SIZE],Out[SIZE],*p1=In,*p2=In,*p3=Out;
//#define getchar() (p1==p2&&(p2=(p1=In)+fread(In,1,SIZE,stdin),p1==p2)?EOF:*p1++)
inline int read()
{
int x=0,f=1;char c=getchar();
while(c<'0'||c>'9'){if(c=='-') f=-1;c=getchar();}
while(c>='0'&&c<='9'){x=x*10+(c-48),c=getchar();}
return x*f;
}
template<typename _T>
inline void write(_T x)
{
static _T sta[35];int top=0;
if(x<0) putchar('-'),x=-x;
do{sta[top++]=x%10,x/=10;}while(x);
while(top) putchar(sta[--top]+'0');
}
const int N=2e5+509,M=1e6+509,mod=1e9+7;
int n,m;
struct matrix
{
int siz=0;
vector< vector<int> > a;
matrix(int SIZ=0,int op=0){siz=SIZ;a.resize(siz);fd(i,0,siz-1) a[i].resize(siz);fd(i,0,siz-1) fd(j,0,siz-1) a[i][i]=0;if(op) fd(i,0,siz-1) a[i][i]=1;}
inline void resize(int siz){a.resize(siz);fd(i,0,siz-1) a[i].resize(siz);}
inline void print(){fd(i,0,siz-1){fd(j,0,siz-1) printf("%lld ",a[i][j]);puts("");}}
friend matrix operator*(matrix x,matrix y){matrix res(x.siz);fd(i,0,x.siz-1) fd(j,0,x.siz-1) fd(k,0,x.siz-1) (res.a[i][j]+=x.a[i][k]*y.a[k][j]%mod)%=mod;return res;};
friend matrix operator+(matrix x,matrix y){matrix res(x.siz);fd(i,0,x.siz-1) fd(j,0,x.siz-1) (res.a[i][j]+=x.a[i][j]+y.a[i][j]%mod)%=mod;return res;};
friend matrix operator%(matrix x,const int y){matrix res(x.siz);fd(i,0,x.siz-1) fd(j,0,x.siz-1) res.a[i][j]=x.a[i][j]%mod;return res;};
friend matrix operator^(matrix x,int y){matrix res(x.siz,1);while(y){if(y&1)res=res*x;x=x*x,y>>=1;}return res;};
inline void operator*=(matrix y){(*this)=(*this)*y;}
inline void operator+=(matrix y){(*this)=(*this)+y;}
inline void operator%=(int y){(*this)=(*this)%y;}
inline void operator^=(int y){(*this)=(*this)^y;}
};
signed main()
{
//#define FJ
#ifdef FJ
freopen(".in","r",stdin);
freopen(".out","w",stdout);
#else
// freopen("A.in","r",stdin);
// freopen("A.out","w",stdout);
#endif
//#define io
#ifdef io
ios::sync_with_stdio(0);
cin.tie(0); cout.tie(0);
#endif
n=read(),m=read();
matrix A(n);
fd(i,0,n-1) fd(j,0,n-1) A.a[i][j]=read();
A^=m;
A.print();
return 0;
}
本文来自博客园，作者：whrwlx，转载请注明原文链接：https://www.cnblogs.com/whrwlx/p/18276108
