题面

题解

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

\(d = 1, 2, 3\)

同复读机。

\(d = 4\)

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

\[\begin{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} = \binom ki\)，\(g_{i, 0} = g_{0, i} = \binom 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\n", 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\n", 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\n", 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;
}