POJ 3233 Matrix Power Series （矩阵快速幂+二分求解）

题意：求S=(A+A^2+A^3+...+A^k)%m的和

方法一：二分求解
S=A+A^2+...+A^k
若k为奇数:
S=(A+A^2+...+A^(k/2))+A^(k/2)*(A+A^2+...+A^(k/2))+A^k
若k为偶数：
S=(A+A^2+...+A^(k/2))+A^(k/2)*(A+A^2+...+A^(k/2))

也可以这么二分(其实和前面的差不多)：
S(2n)=A+A^2+...+A^2n=(1+A^n)*(A+A^2+...+A^n)=(1+A^n)*S(n)
S(2n+1)=A+A^2+...+A^(2n+1)=A(1+A+..+A^2n)=A+(A+A^(n+1))*S(n)

一开始1900+ms，优化了下1500ms...还是太慢了。。。
本来在递归的时候，用快速幂计算A^(k/2)
后来改用直接递归的同时，计算A^(k/2)，直接变成200ms左右。。。瞬间提升了10倍。。。

#include <iostream>
#include <iostream>
#include <cstdio>
#include <string.h>

using namespace std;
const int maxn=65;
int mod;
int n,k,m;
struct Matrix{
int m[maxn][maxn];
void init(){
memset(m,0,sizeof(m));
}
void initE(){
memset(m,0,sizeof(m));
for(int i=0;i<n;i++)
m[i][i]=1;
}
}A;
//重载+运算符
Matrix operator+(Matrix a,Matrix b){
Matrix c;
for(int i=0;i<n;i++){
for(int j=0;j<n;j++)
c.m[i][j]=(a.m[i][j]+b.m[i][j])%mod;
}
return c;
}
//重载*运算符
Matrix operator*(Matrix a,Matrix b){
Matrix c;
for(int i=0;i<n;i++){
for(int j=0;j<n;j++){
c.m[i][j]=0;
for(int k=0;k<n;k++){
c.m[i][j]=(c.m[i][j]+a.m[i][k]*b.m[k][j]%mod)%mod;
}
}
}
return c;
}
//矩阵快速幂
Matrix MquickPow(Matrix A,int b){
Matrix ret;
ret.initE();
while(b){
if(b&1)
ret=ret*A;
A=A*A;
b=b>>1;
}
return ret;
}
int main()
{
Matrix B;
B.init();
scanf("%d%d%d",&n,&k,&m);
mod=m;
for(int i=0;i<n;i++){
for(int j=0;j<n;j++){
scanf("%d",&B.m[i][j]);
}
B.m[i+n][i]=1;
B.m[i+n][i+n]=1;
}
n=n*2;
Matrix ans=MquickPow(B,k+1);
for(int i=n/2;i<n;i++){
for(int j=0;j<n/2;j++){
if(i==j+n/2)
ans.m[i][j]=(ans.m[i][j]-1+mod)%mod;   //对角线还要减去单位矩阵的1
printf("%d ",ans.m[i][j]);
}
printf("\n");
}
return 0;
}

View Code