[學習筆記]矩陣快速冪
阿新 • • 發佈：2018-12-24
入門的題目就不放了……我放一些進階的題目好了
1、P哥破解密碼
比賽的時候還是 \(ljc1301\) 首切了以後再給我們切的 \(\%\%\%\)
沒有連續的三個 \(A\)，矩陣為
\(1,1,1\)
\(1,0,0\)
\(0,1,0\)
\(Code\ Below:\)
#include <cstdio> #define Int long long Int x[999][999]; Int ans[999][999]; Int dx[999][999]; Int n,k; const int p=19260817; void ans_cf() { for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) dx[i][j]=ans[i][j],ans[i][j]=0; for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) for(int k=1;k<=n;k++) ans[i][j]=(ans[i][j]+(x[i][k]*dx[k][j])%p)%p; } void x_cf() { for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) dx[i][j]=x[i][j],x[i][j]=0; for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) for(int k=1;k<=n;k++) x[i][j]=(x[i][j]+(dx[i][k]*dx[k][j])%p)%p; } void fast_pow(Int w) { while(w) { if(w%2==1) ans_cf(); w/=2; x_cf(); } } int main() { int t; scanf("%d",&t); while(t--){ n=3;scanf("%d",&k); ans[1][1]=x[1][1]=1; ans[1][2]=x[1][2]=1; ans[1][3]=x[1][3]=1; ans[2][1]=x[2][1]=1; ans[2][2]=x[2][2]=0; ans[2][3]=x[2][3]=0; ans[3][1]=x[3][1]=0; ans[3][2]=x[3][2]=1; ans[3][3]=x[3][3]=0; fast_pow(k-1); printf("%d\n",(ans[1][1]+ans[2][1]+ans[3][1])%p); } return 0; }
好早以前的碼風了……不要在意
2、[JLOI2015]有意義的字串
思路還是比較好想的，可以得到 \((\frac{b+\sqrt{d}}{2})^n+(\frac{b+\sqrt{d}}{2})^n\) 必為整數。又題目給定 \(b^2\leq d<(b+1)^2\)，所以 \(|(\frac{b-\sqrt{d}}{2})^n|<1\)
現在我們只用考慮 \((\frac{b-\sqrt{d}}{2})^n>0\) 的情況
首先，\(b\not =\sqrt{d}\)，其次，\(n\%2=0\)。所以我們在矩陣快速冪中特判一下這種情況就好了
不過這題還要龜速乘，時間複雜度 \(O(\log^2n)\)
\(Code\ Below:\)
#include <bits/stdc++.h> #define int long long #define ull unsigned long long using namespace std; const int p=7528443412579576937; int n,b,d,ans; inline int add(int x,int y){ return (1ull*x+1ull*y)%p; } int mul(int a,int b){ int ret=0; for(;b;b>>=1,a=add(a,a)) if(b&1) ret=add(ret,a); return ret; } struct Matrix{ int mat[2][2]; Matrix(){ memset(mat,0,sizeof(mat)); } void clear(){ for(int i=0;i<2;i++) mat[i][i]=1; } }; Matrix operator * (const Matrix &a,const Matrix &b){ Matrix c; for(int i=0;i<2;i++) for(int j=0;j<2;j++) for(int k=0;k<2;k++) c.mat[i][j]=add(c.mat[i][j],mul(a.mat[i][k],b.mat[k][j])); return c; } Matrix operator ^ (Matrix a,int b){ Matrix ret;ret.clear(); for(;b;b>>=1,a=a*a) if(b&1) ret=ret*a; return ret; } Matrix x,y; signed main() { scanf("%lld%lld%lld",&b,&d,&n); if(n==0){ printf("1\n"); return 0; } x.mat[0][0]=b;x.mat[0][1]=2; y.mat[0][0]=b;y.mat[0][1]=1; y.mat[1][0]=(d-b*b)/4; y=y^(n-1);x=x*y;ans=x.mat[0][0]; if(b*b!=d&&n%2==0) ans=add(ans-1,p); printf("%lld\n",ans); return 0; }
3、[HNOI2011]數學作業
話說以前 \(noip\) 模擬賽出過這題，只不過我沒做出來。。。
這道題真的太好了，徹底讓我懂了矩陣快速冪。
小矩陣為
\(ans,10^x,1\)
大矩陣為
\(10^{x+1},0,0\)
\(1,1,0\)
\(0,1,1\)
然後就對於每一個 \(10^x\) 都重新矩陣快速冪一下，時間複雜度 \(O(\log^2 n)\)
\(Code\ Below:\)
// luogu-judger-enable-o2
#include <bits/stdc++.h>
#define ll long long
using namespace std;
ll n,p,h[20],sum;
struct matrix{
ll x[4][4];
matrix(){
memset(x,0,sizeof(x));
}
}a,b,c;
matrix operator * (matrix a,matrix b){
matrix c;
for(ll i=1;i<=3;i++)
for(ll j=1;j<=3;j++)
for(ll k=1;k<=3;k++)
c.x[i][j]=(c.x[i][j]+a.x[i][k]%p*b.x[k][j]%p)%p;
return c;
}
matrix fast_pow(matrix a,ll b){
matrix ret=a;b--;
for(;b;b>>=1,a=a*a)
if(b&1) ret=ret*a;
return ret;
}
matrix mul(matrix a,matrix b){
matrix c;
for(ll i=1;i<=3;i++)
for(ll j=1;j<=3;j++)
c.x[1][i]=(c.x[1][i]+a.x[1][j]%p*b.x[j][i]%p)%p;
return c;
}
int main()
{
scanf("%lld%lld",&n,&p);
if(n<=9){
ll ans=0;
for(ll i=1;i<=n;i++)
ans=(ans*10+i)%p;
printf("%lld\n",ans);
return 0;
}
for(ll i=1;i<=9;i++)
sum=(sum*10+i)%p;
h[1]=10;
for(ll i=2;i<=18;i++)
h[i]=h[i-1]*10;
a.x[1][1]=sum;a.x[1][2]=h[1];a.x[1][3]=1;
b.x[1][1]=h[2];b.x[2][1]=b.x[2][2]=b.x[3][2]=b.x[3][3]=1;
for(ll i=2;i<=18;i++){
if(h[i-1]<=n&&n<h[i]){
c=fast_pow(b,n-h[i-1]+1);
a=mul(a,c);
printf("%lld\n",a.x[1][1]%p);
return 0;
}
c=fast_pow(b,h[i]-h[i-1]);
a=mul(a,c);
a.x[1][2]=h[i]%p;a.x[1][3]=1;
b.x[1][1]=b.x[1][1]*10%p;
}
return 0;
}4、[SCOI2009]迷路
\(Code\ Below:\)
#include <bits/stdc++.h>
using namespace std;
const int p=2009;
int n,T;
struct matrix{
int m[110][110];
matrix(){
memset(m,0,sizeof(m));
}
void clear(){
for(int i=1;i<=n;i++)
m[i][i]=1;
}
}f;
matrix operator * (matrix a,matrix b){
matrix c;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
for(int k=1;k<=n;k++)
c.m[i][j]=(c.m[i][j]+a.m[i][k]*b.m[k][j])%p;
return c;
}
matrix operator ^ (matrix a,int b){
matrix ret;ret.clear();
for(;b;b>>=1,a=a*a)
if(b&1) ret=ret*a;
return ret;
}
inline int pos(int i,int j){
return i+j*(n/9);
}
int main()
{
scanf("%d%d",&n,&T);
n*=9;
int x;
for(int i=1;i<=n/9;i++){
for(int j=1;j<=8;j++)
f.m[pos(i,j)][pos(i,j-1)]=1;
for(int j=1;j<=n/9;j++){
scanf("%1d",&x);
f.m[i][pos(j,x-1)]=1;
}
}
f=f^T;
printf("%d\n",f.m[1][n/9]);
return 0;
}5、[HNOI2008]GT考試
用 \(kmp\) 解決矩陣快速冪問題，我也是頭一回遇到
\(Code\ Below:\)
#include <bits/stdc++.h>
using namespace std;
int n,m,mod,fail[30];
char s[30];
struct Matrix{
int mat[30][30];
Matrix(){
memset(mat,0,sizeof(mat));
}
void clear(){
for(int i=0;i<m;i++) mat[i][i]=1;
}
};
Matrix operator * (Matrix a,Matrix b){
Matrix c;
for(int i=0;i<m;i++)
for(int j=0;j<m;j++)
for(int k=0;k<m;k++)
c.mat[i][j]=(c.mat[i][j]+a.mat[i][k]*b.mat[k][j])%mod;
return c;
}
Matrix operator ^ (Matrix a,int b){
Matrix ret;ret.clear();
for(;b;b>>=1,a=a*a)
if(b&1) ret=ret*a;
return ret;
}
Matrix a;
inline void read(int &x){
x=0;int f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
while(isdigit(ch)){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
if(f==-1) x=-x;
}
void print(int x){
if(x<0){putchar('-');x=-x;}
if(x>9) print(x/10);
putchar(x%10+'0');
}
int main()
{
read(n),read(m),read(mod);
scanf("%s",s);
fail[0]=fail[1]=0;
int p=0;
for(int i=1;i<m;i++){
while(p&&s[p]!=s[i]) p=fail[p];
fail[i+1]=(s[p]==s[i])?++p:p;
}
for(int i=0;i<m;i++){
for(int j=0;j<10;j++){
p=i;
while(p&&s[p]-'0'!=j) p=fail[p];
(s[p]-'0'==j)?++p:0;
if(p<m) a.mat[i][p]++;
}
}
a=a^n;
int ans=0;
for(int i=0;i<m;i++)
ans=(ans+a.mat[0][i])%mod;
printf("%d\n",ans);
return 0;
}6、能量採集
公開賽的題，不是莫比烏斯反演
思路比較好想的，接下來可以矩陣快速冪了
考慮優化：二進位制拆分成 \(31\) 份是可以通過此題的
我們算了一下塊取 \(2^{16}\) 空間過不了，不過沒關係，我們將塊設為 \(2^{11}\)
把 \(\%\) 換成 \(-\) 也能提速哦！
\(Code\ Below:\)
#include <bits/stdc++.h>
#define ll long long
#define res register int
using namespace std;
const int p=998244353;
int n,m,q,H[40];
inline void add(int &x,int y){
x=x+y>=p?x+y-p:x+y;
}
struct Matrix{
int n,m,mat[60][60];
Matrix(){
memset(mat,0,sizeof(mat));
}
void clear(){
for(res i=0;i<n;i++) mat[i][i]=1;
}
};
inline Matrix operator * (const Matrix &a,const Matrix &b){
Matrix c;
for(res i=0;i<a.n;i++)
for(res j=0;j<b.m;j++)
for(res k=0;k<a.m;k++) add(c.mat[i][j],1ll*a.mat[i][k]*b.mat[k][j]%p);
c.n=a.n;c.m=b.m;
return c;
}
Matrix f[40],a;
inline void read(int &x){
x=0;int f=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
while(isdigit(ch)){x=(x<<3)+(x<<1)+ch-'0';ch=getchar();}
if(f==-1) x=-x;
}
void print(int x){
if(x<0){putchar('-');x=-x;}
if(x>9) print(x/10);
putchar(x%10+'0');
}
int fast_pow(int a,int b){
res ret=1;
for(;b;b>>=1,a=1ll*a*a%p)
if(b&1) ret=1ll*ret*a%p;
return ret;
}
int main()
{
H[0]=1;
for(res i=1;i<=31;i++) H[i]=H[i-1]<<1;
read(n),read(m),read(q);
a.n=a.m=1;f[0].n=f[0].m=n;
res x,y;f[0].clear();
for(res i=0;i<n;i++) read(a.mat[i][0]);
for(res i=1;i<=m;i++){
read(x),read(y);
f[0].mat[y-1][x-1]++;
}
for(res i=0;i<n;i++){
x=0;
for(res j=0;j<n;j++) add(x,f[0].mat[j][i]);
x=fast_pow(x,p-2);
for(res j=0;j<n;j++) f[0].mat[j][i]=1ll*f[0].mat[j][i]*x%p;
}
for(res i=1;i<=31;i++) f[i]=f[i-1]*f[i-1];
res sum;
while(q--){
read(x);Matrix ans=a;sum=0;
for(res i=0;i<=31;i++)
if(x&H[i]) ans=f[i]*ans;
for(res i=0;i<n;i++) sum^=ans.mat[i][0];
print(sum%p);putchar('\n');
}
return 0;
}7、塊速遞推
聽 \(shadowice1984\) 說最優解用的是生成函式，但是我只會 \(O(1)\) 快速冪
也就是分塊快速冪唄，矩陣快速冪完將 \(1\sim blo\) 的所有矩陣算出來，然後兩個矩陣相乘就好了。因為只用兩個矩陣相乘只詢問一個數，我們大可以不用把所有矩陣都算出來，算一個點就好了。
能不取模就不要取模，這樣快了很多！！！
矩陣：
\(233,1\)
\(666,0\)
考慮如何卡常？
我們發現這個數列迴圈節為 \(p-1\)，所以將 \(n\%(p-1)\)。然後那個塊的大小可以選用 \(2^16\)，這樣位運算會很快，不容易被卡掉。
\(Code\ Below:\)
// luogu-judger-enable-o2
#pragma GCC optimize("Ofast")
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#pragma GCC optimize("inline")
#pragma GCC optimize("-fgcse")
#pragma GCC optimize("-fgcse-lm")
#pragma GCC optimize("-fipa-sra")
#pragma GCC optimize("-ftree-pre")
#pragma GCC optimize("-ftree-vrp")
#pragma GCC optimize("-fpeephole2")
#pragma GCC optimize("-ffast-math")
#pragma GCC optimize("-fsched-spec")
#pragma GCC optimize("unroll-loops")
#pragma GCC optimize("-falign-jumps")
#pragma GCC optimize("-falign-loops")
#pragma GCC optimize("-falign-labels")
#pragma GCC optimize("-fdevirtualize")
#pragma GCC optimize("-fcaller-saves")
#pragma GCC optimize("-fcrossjumping")
#pragma GCC optimize("-fthread-jumps")
#pragma GCC optimize("-funroll-loops")
#pragma GCC optimize("-fwhole-program")
#pragma GCC optimize("-freorder-blocks")
#pragma GCC optimize("-fschedule-insns")
#pragma GCC optimize("inline-functions")
#pragma GCC optimize("-ftree-tail-merge")
#pragma GCC optimize("-fschedule-insns2")
#pragma GCC optimize("-fstrict-aliasing")
#pragma GCC optimize("-fstrict-overflow")
#pragma GCC optimize("-falign-functions")
#pragma GCC optimize("-fcse-skip-blocks")
#pragma GCC optimize("-fcse-follow-jumps")
#pragma GCC optimize("-fsched-interblock")
#pragma GCC optimize("-fpartial-inlining")
#pragma GCC optimize("no-stack-protector")
#pragma GCC optimize("-freorder-functions")
#pragma GCC optimize("-findirect-inlining")
#pragma GCC optimize("-fhoist-adjacent-loads")
#pragma GCC optimize("-frerun-cse-after-loop")
#pragma GCC optimize("inline-small-functions")
#pragma GCC optimize("-finline-small-functions")
#pragma GCC optimize("-ftree-switch-conversion")
#pragma GCC optimize("-foptimize-sibling-calls")
#pragma GCC optimize("-fexpensive-optimizations")
#pragma GCC optimize("-funsafe-loop-optimizations")
#pragma GCC optimize("inline-functions-called-once")
#pragma GCC optimize("-fdelete-null-pointer-checks")
#include <bits/stdc++.h>
#define ll long long
#define ull unsigned long long
using namespace std;
const int maxn=(1ll<<16)+10;
const int p=1e9+7;
const ull base=(1ll<<16)-1;
int T,now;ull n;
inline int add(int x,int y){
x+=y;x>=p?x-=p:0;
return x;
}
struct Matrix{
int mat[2][2];
Matrix(){
memset(mat,0,sizeof(mat));
}
void clear(){
mat[0][0]=mat[1][1]=1;
}
};
inline Matrix operator * (const Matrix &a,const Matrix &b){
Matrix c;
register int i,j,k;
for(i=0;i<2;++i)
for(j=0;j<2;++j)
for(k=0;k<2;++k) c.mat[i][j]=add(c.mat[i][j],1ll*a.mat[i][k]*b.mat[k][j]%p);
return c;
}
inline Matrix operator ^ (Matrix a,ull b){
Matrix ret;ret.clear();
for(;b;b>>=1,a=a*a)
if(b&1) ret=ret*a;
return ret;
}
Matrix a[maxn],b[maxn],x,y;
ull SA,SB,SC;
void init(){scanf("%llu%llu%llu",&SA,&SB,&SC);}
inline ull Rand(){
SA^=SA<<32,SA^=SA>>13,SA^=SA<<1;
ull t=SA;SA=SB,SB=SC,SC^=t^SA;
return SC;
}
int main()
{
register int i;
x.mat[0][0]=233;x.mat[0][1]=1;x.mat[1][0]=666;
y=x^(1ll<<16);
a[0].clear();a[1]=x;
for(i=2;i<=base;++i) a[i]=a[i-1]*x;
b[0].clear();b[1]=y;
for(i=2;i<=base;++i) b[i]=b[i-1]*y;
scanf("%d",&T);
init();
Matrix c,d;
for(i=1;i<=T;++i){
n=(Rand()-1)%(p-1);
const int s=(int)n&base,t=((int)n>>16)&base;
now^=(1ll*a[s].mat[0][0]*b[t].mat[0][0]+1ll*a[s].mat[0][1]*b[t].mat[1][0])%p;
}
printf("%d\n",now%p);
return 0;
}