Matrix Power Series
给出矩阵A，求矩阵\(A+A^2+...+A^k\)各个元素\(mod\ yyb\)的值，\(n\leq 30,k\leq 10^9,yyb\leq 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('\n'))
for(j=0;j<n;++j)
printf("%d ",jz[i][j]);
putchar('\n');
}
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\),故状态矩阵为
\[\begin{bmatrix}A_i&f_{i-1}\end{bmatrix}
\]
转移矩阵为
\[\begin{bmatrix}A&I\\0&I\end{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('\n'))
for(j=0;j<n;++j)
printf("%d ",jz[i][j]);
putchar('\n');
}
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('\n'),++i)
for(j=0;j<n;++j)
printf("%d ",jz[i][j]);
putchar('\n');
}
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;
}
总上所诉，不难得知数列题目的一般求法，矩阵快速幂和分治，而特殊地对于数列前缀和的问题，我们可以利用倍增前缀和优化，但无论如何，数列问题绝不只一条道路通往罗马。
