【算法】矩阵快速求幂

1.算法描述

由于矩阵乘法具有结合律，因此A^4 = A * A * A * A = (A*A) * (A*A) = A^2 * A^2。我们可以得到这样的结论：当n为偶数时，A^n = A^(n/2) * A^(n/2)；当n为奇数时，A^n = A^(n/2) * A^(n/2) * A （其中n/2取整）。这就告诉我们，计算A^n也可以使用二分快速求幂的方法。例如，为了算出A^25的值，我们只需要递归地计算出A^12、A^6、A^3的值即可[1]。

在运算过程中，为了避免高精度运算，可采用取模（具体原因参看[1, 2]）。

2.Referrence

[1] Matrix67,  十个利用矩阵乘法解决的经典题目

[2] Matrix67,  同余运算及其基本性质

[3] shiwei408，矩阵快速幂 poj3070 3233
3735 3150

3.问题

3.1 POJ 3070

求Fibonacci数Fn = Fn −
1 + Fn − 2 for n ≥
2  对10000取余。

题目已经给出了优化解决方案：对矩阵求n次幂，取n次幂的[0][1]元素，即为Fn
。

源代码：

 

3070

Accepted

116K

0MS

C

1178B

2013-08-25 22:27:23

#include "stdio.h"
#include "string.h"

#define Mod 10000

typedef struct
{
int arr[2][2];
}Matrix;

/*multiply the two matrices*/
Matrix mMulti(Matrix a, Matrix b)
{
int i, j, k;
Matrix d;
memset(d.arr,0,sizeof(d.arr));
for(i=0;i<2;i++)
for(j=0;j<2;j++)
{
for(k=0;k<2;k++)
d.arr[i][j]+=a.arr[i][k]*b.arr[k][j];
d.arr[i][j]%=Mod;
}
return d;
}

/*the k-th power of matrix a*/
Matrix mPow(Matrix a, int k)
{
Matrix modd,meven;
modd.arr[0][0]=1; modd.arr[0][1]=0;  //matrix modd is initialized to be a unit matrix
modd.arr[1][0]=0; modd.arr[1][1]=1;
meven=a;                             //matrix meven is initialized

if(k==0)
return modd;
else
{
while(k!=1)
{
if(k&1)
{
k--;
modd=mMulti(modd,meven);
}
else
{
k/=2;
meven=mMulti(meven,meven);
}
}
return mMulti(modd,meven);
}
}

int main()
{
int n;
Matrix fibonacci, goal;
fibonacci.arr[0][0]=1; fibonacci.arr[0][1]=1;
fibonacci.arr[1][0]=1; fibonacci.arr[1][1]=0;

while(scanf("%d",&n)&&n!=-1)
{
goal=mPow(fibonacci,n);
printf("%d\n",goal.arr[0][1]%Mod);
}
return 0;
}

3.2 POJ 3233

求n × n方阵A + A2 + A3 +
… + Ak 

在Discuss中给出了优化解决思路：

| A+A^2+A^3+……A^k |   |A   A| (k-1）次方   | A |
|                   | = |     |              |   |
|     E             |   |0   E|              | E |

一次提交就AC了，不过代码还可以继续优化。

源代码：

3233

Accepted

284K

719MS

C

1451B

2013-08-27 23:04:22

#include "stdio.h"
#include "string.h"

#define MAX 60

typedef struct
{
int arr[MAX][MAX];
}Matrix;

int n,m;

/*multiplication*/
Matrix mMulti(Matrix a, Matrix b)
{
int i, j, k;
Matrix d;
memset(d.arr,0,sizeof(d.arr));
for(i=0;i<2*n;i++)
for(j=0;j<2*n;j++)
for(k=0;k<2*n;k++)
d.arr[i][j]=(d.arr[i][j]+a.arr[i][k]*b.arr[k][j])%m;
return d;
}

/*the k-th power of matrix a*/
Matrix mPow(Matrix a,int k)
{
int i;
Matrix modd,meven;
memset(modd.arr,0,sizeof(modd.arr));  //matrix modd is initialized to be a unit matrix
for(i=0;i<n;i++)
modd.arr[i][i]=1;
meven=a;                              //matrix meven is initialized

if(k==0)
return modd;
else
{
while(k!=1)
{
if(k&1)
{
k--;
modd=mMulti(modd,meven);
}
else
{
k/=2;
meven=mMulti(meven,meven);
}
}
return mMulti(modd,meven);
}
}

int main()
{
int i,j,k;
Matrix A,B,C;
scanf("%d%d%d",&n,&k,&m);
memset(A.arr,0,sizeof(A.arr));
for(i=0;i<n;i++)
for(j=0;j<n;j++)
scanf("%d",&A.arr[i][j]);

memcpy(B.arr,A.arr,sizeof(A.arr));
for(i=0;i<n;i++)
for(j=0;j<n;j++)
A.arr[i][j+n]=A.arr[i][j];
for(i=n;i<2*n;i++)
{
A.arr[i][i]=1;
B.arr[i][i-n]=1;
}

C=mMulti(mPow(A,k-1),B);

/*print the sum*/
for(i=0;i<n;i++)
{
for(j=0;j<n-1;j++)
printf("%d ",C.arr[i][j]);
printf("%d\n",C.arr[i][n-1]);
}
return 0;
}