UVa 11149 Power of Matrix 矩阵快速幂

题意：

给出一个(n imes n)的矩阵(A)，求(A+A^2+A^3+ cdots + A^k)。

分析：

这题是有(k=0)的情况，我们一开始先特判一下，直接输出单位矩阵(E)。
下面讨论(k > 0)的情况：

方法一

设答案为(S_k(k > 0))
把矩阵增广一下
$egin{bmatrix}
A & O
E & E
end{bmatrix}
egin{bmatrix}
A^n
S_{n-1}
end{bmatrix}

egin{bmatrix}
A^{n+1}
S_n
end{bmatrix}( )E(表示单位矩阵，)O(是全为零的矩阵。 ){egin{bmatrix}
A & O
E & E
end{bmatrix}}^k
egin{bmatrix}
A
O
end{bmatrix}

egin{bmatrix}
A^{k+1}
S_k
end{bmatrix}$

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

const int maxn = 40;
const int MOD = 10;

int mul(int a, int b) { return a * b % MOD; }

void add(int& a, int b) { a += b; if(a >= MOD) a -= MOD; }

int n, k, sz;
int a[maxn][maxn], b[maxn][maxn];

struct Matrix
{
	int a[maxn * 2][maxn * 2];

	Matrix() { memset(a, 0, sizeof(a)); }

	Matrix operator * (const Matrix& t) const {
		Matrix ans;
		for(int i = 0; i < sz; i++)
			for(int j = 0; j < sz; j++)
				for(int k = 0; k < sz; k++)
					add(ans.a[i][j], mul(a[i][k], t.a[k][j]));
		return ans;
	}

	void output() {
		for(int i = 0; i < sz; i++) {
			for(int j = 0; j < sz - 1; j++)
				printf("%d ", a[i][j]);
			printf("%d
", a[i][sz - 1]);
		}
		printf("
");
	}
};

Matrix pow_mod(Matrix a, int p) {
	Matrix ans;
	for(int i = 0; i < sz; i++) ans.a[i][i] = 1;
	while(p) {
		if(p & 1) ans = ans * a;
		a = a * a;
		p >>= 1;
	}
	return ans;
}

int main()
{
	while(scanf("%d%d", &n, &k) == 2 && n) {
		Matrix M;
		for(int i = 0; i < n; i++)
			for(int j = 0; j < n; j++) {
				scanf("%d", &a[i][j]);
				a[i][j] %= MOD;
				M.a[i][j] = a[i][j];
			}

		if(!k) {
			for(int i = 0; i < n; i++) {
				for(int j = 0; j < n - 1; j++)
					printf("%d ", i == j ? 1 : 0);
				printf("%d
", i == n - 1 ? 1 : 0);
			}
			printf("
");
			continue;
		}
		
		for(int i = n; i < n * 2; i++)
			M.a[i][i] = M.a[i][i - n] = 1;

		sz = n * 2;
		M = pow_mod(M, k);
		memset(b, 0, sizeof(b));
		for(int i = 0; i < n; i++)
			for(int j = 0; j < n; j++) if(M.a[i + n][j])
				for(int k = 0; k < n; k++)
					add(b[i][k], mul(M.a[i + n][j], a[j][k]));

		for(int i = 0; i < n; i++) {
			for(int j = 0; j < n - 1; j++)
				printf("%d ", b[i][j]);
			printf("%d
", b[i][n - 1]);
		}
		printf("
");
	}

	return 0;
}

方法二

有如下递归式：

• (S_k=(E+A^{frac{k}{2}})S^{frac{k}{2}})，k是偶数
• (S_k=(E+A^{frac{k}{2}})S^{frac{k}{2}}+A^k)，k是奇数

所以也可以直接递归求解。

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;

const int maxn = 40;
const int MOD = 10;

int mul(int a, int b) { return a * b % MOD; }

void add(int& a, int b) { a += b; if(a >= MOD) a -= MOD; }

int n, k;

struct Matrix
{
	int a[maxn][maxn];

	Matrix() { memset(a, 0, sizeof(a)); }
	
	void E() {
		memset(a, 0, sizeof(a));
		for(int i = 0; i < maxn; i++) a[i][i] = 1;
	}
	
	Matrix operator + (const Matrix& t) const {
		Matrix ans;
		for(int i = 0; i < n; i++)
			for(int j = 0; j < n; j++) {
				add(ans.a[i][j], a[i][j]);
				add(ans.a[i][j], t.a[i][j]);
			}
		return ans;
	}
	
	Matrix operator * (const Matrix& t) const {
		Matrix ans;
		for(int i = 0; i < n; i++)
			for(int j = 0; j < n; j++) if(a[i][j])
				for(int k = 0; k < n; k++)
					add(ans.a[i][k], mul(a[i][j], t.a[j][k]));
		return ans;
	}

	void output() {
		for(int i = 0; i < n; i++) {
			for(int j = 0; j < n - 1; j++)
				printf("%d ", a[i][j]);
			printf("%d
", a[i][n - 1]);
		}
		printf("
");
	}
};

Matrix pow_mod(Matrix a, int p) {
	Matrix ans;
	ans.E();
	while(p) {
		if(p & 1) ans = ans * a;
		a = a * a;
		p >>= 1;
	}
	return ans;
}

Matrix E;

Matrix sum(Matrix a, int p) {
	if(p == 1) return a;
	Matrix ans;
	ans = (E + pow_mod(a, p / 2)) * sum(a, p / 2);
	if(p & 1) ans = ans + pow_mod(a, p);
	return ans;
}

int main()
{
	E.E();
	while(scanf("%d%d", &n, &k) == 2) {
		if(n == 0 && k == 0) break;
		Matrix a;
		for(int i = 0; i < n; i++)
			for(int j = 0; j < n; j++) {
				scanf("%d", &a.a[i][j]);
				a.a[i][j] %= 10;
			}

		if(k == 0) {
			E.output();
			continue;
		}

		a = sum(a, k);
		a.output();
	}

	return 0;
}