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 <cstdio>
#include <string.h>

using namespace std;
const int maxn=31;
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;
}
Matrix p;
Matrix dfs(Matrix A,int k){
    if(k==1){
        p=A;
        return A;
    }
    Matrix ret,ans;
    ret=dfs(A,k/2);
    //Matrix p=MquickPow(A,k/2);如果用快速幂计算p=A^(k/2)，则要1500ms，而直接在递归的时候同时计算p，则只要188ms。
    if(k&1){
        //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=0;i<n;i++){
        for(int j=0;j<n;j++){
            scanf("%d",&A.m[i][j]);
        }
    }
    Matrix ans;
    ans=dfs(A,k);
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            printf("%d ",ans.m[i][j]);
        }
        printf("
");
    }
    return 0;
}

View Code

方法二：

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=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("
");
    }
    return 0;
}

View Code