矩阵的快速幂

原博主：http://www.cnblogs.com/vongang/archive/2012/04/01/2429015.html

据说，矩阵快速幂在递推式优化上相当神奇，而且效率很高。。。

　　两矩阵相乘，朴素算法的复杂度是O(N^3)。如果求一次矩阵的M次幂，按朴素的写法就是O(N^3*M)。既然是求幂，不免想到快速幂取模的算法，这里有快速幂取模的介绍，a^b %m 的复杂度可以降到O(logb)。如果矩阵相乘是不是也可以实现O(N^3 * logM)的时间复杂度呢？答案是肯定的。

　　先定义矩阵数据结构：　　

struct Mat {
    double mat[N][N];
};

　　O(N^3)实现一次矩阵乘法

Mat operator * (Mat a, Mat b) {
    Mat c;
    memset(c.mat, 0, sizeof(c.mat));
    int i, j, k;
    for(k = 0; k < n; ++k) {
        for(i = 0; i < n; ++i) {
            if(a.mat[i][k] <= 0)  continue;   //(针对ZOJ2853)剪枝，cpu运算乘法的效率并不是想像的那么理想（加法的运算效率高于乘法，比如Strassen矩阵乘法）
            for(j = 0; j < n; ++j) {
                if(b.mat[k][j] <= 0)    continue;    //剪枝
                c.mat[i][j] += a.mat[i][k] * b.mat[k][j];
            }
        }
    }
    return c;
}

自己写了一下。
题目链接： http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1113

代码:

#include<iostream>
#include<cstdio>
#include<cstring>
#define Mod 1000000007
typedef long long int ll;
using namespace std;
ll n,m;
struct Martix
{
    ll a[110][110];
};
Martix mul(Martix x,Martix y)
{
    Martix ans;
    memset(ans.a,0,sizeof(ans.a));
    ll i,j,k;
    for(i=1;i<=n;i++)
    {
        for(j=1;j<=n;j++)
        {
            for(k=1;k<=n;k++)
            {
                ans.a[i][j]+=((x.a[i][k]%Mod)*(y.a[k][j]%Mod))%Mod;
            }
        }
    }
    return ans;
}
Martix quickpow(Martix m,ll k)
{
    ll i,j;
    Martix fuck;
    for(i=1;i<=n;i++)
    {
        for(j=1;j<=n;j++)
        {
            if(i==j)
                fuck.a[i][j]=1;
            else
                fuck.a[i][j]=0;
        }
    }
    while(k>0)
    {
        if(k&1)
            fuck=mul(fuck,m);
        m=mul(m,m);
        k=k>>1;
    }
    return fuck;
}
int main()
{
    int i,j;
    Martix start,ans;
    scanf("%I64d%I64d",&n,&m);
    for(i=1;i<=n;i++)
    {
        for(j=1;j<=n;j++)
        {
            scanf("%I64d",&start.a[i][j]);
        }
    }
    ans =quickpow(start,m);
    for(i=1;i<=n;i++)
    {
        for(j=1;j<=n;j++)
        {
            if(j==n)
                printf("%I64d\n",ans.a[i][j]%Mod);
            else
                printf("%I64d ",ans.a[i][j]%Mod);
        }
    }
    return 0;
}