洛谷P5548 [BJ United Round #3] 押韵

题面

题解

容易发现答案就是 (displaystyle [x^{nd}] left(sum_{i geq 0} [d|i] frac {x^i}{i!} ight)^k)，对其进行单位根反演就是 (displaystyle [x^{nd}] left(frac 1d sum_{j=0}^{d-1} exp omega_d^i x ight)^k)。

(d = 1, 2, 3)

同复读机。

(d = 4)

通过观察可以发现其实是对于所有 (exp (a + omega b)) 求出它的系数，相当于你现在在 ((0, 0))，可以上下左右走，走 (k) 步恰好到 ((a, b)) 的方案数。将坐标系旋转一下可以发现两维就独立了，方案数就是 (displaystyle inom k {frac {k + |a + b|} 2} inom k {frac {k + |a - b|} 2})。

当然也可以像 (d = 6) 一样推。

(d = 6)

此时可以发现任意 (omega^k) 都可以用 (1) 和 (omega) 表示出来。

那么问题也就变成了从 ((0, 0)) 开始走 (k) 步，可以走 ((1, 0), (0, 1), (-1, 1), (-1, 0), (0, -1), (1, -1))，走到 ((a, b)) 的方案数。

考虑到负数的下标不好处理，于是整体加一个向量 ((1, 1))。

令 (F(x, y) = x + x^2 + x^2y + y + y^2 + y^2x)

这样现在的问题就是求出 (G(x, y) = F(x, y)^k) 的系数。

对 (x) 求导，有 (displaystyle frac {partial G} {partial x} = k F^{k - 1} frac {partial F} {partial x})，令两边系数对应相等，有：

[egin{aligned} & ng_{n, m} + (n - 1) g_{n - 1, m} + (n - 1) g_{n - 1, m - 1} + (n + 1) g_{n + 1, m - 1} + (n + 1)g_{n + 1, m - 2} + n g_{n, m - 2}\ =& k(g_{n, m} + 2g_{n - 1, m} + 2g_{n - 1, m - 1} + g_{n, m - 2}) end{aligned} ]

将 (g_{n + 1, m - 1}) 移到一边然后令 (n gets n - 1, m gets m + 1) 可得：

[n g_{n, m} = (k - n + 1) (g_{n - 1, m + 1} + g_{n - 1, m - 1}) + (2k - n + 2) (g_{n - 2, m + 1} + g_{n - 2, m}) - n g_{n, m - 1} ]

递推的边界是 (g_{i, k - i} = inom ki)，(g_{i, 0} = g_{0, i} = inom k{i - k})，这样就做完了。

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
#define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout)

inline int read()
{
	int data = 0, w = 1; char ch = getchar();
	while (ch != '-' && (ch < '0' || ch > '9')) ch = getchar();
	if (ch == '-') w = -1, ch = getchar();
	while (ch >= '0' && ch <= '9') data = data * 10 + (ch ^ 48), ch = getchar();
	return data * w;
}

const int N(2010), Mod(1049874433);
inline int upd(const int &x) { return x + (x >> 31 & Mod); }
int fastpow(int x, int y)
{
	int ans = 1;
	for (; y; y >>= 1, x = 1ll * x * x % Mod)
		if (y & 1) ans = 1ll * ans * x % Mod;
	return ans;
}

int n, K, D, w, pw[10], fac[N], inv[N], dK;
int C(int n, int m) { return 1ll * fac[n] * inv[m] % Mod * inv[n - m] % Mod; }

namespace Subtask1
{
	int a[10], ans, res;
	void calc()
	{
		for (int i = res = 0; i < D; i++) res = (res + 1ll * a[i] * pw[i]) % Mod;
		res = fastpow(res, n);
		for (int i = 0; i < D; i++) res = 1ll * res * inv[a[i]] % Mod;
		ans = (ans + 1ll * res * fac[K]) % Mod;
	}

	void dfs(int x, int k)
	{
		if (x == D - 1) return a[x] = k, calc();
		for (int i = 0; i <= k; i++) a[x] = i, dfs(x + 1, k - i), a[x] = 0;
	}

	void main()
	{
		dfs(0, K);
		printf("%lld
", 1ll * ans * dK % Mod);
	}
}

namespace Subtask2
{
	void main()
	{
		int ans = 0;
		for (int i = 1; i <= K; i++)
			for (int j = 0; i + j <= K; j++)
			{
				if ((i + j + K) & 1) continue;
				int cnt = 1ll * C(K, (K + i + j) >> 1) * C(K, (K + std::abs(i - j)) >> 1) % Mod;
				ans = (ans + 1ll * fastpow((i + 1ll * j * w) % Mod, n) * cnt) % Mod;
			}
		printf("%lld
", 4ll * ans * dK % Mod);
	}
}

namespace Subtask3
{
	int g[N << 1][N << 1], m, ans, rev[N << 1];

	void main()
	{
		m = K << 1, rev[1] = 1;
		for (int i = 2; i <= m; i++) rev[i] = 1ll * (Mod - Mod / i) * rev[Mod % i] % Mod;
		for (int i = 0; i <= K; i++) g[i][K - i] = C(K, i);
		for (int i = K + 1; i <= m; i++) g[i][0] = g[0][i] = C(K, i - K);
		for (int i = 1; i <= m; i++)
			for (int j = std::max(K - i + 1, 1); j <= m; j++)
			{
				g[i][j] = 1ll * (K - i + 1) * (g[i - 1][j + 1] + g[i - 1][j - 1]) % Mod;
				if (i > 1) g[i][j] = (g[i][j] + 1ll * (((K + 1) << 1) - i) * (g[i - 2][j + 1] + g[i - 2][j])) % Mod;
				g[i][j] = (g[i][j] - 1ll * i * g[i][j - 1]) % Mod;
				g[i][j] = 1ll * g[i][j] * rev[i] % Mod;
			}
		for (int i = 0; i <= m; i++) for (int j = 0; j <= m; j++) if (g[i][j])
			ans = (ans + 1ll * g[i][j] * fastpow((i - K + 1ll * (j - K) * w) % Mod, n)) % Mod;
		printf("%lld
", 1ll * upd(ans) * dK % Mod);
	}
}

int main()
{
#ifndef ONLINE_JUDGE
	file(cpp);
#endif
	n = read(), K = read(), D = read(), n = 1ll * n * D % (Mod - 1);
	w = fastpow(7, (Mod - 1) / D), fac[0] = pw[0] = 1;
	for (int i = 1; i < D; i++) pw[i] = 1ll * pw[i - 1] * w % Mod;
	for (int i = 1; i <= K; i++) fac[i] = 1ll * fac[i - 1] * i % Mod;
	inv[K] = fastpow(fac[K], Mod - 2), dK = fastpow(fastpow(D, K), Mod - 2);
	for (int i = K; i; i--) inv[i - 1] = 1ll * inv[i] * i % Mod;
	if (D <= 3) Subtask1::main();
	if (D == 4) Subtask2::main();
	if (D == 6) Subtask3::main();
	return 0;
}