矩阵快速幂可认为是快速幂算法的扩展。
有一种借助于 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行移到参数里。