POJ 3233 Matrix Power Series 矩阵快速幂

最新推荐文章于 2022-05-05 20:34:22 发布

原创最新推荐文章于 2022-05-05 20:34:22 发布 · 241 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#poj #矩阵快速幂

挑战程序设计同时被 2 个专栏收录

70 篇文章

订阅专栏

数论

26 篇文章

订阅专栏

本文介绍了一种解决矩阵幂级数求和问题的有效方法，通过构造辅助矩阵并利用矩阵快速幂运算，实现了对给定矩阵A从A到A^k的所有幂次的和的高效计算。

题目链接

Matrix Power Series
Time Limit: 3000MS Memory Limit: 131072K
Description
Given a n × n matrix A and a positive integer k, find the sum S = A + A^2 + A^3 + … + A^k.
Input
The input contains exactly one test case. The first line of input contains three positive integers n (n ≤ 30), k (k ≤ 109) and m (m < 104). Then follow n lines each containing n nonnegative integers below 32,768, giving A’s elements in row-major order.
Output
Output the elements of S modulo m in the same way as A is given.
Sample Input
2 2 4
0 1
1 1
Sample Output
1 2
2 3

题目大意：
求矩阵A的幂的和。

解题思路：
S_k = A + A^2 + A^3 + … + A^k
S_k 与A^k 一同变化，那么我们便可构造如下辅助矩阵( I 代表单位矩阵)：
$\left[ \begin{matrix} A^k\\ S_k \end{matrix} \right] =\left[ \begin{matrix} A & 0\\ A & I \end{matrix} \right] *\left[ \begin{matrix} A^{k-1}\\ S_{k-1} \end{matrix} \right] =\left[ \begin{matrix} A & 0\\ A & I \end{matrix} \right]^{k-1} *\left[ \begin{matrix} A^1\\ S_1 \end{matrix} \right]$
又有：
$\left[ \begin{matrix} A^1\\ S_1 \end{matrix} \right] =\left[ \begin{matrix} A\\ A \end{matrix} \right]=\left[ \begin{matrix} A & 0\\ A & I \end{matrix} \right] *\left[ \begin{matrix} I\\ 0 \end{matrix} \right]$
最后可得：
$\left[ \begin{matrix} A^k\\ S_k \end{matrix} \right] =\left[ \begin{matrix} A & 0\\ A & I \end{matrix} \right]^{k} *\left[ \begin{matrix} I\\ 0 \end{matrix} \right]$

AC代码：

#include<iostream>
#include<vector>
#include<stdio.h>
using namespace std;
typedef long long ll;
typedef vector<ll> vec;
typedef vector<vec> mat;
int M;
mat mul(mat a, mat b) {
	mat c(a.size(), vec(b[0].size()));
	for (int i = 0; i < a.size(); i++) {
		for (int k = 0; k < b.size(); k++) {
			for (int j = 0; j < b[0].size(); j++) {
				c[i][j] = (c[i][j] + a[i][k] * b[k][j]) % M;
			}
		}
	}
	return c;
}
mat mod_pow(mat a, ll n) {
	mat b(a.size(), vec(a[0].size()));
	for (int i = 0; i < a.size(); i++) {
		b[i][i] = 1;
	}
	while (n > 0) {
		if (n & 1) b = mul(a, b);
		a = mul(a, a);
		n >>= 1;
	}
	return b;
}
int main() {
	int n; ll k;
	scanf("%d%lld%d", &n, &k, &M);
	mat a(2*n, vec(2*n));//构造辅助矩阵
	for (int i = 0; i < n; i++) {
		for (int j = 0; j < n; j++) {
			scanf("%lld", &a[i][j]);
			a[i + n][j] = a[i][j];
		}
		a[i + n][i + n] = 1;
	}
	a = mod_pow(a, k);
	//此时矩阵中的左下角即为所求
	for (int i = 0; i < n; i++) {
		printf("%lld", a[i + n][0]);
		for (int j = 1; j < n; j++) {
			printf(" %lld", a[i + n][j]);
		}
		printf("\n");
	}
	//system("pause");
}