一般意义上,求某一个矩阵的幂次,我们可以通过不断的迭代进行,但是要注意每一次运算都是O(n^3)的复杂度,如果幂次比较大,对时间要求紧一点,直接暴力明显就不行了。那应该怎么办呢?答案就是快速幂算法。这个算法的实现是依据矩阵乘法的结合率以及二分法。我们知道矩阵乘法具有结合律,如A^19 = A*A*...*A*A = (A*...*A) * (A*A)*A = A^16 * A^2 * A,走到这里大概就能知道快速幂的思想了,即二分因子,比如我们可以先算出因子A^2,然后由A^2推出A^4,再由A^4推出A^16,这样一来复杂度就从原来的O(n)降到了O(log(n))。
规律倒是出来了,那具体应该怎样实现呢?这里需要借助于二进制,我们知道任何整数都可以写成二进制的形式,如果将二进制看成离散的位的话,并且和前面所说的因子有联系的话,那么就可以通过这种内在的关系写出快速幂算法。可能大家已经发现了,在等式A^19 = A^16 * A^2 * A中,16,2,1不就是19对应二进制(10011)中值为1的位的权值么,所以我们就可以利用这个方法去算其中的某一个因子,进而求出答案了,下面给出核心算法:
Mat pow(Mat a, int k)
{
Mat c;
int i, j;
c = Mat();
c.n = a.n;
c.m = a.n;
for(i=1; i<=a.n; ++i) // 初始化为单位矩阵
c.mat[i][i] = 1;
while(k)
{
if(k & 1) // 二进制位为1则执行
c = c*a;
a = a*a;
k >>= 1;
}
return c;
}
上述代码中的Mat是上一次我们谈到的矩阵结构体,c矩阵需要初始化为单位矩阵,矩阵乘法也是前面已经重载过的。上述代码中,需要注意的是二进制的用法,即while循环,其实还有其它许多类算法都巧妙的用到了二进制,如树状数组、状压dp等,这些需要慢慢体会了。其实,上面的思想也适用于普通整数的大幂次求解,只需要将矩阵看成普通的整数即可,详见这里,传送门(点击打开链接)。
好了,有了快速幂算法,我们就能比较快速的算出矩阵的幂次了,下面就练习一道题,HDOJ:1575,时空转移(点击打开链接),题目如下:
Tr A
Time Limit: 1000/1000 MS (Java/Others) Memory Limit: 32768/32768 K (Java/Others)Total Submission(s): 2953 Accepted Submission(s): 2200
每组数据的第一行有n(2 <= n <= 10)和k(2 <= k < 10^9)两个数据。接下来有n行,每行有n个数据,每个数据的范围是[0,9],表示方阵A的内容。
2 2 2 1 0 0 1 3 99999999 1 2 3 4 5 6 7 8 9
2 2686
赤裸裸的矩阵幂运算。
分析:
一看到幂次k的范围就知道不能暴力,必须要用快速幂了,最后再将对角线上的元素加起来即可。
源代码:
#include <cstdio>
#include <cstring>
const int MAXN = 15;
const int MOD = 9973;
struct Mat
{
int n, m;
int mat[MAXN][MAXN];
Mat()
{
memset(mat, 0, sizeof(mat));
n = m = MAXN;
};
Mat operator * (Mat b)
{
Mat c;
c = Mat();
c.n = n;
c.m = b.m;
for(int i=1; i<=n; ++i) // 注意这里从1开始
for(int j=1; j<=b.m; ++j)
{
//if(mat[j][i] <= 0) continue; // 剪枝
for(int k=1; k<=m; ++k)
{
//if(b.mat[i][k] <= 0) continue; // 剪枝
c.mat[i][j] += (mat[i][k]*b.mat[k][j]) % MOD;
c.mat[i][j] %= MOD;
}
}
return c;
}
Mat pow(Mat a, int k)
{
Mat c;
c.n = a.n;
c.m = a.n;
for(int i=1; i<=a.n; ++i) // 初始化为单位矩阵
c.mat[i][i] = 1;
while(k)
{
if(k & 1)
c = c*a;
a = a*a;
k >>= 1;
}
return c;
}
void in(int x, int y)
{
n = x;
m = y;
for(int i=1; i<=x; ++i)
for(int j=1; j<=y; ++j)
scanf("%d", &mat[i][j]);
}
};
int main()
{//freopen("sample.txt", "r", stdin);
int cas;
scanf("%d", &cas);
while(cas--)
{
int n, k;
Mat data;
scanf("%d%d", &n, &k);
data.in(n, n);
data = data.pow(data, k);
int ans = 0;
for(int i=1; i<=n; ++i)
ans = (ans+data.mat[i][i]) % MOD;
printf("%d\n", ans);
}
return 0;
}
好了,到这里快速幂就告一段落了,不过还有一种快速幂的经典变形不得不说,那就是有矩阵A,求A + A^2 + A^3 + ... + A^k,当k很大时,就算用快速幂也只能快速的求出单个幂次,所以它的经典性就出来了,即两次二分。用快速幂二分求出A^i,然后再二分求和,即A + A^2 + A^3 + A^4 + A^5 + A^6 = (A + A^2 + A^3) + A^3*(A + A^2 + A^3)。如此一来,就不会超时啦!有一道类似的题目,POJ:3233,时空转移(点击打开链接),下面也给出一份源代码吧:
#include <cstdio>
#include <cstring>
const int MAXN = 35;
int MOD;
struct Mat
{
int n, m;
int mat[MAXN][MAXN];
Mat()
{
memset(mat, 0, sizeof(mat));
n = m = MAXN;
};
Mat operator * (Mat b) // 重载乘法
{
Mat c;
c = Mat();
c.n = n;
c.m = b.m;
for(int i=1; i<=n; ++i) // 注意这里从1开始
for(int j=1; j<=b.m; ++j)
{
//if(mat[j][i] <= 0) continue; // 剪枝
for(int k=1; k<=m; ++k)
{
//if(b.mat[i][k] <= 0) continue; // 剪枝
c.mat[i][j] += (mat[i][k]*b.mat[k][j]) % MOD;
c.mat[i][j] %= MOD;
}
}
return c;
}
Mat operator + (Mat b) // 重载加法
{
Mat c;
c.n = n;
c.m = m;
for(int i=1; i<=n; ++i)
for(int j=1; j<=m; ++j)
c.mat[i][j] = (mat[i][j]+b.mat[i][j]) % MOD;
return c;
}
Mat pow(Mat a, int k) // 快速幂
{
Mat c;
c.n = a.n;
c.m = a.n;
for(int i=1; i<=a.n; ++i) // 初始化为单位矩阵
c.mat[i][i] = 1;
while(k)
{
if(k & 1)
c = c*a;
a = a*a;
k >>= 1;
}
return c;
}
void in(int x, int y) // 输入矩阵
{
n = x;
m = y;
for(int i=1; i<=x; ++i)
for(int j=1; j<=y; ++j)
{
scanf("%d", &mat[i][j]);
mat[i][j] %= MOD;
}
}
void out() // 输出矩阵
{
for(int i=1; i<=n; ++i)
{
printf("%d ", mat[i][1]);
for(int j=2; j<=m; ++j)
printf("%d ", mat[i][j]);
puts("");
}
}
};
Mat tmp;
Mat slove(Mat init, int k)
{
if(k == 1)
return init;
tmp = slove(init, k>>1); // 递归二分求解
tmp = tmp + tmp * init.pow(init,k>>1);
if(k & 1)
return tmp + init.pow(init,k);
else
return tmp;
}
int main()
{//freopen("sample.txt", "r", stdin);
int n, k;
while(~scanf("%d%d%d", &n, &k, &MOD))
{
Mat data;
tmp.n = tmp.m = n;
data.in(n, n);
data = slove(data, k);
data.out();
}
return 0;
}
这份代码重用了几个以前写过的函数,所以看起来和前面的有些相似,但是原理需要注意,那就是二进制与二分的思想。下一文章我们继续讨论矩阵的另一个经典应用,递推式的求解,传送门( 点击打开链接)。
好了,快速幂就说到这里了,需要熟练掌握的话就只有多多练习了,HDOJ:2157。