POJ - 3233 - Matrix Power Series（矩阵分块 or 分治 + 矩阵快速幂）

本文主要是介绍POJ - 3233 - Matrix Power Series（矩阵分块 or 分治 + 矩阵快速幂），希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

题目链接：http://poj.org/problem?id=3233

题意：已知一个n*n的矩阵A，和一个正整数k，求 $S = A + A^2 + A^3 +...+ A^k$ 。

分治思路：首先我们知道 $A^x$ 可以用矩阵快速幂求出来。其次可以对 $k$ 进行分治，每次将规模减半，如下：

$k = 10$ ： $S(9) = ( A^1+A^2+A^3+A^4+ A^5 ) + A^5 * ( A^1+A^2+A^3+A^4+A^5 )$ $= S(5) + A^5 * S(5)$ 。

$k = 5$ ： $S(5) = ( A^1+A^2 ) + A^3 + A^3 * ( A^1+A^2 ) = S(2) + A^3 + A^3 * S(2)$ 。

$k = 2$ : $S(2) = A^1 + A^2 = S(1) + A^1 * S(1)$ 。

从上面几个式子可以发现，当k为奇数或者偶数的区别。

对于一个 $k$ 是偶数则： $S(k) = A^{k/2} *S(k/2)+S(k/2)$ 。

如果 $k$ 为奇数的话需要加上 $A^{k/2 + 1}$ 也就是： $S(k) = A^{k/2+1}*S(k/2+1)+S(k/2+1)+A^{k/2 + 1}$ 。

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <algorithm>
using namespace std;
typedef long long ll;
const int MAXN = 2;
int n, k, MOD;
struct Matrix
{int m[50][50];void init(int x){memset(m, 0, sizeof(m));for(int i = 0; i < n; i++){m[i][i] = x;}}Matrix operator * (const Matrix &a) const{Matrix res;res.init(0);for(int i = 0; i < n; i++)for(int j = 0; j < n; j++)for(int k = 0; k < n; k++)res.m[i][j] = (res.m[i][j] + (m[i][k] * a.m[k][j]) % MOD) % MOD;return res;}Matrix operator + (const Matrix &a) const{Matrix ans;ans.init(0);for(int i = 0; i < n; i++)for(int j = 0; j < n; j++)ans.m[i][j] = (ans.m[i][j] + (m[i][j] + a.m[i][j]) % MOD) % MOD;return ans;}
};Matrix quickPow(Matrix x, int p)
{Matrix res;res.init(1);x = x * res;while(p){if(p & 1) res = res * x;p >>= 1;x = x * x;}return res;
}Matrix a, x;void init()
{for(int i = 0; i < n; i++){for(int j = 0; j < n; j++){scanf("%d", &a.m[i][j]);a.m[i][j] %= MOD;}}
}
Matrix dfs(Matrix c, int k)
{if(k == 1) return c;Matrix dx, dy;dx = dfs(c, k / 2);if(k & 1){dy = quickPow(c, k / 2 + 1);return dx + dy + (dx * dy);}else{dy = quickPow(c, k/2);return dx + (dx * dy);}
}
int main()
{scanf("%d%d%d", &n, &k, &MOD);init();x = dfs(a, k);for(int i = 0; i < n; i++)for(int j = 0; j < n; j++)printf("%d%c",x.m[i][j], j == n - 1 ? '\n' : ' ');return 0;
}
/*
2 2 4
0 1
1 1
------
1 2
2 3
*/

矩阵分块：构造一个 $2n*2n$ 的矩阵： $\begin{bmatrix} A&E\\ zero&E \end{bmatrix}$ 其中 $E$ 为单位矩阵， $zero$ 为 $0$ 矩阵，矩阵的每个元素都是一个矩阵，我们可以把分块矩阵展开，它的乘法和普通矩阵的乘法是一样的。 $\begin{bmatrix} A&E\\ zero&E \end{bmatrix}^{k+1}=\begin{bmatrix} A ^{k+1}&A^{k}+A^{k-1}+...+A+E\\ zero&E \end{bmatrix}^{k+1}$ 取右上角的矩阵然后减去一个单位矩阵 $E$ 就是答案。

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <algorithm>
using namespace std;
typedef long long ll;
const int MAXN = 2;
ll n, k, MOD;
struct Matrix
{ll m[80][80];void init(int x){memset(m, 0, sizeof(m));for(int i = 0; i < n * 2; i++){m[i][i] = x;}}Matrix operator * (const Matrix &a) const{Matrix res;res.init(0);for(int i = 0; i < n * 2; i++)for(int j = 0; j < n* 2; j++)for(int k = 0; k < n * 2; k++)res.m[i][j] = (res.m[i][j] + (m[i][k] * a.m[k][j]) % MOD) % MOD;return res;}Matrix operator + (const Matrix &a) const{Matrix ans;ans.init(0);for(int i = 0; i < MAXN; i++)for(int j = 0; j < MAXN; j++)ans.m[i][j] = (ans.m[i][j] + (m[i][j] + a.m[i][j]) % MOD) % MOD;return ans;}
};Matrix quickPow(Matrix x, ll p)
{Matrix res;res.init(1);x = x * res;while(p){if(p & 1) res = res * x;p >>= 1;x = x * x;}return res;
}Matrix a, x;void init()
{for(int i = 0; i < n; i++){for(int j = 0; j < n; j++){scanf("%lld", &a.m[i][j]);a.m[i][j] %= MOD;}a.m[i][i+n] = 1; //构造右上角a.m[i+n][i+n] = 1; //构造左下角}
}
int main()
{while(~scanf("%lld%lld%lld", &n, &k, &MOD)){init();a = quickPow(a, k + 1);for(int i = 0; i < n; i++){for(int j = 0; j < n; j++){x.m[i][j] = a.m[i][j+n];if(i == j)//减去单位矩阵Ex.m[i][j] = (x.m[i][j] - 1 + MOD) % MOD;}}for(int i = 0; i < n; i++)for(int j = 0; j < n; j++)printf("%lld%c",x.m[i][j], j == n - 1 ? '\n' : ' ');}return 0;
}
/*
2 2 4
0 1
1 1
------
1 2
2 3
*/

这篇关于POJ - 3233 - Matrix Power Series（矩阵分块 or 分治 + 矩阵快速幂）的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！