[POJ 3233] Matrix Power Series【分治 / 矩乘优化】

  • 2018-01-07
  • 0
  • 0

Problem:

Time Limit: 3000MS Memory Limit: 131072K

Description

Given a n × n matrix A and a positive integer k, find the sum S = A + A2 + A3 + … + Ak.

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

Source

POJ Monthly--2007.06.03, Huang, Jinsong

Solution:

本题是一道很有思维含量的题。

首先,我们可以明确暴力模拟的复杂度高达 O(kn3logk),TLE 炸上天。。

于是我们想到可以用分治来降低复杂度,即

  • 若 k 为偶数,则 S = A + A2 + A3 + … + Ak = (I + Ak/2) × (A + A2 + A3 + … + Ak/2)
  • 若 k 为奇数,则 S = A + A2 + A3 + … + Ak = (I + A[k/2]) × (A + A2 + A3 + … + A[k/2]) + Ak

这种方法与快速幂的思想是一致的,单层复杂度也为 O(logk)。

所以分治的复杂度为 O(n3log2k),可以 AC。详见 Code #1

但是分治并不是该题的最优算法,矩乘优化可以将复杂度进一步降低。

矩乘优化?这个 A 不是本身就是矩阵了吗?

我们可以矩阵套矩阵啊。。

构造矩阵 B,对其进行快速幂运算后,取出右上角的小矩阵即为答案。

所以矩乘优化的复杂度为 O(n3logk)。详见 Code #2

Code #1: O(n3log2k) [1000K, 1735MS]

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<cassert>
#include<iostream>
#include<algorithm>
using namespace std;

int k, m;

struct Matrix{
	int ele[32][32], size;
	
	inline void setIdentity(int idensize){
		memset(ele, 0, sizeof(ele));
		size = idensize;
		for(register int i = 1; i <= size; i++) ele[i][i] = 1;
	}
	
	inline Matrix operator + (const Matrix &mat2) const {
	//	assert(size == mat2.size);
		Matrix res; res.size = size;
		for(register int i = 1; i <= size; i++)
			for(register int j = 1; j <= size; j++)
				res.ele[i][j] = (ele[i][j] + mat2.ele[i][j]) % m;
		return res;
	}
	
	inline Matrix operator * (const Matrix &mat2) const {
	//	assert(size == mat2.size);
		Matrix res; res.size = size;
		for(register int i = 1; i <= size; i++)
			for(register int j = 1; j <= mat2.size; j++){
				res.ele[i][j] = 0;
				for(register int k = 1; k <= size; k++)
					res.ele[i][j] = (res.ele[i][j] + ele[i][k] * mat2.ele[k][j]) % m;
			}
		return res;
	}
	
	inline Matrix pow(int ex){
		Matrix bas = *this, res; res.setIdentity(size);
		while(ex){
			if(ex & 1) res = res * bas;
			ex >>= 1, bas = bas * bas;
		}
		return res;
	}
} A, I;

inline Matrix sum(int ex){  // Return the sum of the geometric sequence
	if(ex == 1) return A;
	if(ex & 1) return (I + A.pow(ex >> 1)) * sum(ex >> 1) + A.pow(ex);
	return (I + A.pow(ex >> 1)) * sum(ex >> 1);
}

int main(){
	scanf("%d%d%d", &A.size, &k, &m);
	I.setIdentity(A.size);
	for(register int i = 1; i <= A.size; i++)
		for(register int j = 1; j <= A.size; j++)
			scanf("%d", &A.ele[i][j]);
	A = sum(k);
	for(register int i = 1; i <= A.size; i++){
		for(register int j = 1; j < A.size; j++)
			printf("%d ", A.ele[i][j]);
		printf("%d\n", A.ele[i][A.size]);
	}
	return 0;
} 

Code #2: O(n3logk) [372K, 719MS]

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<cassert>
#include<iostream>
#include<algorithm>
using namespace std;

int k, m;

struct Matrix{
	int ele[32][32], size;
	
	inline void setIdentity(int idsize, int x){
		size = idsize;
		for(register int i = 1; i <= size; i++)
			for(register int j = 1; j <= size; j++)
				if(i == j) ele[i][j] = x;
				else ele[i][j] = 0;
	}
	
	inline Matrix operator + (const Matrix &mat2) const {
	//	assert(size == mat2.size);
		Matrix res; res.size = size;
		for(register int i = 1; i <= size; i++)
			for(register int j = 1; j <= size; j++)
				res.ele[i][j] = (ele[i][j] + mat2.ele[i][j]) % m;
		return res;
	}
	
	inline Matrix operator * (const Matrix &mat2) const {
	//	assert(size == mat2.size);
		Matrix res; res.size = size;
		for(register int i = 1; i <= size; i++)
			for(register int j = 1; j <= mat2.size; j++){
				res.ele[i][j] = 0;
				for(register int k = 1; k <= size; k++)
					res.ele[i][j] = (res.ele[i][j] + ele[i][k] * mat2.ele[k][j]) % m;
			}
		return res;
	}
} A, I, O;

struct matMatrix{
	Matrix ele[3][3];
	int size;
	
	inline matMatrix operator * (const matMatrix &mat2) const {
	//	assert(size == mat2.size);
		matMatrix res; res.size = size;
		for(register int i = 1; i <= size; i++)
			for(register int j = 1; j <= mat2.size; j++){
				res.ele[i][j].setIdentity(A.size, 0);
				for(register int k = 1; k <= size; k++)
					res.ele[i][j] = res.ele[i][j] + ele[i][k] * mat2.ele[k][j];
			}
		return res;
	}
	
	inline matMatrix pow(int ex){	
		matMatrix bas = *this, res = *this; ex--;
		while(ex){
			if(ex & 1) res = res * bas;
			ex >>= 1, bas = bas * bas;
		}
		return res;
	}
} B;

int main(){
	scanf("%d%d%d", &A.size, &k, &m);
	I.setIdentity(A.size, 1), O.setIdentity(A.size, 0);
	for(register int i = 1; i <= A.size; i++)
		for(register int j = 1; j <= A.size; j++)
			scanf("%d", &A.ele[i][j]);
	B.size = 2;
	B.ele[1][1] = A, B.ele[1][2] = A;
	B.ele[2][1] = O, B.ele[2][2] = I;
	B = B.pow(k), A = B.ele[1][2];
	for(register int i = 1; i <= A.size; i++){
		for(register int j = 1; j < A.size; j++)
			printf("%d ", A.ele[i][j]);
		printf("%d\n", A.ele[i][A.size]);
	}
	return 0;
}

评论

还没有任何评论,你来说两句吧



常年不在线的QQ:
49750

不定期更新的GitHub:
https://github.com/Darkleafin


OPEN AT 2017.12.10

如遇到代码不能正常显示的情况,请刷新页面。
If the code cannot be displayed normally, please refresh the page.


发现一个优美的网站:
https://visualgo.net/en
















- Theme by Qzhai