矩阵（matrix）应用大总结（一）WOJ 642 Lost In WHU + POJ 3233

        矩阵、矩阵乘法是线性代数中最基本的东西，然而却在ACM竞赛以及各个领域发挥着不可替代的作用。

        好了不装逼了。矩阵可以说在很早的时候就被动的遇到过了，但是一直没有做一个系统的总结。只是简单的知道有矩阵快速幂，然后线性递推可以利用矩阵乘法加速。但是这些的具体原理也不是非常的理解。直到最近几个月，实在遇到了很多用矩阵的题目，不得不系统的学习一下。按照惯例，先膜一下鼻祖matrix67——>《十个利用矩阵乘法解决的经典题目》。果然是人如其id，对矩阵就是很了解orz……

        另外，我在写武汉大学校赛现场赛的一题的题解的时候也说了，要补线性递推的构造方法，现在就来填坑了。

        这篇绝对是干货！

        从最基础的开始（矩阵乘法和快速幂我就不说了），线性递推。这个的话，最原始的问题就是FIbonacci数列了，题面很简单，求Fibonacci数列的第n项mod一个数k。n可以很大。暴力做的话不用我说，直接tle，所以得用矩阵优化。第一次听到用矩阵优化递推，我（当时才是高中生）整个人都斯巴达了。记得应该是学习了这个文章，后来才发现作者是著名的叉姐orz……学了之后几年没有碰然后又全忘了……

        主要思想就是，对于形如f(n)=a1*f(n-1)+a2*f(n-2)+……+ak-1*f(n-k+1)+ak*f(n-k)的k阶线性递推式，构造矩阵，利用矩阵乘法的结合律用快速幂求数列第n项。具体构造方法的话，设构造矩阵为x（在递推式没有常数项的时候，矩阵大小为k*k），首先在x[i+1][i]上都放上1，然后再在x[1][1..k]上放对应的系数a1,a2,...,ak。对应的还会有一个向量，储存初始值，自上而下依次为f(k)到f(1)。如果要得到第n（n>k）项的结果，只需要让x的n-k次方乘初始向量，然后向量的最上面的数字即为f(n)。很自然的会问，这样的原理是什么，下面就来解释。首先，在对角线往下一条斜线全部赋值1，是在维护最后向量的f(n-1)到f(n-k)项，在矩阵上体现就是一个初等变化，使得向量的坐标整体下移，原来的f(n)，成了一次之后就到了f(n-1)的位置。一方面维护最后的结果向量，另一方面与系数对应上。其次，第一行的系数则是巧妙的对上了，矩阵乘法恰好对应每一项。如此一来巧妙配合上快速幂巧妙的加速。对于k阶的递推，复杂度为O(k^3logN)。

        这个比较简单这里就不给例子了。

        然后就是Lost In WHU了。我在之前写这道题的文章里也讲的比较清楚了，邻接矩阵作为一种特殊的矩阵，其N次方对应位置的数值，代表经过N步从i到达j的路径条数。但是这题给出了另一个问题，要求的是小于等与次数T的所有方案数之和。我们设邻接矩阵为A，相当于求A+A^2+A^3+...+A^N。虽然说我们有快速幂，但是当N很大的时候仍然会tle。解决这个问题的其中一种取巧的方法我之前已经说过了，这里不在赘述。我来说说另外两种相对“暴力”的方法。我们观察到，原式=f(N)=(A+A^2+...+A^N/2)+(A+A^2+...+A^N/2)*A^N/2(N为偶数)，当N为奇数的时候，在后面加一项即可。利用这个性质我们又可以二分了，和快速幂是一个道理，对半的去求解。具体代码如下：

#include<bits/stdc++.h>#define mod 1000000007#define N 110using namespace std;int n,m,k;struct matrix{long long a[N][N];void init(){for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)a[i][j]=0;}} x;matrix matmulti(matrix x,matrix y){matrix ans; ans.init();for(int i=1;i<=n;i++)for(int j=1;j<=n;j++){for(int k=1;k<=n;k++)        ans.a[i][j]+=(x.a[i][k]*y.a[k][j])%mod;    ans.a[i][j]%=mod;}    return ans;}matrix matadd(matrix x,matrix y){for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)    x.a[i][j]=(x.a[i][j]+y.a[i][j])%mod;return x;}matrix matpow(matrix x,int y){ while ((y&1)==0){y>>=1;x=matmulti(x,x);}matrix ans=x; y>>=1;while (y!=0){x=matmulti(x,x);if ((y&1)!=0) ans=matmulti(ans,x);y>>=1;}return ans;}matrix matmultisum(matrix x,int y){if (y==1) return x;matrix ans;ans=matmultisum(x,y>>1);matrix z=matpow(x,y>>1);ans=matadd(ans,matmulti(ans,z));if ((y&1)!=0) ans=matadd(ans,matpow(x,y));return ans;}int main(){scanf("%d%d",&n,&m);for(int i=1;i<=m;i++){int u,v; scanf("%d%d",&u,&v);if (v==n) x.a[u][v]=1;else if (u==n) x.a[v][u]=1;else x.a[u][v]=x.a[v][u]=1;  }scanf("%d",&k);matrix ans=matmultisum(x,k);printf("%lld\n",ans.a[1][n]);return 0;}

        这个还是很好理解的，我就不多解释了，当时就觉得非常的厉害orz……

        但是故事并没有结束，对于刚刚那个优化的式子，我们注意到可以写成这样的形式f(n)=(1+f(n/2))*f(n/2)，这又是一个递推，只不过问题出现了，这个并不是线性的，但不用怕，用前人的经验，我们也能发现，同样也可以构造矩阵解决它。注意到对于矩阵{{k,1}{0,1}}k为常数，它n次方后的结果是{{k^n,1+k+k^2+...+k^n-1}{0,1}}，注意到它的第一行第二列，正好是我们想要的形式。当k从常数变成一个矩阵A之后，相应的1也应该变成等尺寸的单位矩阵（线代中也经常这么推广）。如此对矩阵求k+1次方，它对应位置的矩阵减去一个等尺寸的单位矩阵即可。这里对于非线性的递推的矩阵应该也是有规律的，但是这次就先不总结了，还是要靠自己体会吧~具体代码如下：

#include<iostream>#include<cstdio>#include<cstdlib>#include<cstdio>#include<cstring>#include<algorithm>#include<cmath>#define LL long long#define mod 1000000007#define N 210//千万注意用这种方法的话矩阵大小要增加一倍using namespace std;int n,m,k;struct matrix{LL a[N][N];void init(){for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)a[i][j]=0;}} x;matrix matmulti(matrix x,matrix y){matrix ans; ans.init();for(int i=1;i<=n;i++)for(int j=1;j<=n;j++){for(int k=1;k<=n;k++)        ans.a[i][j]+=(x.a[i][k]*y.a[k][j])%mod;    ans.a[i][j]%=mod;}    return ans;}matrix matadd(matrix x,matrix y){for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)    x.a[i][j]=(x.a[i][j]+y.a[i][j])%mod;return x;}matrix matpow(matrix x,int y){ while ((y&1)==0){y>>=1;x=matmulti(x,x);}matrix ans=x; y>>=1;while (y!=0){x=matmulti(x,x);if ((y&1)!=0) ans=matmulti(ans,x);y>>=1;}return ans;}int main(){scanf("%d%d",&n,&m);for(int i=1;i<=m;i++){int u,v; scanf("%d%d",&u,&v);if (v==n) x.a[u][v]=1;else if (u==n) x.a[v][u]=1;else x.a[u][v]=x.a[v][u]=1;  }scanf("%d",&k);for(int i=n+1;i<=n+n;i++)x.a[i-n][i]=x.a[i][i]=1;//构造两个单位矩阵n*=2;//对应尺寸也要扩大matrix ans=matpow(x,k+1);printf("%lld\n",ans.a[1][n]);return 0;}

        所以说，我用三种不同的姿势把Lost In WHU给A了，现在我们来看看具体效果吧：

        自上而下依次是标准解法、本文第二种解法和本文第一种解法。可以看出本文第一种解法速度较慢，原因的话是因为那个二分的过程是用递归的，而且在递归过程中每次都新生成了矩阵，导致内存也较大（这个可以避免）；本文第二种解法缺点之处则在于把矩阵扩大了，根据矩阵乘法O(k^3*logN)的复杂度，比没有扩大的时候是会慢挺多的。所以我才说这两种方法相对暴力，但是却体现了矩阵本身的优美。

        然后POJ 3233的那题实际上比Lost In WHU还要裸我就不多说了。

        最后，矩阵还可以用于一些图形上的变化，最著名的莫过于转轴公式了。当然除了这个，什么平移旋转放缩它都是可以的。再留一点尾巴，对于矩阵乘法目前的复杂度是O(N^3)，但是实际上在某些情况下貌似可以优化到O(N^2)。真是快上加快啊，这个的话就下次再说了。

        最后的最后再来小结一下吧，通过这些介绍，我们发现矩阵可以用在很多地方。线性递推、甚至是所有的递推，然后不止于递推，优化dp也是可以的。另一方面，运用与图论的邻接矩阵，解决规定步骤数的方案数问题。然后对于一些特殊的求和，也可以构造矩阵求解。关于构造的话，有很深的学问，彻底理解透了就不只是简单的线性那么简单，你可以在初始向量里放很多你想要的东西，通过一个矩阵解决多个问题也不是梦想，放各种状态也是种不错的选择。这一方面的话，我们人也还要再好好研习，但是矩阵最常用的部分几乎都已经在这篇文章中呈现了出来。

        好了，我决定停手睡觉了~