题解：P10502 Matrix Power Series

解法

提供一种较新颖的做法，首先由题会先想到一个经典的问题，求斐波那契数列的第 \(n\) 项前缀和，即 \(\sum_{i=1}^{n}f_{i}\)。这道题的解法即为维护总和与目前的 \(f_{i}\)，构造出转移矩阵后矩阵快速幂解决即可。

考虑 \(A_{i}\) 表示 \(A^i\)，那么此题就是求 \(\sum_{i=1}^{n}A_{i}\)，与上段介绍的问题极其相似，于是考虑上问的类似解法。

但是我们想知道的答案是一个矩阵怎么办呢？创新点！由于矩阵加满足交换律，并且满足分配律，所以就可以推出矩阵套矩阵满足结合律，那我们就可以矩阵套矩阵快速幂。

维护由矩阵构成的答案，然后设计一个由矩阵构成的状态转移矩阵！

状态矩阵套矩阵为：

\[\huge \begin{bmatrix} S & A_{i} \end{bmatrix} \]

其中，\(S\) 表示目前的前缀和，\(A_{i}\) 表示 \(A^{i}\)。

设全零矩阵为 \(E\)，单位矩阵为 \(R\)，\(A\) 表示给定矩阵，那么显然转移矩阵如下：

\[\huge \begin{bmatrix} R & E\\ A & A\\ \end{bmatrix} \]

矩阵快速幂转移即可，时间复杂度 \(O(S^3n^3 \log k)\)，其中 \(S\) 为转移矩阵套矩阵行数与列数，在刚才的推导中 \(S=2\)。

Code

#include <bits/stdc++.h>
using namespace std;
const int N=35;
struct matrix {
	int n,m,a[N][N];
	matrix(){memset(a,0,sizeof a);}
}E,R,A;
struct Smatrix {
	int n,m;
	matrix a[4][4];
};
int p;
matrix operator*(matrix a,matrix b) {
	if(a.n>b.n) swap(a,b);
	matrix c;c.n=a.n,c.m=a.m;
	for(int i=1;i<=a.n;i++)
		for(int j=1;j<=a.m;j++)
			for(int k=1;k<=a.m;k++)
				c.a[i][j]=(c.a[i][j]+a.a[i][k]*b.a[k][j]%p)%p;
	return c;
}
matrix operator+(matrix a,matrix b) {
	matrix c;c.n=max(a.n,b.n),c.m=max(a.m,b.m);
	for(int i=1;i<=max(a.n,b.n);i++)
		for(int j=1;j<=max(a.m,b.m);j++)
			c.a[i][j]=(a.a[i][j]+b.a[i][j])%p;
	return c;
}
Smatrix operator*(Smatrix a,Smatrix b) {
	if(a.n>b.n) swap(a,b);
	Smatrix c;c.n=a.n,c.m=a.m;
	for(int i=1;i<=a.n;i++)
		for(int j=1;j<=a.m;j++)
			c.a[i][j].n=a.a[i][j].n,c.a[i][j].m=a.a[i][j].m;
	for(int i=1;i<=a.n;i++)
		for(int j=1;j<=a.m;j++)
			for(int k=1;k<=a.m;k++)
				c.a[i][j]=c.a[i][j]+a.a[i][k]*b.a[k][j];
	return c;
}
Smatrix qmi(Smatrix a,int b) {
	b--;
	Smatrix res;res.n=a.n,res.m=a.m;
	for(int i=1;i<=a.n;i++)
		for(int j=1;j<=a.m;j++) res.a[i][j]=a.a[i][j];
	while(b) {
		if(b&1) res=res*a;
		a=a*a;
		b>>=1;
	}
	return res;
}
int n,k;
signed main() {
	cin>>n>>k>>p;
	A.n=n,A.m=n;R.n=R.m=n;E.n=E.m=n;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++) cin>>A.a[i][j],R.a[i][j]=(i==j);
	Smatrix O,G;
	O.n=1,O.m=2,G.n=2,G.m=2;
	O.a[1][1]=A,O.a[1][2]=A;
	G.a[1][1]=R,G.a[2][1]=A;
	G.a[1][2]=E,G.a[2][2]=A;
	if(k-1>0) {
		G=qmi(G,k-1);
		O=O*G;
	}
	for(int i=1;i<=n;i++,puts(""))
		for(int j=1;j<=n;j++)
			cout<<O.a[1][1].a[i][j]<<' ';
	return 0;
}