矩阵快速幂优化的动态规划

因为最近写矩阵快速幂总是搞反相乘的顺序所以来写一发博客

不过突然写这么简单的东西似乎会被鄙视？

矩阵和矩阵乘法

（没学过矩乘的同学最好不要通过此博客学习，它的目的不在于此）

矩阵是一个由数组成的矩形，特殊地，行数为1的矩阵被称作行向量，列数为1的矩阵被称作列向量。

矩阵A（m*n）和B（n*p）可以进行乘法得到矩阵C（m*p），也就是说当前者的列数等于后者行数的时候。乘法的定义是

cij=∑mk=1Aik∗Bkj

按照定义，可以写出矩阵和矩乘的代码，其中矩乘的效率是O(n*m*p)

struct matrix{int m[maxm][maxm];int x,y;inline int *operator [](int x){return m[x];}matrix(int x,int y):x(x),y(y){memset(m,0,sizeof(m));}matrix operator *(matrix &b){matrix ans(x,b.y);for(int i=0;i<x;i++)for(int j=0;j<b.y;j++)for(int k=0;k<y;k++)ans[i][j]+=m[i][k]*b[k][j];//一般的题目都会在此处要求取模，请注意相乘时int溢出return ans;}};

矩阵乘法是满足结合律的，也就是说，(A*B)*C=A*(B*C)

矩阵乘法与DP

对于满足DP[i][j]需要从若干个DP[i-1][k]处取值乘上常数求和的，且每一轮从i-1到i转移方式固定的动态规划（如果上面的解释很难懂可以看公式）

DP[i][j]=∑mk=1DP[i−1][k]∗A[j][k]

可以看做是每次把初始的DP[i-1]向量乘上一个行数和列数相等的转移矩阵A，变换成一个长度相等的DP[i]向量。若使得DP数组对应向量是行向量，那么A[j][k]表示从DP[i-1][j]转移到DP[i][k]时需要乘的数，需要将行向量乘以转移矩阵（顺序很重要！不满足交换律）。若DP数组对应向量是列向量，那么A[j][k]表示从DP[i-1][k]转移到DP[i][j]次需要乘以的数，需要将转移矩阵乘以向量

如果转移次数为n，数组长度为m，那么普通DP的复杂度为O(n*m^2)。当n很大（1e18），m很小（<=100）的时候，这样的转移显然是耗时的。这时我们便可以用矩阵快速幂优化。

矩阵快速幂及其优化

前面提到矩乘满足结合律，那么设DP[0]对应的向量是B，转移矩阵为A，B*A*A*...*A=B*(A*A*...*A)=B*(A^n)

其中A^n可以在O(m^3*log n)时间内完成，其中m^3是矩乘的时间

代码如下

matrix operator ^(matrix b,long t){matrix ans;//此处ans应初始化为单位矩阵while(t){if(t&1) ans=ans*b;b=b*b;t>>=1;}return ans;}

但是上文中的乘法，传了一整个matrix结构体作为参数，又返回了一个结构体，这使得赋值占用了大量的时间，可能会被卡常数。那么能不能通过传指针的方法避免这部分常数呢？答案是肯定的

void quickmul(matrix *a,matrix *b,matrix *ans){memset(ans->m,0,sizeof(ans->m));for(int i=0;i<a->x;i++)for(int j=0;j<b->y;j++)for(int k=0;k<a->y;k++)ans->m[i][j]+=a->m[i][k]*b->m[k][j];}matrix operator ^(matrix in,int t){matrix tmpbase(in.x,in.y);matrix tmp1(in.x),tmp2(in.x);matrix *ans=&tmp1,*swapans=&tmp2;matrix *b=&in,*swapb=&tmpbase;while(t){if(t&1){quickmul(ans,b,swapans);swap(ans,swapans);//交换指针地址，不改变其中的值}quickmul(b,b,swapb);swap(b,swapb);t>>=1;}return *ans;}