本文主要是介绍POJ - 3233 - Matrix Power Series(矩阵分块 or 分治 + 矩阵快速幂),希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!
题目链接:http://poj.org/problem?id=3233
题意:已知一个n*n的矩阵A,和一个正整数k,求。
分治思路:首先我们知道 可以用矩阵快速幂求出来。其次可以对 进行分治,每次将规模减半,如下:
: 。
: 。
: 。
从上面几个式子可以发现,当k为奇数或者偶数的区别。
对于一个 是偶数则:。
如果 为奇数的话需要加上也就是:。
#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
*/
矩阵分块:构造一个的矩阵:其中为单位矩阵,为矩阵,矩阵的每个元素都是一个矩阵,我们可以把分块矩阵展开,它的乘法和普通矩阵的乘法是一样的。取右上角的矩阵然后减去一个单位矩阵就是答案。
#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 分治 + 矩阵快速幂)的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!