[学习笔记]矩阵快速幂

入门的题目就不放了……我放一些进阶的题目好了

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
",(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^2leq d<(b+1)^2)，所以 (|(frac{b-sqrt{d}}{2})^n|<1)

现在我们只用考虑 ((frac{b-sqrt{d}}{2})^n>0) 的情况

首先，(b ot =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
");
		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
",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
",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
",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
",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
",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('
');
	}
	return 0;
}

7、块速递推

听 (shadowice1984) 说最优解用的是生成函数，但是我只会 (O(1)) 快速幂

也就是分块快速幂呗，矩阵快速幂完将 (1sim 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
",now%p);
    return 0;
}