【算法】矩阵快速求幂

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数F_n=F_{n− 1}+F_{n− 2}forn≥ 2 对10000取余。

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

_源代码：

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+A²+A³+ … +A^k

^{在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;
}