给出矩阵A，求矩阵(A+A^2+...+A^k)各个元素(mod yyb)的值，(nleq 30,kleq 10^9,yybleq 10^4)。
解
法一：分治
显然是数列题，故数列最浅显的减法是分治，寻找其中重复计算的部分，故可以
[A+A^2+...+A^{k/2}+A^{k/2+1}+...+A^k=
]
[A+A^2+...+A^{k/2}+A^{k/2}(A+...+A^{k/2})+(k&1)A^k=
]
[(A^{k/2}+1)(A+...+A^{k/2})+(k&1)A^k
]
对式子主体进行分治，其他的部分快速幂即可。
参考代码：
#include <iostream>
#include <cstdio>
#include <cstring>
#define il inline
#define ri register
using namespace std;
int n,yyb;
struct matrix{
int jz[30][30];
il void clear(){
memset(jz,0,sizeof(jz));
}
il void unit(){
clear();ri int i;
for(i=0;i<n;++i)jz[i][i]|=true;
}
il void read(){
ri int i,j;
for(i=0;i<n;++i)
for(j=0;j<n;++j)
scanf("%d",&jz[i][j]);
}
il void print(){
ri int i,j;
for(i=0;i<n;++i,putchar('
'))
for(j=0;j<n;++j)
printf("%d ",jz[i][j]);
putchar('
');
}
il matrix operator*(matrix x){
matrix y;y.clear();ri int i,j,k;
for(i=0;i<n;++i)
for(j=0;j<n;y.jz[i][j]%=yyb,++j)
for(k=0;k<n;++k)
y.jz[i][j]+=jz[i][k]*x.jz[k][j]%yyb;
return y;
}
il matrix operator+(matrix x){
matrix y;ri int i,j;
for(i=0;i<n;++i)
for(j=0;j<n;++j)
y.jz[i][j]=(jz[i][j]+x.jz[i][j])%yyb;
return y;
}template<class free>
il matrix operator^(free y){
matrix x(*this),ans;ans.unit();
while(y){
if(y&1)ans=ans*x;
x=x*x,y>>=1;
}return ans;
}
}state,unit,zero;
il matrix efen(int);
int main(){
int k,i;
scanf("%d%d%d",&n,&k,&yyb);
unit.unit(),state.read();
efen(k).print();
return 0;
}
il matrix efen(int y){
if(!y)return unit;
if(y==1)return state;
return ((state^(y>>1))+unit)*efen(y>>1)
+((y&1)?(state^y):zero);
}
法二：矩阵快速幂
显然数列的题目，经常会存在递推方程，于是矩阵快速幂会在其中大有用武之地，于是设(f_i=A+A^2+...+A^i),不难有递推方程(f_i=f_{i-1}+A^i)，于是我们可以同时转移(f_i,A^i),故状态矩阵为
[egin{bmatrix}A_i&f_{i-1}end{bmatrix}
]
转移矩阵为
[egin{bmatrix}A&I\0&Iend{bmatrix}
]
以此矩阵中套矩阵仿照套路即可。
参考代码：
#include <iostream>
#include <cstdio>
#include <cstring>
#define il inline
#define ri register
using namespace std;
int n,yyb;
struct matrix1{
int jz[30][30];
il void read(){
ri int i,j;
for(i=0;i<n;++i)
for(j=0;j<n;++j)
scanf("%d",&jz[i][j]);
}
il void print(){
ri int i,j;
for(i=0;i<n;++i,putchar('
'))
for(j=0;j<n;++j)
printf("%d ",jz[i][j]);
putchar('
');
}
il void clear(){
memset(jz,0,sizeof(jz));
}
il void unit(){
clear();ri int i;
for(i=0;i<n;++i)jz[i][i]|=true;
}
il matrix1 operator*(matrix1 x){
matrix1 y;y.clear();
ri int i,j,k;
for(i=0;i<n;++i)
for(j=0;j<n;y.jz[i][j]%=yyb,++j)
for(k=0;k<n;++k)
y.jz[i][j]+=jz[i][k]*x.jz[k][j]%yyb;
return y;
}
il matrix1 operator+(matrix1 x){
matrix1 y;ri int i,j;
for(i=0;i<n;++i)
for(j=0;j<n;++j)
y.jz[i][j]=(jz[i][j]+x.jz[i][j])%yyb;
return y;
}template<class free>
il matrix1 operator^(free y){
matrix1 ans,x(*this);ans.unit();
while(y){
if(y&1)ans=ans*x;
x=x*x,y>>=1;
}return ans;
}
}A;
struct matrix2{
matrix1 jz[2][2];
il void clear(){
memset(jz,0,sizeof(jz));
}
il void unit(){
clear();ri int i;
for(i=0;i<2;++i)jz[i][i].unit();
}
il matrix2 operator*(matrix2 x){
matrix2 y;y.clear();
ri int i,j,k;
for(i=0;i<2;++i)
for(j=0;j<2;++j)
for(k=0;k<2;++k)
y.jz[i][j]=y.jz[i][j]+jz[i][k]*x.jz[k][j];
return y;
}template<class free>
il matrix2 operator^(free y){
matrix2 ans,x(*this);ans.unit();
while(y){
if(y&1)ans=ans*x;
x=x*x,y>>=1;
}return ans;
}
}state,tran;
int main(){
int k;
scanf("%d%d%d",&n,&k,&yyb);
A.read(),state.jz[0][0]=A;
tran.jz[0][0]=A,tran.jz[0][1].unit();
tran.jz[1][1].unit(),state=state*(tran^k);
state.jz[0][1].print();
return 0;
}
法三：倍增前缀和优化
数列求和可以利用倍增优化，主要思想是维护(A^{2^i})，这个显然可通过(A^{2^{i+1}}=(A^{2^i})^2)来暴力递推，再维护一个倍增的数列，(f_i=A^1+A^2+...+A^{2^i})，不难得知其转移方程为(f_i=f_{i-1}A^{2^{i-1}}+f_{i-1})，而这个的维护也可以暴力维护，于是对于我们的求和，以(A+A^2+...+A^{15})为例，有
[A+A^2+...+A^{15}=f_3+A^9+...+A^{15}=
]
[f_3+A^8(A+...+A^7)=f_3+A^8(f_2+A^5+A^6+A^7)=
]
[f_3+A^8[f_2+A^4(A+A^2+A^3)]=
]
[f_3+A^8[f_2+A^4(f_1+A^2A^1)]=f_3+A^8[f_2+A^4(f_1+A^2f_0)]
]
于是按照这个先维护好一段倍增前缀和，再对后式提公因式，继续利用倍增前缀和优化，不停地继续，即可得到答案，当然你也可以递归处理，以下参考程序用的是非递归的方式。
参考代码：
#include <iostream>
#include <cstdio>
#include <cstring>
#define il inline
#define ri register
using namespace std;
int n,yyb;
struct matrix{
int jz[30][30];
il void clear(){
memset(jz,0,sizeof(jz));
}
il void unit(){
clear();ri int i;
for(i=0;i<n;++i)jz[i][i]|=true;
}
il void read(){
ri int i,j;
for(i=0;i<n;++i)
for(j=0;j<n;++j)
scanf("%d",&jz[i][j]);
}
il void print(){
ri int i,j;
for(i=0;i<n;putchar('
'),++i)
for(j=0;j<n;++j)
printf("%d ",jz[i][j]);
putchar('
');
}
il matrix operator*(matrix x){
matrix y;ri int i,j,k;y.clear();
for(i=0;i<n;++i)
for(j=0;j<n;y.jz[i][j]%=yyb,++j)
for(k=0;k<n;++k)
y.jz[i][j]+=jz[i][k]*x.jz[k][j]%yyb;
return y;
}
il matrix operator+(matrix x){
matrix y;ri int i,j;
for(i=0;i<n;++i)
for(j=0;j<n;++j)
y.jz[i][j]=(jz[i][j]+x.jz[i][j])%yyb;
return y;
}
}A,p[31],sum[31],ans,jilu;
int main(){
int k,i;
scanf("%d%d%d",&n,&k,&yyb);
A.read(),p[0]=A,jilu.unit(),sum[0]=A;
for(i=1;i<=30;++i)p[i]=p[i-1]*p[i-1];
for(i=1;i<=30;++i)sum[i]=sum[i-1]*p[i-1]+sum[i-1];
for(i=30;i>=0;--i)
if(k>>i&1){
ans=ans+sum[i]*jilu;
jilu=jilu*p[i];
}ans.print();
return 0;
}
总上所诉，不难得知数列题目的一般求法，矩阵快速幂和分治，而特殊地对于数列前缀和的问题，我们可以利用倍增前缀和优化，但无论如何，数列问题绝不只一条道路通往罗马。