UVA1386 Cellular Automaton

题目链接：https://www.luogu.com.cn/problem/UVA1386

看到题目，每次操作都会影响每个格子，很容易列出操作矩阵。
但是直接矩阵快速幂会T，必须优化才行
观察整个矩阵，发现由于距离这一因素，这个操作矩阵是循环矩阵（从第二行开始每一行都是上一行的循环右移）
可以证明，两个循环矩阵的乘积的矩阵还是循环矩阵

证：令$A,B$为两个循环矩阵，$C=A\ast B$
对于$C_{i,j}$，有$$C_{i,j}=\sum _{k=1}^r{A}_{i,k}\ast B_{k,j}$$
$$C_{i-1,j-1}=\sum _{k=2}^{r+1}{A}_{i-1,k-1}\ast B_{i-1,k-1}$$
由循环矩阵的性质： $A_{i,k}=A_{i-1,k-1}，B$同理
可得：$C_{i,j}=C_{i-1,j-1}$（每次累加可以一一对应）
证毕.

由此，我们每次只需要把第一行乘出来即可，单次乘法时间复杂度变为了$O(n^2)$，再矩阵快速幂即可

$Code$

#include <iostream>
#include <cstdlib>
#include <cstdio>
#include <cstring>
using namespace std;
#define int long long
const int MAXN = 500;
int Mod, a[MAXN + 10], g[MAXN + 10][MAXN + 10], f[MAXN + 10][MAXN + 10], ans[MAXN + 10];
void fastpow(int, int);
inline int read();
void MUL(int, int);

signed main(){
	//freopen ("std.in","r",stdin);
	//freopen ("std.out","w",stdout);
	int n, m, d, k;
	while (scanf("%lld%lld%lld%lld", &n, &m, &d, &k) == 4){
		memset(g, 0, sizeof(g));
		Mod = m;
		for (register int i = 1; i <= n; ++i)	a[i] = read();
		for (register int i = 1; i <= n; ++i)
			if (i - 1 <= d || n - (i - 1) <= d)	g[1][i] = 1;
		for (register int i = 2; i <= n; ++i){
			g[i][1] = g[i - 1][n];
			for (register int j = 2; j <= n; ++j)
				g[i][j] = g[i - 1][j - 1];
		}
		fastpow(n, k);
		for (register int i = 1; i <= n; ++i){
			ans[i] = 0;
			for (register int j = 1; j <= n; ++j)
				ans[i] = (ans[i] + f[i][j] * a[j]) % Mod;
			}
		for (register int i = 1; i < n; ++i) cout << ans[i] << " ";
		cout << ans[n] << endl;
	}
	return 0;
}


inline int read(){
	int x = 0;
	char c = getchar();
	while (!isdigit(c)) c = getchar();
	while (isdigit(c))x = (x << 1) + (x << 3) + (c & 15), c = getchar();
	return x;
}

void fastpow(int n, int p){
	for (register int i = 1; i <= n; ++i)
		for (register int j = 1; j <= n; ++j)
			f[i][j] = g[i][j];
	--p;
	while (p){
		if (p & 1)	MUL(0, n);
		MUL(1, n);
		p >>= 1;
	}
}

void MUL(int op, int n){
	static int b[MAXN + 10][MAXN + 10], a[MAXN + 10][MAXN + 10], c[MAXN + 10][MAXN + 10];
	for (register int i = 1; i <= n; ++i)
		for (register int j = 1; j <= n; ++j){
			if (op)	a[i][j] = g[i][j];
			else    a[i][j] = f[i][j];
			b[i][j] = g[i][j];	c[i][j] = 0;
		}
	for (register int i = 1; i <= n; ++i)
		for (register int j = 1; j <= n; ++j)
			c[1][i] = (c[1][i] + a[1][j] * b[j][i]) % Mod;
	
	for (register int i = 2; i <= n; ++i){
		c[i][1] = c[i - 1][n];
		for (register int j = 2; j <= n; ++j)
			c[i][j] = c[i - 1][j - 1];
	}
	for (register int i = 1; i <= n; ++i)
		for (register int j = 1; j <= n; ++j){
			if (op)	g[i][j] = c[i][j];
			else    f[i][j] = c[i][j];
		}
}