题意:已知一个n*n的矩阵A,和一个正整数k,求S
思路:矩阵快速幂。首先我们知道 A^x 可以用矩阵快速幂求出来(具体上篇博文)。其次可以对k进行二分,每次将规模减半,分k为奇偶两种情况,如当k = 6和k = 7时有:
题意:已知一个n*n的矩阵A,和一个正整数k,求S = A + A2 + A3 + … + Ak。
思路:矩阵快速幂。首先我们知道 A^x 可以用矩阵快速幂求出来(具体可见poj 3070)。其次可以对k进行二分,每次将规模减半,分k为奇偶两种情况,如当k = 6和k = 7时有:
k = 6 有: S(6) = (1 + A^3) * (A + A^2 + A^3) = (1 + A^3) * S(3)。
k = 7 有: S(7) = A + (A + A^4) * (A + A^2 + A^3) = A + (A + A^4) * S(3)。
ps:对矩阵定义成结构体Matrix,求S时用递归,程序会比较直观,好写一点。当然定义成数组,然后再进行一些预处理,效率会更高些。
源代码:(880K,1188MS)
#include<iostream>
using namespace std;
const int MAX = 32;
struct Matrix
{
int v[MAX][MAX];
};
int n, k, M;
Matrix mtAdd(Matrix A, Matrix B) // 求矩阵 A + B
{
int i, j;
Matrix C;
for(i = 0; i < n; i ++)
for(j = 0; j < n; j ++)
{
C.v[i][j] = (A.v[i][j] + B.v[i][j]) % M;
}
return C;
}
Matrix mtMul(Matrix A, Matrix B) // 求矩阵 A * B
{
int i, j, k;
Matrix C;
for(i = 0; i < n; i ++)
for(j = 0; j < n; j ++)
{
C.v[i][j] = 0;
for(k = 0; k < n; k ++)
{
C.v[i][j] = (A.v[i][k] * B.v[k][j] + C.v[i][j]) % M;
}
}
return C;
}
Matrix mtPow(Matrix A, int k) // 求矩阵 A ^ k
{
if(k == 0) {
memset(A.v, 0, sizeof(A.v));
for(int i = 0; i < n; i ++)
{
A.v[i][i] = 1;
}
return A;
}
if(k == 1) return A;
Matrix C = mtPow(A, k / 2);
if(k % 2 == 0) {
return mtMul(C, C);
} else {
return mtMul(mtMul(C, C), A);
}
}
Matrix mtCal(Matrix A, int k) // 求S (k) = A + A2 + A3 + … + Ak
{
if(k == 1) return A;
Matrix B = mtPow(A, (k+1) / 2);
Matrix C = mtCal(A, k / 2);
if(k % 2 == 0) {
return mtMul(mtAdd(mtPow(A, 0), B), C); // 如S(6) = (1 + A^3) * S(3)。
} else {
return mtAdd(A, mtMul(mtAdd(A, B), C)); // 如S(7) = A + (A + A^4) * S(3)
}
}
int main()
{
int i, j;
Matrix A;
cin >> n >> k >> M;
for(i = 0; i < n; i ++)
for(j = 0; j < n; j ++)
{
cin >> A.v[i][j];
}
A = mtCal(A, k);
for(i = 0; i < n; i ++)
{
for(j = 0; j < n; j ++)
{
cout << A.v[i][j] << ' ';
}
cout << endl;
}
return 0;
}