【LOJ2325】「清华集训 2017」小 Y 和恐怖的奴隶主

【题目链接】

【思路要点】

  • 由于场上的每一个位置只有可能包含一个0,1,2或3血的奴隶主,因此本质不同的状态数只有\(\binom{11}{3}=165\)种,矩阵乘法+快速幂即可,询问时只需要矩阵乘向量。
  • 时间复杂度\(O(C^3LogN+TC^2LogN)\),其中\(C=165\)。
  • 上述做法可能会因为常数问题被卡到90~95分。
  • 考虑对于已知的\(M\)和\(K\),DP得到\(N≤2*D\)时的答案,然后运行BM算法。
  • 我们发现在题目限制下,答案是一个长度至多为146的常系数齐次线性递推数列。
  • 用Cayley-Hamilton定理优化递推即可。
  • 时间复杂度\(O(DC^2+D^2+TD^2LogN)\),其中\(C=165,D=146\)。

【代码】

#include<bits/stdc++.h>
using namespace std;
const int MAXN = 205;
const int MAXT = 1005;
const int MAXLOG = 62;
const int P = 998244353;
template <typename T> void read(T &x) {
	int f = 1; x = 0;
	char c = getchar();
	for (; !isdigit(c); c = getchar()) if (c == '-') f = -f;
	for (; isdigit(c); c = getchar()) x = x * 10 + c - '0';
	x *= f;
}
template <typename T> void write(T x) {
	if (x < 0) x = -x, putchar('-');
	if (x > 9) write(x / 10);
	putchar(x % 10 + '0');
}
template <typename T> void writeln(T x) {
	write(x);
	puts("");
}
namespace CayleyHamilton {
	const int MAXN = 305;
	int n, k;
	int a[MAXN], now[MAXN];
	int res[MAXT][MAXN], h[MAXN];
	void times(int *x, int *y) {
		static int tmp[MAXN];
		memset(tmp, 0, sizeof(tmp));
		for (int i = 0; i <= k - 1; i++)
		for (int j = 0; j <= k - 1; j++)
			tmp[i + j] = (tmp[i + j] + 1ll * x[i] * y[j]) % P;
		for (int i = 2 * k - 2; i >= k; i--) {
			for (int j = 1; j <= k; j++)
				tmp[i - j] = (tmp[i - j] + 1ll * tmp[i] * a[j]) % P;
			tmp[i] = 0;
		}
		for (int i = 0; i <= k - 1; i++)
			x[i] = tmp[i];
	}
	void init(vector <int> &res, int *val) {
		k = res.size();
		for (int i = 1; i <= k; i++)
			a[i] = res[i - 1];
		for (int i = 1; i <= k * 2; i++)
			h[i] = val[i];
	}
	void query(int cnt, long long *n, long long *ans) {
		for (int i = 1; i <= cnt; i++) {
			if (n[i] <= 2 * k) {
				ans[i] = h[n[i]];
				n[i] = 0;
			} else {
				n[i] -= k;
				res[i][0] = 1;
			}
		}
		now[1] = 1;
		for (int i = 0; i <= MAXLOG; i++) {
			long long tmp = 1ll << i;
			for (int j = 1; j <= cnt; j++)
				if (n[j] & tmp) times(res[j], now);
			times(now, now);
		}
		for (int i = 1; i <= cnt; i++) {
			if (n[i] == 0) continue;
			for (int j = 0; j <= k - 1; j++)
				ans[i] = (ans[i] + 1ll * h[k + j] * res[i][j]) % P;
		}
	}
}
namespace BerlekampMassey {
	const int MAXN = 305;
	int n, val[MAXN], cnt, fail[MAXN], delta[MAXN];
	vector <int> ans[MAXN];
	int power(int x, int y) {
		if (y == 0) return 1;
		int tmp = power(x, y / 2);
		if (y % 2 == 0) return 1ll * tmp * tmp % P;
		else return 1ll * tmp * tmp % P * x % P;
	}
	void work() {
		ans[cnt = 0].clear();
		for (int i = 1; i <= n; i++) {
			int now = val[i];
			for (unsigned j = 0; j < ans[cnt].size(); j++)
				now = (now - 1ll * val[i - j - 1] * ans[cnt][j] % P + P) % P;
			delta[i] = now; if (now == 0) continue;
			fail[cnt] = i;
			if (cnt == 0) {
				ans[++cnt].clear();
				ans[cnt].resize(i);
				continue;
			}
			ans[++cnt].clear();
			ans[cnt].resize(i - fail[cnt - 2] - 1);
			int mul = 1ll * now * power(delta[fail[cnt - 2]], P - 2) % P;
			ans[cnt].push_back(mul);
			for (unsigned j = 0; j < ans[cnt - 2].size(); j++)
				ans[cnt].push_back(1ll * ans[cnt - 2][j] * (P - mul) % P);
			if (ans[cnt].size() < ans[cnt - 1].size()) ans[cnt].resize(ans[cnt - 1].size());
			for (unsigned j = 0; j < ans[cnt - 1].size(); j++)
				ans[cnt][j] = (ans[cnt][j] + ans[cnt - 1][j]) % P;
		}
	}
}
int T, m, k, tot;
int cnt[MAXN][4], sum[MAXN], res[4];
int inv[MAXN], vec[1][MAXN], tmp[1][MAXN];
int matrix[MAXN][MAXN];
long long x[MAXT], ans[MAXT];
int power(int x, int y) {
	if (y == 0) return 1;
	int tmp = power(x, y / 2);
	if (y % 2 == 0) return 1ll * tmp * tmp % P;
	else return 1ll * tmp * tmp % P * x % P;
}
int find(int *res) {
	for (int i = 1; i <= tot; i++) {
		bool flg = true;
		for (int j = 0; j <= 3; j++)
			if (res[j] != cnt[i][j]) flg = false;
		if (flg) return i;
	}
	return -1;
}
int main() {
	read(T), read(m), read(k);
	for (int i = 0; i <= k; i++) {
		if (m == 1) {
			tot++;
			cnt[tot][0] = i;
			cnt[tot][1] = k - i;
			sum[tot] = k + 1 - i;
			continue;
		}
		for (int j = 0; j <= k - i; j++) {
			if (m == 2) {
				tot++;
				cnt[tot][0] = i;
				cnt[tot][1] = j;
				cnt[tot][2] = k - i - j;
				sum[tot] = k + 1 - i;
				continue;
			}
			for (int l = 0; l <= k - i - j; l++) {
				tot++;
				cnt[tot][0] = i;
				cnt[tot][1] = j;
				cnt[tot][2] = l;
				cnt[tot][3] = k - i - j - l;
				sum[tot] = k + 1 - i;
			}
		}
	}
	for (int i = 1; i <= k + 1; i++)
		inv[i] = power(i, P - 2);
	for (int i = 1; i <= tot; i++) {
		matrix[i][tot + 1] = inv[sum[i]];
		matrix[i][i] = inv[sum[i]];
		memcpy(res, cnt[i], sizeof(cnt[i]));
		if (res[1]) {
			int tmp = res[1];
			res[1]--, res[0]++;
			matrix[i][find(res)] = 1ll * inv[sum[i]] * tmp % P;
		}
		memcpy(res, cnt[i], sizeof(cnt[i]));
		if (res[2]) {
			int tmp = res[2];
			res[2]--, res[1]++;
			if (res[0]) res[0]--, res[m]++;
			matrix[i][find(res)] = 1ll * inv[sum[i]] * tmp % P;
		}
		memcpy(res, cnt[i], sizeof(cnt[i]));
		if (res[3]) {
			int tmp = res[3];
			res[3]--, res[2]++;
			if (res[0]) res[0]--, res[m]++;
			matrix[i][find(res)] = 1ll * inv[sum[i]] * tmp % P;
		}
	}
	tot++; matrix[tot][tot] = 1;
	memset(res, 0, sizeof(res));
	res[m] = 1, res[0] = k - 1;
	memset(vec, 0, sizeof(vec));
	vec[0][find(res)] = 1;
	for (int q = 1; q <= 300; q++) {
		memset(tmp, 0, sizeof(tmp));
		for (int i = 0; i <= 0; i++)
		for (int j = 1; j <= tot; j++)
		for (int k = 1; k <= tot; k++)
			tmp[i][j] = (tmp[i][j] + 1ll * vec[i][k] * matrix[k][j]) % P;
		for (int i = 1; i <= tot; i++)
			vec[0][i] = tmp[0][i];
		BerlekampMassey :: val[q] = vec[0][tot];
	}
	BerlekampMassey :: n = 300;
	BerlekampMassey :: work();
	CayleyHamilton :: init(BerlekampMassey :: ans[BerlekampMassey :: cnt], BerlekampMassey :: val);
	for (int i = 1; i <= T; i++)
		read(x[i]);
	CayleyHamilton :: query(T, x, ans);
	for (int i = 1; i <= T; i++)
		writeln(ans[i]);
	return 0;
}

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值