首页 > 代码库 > 【矩阵加速】 矩阵 快速幂

【矩阵加速】 矩阵 快速幂

矩阵的快速幂是用来高效地计算矩阵的高次方的。将朴素的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;}

 

【矩阵加速】 矩阵 快速幂