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

• 2018-01-07
• 0
• 0

Problem:

Description

Input

Output

Sample Input

2 2 4
0 1
1 1

Sample Output

1 2
2 3

Source

Solution:

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

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

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

• 若 k 为偶数，则 S = A + A² + A³ + … + A^k = (I + A^k/2) × (A + A² + A³ + … + A^k/2)
• 若 k 为奇数，则 S = A + A² + A³ + … + A^k = (I + A^[k/2]) × (A + A² + A³ + … + A^[k/2]) + A^k

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

所以分治的复杂度为 O(n³log²k)，可以 AC。详见 Code #1。

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

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

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

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

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

Code #1: O(n³log²k) [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(n³logk) [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;
}