/*
* 思路:这题数据k<10^9,很显然需要O(log n)的算法,自然就是矩阵
* 快速幂.矩阵快速幂和快速幂差不多就是求矩阵A的k次方时,每次折半
* 类似于二分的思想,比如要求A的100次方,我们先求出A的50次方,
* 然后用得到的结果平方就好了.试想如果要算16次方,每次折半,我们
* 只需要做4次矩阵乘法,但是如果采用朴素的连乘则需要15次,效率对比
* 可想而知.
* 但是还有一个问题,我们用矩阵快速幂可以求得一个矩阵的n次方.但是
* 这题要求的是一系列的矩阵幂之和,就算每个矩阵幂都用上上面的快速
* 幂也会超时啊!表示俺线性代数学得太糟糕了,看了网友的构造方法,
* 才学会转换.如下:
* 构造一个矩阵
* B = |A En|
* |0 En|
* 那么B^(k+1) = |A^(k+1) En+A+A^2+A^3+...+A^k|
* | 0 En |
* 其中En表示n阶单位矩阵.
* 所以要求A+A^2+...+A^k,只需要求出B^(k+1),然后右上角那个分
* 块矩阵减单位矩阵就可以了.方法真是看似简单,但是炒鸡厉害...
* 计算思想乳此感人,代码就算了,各有各的风格.如果你看到了本文,
* 程序千万要自己实现,别直接模仿或是抄写别人的代码.
*/
#include <cstdio>
#include <cstring>
#include <vector>
#include <iostream>
using namespace std;
const int MAX = 65;
typedef long long ll;
int n, k, m;
int SIZE;
struct mat {
int a[MAX][MAX];
mat(){}
mat(int b) {
memset(a, 0, sizeof(a));
for (int i = 0; i < MAX; ++i) {
a[i][i] = b;
}
}//init
mat operator*(mat b) {
mat tmp(0);
int sum;
/************矩阵乘法***************/
/*for (int i = 0; i < MAX; ++i) {
for (int j = 0; j < MAX; ++j) {
sum = 0;
for (int k = 0; k < MAX; ++k) {
sum += a[i][k]*b.a[k][j]; sum %= m;
}
tmp.a[i][j] = sum;
}
}*/
/**********优化版,原理以后学了计组再说*********/
for (int i = 0; i < MAX; ++i) {
for (int k = 0; k < MAX; ++k) {
for (int j = 0; j < MAX; ++j) {
tmp.a[i][j] += a[i][k]*b.a[k][j]; tmp.a[i][j] %= m;
}
}
}
return tmp;
}
mat operator=(mat b) {
memcpy(a, b.a, sizeof(b.a));
return *this;
}
};
mat YY(1);
mat fast_pow(mat a, int x) {
mat res(1);
while (x) {
if (x&1) res = res*a;
a = a*a;
x >>= 1;
}
return res;
}
int main() {
mat M(MAX);
while (~scanf(" %d %d %d", &n, &k, &m)) {
SIZE = n << 1;
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
scanf(" %d", &M.a[i][j]);
}
}
/* 我们要的是构造后的矩阵B,上面输入存入的只是B左上角的
* 分块矩阵,另外三块还要手动搞
*/
/* 右上角n阶单位矩阵 */
for (int i = 0; i < n; ++i) {
for (int j = n; j < SIZE; ++j) {
if (j-n == i) M.a[i][j] = 1;
else M.a[i][j] = 0;
}
}
/* 左下角0矩阵 */
for (int i = n; i < SIZE; ++i) {
for (int j = 0; j < n; ++j) {
M.a[i][j] = 0;
}
}
/* 右下角n阶单位矩阵 */
for (int i = n; i < SIZE; ++i) {
for (int j = n; j < SIZE; ++j) {
M.a[i][j] = (i == j ? 1 : 0);
}
}
mat res = fast_pow(M, k+1);
/* 减去单位矩阵 */
for (int i = 0; i < n; ++i) {
res.a[i][i+n]--;
if (res.a[i][i+n] < 0) res.a[i][i+n] += m;
}
/* 打印 */
for (int i = 0; i < n; ++i) {
for (int j = n; j < SIZE; ++j) {
printf("%d%c",res.a[i][j], (j^(SIZE-1) ? ' ' : '\n'));
}
}
}
return 0;
}
poj3233Matrix Power Series矩阵快速幂
最新推荐文章于 2022-05-05 20:34:22 发布