【题解】Luogu P10502 Matrix Power Series

题意分析

给定一个 \(n \times n\) 的矩阵 \(A\) 和正整数 \(k\)，求 \(S=A^1+A^2+\cdots+A^k\)。

解题思路

求 \(A^n\) 要用到矩阵快速幂。但是 \(k \le 10^9\)，求 \(k\) 个幂会超时，所以需要用到分治的策略。

我们知道，矩阵乘法满足乘法分配律，即 \(A \times C + B \times C = (A + B) \times C\)。那么我们就可以二分 \(k\)，\(A^{mid+1}\) 加到 \(A^r\) 就等于 \(A^l\) 加到 \(A^{mid}\) 乘 \(A^{\frac{r-l+1}{2}}\)。当 \((r-l+1)\) 是偶数即 \(l\) 到 \(r\) 之间的项数为偶数时，求 \(A^l\) 加到 \(A^{mid}\) 乘 \(A^{\frac{r-l+1}{2}}\) 再加 \(A^l\) 加到 \(A^{mid}\)，而求 \(A^l\) 加到 \(A^{mid}\) 时可以继续二分，直到 \(l=r\) 为止。当 \((r-l+1)\) 是奇数时，我们把多出来的一项单独求（这里选择第 \(l\) 项），剩下的用二分求。如此一来时间复杂度就会压缩到 \(O(\log k)\)。

写一个简单的递归即可。

代码实现

#include<iostream>
#include<cstring>
#define int long long 
using namespace std;
struct Matrix{
	int mx[110][110];
}a,s,c;
int n,m,k;
Matrix operator *(const Matrix &a,const Matrix &b){//运算符重载为矩阵乘法 
	memset(c.mx,0,sizeof(c.mx));
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			for(int d=1;d<=n;d++){
				c.mx[i][j]=(c.mx[i][j]+(a.mx[i][d]*b.mx[d][j])%m)%m;
			}
		}
	}
	return c;
}
Matrix operator +(const Matrix &a,const Matrix &b){//运算符重载为矩阵加法 
	memset(c.mx,0,sizeof(c.mx));
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			c.mx[i][j]=(a.mx[i][j]+b.mx[i][j])%m;
		}
	}
	return c;
}
Matrix ksm(Matrix a,int k){//快速幂 
	Matrix d,e=a;
	memset(d.mx,0,sizeof(d.mx));
	for(int i=1;i<=n;i++) d.mx[i][i]=1;
	while(k){
		if(k&1) d=d*e;
		e=e*e;
		k>>=1;
	}
	return d;
}
Matrix sum(int l,int r){//二分 
	if(l==r){
		return ksm(a,l);
	}
	if((r-l+1)%2==0){
		int mid=(l+r)/2;
		Matrix x=ksm(a,mid+1-l),y=sum(l,mid);
		return y+y*x;
	}else{
		return sum(l,l)+sum(l+1,r);
	}
}
signed main(){
	cin>>n>>k>>m;
	for(int i=1;i<=n;i++) s.mx[i][i]=1;
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			cin>>a.mx[i][j];
			a.mx[i][j]%=m;
		}
	}
	s=sum(1,k);
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++) cout<<s.mx[i][j]<<" ";
		cout<<endl;
	}
	return 0;
}

注意事项

记得随时取余 \(m\)。