传送门
题意:给定一个(n*n)的矩阵A,求(sum_{i=1}^{k}A^i),运算值对m取模.((n<=30,k<=10^9))
分析:如果题目是求(A^k),就是矩阵快速幂的模板题了.但现在我们要求的是(sum_{i=1}^{k}A^i),如果直接求时间复杂度是(O(kn^3)),T飞.因为本题最基础的肯定是矩阵乘法,n又这么小,所以时间复杂度(n^3)是不可避免的,所以我们要从k下手,因为(sqrt{k})都会超时,所以肯定是要想方设法优化到(logk).
对于(sum_{i=1}^{k}A^i)这种式子,它一定是可以由前面的(k比较小的)推出后面的(k比较大的),例如,假设我们已知原式(A^1+A^2+A^3),要求(A^4+A^5+A^6),则可以让原式乘上(A^3),不难发现如果再加上一个原式,就得到(A^1+A^2+A^3+A^4+A^5+A^6),这不就是我们想要的吗?
为了方便,我们设(B_x=A^1+A^2+...+A^x),即(B_3=A^1+A^2+A^3),则(B_6=B_3*A^3+B_3),我们要求的答案就是(B_k).
而且我们现在要把k优化到(logk),就要想到可以二进制优化.我们预处理求出(A^{2^0},A^{2^1},A^{2^2},A^{2^3},...,A^{2^{31}}),这同时也可以得到(B_{2^0},B_{2^1},B_{2^2},B_{2^3},...,B_{2^{31}}).当然了,数组不能开(2^{31})这么大,所以我们直接用A[i]表示(A^{2^i}),B[i]表示(B^{2^i}).
我们预处理完A,B数组后,考虑我们的答案(B_k),首先k一定能分解成几个2的整数次幂相加的形式,例如(k=7=2^0+2^1+2^2),这时我们想,要是(B_7=B_{1+2+4}=B_1+B_2+B_4)那该多好啊!
虽然这是不合法的,但它的确为我们提供了一种思路.结合上面那个(B_6=B_3*A^3+B_3)的例子,你能否推出(B_{x+y}=B_x+A^x*B_y)这个式子呢?(其中x和y都是2的整数次幂)
int n,m,K;
struct mat{
int jz[35][35];
void in(){
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
jz[i][j]=read()%m;
}//读入
void out(){
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++)
printf("%d ",jz[i][j]%m);
printf("
");
}
}//输出
//直接重载运算符'*'为矩阵乘法
mat operator *(mat x){
mat y;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
y.jz[i][j]=0;//记得初始化为0
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
for(int k=1;k<=n;k++)
y.jz[i][j]=y.jz[i][j]%m+(jz[i][k]*x.jz[k][j])%m;
return y;
}
//直接重载运算符'*'为矩阵加法
mat operator +(mat x){
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
x.jz[i][j]=(jz[i][j]+x.jz[i][j])%m;
return x;
}
}ans,a[32],b[32];
int main(){
n=read();K=read();m=read();
a[1].in();
//上面的分析中,说a[i]表示a[1<<i],但为了能够从
//下标1开始操作,代码中a[i]表示a[1<<(i-1)],b同理
//所以a[1]就是A^1
b[1]=a[1];
for(int i=1;i<=30;i++){
a[i+1]=a[i]*a[i];//我们重载了乘号和加号
b[i+1]=b[i]*a[i]+b[i];
}//预处理a,b数组
for(int i=1;i<=31;i++){//二进制分解K
if(K&1)ans=ans*a[i]+b[i];
//K&1表示K分解成二进制数后的第i位是1,
//这里的ans可以理解为b[ans]
K>>=1;
}
ans.out();
return 0;
}