矩阵快速幂

矩阵快速幂可认为是快速幂算法的扩展。

有一种借助于 typedef 的较为方便的写法。

#define MAX_D 100
#define mod 1000000007

typedef long long ll;
typedef ll mat[MAX_D][MAX_D];
typedef ll vec[MAX_D];

采用这种方式定义矩阵，既简洁又可以很方便地调用 memset, memcpy 等函数。

矩阵（方阵）乘法，推荐写成参数较少的形式

mat_mul(mat a, mat b, int dim)

其中 a 相当于目的操作数 (destination)，b 相当于源操作数 (source)，结果 (result) 存在 a 中。

void mat_mul(mat a, mat b, int dim){
    mat res;
    memset(res, 0, sizeof(res));
    for(int i=0; i<dim; i++)
        for(int j=0; j<dim; j++)
            if(a[i][j])　　
                for(int k=0; k<dim; k++)
                    res[i][k]+=a[i][j]*b[j][k],
                    res[i][k]%=mod;
    memcpy(a, res, sizeof(res));
}

这种写法针对稀疏矩阵进行了优化。

还可以重载 (overload)  一个行矩阵（向量）乘方阵的版本。

注意：数学上一般将向量视做列矩阵，但在某些涉及状态转移矩阵的计数问题中，将“向量” （其实这里是否称作"向量"也无关紧要）视作行矩阵更方便。

void mat_mul(vec a, mat b, int dim){
    vec res;
    memset(res, 0, sizeof(res));
    for(int i=0; i<dim; i++)
        for(int j=0; j<dim; j++)
            res[j]+=a[i]*b[i][j],
            res[j]%=mod;
    memcpy(a, res, sizeof(res));
}

快速幂采用的是非递归的写法

void mat_pow(mat a, int dim, int n){
    mat res;
    memset(res, 0, sizeof(res));
    for(int i=0; i<dim; i++)
        res[i][i]=1;
    mat tmp;
    memcpy(tmp, a, sizeof(mat));
    while(n){
        if(n&1){
            mat_mul(res, tmp, dim);
            n--;
        }
        else{
            mat_mul(tmp, tmp, dim);
            n>>=1;
        }
    }
    memcpy(a, res, sizeof(mat));
}

这里仍采用参数较少的写法，如需保留矩阵 a （a 往往是状态转移矩阵），就把第2行移到参数里。