矩阵乘法在很多题目中都有显著的应用,在一些递推、计算几何和其他可以用矩阵表示元素操作的题目中,可以借助快速幂的思想优化不少时间。
1.加速递推
最著名的递推当属Fibonacci数列了。
\(
f[0]=f[1]=1
f[n]=f[n-1]+f[n-2]
\)
朴素的递推需要O(n)的时间(应该没人写递归吧……
如果采用矩阵快速幂
\(
\begin{bmatrix}
f(n+1) \\
f(n)
\end{bmatrix}
=
(\begin{bmatrix}
1 & 1\\
1 & 0
\end{bmatrix})^{n}
\cdot
\begin{bmatrix}
1 \\
1
\end{bmatrix}
\)
其中矩阵乘法可以采用快速幂的思想加速运算(矩阵乘法满足结合律)
稍稍复杂的例子还有
http://wikioi.com/problem/1281/
http://acm.buaa.edu.cn/problem/653/
……等等很多
2.计算几何
example:via http://www.matrix67.com/blog/archives/276/
给定n个点,m个操作,构造O(m+n)的算法输出m个操作后各点的位置。操作有平移、缩放、翻转和旋转
这里的操作是对所有点同时进行的。其中翻转是以坐标轴为对称轴进行翻转(两种情况),旋转则以原点为中心。如果对每个点分别进行模拟,那么m个操作总共耗时O(mn)。利用矩阵乘法可以在O(m)的时间里把所有操作合并为一个矩阵,然后每个点与该矩阵相乘即可直接得出最终该点的位置,总共耗时O(m+n)。假设初始时某个点的坐标为x和y,下面5个矩阵可以分别对其进行平移、旋转、翻转和旋转操作。预先把所有m个操作所对应的矩阵全部乘起来,再乘以(x,y,1),即可一步得出最终点的位置。
3.……以及其他
example:黑书p205细菌、VOJ1049序列置换、有向图方案数……
顺次给出m个置换,反复使用这m个置换对初始序列进行操作,问k次置换后的序列。m<=10, k<2^31。
首先将这m个置换“合并”起来(算出这m个置换的乘积),然后接下来我们需要执行这个置换k/m次(取整,若有余数则剩下几步模拟即可)。注意任意一个置换都可以表示成矩阵的形式。例如,将1 2 3 4置换为3 1 2 4,相当于下面的矩阵乘法:
置换k/m次就相当于在前面乘以k/m个这样的矩阵。我们可以二分计算出该矩阵的k/m次方,再乘以初始序列即可。做出来了别忙着高兴,得意之时就是你灭亡之日,别忘了最后可能还有几个置换需要模拟。
给定一个有向图,问从A点恰好走k步(允许重复经过边)到达B点的方案数mod p的值
把给定的图转为邻接矩阵,即A(i,j)=1当且仅当存在一条边i->j。令C=A*A,那么C(i,j)=ΣA(i,k)*A(k,j),实际上就等于从点i到点j恰好经过2条边的路径数(枚举k为中转点)。类似地,C*A的第i行第j列就表示从i到j经过3条边的路径数。同理,如果要求经过k步的路径数,我们只需要二分求出A^k即可。
4.模版
模版代码来源于http://blog.csdn.net/u010153200/article/details/8974369
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 |
template <int A, int B> struct matrix { long long v[A][B]; matrix(long long r = 0) { memset(v, 0, sizeof(v)); for (int i = 0; i < A && i < B; i++) v[i][i] = r; } matrix(long long a[A][B]) { memcpy(v, a, sizeof(a)); } }; inline long long quickadd(long long a, long long b) { long long ans = 0; for (; b; b >>= 1) { if (b & 1) ans = (ans + a) % m; a = (a + a) % m; } return ans; } template <int A, int B, int C> matrix <A, C> operator * (const matrix <A, B> &a, const matrix <B, C> &b) { matrix <A, C> ans; for (int i = 0; i < A; i++) for (int j = 0; j < C; j++) { long long v = 0; for (int k = 0; k < B; k++) v += quickadd(a.v[i][k], b.v[k][j]); ans.v[i][j] = (v % m + m) % m; } return ans; } template <int A> matrix <A, A> power(const matrix<A, A> &a, long long b) { matrix <A, A> ans = 1; for (matrix<A, A> i = a; b; b >>= 1) { if (b & 1) ans = ans * i; i = i * i; } return ans; } |
[Dr.Lib]Note:Algorithms – 矩阵乘法 by Liqueur Librazy is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.