乘法 （矩阵快速幂 矩阵套矩阵？）[2020.5.2]

乘法

题目描述

输入一个 $n ∗ n$ 的矩阵 $A$，请求出 $𝐴 + 𝐴^2+… + 𝐴^k$ 对 m 取模的结果。

输入格式

第一行为三个正整数 n, k, m，含义如上所述； 接下来 n 行，每行输入 n 个非负整数，用于描述矩阵 $A$。

输出格式

输出 n 行，每行 n 个整数，表示答案矩阵。

样例输入

2 2 5
2 1
0 3

样例输出

1 1
0 2

数据范围与约定

对于前 30% 的数据，k<=30；
对于另外 30% 的数据，n=1；
对于 100％的数据，n<=$30$, k<=$10^9$, m<=$10^4$，输入矩阵的数在 0 ~ m-1 范围内。

解题思路

考场上直接想到的是用等比数列的求和直接搞
但是后面就会因为分母要除一个式子而要用高斯消元
而且列出来后系数会有很多0，而且因为是整数，消元的话还得求个公倍数
于是就有了代码长度4600的超长代码（空格码风）
结果还没改出来果断放弃改为暴力（泪）

还是来康正解吧
首先我们令 \(f_i = 𝐴 + 𝐴^2+ … + 𝐴^i\), 还有矩阵\(\begin{bmatrix}f_i\\A^i\\\end{bmatrix}\)
注意该矩阵中的 \(f_i\) 和 \(A^i\) 本身都是矩阵
接下来很容易就可以推导出转移矩阵\(\begin{bmatrix}I_n&A\\0&A\end{bmatrix}\)
即\(\begin{bmatrix}I_n&A\\0&A\end{bmatrix} * \begin{bmatrix}f_i\\A^i\\\end{bmatrix} = \begin{bmatrix}f_{i+1}\\A^{i+1}\\\end{bmatrix}\)
其中 \(I_n\) 表示的是一个n阶的单位矩阵
然后跑个快速幂就出来啦！
不可局限于某个特定的思维，将所学进行拓展，感叹收获颇丰

Code:

Click for the code.

#include<bits/stdc++.h>
using namespace std;

inline int read() {//朴实无华的快读
	int res = 0, f = 1;
	char ch = getchar();
	while (ch < '0' || ch > '9') {
		if (ch == '-') f = -1;
		ch = getchar();
	}
	while (ch >= '0' && ch <= '9') {
		res = (res << 3) + (res << 1) + (ch ^ 48);
		ch = getchar(); 
	}
	return res * f;
}

int n, k, m;

struct arr {//常规的数组
	int a[35][35];
	void st() {//设置为单位矩阵
		for (int i = 1; i <= n; i++)
			a[i][i] = 1;
		return ;
	}
	void init() {
		for (int i = 1; i <= n; i++) 
		for (int j = 1; j <= n; j++)
			a[i][j] = read();
	}
	void clear() {//清零
		for (int i = 1; i <= n; i++) 
		for (int j = 1; j <= n; j++) 
			a[i][j] = 0;
	}
	arr operator * (const arr &x) const {
		arr res; res.clear();
//		for (int i = 1; i <= n; i++) 
//			for (int j = 1; j <= n; j++) res.a[i][j] = a[i][j];
		
		for (int i = 1; i <= n; i++) 
			for (int j = 1; j <= n; j++) 
				for (int k = 1; k <= n; k++) 
					res.a[i][j] = (res.a[i][j] + a[i][k] * x.a[k][j]) % m;
		return res;
	}
	arr operator + (const arr &x) const {
		arr res; res.clear();
		for (int i = 1; i <= n; i++) 
			for (int j = 1; j <= n; j++)
				res.a[i][j] = (a[i][j] + x.a[i][j]) % m;
		return res;
	}
};

struct arrplus {//本题需要的矩阵套矩阵
	arr a[3][3];
	void clear() {//清零
		for (int i = 1; i <= 2; i++) 
			for (int j = 1; j <= 2; j++) a[i][j].clear();
	}
	void st() {//设置单位矩阵
		for (int i = 1; i <= 2; i++) a[i][i].st();
	}
	arrplus operator * (const arrplus &x) const {
		arrplus res; res.clear();
		for (int i = 1; i <= 2; i++) 
			for (int j = 1; j <= 2; j++) 
				for (int k = 1; k <= 2; k++) 
					res.a[i][j] = res.a[i][j] + a[i][k] * x.a[k][j];
		return res;
	}
};

void pt(arr x) {
	for (int i = 1; i <= n; i++) {
		for (int j = 1; j <= n; j++)
			printf("%d ", x.a[i][j]);
		printf("\n");
	}	
}

arrplus fast_pow(arrplus x, int ti) {
	arrplus res; res.clear();
	res.st();
	while (ti) {
		if (ti & 1) res = res * x;
		x = x * x;
		ti >>= 1;
	}
	return res;
}

arr A;
arrplus turn, yuan;


int main() {
	// freopen("mul.in", "r", stdin);
	// freopen("mul.out", "w", stdout);
	n = read(); k = read(); m = read();
	A.init();
	turn.a[1][2] = A; turn.a[2][2] = A;
	turn.a[1][1].st(); turn.a[2][1].clear();
	turn = fast_pow(turn, k - 1);
	yuan.a[1][1] = yuan.a[2][1] = A;
	yuan = turn * yuan;
	
	pt(yuan.a[1][1]);
}