[HNOI2019]白兔之舞

wyy603

于 2020-07-08 12:23:28 发布

阅读量218

点赞数

本文链接：https://blog.csdn.net/wyy603/article/details/107202227

版权

题解

我们枚举 $y$ 坐标序列，并使用组合数计算 $x$ 坐标序列的个数，所以我们要对 $\forall 0\le t<k$ 将

$\sum_{m\bmod k=t} \sum_{a_0=x,a_1,a_2,\cdots,a_m=y} \binom{L}{m} \prod_{i=1}^{m-1} w(a_i,a_{i+1})$

求出，由于 $n$ 很小容易联想到矩阵乘法，令 $G$ 为给定矩阵，则将上式改写为

$[G_{x,y}] \sum_{m\bmod k=t} \binom{L}{m} G^m$

（dbq， $G_{x,y}]$ 这个记号是我瞎掰的）

令 $\omega$ 为 $k$ 次单位根，作用单位根反演：

$\begin{aligned} & [G_{x,y}] \sum_{m\ge 0} \binom{L}{m} G^m \frac{1}{k} \sum_{j=0}^{k-1} \omega^{(m-t)j} \\ &= [G_{x,y}] \frac{1}{k} \sum_{j=0}^{k-1} \omega^{-tj} (G\omega^{j}+I)^L \end{aligned}$

令 $t_j=[G_{x,y}](G\omega^j+1)^L$ ，使用技巧 $ij=\binom{i+j}{2}-\binom{i}{2}-\binom{j}{2}$ ，于是上式化为

$\begin{aligned} & \frac{1}{k} \sum_{j=0}^{k-1} \omega^{-\binom{t+j}{2}+\binom{t}{2}+\binom{j}{2}}t_j \\ &= \frac{1}{k} \omega^{\binom{t}{2}} \sum_{j=0}^{k-1} (t_j\omega^{\binom{j}{2}})\omega^{-\binom{t+j}{2}} \end{aligned}$

将其中一个序列反转，比如令 $f_j=t_j\omega^{\binom{j}{2}}=[G_{x,y}](G\omega^{j}+I)^L\omega^{\binom{j}{2}} (0\le j<k), g_{2k-2-j}=\omega^{-\binom{j}{2}}$ ，则上式将化为卷积：

$\frac{1}{k} \omega^{\binom{t}{2}} \sum_{j=0}^{k-1} f_jg_{2k-2-t-j}$

由于我们只用到了单位根反演中的性质，且 $k ∣ (p - 1)$ ，所以我们可以使用模 $p$ 意义下的原根 $g$ 代替单位根。

接下来我们需要进行模 $p$ 意义下的卷积，使用 MTT 或者三模 NTT 都可以解决。 $O(k\log k)$ 。

代码

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

typedef long long ll;
typedef long double lb;
const int MAXN = 6e5 + 5;
const int MOD = 998244353;
const lb pi = acos(-1);

ll q_pow(ll a, ll b, ll p) {
	ll ret = 1;
	for (; b; a = a * a % p, b >>= 1) if (b & 1) ret = ret * a % p;
	return ret;
}
ll inv(ll x, ll p) { return q_pow(x, p - 2, p); }
inline void ADD(int& x, int y, int p) { x += y; if (x >= p) x -= p; }

int getGen(int p) {
	int g;
	for (g = 2; g <= p; ++g) {
		int tmp = p - 1; bool flag = 1;
		for (int i = 2; i * i <= tmp; ++i) if (tmp % i == 0) {
			while (tmp % i == 0) tmp /= i;
			if (q_pow(g, (p - 1) / i, p) == 1) { flag = 0; break; }
		}
		if (flag) break;
	}
	return g;
}

int N, K, L, X, Y, P, Omega;
struct matrix {
	int a[3][3];
	matrix () { for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) a[i][j] = 0; }
	matrix operator * (int k) {
		matrix C;
		for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) C.a[i][j] = 1ll * a[i][j] * k % P;
		return C;
	}
	friend matrix operator* (const matrix &A, const matrix& B) {
		matrix C;
		for (int i = 0; i < N; ++i) for (int k = 0; k < N; ++k) for (int j = 0; j < N; ++j)
			ADD(C.a[i][j], 1ll * A.a[i][k] * B.a[k][j] % P, P);
		return C;
	}
	matrix operator ^ (int k) {
		matrix RET, A = *this;
		for (int i = 0; i < N; ++i) RET.a[i][i] = 1;
		for (; k; A = A * A, k >>= 1) if (k & 1) RET = RET * A;
		return RET;
	}
} mat;

struct comp {
	lb x, y;
	comp (lb x = 0, lb y = 0) : x(x), y(y) {}
	friend comp operator+ (const comp& l, const comp& r) { return comp(l.x + r.x, l.y + r.y); }
	friend comp operator- (const comp& l, const comp& r) { return comp(l.x - r.x, l.y - r.y); }
	friend comp operator* (const comp& l, const comp& r) { return comp(l.x * r.x - l.y * r.y, l.x * r.y + l.y * r.x); }
};

int rev[MAXN];
void FFT(comp* A, int LIM, int op) {
	for (int i = 0; i < LIM; ++i) if (rev[i] < i) swap(A[rev[i]], A[i]);
	for (int l = 2; l <= LIM; l <<= 1) {
		comp g(cos(2.0 * pi / l), op * sin(2.0 * pi / l));
		for (int i = 0; i < LIM; i += l) {
			comp wn(1, 0);
			for (int j = i; j < i + (l >> 1); ++j, wn = wn * g) {
				comp x = A[j], y = wn * A[j + (l >> 1)];
				A[j] = x + y, A[j + (l >> 1)] = x - y;
			}
		}
	}
	if (op == -1) {
		for (int i = 0; i < LIM; ++i) A[i].x /= LIM, A[i].y /= LIM;
	}
}
comp A0[MAXN], A1[MAXN], B0[MAXN], B1[MAXN];
void MTT(int* A, int* B, int DEG) {
	int LIM = 1, L = 0;
	while (LIM <= DEG * 2) LIM <<= 1, ++L;
	for (int i = 0; i < LIM; ++i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (L - 1));
	for (int i = 0; i < LIM; ++i) {
		A0[i] = A[i] >> 15, A1[i] = A[i] & ((1 << 15) - 1);
		B0[i] = B[i] >> 15, B1[i] = B[i] & ((1 << 15) - 1);
	}
	FFT(A0, LIM, 1), FFT(A1, LIM, 1), FFT(B0, LIM, 1), FFT(B1, LIM, 1);
	for (int i = 0; i < LIM; ++i) {
		comp a0 = A0[i], a1 = A1[i], b0 = B0[i], b1 = B1[i];
		A0[i] = a0 * b0, A1[i] = a0 * b1, B0[i] = a1 * b0, B1[i] = a1 * b1;
	}
	FFT(A0, LIM, -1), FFT(A1, LIM, -1), FFT(B0, LIM, -1), FFT(B1, LIM, -1);
	for (int i = 0; i < LIM; ++i) {
		ll t1 = (ll)(A0[i].x + 0.5) % P, t2 = (ll)(A1[i].x + 0.5) % P;
		ll t3 = (ll)(B0[i].x + 0.5) % P, t4 = (ll)(B1[i].x + 0.5) % P;
		A[i] = (((t1 << 30) + ((t2 + t3) << 15) + t4) % P + P) % P;
	}
}

int F[MAXN], G[MAXN];
int main() {
	scanf("%d%d%d%d%d%d", &N, &K, &L, &X, &Y, &P), --X, --Y;
	int Gp = getGen(P); Omega = q_pow(Gp, (P - 1) / K, P);
	for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) scanf("%d", &mat.a[i][j]);
	for (int i = 0; i < K; ++i) {
		matrix tmp = mat; int pw = q_pow(Omega, i, P);
		for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) tmp.a[i][j] = 1ll * tmp.a[i][j] * pw % P;
		for (int i = 0; i < N; ++i) ADD(tmp.a[i][i], 1, P);
		tmp = tmp ^ L;
		F[i] = 1ll * tmp.a[X][Y] * q_pow(Omega, 1ll * i * (i - 1) / 2, P) % P;
	}
	for (int i = 0; i <= 2 * K - 2; ++i) G[2 * K - 2 - i] = q_pow(inv(Omega, P), 1ll * i * (i - 1) / 2, P);
	MTT(F, G, 2 * K);
	int invK = inv(K, P);
	for (int i = 0; i < K; ++i) {
		printf("%d\n", 1ll * invK * q_pow(Omega, 1ll * i * (i - 1) / 2, P) % P * F[2 * K - 2 - i] % P);
	}
	return 0;
}

wyy603

关注

0
点赞
踩
0

收藏

觉得还不错? 一键收藏
0
评论
[HNOI2019]白兔之舞

题解我们枚举 yyy 坐标序列，并使用组合数计算 xxx 坐标序列的个数，所以我们要对 ∀0≤t<k\forall 0\le t<k∀0≤t<k 将∑m mod k=t∑a0=x,a1,a2,⋯ ,am=y(Lm)∏i=1m−1w(ai,ai+1) \sum_{m\bmod k=t} \sum_{a_0=x,a_1,a_2,\cdots,a_m=y} \binom{L}{m} \prod_{i=1}^{m-1} w(a_i,a_{i+1}) mmodk=t∑a0=x,a1,a2
复制链接

扫一扫