POJ-3233 Matrix Power Series (矩阵快速幂 + 二分 或 构造矩阵)

1.题面

http://poj.org/problem?id=3233

2.题意

给你一个n*n的矩阵A，让你计算A^1 + A^2 + A^3 + ……+A^k，最终答案需要对k取模。

3.思路

这道题目有两种解题思路，都需要用到矩阵快速幂

1.二分递归的方法

我们首先记S = A^1 + A^2 + A^3 + ……+A^k；

如果k是偶数，那么，

S = （A^1 + A^2 + A^3 + ……+A^(k/2)) + A^(k/2)*[A^1 + A^2 + A^3 + ……+A^(k/2)]

然后递归地求解 （A^1 + A^2 + A^3 + ……+A^(k/2))。

如果k是奇数，我们先计算A^k，然后将k-1化为偶数

2.构造矩阵的方法

我使用第一种方法做的，第二种方法就要更巧妙一些。

我们构造上图中的矩阵，最终答案是矩阵的n+1次方的右上角部分。

果然数学才是第一生产力啊。

4.代码

代码贼丑，瞎看看吧。

第一种方法：

template <class RandomAccessIterator>
void Print(RandomAccessIterator first, RandomAccessIterator last, char endl1='\n', char endl2=' '){
	while (first != last)
		cout << (*first++) << (last-1==first?endl1:endl2);
	}

const int debug = 1;
const int size  = 5000 + 10; 
const int INF = INT_MAX>>1;
typedef long long ll;

const int matsize = 31;
int MOD;

int n, k, m;

struct _matrix{
	int m[matsize][matsize];
};

_matrix operator * (const _matrix& a, const _matrix& b){
	_matrix ret;
	memset(ret.m, 0, sizeof(ret.m));
	for (int i = 0; i < n; i++){
		for (int j = 0; j < n; j++){
			for (int k = 0; k < n; k++){
				ret.m[i][j] += a.m[i][k]*b.m[k][j];
				ret.m[i][j] %= MOD;
			}
		}
	}
	return ret;
}

_matrix operator ^ (_matrix a, int p){
	_matrix ret;
	memset(ret.m, 0, sizeof(ret.m));
	for (int i = 0; i < n; i++)
		ret.m[i][i] = 1;
	while (p > 0){
		if (p&1){
			ret = ret*a;
		}
		p >>= 1;
		a = a*a;
	}
	return ret;
}


_matrix operator + (const _matrix& a, const _matrix& b){
	_matrix ret;
	memset(ret.m, 0, sizeof(ret.m));
	for (int i = 0; i < n; i++){
		for (int j = 0; j < n; j++){
			ret.m[i][j] = (a.m[i][j] + b.m[i][j])%MOD;
		}
	}
	return ret;
}

_matrix SUM(_matrix m, int k){
	if (k == 1) return m;

	int mid = k>>1;
	_matrix ret;
	_matrix tmp = SUM(m, mid);
	if (k%2 == 1){
		ret = tmp + (m^mid)*tmp + (m^k);
	}else {
		ret = tmp + (m^mid)*tmp;
	}
	return ret;
}

int main(){
	// std::ios::sync_with_stdio(false);cin.tie(0);
	while (~scanf("%d%d%d", &n, &k, &MOD)){
		_matrix m;
		for (int i = 0; i < n; i++){
			for (int j = 0; j < n; j++){
				scanf("%d", &m.m[i][j]);
			}
		}
		_matrix ans = SUM(m, k);
		for (int i = 0; i < n; i++){
			Print(ans.m[i], ans.m[i]+n, '\n', ' ');
		}

	}
	return 0;
}

第二种方法

template <class RandomAccessIterator>
void Print(RandomAccessIterator first, RandomAccessIterator last, char endl1='\n', char endl2=' '){
	while (first != last)
		cout << (*first++) << (last-1==first?endl1:endl2);
}


const int debug = 1;
const int size  = 5000 + 10; 
const int INF = INT_MAX>>1;
typedef long long ll;

const int matsize = 70;
int MOD;

int n, k, m;

struct _matrix{
	int m[matsize][matsize];
};

_matrix operator * (const _matrix& a, const _matrix& b){
	_matrix ret;
	memset(ret.m, 0, sizeof(ret.m));
	for (int i = 0; i < n; i++){
		for (int j = 0; j < n; j++){
			for (int k = 0; k < n; k++){
				ret.m[i][j] += a.m[i][k]*b.m[k][j];
				ret.m[i][j] %= MOD;
			}
		}
	}
	return ret;
}

_matrix operator ^ (_matrix a, int p){
	_matrix ret;
	memset(ret.m, 0, sizeof(ret.m));
	for (int i = 0; i < n; i++)
		ret.m[i][i] = 1;
	while (p > 0){
		if (p&1){
			ret = ret*a;
		}
		p >>= 1;
		a = a*a;
	}
	return ret;
}


_matrix operator + (const _matrix& a, const _matrix& b){
	_matrix ret;
	memset(ret.m, 0, sizeof(ret.m));
	for (int i = 0; i < n; i++){
		for (int j = 0; j < n; j++){
			ret.m[i][j] = (a.m[i][j] + b.m[i][j])%MOD;
		}
	}
	return ret;
}


int main(){
	// std::ios::sync_with_stdio(false);cin.tie(0);
	while (~scanf("%d%d%d", &n, &k, &MOD)){
		_matrix m;
		memset(m.m, 0, sizeof(m.m));
		for (int i = 0; i < n; i++){
			for (int j = 0; j < n; j++){
				scanf("%d", &m.m[i][j]);
			}
		}
		for (int i = 0; i < n; i++){
			m.m[i][i+n] = m.m[i+n][i+n] = 1;
		}
		n <<= 1;

		m = m^(k+1);
		for (int i = 0; i < n/2; i++){
			for (int j = n/2; j < n; j++){
				if (j == i+n/2)
					printf("%d", (m.m[i][j]-1+MOD)%MOD);
				else 
					printf("%d", m.m[i][j]);
				if (j == n-1) printf("\n");
				else printf(" ");
			}
		}

	}
	return 0;
}