首页 > 代码库 > 【矩阵加速】 矩阵 快速幂
【矩阵加速】 矩阵 快速幂
矩阵的快速幂是用来高效地计算矩阵的高次方的。将朴素的o(n)的时间复杂度,降到log(n)。
这里先对原理(主要运用了矩阵乘法的结合律)做下简单形象的介绍:
一般一个矩阵的n次方,我们会通过连乘n-1次来得到它的n次幂。
但做下简单的改进就能减少连乘的次数,方法如下:
把n个矩阵进行两两分组,比如:A*A*A*A*A*A => (A*A)*(A*A)*(A*A)
这样变的好处是,你只需要计算一次A*A,然后将结果(A*A)连乘自己两次就能得到A^6,即(A*A)^3=A^6。算一下发现这次一共乘了3次,少于原来的5次。
其实大家还可以取A^3作为一个基本单位。原理都一样:利用矩阵乘法的结合律,来减少重复计算的次数。
回头看看矩阵的快速幂问题,我们是不是也能把它离散化呢?比如A^19 => (A^16)*(A^2)*(A^1),显然采取这样的方式计算时因子数将是log(n)级别的(原来的因子数是n),不仅这样,因子间也是存在某种联系的,比如A^4能通过(A^2)*(A^2)得到,A^8又能通过(A^4)*(A^4)得到,这点也充分利用了现有的结果作为有利条件。下面举个例子进行说明:
现在要求A^156,而156(10)=10011100(2)
也就有A^156=>(A^4)*(A^8)*(A^16)*(A^128) 考虑到因子间的联系,我们从二进制10011100中的最右端开始计算到最左端。细节就说到这,下面给核心代码:
1 while(N)
2 {
3 if(N&1)
4 res=res*A;
5 n>>=1;
6 A=A*A;
7 }
里面的乘号,是矩阵乘的运算,res是结果矩阵。
第3行代码每进行一次,二进制数就少了最后面的一个1。二进制数有多少个1就第3行代码就执行多少次。
poj 3070 矩阵加速求斐波那契数列代码:
#include <cstdio>#include <iostream>#include <algorithm>#include <memory.h>using namespace std;int n = 2, m;int s[2][2]={{1, 1}, {1, 0}};struct matrix{ int r, c; int mat[3][3];} res, origin;void init(){ memset(res.mat, 0, sizeof(res.mat)); res.r = n; res.c = n; for(int i = 1; i <= n; i++) { res.mat[i][i] = 1; } origin.c = n; origin.r = n; for(int i = 1; i <= n; i++) { for(int j = 1; j <= n; j++) origin.mat[i][j] = s[i - 1][j - 1]; }}matrix multi(matrix x, matrix y){ matrix t; int i, j, k; memset(t.mat, 0, sizeof(t.mat)); t.r = x.r; t.c = y.c; for(i = 1; i <= x.r; i++) { for(k = 1; k <= x.c; k++) if(x.mat[i][k]) { for(j = 1; j <= y.c; j++) { t.mat[i][j] += (x.mat[i][k] * y.mat[k][j]) % 10000; t.mat[i][j] %= 10000; } } }// for(i=1;i<=2;i++)// {// for(j=1;j<=2;j++)// {// for(k=1;j<=2;k++)// {// t.mat[i][j]+=(x.mat[i][k]*y.mat[k][j])000;// t.mat[i][j]%=10000;// }// }// } return t;}int main(){ freopen("in.txt","r",stdin); int m; while(scanf("%d", &m) != EOF) { if(m == -1) break; init(); while(m) { if(m & 1) { res = multi(origin, res); } origin = multi(origin, origin); m = m >> 1; } printf("%d\n", res.mat[1][2] % 10000); } return 0;}
【矩阵加速】 矩阵 快速幂