题意：求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 <cstdio>
#include <string.h> using namespace std;
const int maxn=;
int mod;
int n,k,m;
struct Matrix{
int m[maxn][maxn];
void init(){
memset(m,,sizeof(m));
}
void initE(){
memset(m,,sizeof(m));
for(int i=;i<n;i++)
m[i][i]=;
}
}A;
//重载+运算符
Matrix operator+(Matrix a,Matrix b){
Matrix c;
for(int i=;i<n;i++){
for(int j=;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=;i<n;i++){
for(int j=;j<n;j++){
c.m[i][j]=;
for(int k=;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&)
ret=ret*A;
A=A*A;
b=b>>;
}
return ret;
}
Matrix p;
Matrix dfs(Matrix A,int k){
if(k==){
p=A;
return A;
}
Matrix ret,ans;
ret=dfs(A,k/);
//Matrix p=MquickPow(A,k/2);如果用快速幂计算p=A^(k/2)，则要1500ms，而直接在递归的时候同时计算p，则只要188ms。
if(k&){
//return ret+ret*p+p*p*A;
ans=ret+ret*p+p*p*A;
p=p*p*A;
}
else{
//return ret+ret*p;
ans=ret+ret*p;
p=p*p;
}
return ans;
}
int main()
{
scanf("%d%d%d",&n,&k,&m);
mod=m;
for(int i=;i<n;i++){
for(int j=;j<n;j++){
scanf("%d",&A.m[i][j]);
}
}
Matrix ans;
ans=dfs(A,k);
for(int i=;i<n;i++){
for(int j=;j<n;j++){
printf("%d ",ans.m[i][j]);
}
printf("\n");
}
return ;
}
方法二：
688ms
http://blog.sina.com.cn/s/blog_9d5278a301015mbd.html
http://blog.csdn.net/wangjian8006/article/details/7868864
S=A(E+A(E+A(E+...A(E+A))))
这样就可以想，不妨构造一个矩阵T使得T*T，这样乘下去每次可以得到A*(A+E)+E
A 0 A^2 0 A^3 0
B= B^2= B^3=
E E A+E E A^2+A+E E
不难得出： A^(k+1) 0
B^(k+1)=
S(k) E
#include <iostream>
#include <iostream>
#include <cstdio>
#include <string.h> using namespace std;
const int maxn=;
int mod;
int n,k,m;
struct Matrix{
int m[maxn][maxn];
void init(){
memset(m,,sizeof(m));
}
void initE(){
memset(m,,sizeof(m));
for(int i=;i<n;i++)
m[i][i]=;
}
}A;
//重载+运算符
Matrix operator+(Matrix a,Matrix b){
Matrix c;
for(int i=;i<n;i++){
for(int j=;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=;i<n;i++){
for(int j=;j<n;j++){
c.m[i][j]=;
for(int k=;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&)
ret=ret*A;
A=A*A;
b=b>>;
}
return ret;
}
int main()
{
Matrix B;
B.init();
scanf("%d%d%d",&n,&k,&m);
mod=m;
for(int i=;i<n;i++){
for(int j=;j<n;j++){
scanf("%d",&B.m[i][j]);
}
B.m[i+n][i]=;
B.m[i+n][i+n]=;
}
n=n*;
Matrix ans=MquickPow(B,k+);
for(int i=n/;i<n;i++){
for(int j=;j<n/;j++){
if(i==j+n/)
ans.m[i][j]=(ans.m[i][j]-+mod)%mod; //对角线还要减去单位矩阵的1
printf("%d ",ans.m[i][j]);
}
printf("\n");
}
return ;
}