矩阵快速幂

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

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

1 #define MAX_D 100
2 #define mod 1000000007
3 
4 typedef long long ll;
5 typedef ll mat[MAX_D][MAX_D];
6 typedef ll vec[MAX_D];

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

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

mat_mul(mat a, mat b, int dim)

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

 1 void mat_mul(mat a, mat b, int dim){
 2     mat res;
 3     memset(res, 0, sizeof(res));
 4     for(int i=0; i<dim; i++)
 5         for(int j=0; j<dim; j++)
 6             if(a[i][j])　　
 7                 for(int k=0; k<dim; k++)
 8                     res[i][k]+=a[i][j]*b[j][k],
 9                     res[i][k]%=mod;
10     memcpy(a, res, sizeof(res));
11 }

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

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

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

1 void mat_mul(vec a, mat b, int dim){
2     vec res;
3     memset(res, 0, sizeof(res));
4     for(int i=0; i<dim; i++)
5         for(int j=0; j<dim; j++)
6             res[j]+=a[i]*b[i][j],
7             res[j]%=mod;
8     memcpy(a, res, sizeof(res));
9 }

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

 1 void mat_pow(mat a, int dim, int n){
 2     mat res;
 3     memset(res, 0, sizeof(res));
 4     for(int i=0; i<dim; i++)
 5         res[i][i]=1;
 6     mat tmp;
 7     memcpy(tmp, a, sizeof(mat));
 8     while(n){
 9         if(n&1){
10             mat_mul(res, tmp, dim);
11             n--;
12         }
13         else{
14             mat_mul(tmp, tmp, dim);
15             n>>=1;
16         }
17     }
18     memcpy(a, res, sizeof(mat));
19 }

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