在应用矩阵的快速幂运算解决一些递推问题时,由于递推式不是一个直接的线性关系,这样不能直接简单地得到用于运算的矩阵,需要进行适当的构造。下面先看一道POJ 上的经典题目。
【例1】Matrix Power Series (POJ 3233)
Description
Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak.
Input
The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 109) and m (m < 104). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order.
Output
Output the elements of S modulo m in the same way as A is given.
Sample Input
2 2 4
0 1
1 1
Sample Output
1 2
2 3
(1)编程思路。
题目的意思是:输入n*n矩阵A,求S=A + A^2 + A^3 + … + A^k的结果(两个矩阵相加就是对应位置分别相加),输出的数据mod m。
显然不能将各个矩阵的2~k次幂均计算出来,会超时的。
可以构造出一个2n*2n的矩阵P,采用分块的方法表示如下,,其中 A 为原矩阵,I 为单位矩阵,O 为0矩阵。
我们发现 Pk+1 右上角那一个分块n*n矩阵正是要求的 A+A2+...+Ak 。
于是我们构造出 P矩阵,然后对它求矩阵快速幂,最后减去一个单位阵即可得到结果。
(2)源程序。
#include <stdio.h>
#include <string.h>
struct Matrix
{
int mat[61][61]; // 存储矩阵中各元素
};
Matrix matMul(Matrix a ,Matrix b,int n,int m)
{
Matrix c;
memset(c.mat,0,sizeof(c.mat));
int i,j,k;
for (i = 1; i<=n ; i++)
for (j=1 ;j<=n ; j++)
for (k = 1 ;k<=n ;k++)
{
c.mat[i][j]=(c.mat[i][j]+a.mat[i][k] * b.mat[k][j]) % m;
}
return c;
}
Matrix quickMatPow(Matrix a ,int n,int b,int m) // n阶矩阵a快速b次幂
{
Matrix c;
memset(c.mat ,0 ,sizeof(c.mat));
int i;
for (i = 1 ;i <= n ;i++)
c.mat[i][i] = 1;
while (b!=0)
{
if (b & 1)
c = matMul(c ,a ,n,m); // c=c*a;
a = matMul(a ,a ,n,m); // a=a*a
b /= 2;
}
return c;
}
int main()
{
int n,k,m,i,j;
Matrix p ;
scanf("%d%d%d",&n,&k,&m);
memset(p.mat,0,sizeof(p.mat));
for (i=1;i<=n;i++)
for (j=1;j<=n;j++)
scanf("%d",&p.mat[i][j]);
for (i=1;i<=n;i++)
{
p.mat[i][n+i]=1; // 右上角的单位矩阵
p.mat[n+i][n+i]=1; // 右下角的单位矩阵
}
p = quickMatPow(p,2*n,k+1,m);
for (i=1;i<=n;i++)
p.mat[i][i+n]--; // 减去单位矩阵
for (i=1;i<=n;i++)
{
for (j=n+1;j<