矩阵快速幂

http://www.cnblogs.com/vongang/archive/2012/04/01/2429015.html
两矩阵相乘，朴素算法的复杂度是O(N^3)。如果求一次矩阵的M次幂，按朴素的写法就是O(N^3*M)。既然是求幂，不免想到快速幂取模的算法，这里有快速幂取模的介绍，a^b %m 的复杂度可以降到O(logb)。
先定义矩阵数据结构：　　

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

任何矩阵乘以单位矩阵，其值不改变。

Mat operator ^ (Mat a, int k) {
    Mat c;
    int i, j;
    for(i = 0; i < n; ++i)
        for(j = 0; j < n; ++j)
            c.mat[i][j] = (i == j);    //初始化为单位矩阵

    for(; k; k >>= 1) {
        if(k&1) c = c*a;

        a = a*a;
    }
    return c;
}

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;   //不要小看这里的剪枝，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;
}