Luogu LOJ orz laofu！!1 分析 $op = 0$ 只需要统计两棵树有多少条重边即可。假设有 $c$ 条，答案即为 $y^{n - c}$。 $op = 1$ 先特判掉 $y = 1$ 的情况，答案为 $n^{n - 2}$。 我们要求的是 $$ \sum_{T_2} y^{n - |T_1 \cap T_2|} $$ 交集不好处理，我们考虑这样一个式子 $$ f(S) = \sum_{T \subseteq S} \sum_{P \subseteq T} (-1)^{|T| - |P|} f(P) \tag{1} \label{eq1} $$ [collapse title="证明（可能是伪的）" status="false"] 规定 $0 ^ 0 = 1$，可以证明在下面的推导中这是成立的。 $$ \begin{align} \text{右边} = & \sum_{P \subseteq S} f(P) \sum_{j = 0}^{|S| - |P|} (-1)^j {|S| - |P| \choose j} \\ = & \sum_{P \subseteq S} f(P) (1 - 1)^{|S| - |P|} \\ = & f(S) = \text{左边} \end{align} $$ [/collapse] 把 $\eqref{eq1}$ 套进去得到 $$ \sum_{T_2} \sum_{S \subseteq T_1 \cap T_2} \sum_{T \subseteq S} (-1) ^ {|S| - |T|} y^{n - |T|} $$ 设 $g(S)$ 为包含 $S$ 中边的生成树个数，上式可以变为 $$ \begin{align} & \sum_{S \subseteq T_1} \left(\sum_{T \subseteq S} (-1)^{|S| - |T|} y^{n - |T|}\right) g(S) \\ = & \sum_{S \subseteq T_1} y^{n - |S|} \left(\sum_{T \subseteq S} (-y)^{n - |T|}\right) g(S) \\ = & \sum_{S \subseteq T_1} y^{n - |S|} \left(\sum_{i = 0}^{|S|} {|S| \choose i} (-y)^{n - i}\right) g(S) \\ = & \sum_{S \subseteq T_1} y^{n - |S|} (1 - y)^{|S|} g(S) \\ \end{align} $$ 考虑 $g(S)$ 等于什么。我们把 $S$ 中的边加入，那么原图会构成 $k$ 个大小为 $a_1, a_2, \cdots, a_k$ 的连通块，将它们连成一棵树的方案数为 $$ g(S) = n^{k - 2} \prod_{i = 1}^k a_i $$ [collapse title="证明" status="false"] link [/collapse] 于是上式变为 $$ \begin{align} & \sum_{S \subseteq T_1} y^k (1 - y)^{n - k} n^{k - 2} \prod_{i = 1}^k a_i \\ = & \frac{(1 - y)^n}{n^2} \sum_{S \subseteq T_1} \left(\frac{ny}{1 - y}\right)^k \prod_{i = 1}^k a_i \\ = & \frac{(1 - y)^n}{n^2} \sum_{S \subseteq T_1} \prod_{i = 1}^k \frac{ny}{1 - y} a_i\\ \end{align} $$ 不妨设 $k = \frac{ny}{1 - y}$，上式的实际意义相当于每个大小为 $a$ 的连通块产生 $ka$ 的乘积贡献。 这样子就可以想到一个 DP：设 $f_{i, j}$ 表示以 $i$ 为根的子树、$i$ 所在的连通块大小为 $j$ 的答案，转移为树形背包。然而这样子是 $\mathcal{O}(n^2)$ 的。 考虑把 $a_i$ 写成 $a_i \choose 1$，相当于在每个连通块内选一个点产生 $ka_i$ 的乘积贡献。这样子就有一个更优的 DP：设 $f_{i, 0/1}$ 表示以 $i$ 为根的子树、$i$ 所在连通块是否已经选了一个点产生贡献，转移讨论一下即可。这样子即可做到 $\mathcal{O}(n)$。 $op = 2$ 先判掉 $y = 1$ 的情况，答案为 $n^{2(n - 2)}$。 我们要求的是 $$ \sum_{T_1} \sum_{T_2} y^{n - |T_1 \cap T_2|} $$ 同样把 $\eqref{eq1}$ 套进去得到 $$ \sum_{T_1} \sum_{T_2} \sum_{S\subseteq T_1 \cap T_2} \sum_{T\subseteq S} (-1)^{|S| - |T|} y^{n - |T|} $$ 设 $g(S)$ 表示包含 $S$ 中边的生成树个数，上式可以变为（过程和上面类似就不写了） $$ \sum_{S} y^{n - |S|} (1 - y)^{|S|} g(S)^2 $$ 而 $$ g(S)^2 = n^{2k - 4} \prod_{i = 1}^k a_i^2 $$ 代回去得到 $$ \begin{align} & \sum_{S} y^k (1 - y)^{n - k} n^{2k - 4} \prod_{i = 1}^k a_i^2 \\ = & \frac{(1 - y)^n}{n^4}\sum_S\prod_{i = 1}^k \frac{n^2 y}{1 - y} a_i^2 \end{align} $$ 考察每个连通分量。一个大小为 $a$ 的连通分量会产生 $\frac{n^2 y}{1 - y} a^2$ 的乘积贡献，而 $a$ 个点的生成树有 $a^{a-2}$ 棵，所以其总贡献为 $\frac{n^2 y}{1 - y} a^a$。 而整张图是由若干个连通分量组合而成的，这让我们想到 EGF 中的 $\exp$。 于是构造 EGF $F(x) = \sum_{i \geq 1} \frac{n^2 y}{(1-y) i!} i^i x^i$，然后求 $G(x) = \exp(F(x))$，答案即为 $\frac{(1 - y)^n n!}{n^4}[x^n] G(x)$。 总结：这是一道好题，可是我什么都不会。 代码 // ==================================== // author: M_sea // website: https://m-sea-blog.com/ // ==================================== #include <bits/stdc++.h> #define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout) #define debug(...) fprintf(stderr, __VA_ARGS__) typedef long long ll; int read() { int X = 0, w = 1; char c = getchar(); while (c < '0' || c > '9') { if (c == '-') w = -1; c = getchar(); } while (c >= '0' && c <= '9') X = X * 10 + c - '0', c = getchar(); return X * w; } const int mod = 998244353; int qpow(int a, int b) { int c = 1; for (; b; b >>= 1, a = 1ll * a * a % mod) if (b & 1) c = 1ll * c * a % mod; return c; } int n, y; namespace Subtask0 { std::map<std::pair<int,int>, bool> M; void main() { for (int i = 1; i < n; ++i) { int u = read(), v = read(); M[std::make_pair(u, v)] = M[std::make_pair(v, u)] = 1; } int c = n; for (int i = 1; i < n; ++i) { int u = read(), v = read(); if (M.count(std::make_pair(u, v))) --c; } printf("%d\n", qpow(y, c)); } } namespace Subtask1 { const int N = 100000 + 10; std::vector<int> E[N]; int k, f[N][2]; void dfs(int u, int fa) { f[u][0] = 1, f[u][1] = k; for (int v : E[u]) { if (v == fa) continue; dfs(v, u); int f0 = 1ll * f[u][0] * (f[v][0] + f[v][1]) % mod; int f1 = (1ll * f[u][0] * f[v][1] + 1ll * f[u][1] * f[v][0] + 1ll * f[u][1] * f[v][1]) % mod; f[u][0] = f0, f[u][1] = f1; } } void main() { if (y == 1) { printf("%d\n", qpow(n, n - 2)); return; } for (int i = 1; i < n; ++i) { int u = read(), v = read(); E[u].emplace_back(v), E[v].emplace_back(u); } k = 1ll * n * y % mod * qpow(1 - y + mod, mod - 2) % mod; dfs(1, 0); int ans = 1ll * f[1][1] * qpow(1 - y + mod, n) % mod * qpow(n, mod - 3) % mod; printf("%d\n", ans); } } namespace Subtask2 { const int N = 262144 + 10; int r[N]; void NTT(int *A, int n, int op) { for (int i = 0; i < n; ++i) if (i < r[i]) std::swap(A[i], A[r[i]]); for (int i = 1; i < n; i <<= 1) { int rot = qpow(op == 1 ? 3 : 332748118, (mod - 1) / (i << 1)); for (int j = 0; j < n; j += i << 1) for (int k = 0, w = 1; k < i; ++k, w = 1ll * w * rot % mod) { int x = A[j + k], y = 1ll * w * A[j + k + i] % mod; A[j + k] = (x + y) % mod, A[j + k + i] = (x - y + mod) % mod; } } if (op == -1) { int inv = qpow(n, mod - 2); for (int i = 0; i < n; ++i) A[i] = 1ll * A[i] * inv % mod; } } int NTT_init(int n) { int lim = 1, l = 0; for (; lim < n; lim <<= 1, ++l); for (int i = 0; i < lim; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); return lim; } void PolyInv(int *F, int *G, int n) { static int A[N], B[N]; if (n == 1) { G[0] = qpow(F[0], mod - 2); return; } PolyInv(F, G, n >> 1); for (int i = 0; i < n; ++i) A[i] = F[i], B[i] = G[i]; int lim = NTT_init(n << 1); NTT(A, lim, 1), NTT(B, lim, 1); for (int i = 0; i < lim; ++i) A[i] = 1ll * A[i] * B[i] % mod * B[i] % mod; NTT(A, lim, -1); for (int i = 0; i < n; ++i) G[i] = (2ll * G[i] - A[i] + mod) % mod; for (int i = 0; i < lim; ++i) A[i] = B[i] = 0; } void PolyDeri(int *F, int *G, int n) { for (int i = 1; i < n; ++i) G[i - 1] = 1ll * F[i] * i %mod; G[n - 1] = 0; } void PolyInte(int *F, int *G, int n) { for (int i = 1; i < n; ++i) G[i] = 1ll * F[i-1] * qpow(i, mod - 2) % mod; G[0] = 0; } void PolyLn(int *F, int *G, int n) { static int A[N], B[N]; PolyDeri(F, A, n), PolyInv(F, B, n); int lim = NTT_init(n << 1); NTT(A, lim, 1), NTT(B, lim, 1); for (int i = 0; i < lim; ++i) A[i] = 1ll * A[i] * B[i] % mod; NTT(A, lim, -1); PolyInte(A, G, n); for (int i = 0; i < lim; ++i) A[i] = B[i] = 0; } void PolyExp(int *F, int *G, int n) { static int A[N], B[N]; if (n == 1) { G[0] = 1; return; } PolyExp(F, G, n >> 1); for (int i = 0; i < n; ++i) A[i] = G[i]; PolyLn(G, B, n); for (int i = 0; i < n; ++i) B[i] = (mod - B[i] + F[i]) % mod; B[0] = (B[0] + 1) % mod; int lim = NTT_init(n << 1); NTT(A, lim, 1), NTT(B, lim, 1); for (int i = 0; i < lim; ++i) A[i] = 1ll * A[i] * B[i] % mod; NTT(A, lim, -1); for (int i = 0; i < n; ++i) G[i] = A[i]; for (int i = 0; i < lim; ++i) A[i] = B[i] = 0; } int fac[N], ifac[N]; int F[N], G[N]; void main() { if (y == 1) { printf("%d\n", qpow(n, 2 * (n - 2))); return; } fac[0] = 1; for (int i = 1; i <= n; ++i) fac[i] = 1ll * fac[i - 1] * i % mod; ifac[n] = qpow(fac[n], mod - 2); for (int i = n; i; --i) ifac[i - 1] = 1ll * ifac[i] * i % mod; int k = 1ll * n * n % mod * y % mod * qpow(1 - y + mod, mod - 2) % mod; for (int i = 1; i <= n; ++i) F[i] = 1ll * k * qpow(i, i) % mod * ifac[i] % mod; int lim = 1; for (; lim <= n; lim <<= 1); PolyExp(F, G, lim); int ans = 1ll * G[n] * fac[n] % mod * qpow(1 - y + mod, n) % mod * qpow(n, mod - 5) % mod; printf("%d\n", ans); } } int main() { n = read(), y = read(); int op = read(); switch (op) { case 0: Subtask0::main(); break; case 1: Subtask1::main(); break; case 2: Subtask2::main(); break; } return 0; }
Luogu 分析 为了方便我们先令 $y$ 等于 $x + 1$，则 $S_k(x) = \sum_{i = 0}^{y - 1} i ^ k$。 根据伯努利数的结论有 $$ S_k(y) = \frac{1}{k + 1} \sum_{i = 0}^k {k + 1 \choose i} B_i y^{k + 1 - i} $$ 则答案为 $$ \sum_{k = 0}^n \frac{a_k}{k + 1} \sum_{i = 0}^k {k + 1 \choose i} B_i y ^ {k + 1 - i} $$ 稍作变形得到 $$ \sum_{k = 0}^n \frac{a_k}{k + 1} \sum_{i = 0}^k {k + 1 \choose i + 1} B_{k - i} y ^ {i + 1} $$ 交换求和号并拆掉组合数 $$ \sum_{i = 0}^k \frac{y ^ {i + 1}}{(i + 1)!} \sum_{k = i}^n a_k k! \frac{B_{k - i}}{(k - i)!} $$ 后面是差卷积的形式，翻转后一遍 NTT 即可。 然而 $y = x + 1$，考虑把求出的多项式改写为关于 $x$ 的形式。这相当于已知 $\sum_{i \geq 0} f_i x ^ i$，求 $\sum_{i \geq 0} f_i (x + 1) ^ i$。 用二项式定理拆开 $$ \sum_{i = 0}^{n + 1} f_i \sum_{j = 0}^i {i \choose j} x^j $$ 交换求和号并拆掉组合数 $$ \sum_{j = 0}^{n + 1} \frac{x ^ j}{j!} \sum_{i = j}^n i! f_i \frac{1}{(i - j)!} $$ 后面同样是差卷积的形式，翻转后一遍 NTT 即可。 代码 // ==================================== // author: M_sea // website: https://m-sea-blog.com/ // ==================================== #include <bits/stdc++.h> #define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout) #define debug(...) fprintf(stderr, __VA_ARGS__) using namespace std; typedef long long ll; int read() { int X = 0, w = 1; char c = getchar(); while (c < '0' || c > '9') { if (c == '-') w = -1; c = getchar(); } while (c >= '0' && c <= '9') X = X * 10 + c - '0', c = getchar(); return X * w; } const int N = 524288 + 10; const int mod = 998244353; int qpow(int a, int b) { int c = 1; for (; b; b >>= 1, a = 1ll * a * a % mod) if (b & 1) c = 1ll * c * a % mod; return c; } int r[N]; void NTT(int *A, int n, int op) { for (int i = 0; i < n; ++i) if (i < r[i]) swap(A[i], A[r[i]]); for (int i = 1; i < n; i <<= 1) { int rot = qpow(op == 1 ? 3 : 332748118, (mod - 1) / (i << 1)); for (int j = 0; j < n; j += i << 1) for (int k = 0, w = 1; k < i; ++k, w = 1ll * w * rot % mod) { int x = A[j + k], y = 1ll * w * A[j + k + i] % mod; A[j + k] = (x + y) % mod, A[j + k + i] = (x - y + mod) % mod; } } if (op == -1) { int inv = qpow(n, mod - 2); for (int i = 0; i < n; ++i) A[i] = 1ll * A[i] * inv % mod; } } int NTT_init(int n) { int lim = 1, l = 0; for (; lim < n; lim <<= 1, ++l); for (int i = 0; i < lim; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); return lim; } void PolyInv(int *F, int *G, int n) { static int A[N], B[N]; if (n == 1) { G[0] = qpow(F[0], mod - 2); return; } PolyInv(F, G, n >> 1); for (int i = 0; i < n; ++i) A[i] = F[i], B[i] = G[i]; int lim = NTT_init(n << 1); NTT(A, lim, 1), NTT(B, lim, 1); for (int i = 0; i < lim; ++i) A[i] = 1ll * A[i] * B[i] % mod * B[i] % mod; NTT(A, lim, -1); for (int i = 0; i < n; ++i) G[i] = (2ll * G[i] - A[i] + mod) % mod; for (int i = 0; i < lim; ++i) A[i] = B[i] = 0; } int n, a[N]; int fac[N], ifac[N]; int B[N], F[N], G[N]; void init(int n) { fac[0] = 1; for (int i = 1; i <= n; ++i) fac[i] = 1ll * fac[i - 1] * i % mod; ifac[n] = qpow(fac[n], mod - 2); for (int i = n; i; --i) ifac[i - 1] = 1ll * ifac[i] * i % mod; } int main() { n = read(); init(n + 1); for (int i = 0; i <= n; ++i) a[i] = read(); for (int i = 0; i <= n; ++i) F[i] = ifac[i + 1]; int lim = 1; for (; lim <= n; lim <<= 1); PolyInv(F, B, lim); for (int i = 0; i <= n; ++i) F[i] = 1ll * a[i] * fac[i] % mod; for (int i = 0; i <= n; ++i) G[i] = B[i]; reverse(G, G + n + 1); lim = NTT_init(n << 1 | 1); NTT(F, lim, 1), NTT(G, lim, 1); for (int i = 0; i < lim; ++i) G[i] = 1ll * F[i] * G[i] % mod; NTT(G, lim, -1); for (int i = 0; i < lim; ++i) F[i] = 0; for (int i = 0; i <= n; ++i) F[i + 1] = 1ll * ifac[i + 1] * G[n + i] % mod; lim = NTT_init((n + 1) << 1 | 1); for (int i = 0; i < lim; ++i) G[i] = 0; for (int i = 0; i <= n + 1; ++i) F[i] = 1ll * F[i] * fac[i] % mod; for (int i = 0; i <= n + 1; ++i) G[i] = ifac[i]; reverse(G, G + n + 2); NTT(F, lim, 1), NTT(G, lim, 1); for (int i = 0; i < lim; ++i) F[i] = 1ll * F[i] * G[i] % mod; NTT(F, lim, -1); for (int i = 0; i <= n + 1; ++i) F[i] = 1ll * F[n + 1 + i] * ifac[i] % mod; for (int i = 0; i <= n + 1; ++i) printf("%d%c", F[i], " \n"[i == n + 1]); return 0; }
Luogu 分析 首先求出烷基的 OGF $F(x)$，具体求解过程可以看这里。 考虑一个烯烃断掉双键后会构成一棵根的儿子数 $\leq 2$、其余节点的儿子数 $\leq 3$ 的树，我们考虑求出其 OGF $G(x)$。 根据 Burnside 引理可以得到 $$ G(x) = 1 + x \frac{F(x) ^ 2 + F(x ^ 2)}{2} $$ 然后求出烯烃的 OGF $H(x)$。同样根据 Burnside 引理可以得到 $$ H(x) = \frac{(G(x) - 1) ^ 2 + G(x ^ 2) - 1}{2} $$ 这里减 $1$ 是因为双键两端都必须有碳原子（双键显然不能连在氢原子上）。 依据上式直接计算即可。时间复杂度 $\mathcal{O}(n \log n)$。 代码 // ==================================== // author: M_sea // website: https://m-sea-blog.com/ // ==================================== #include <bits/stdc++.h> #define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout) #define debug(...) fprintf(stderr, __VA_ARGS__) using namespace std; typedef long long ll; int read() { int X = 0, w = 1; char c = getchar(); while (c < '0' || c > '9') { if (c == '-') w = -1; c = getchar(); } while (c >= '0' && c <= '9') X = X * 10 + c - '0', c = getchar(); return X * w; } const int N = 262144 + 10; const int mod = 998244353, inv2 = 499122177; int qpow(int a, int b) { int c = 1; for (; b; b >>= 1, a = 1ll * a * a % mod) if (b & 1) c = 1ll * c * a % mod; return c; } namespace Poly { int r[N]; void NTT(int *A, int n, int op) { for (int i = 0; i < n; ++i) if (i < r[i]) swap(A[i], A[r[i]]); for (int i = 1; i < n; i <<= 1) { int rot = qpow(op == 1 ? 3 : 332748118, (mod - 1) / (i << 1)); for (int j = 0; j < n; j += i << 1) for (int k = 0, w = 1; k < i; ++k, w = 1ll * w * rot % mod) { int x = A[j + k], y = 1ll * w * A[j + k + i] % mod; A[j + k] = (x + y) % mod, A[j + k + i] = (x - y + mod) % mod; } } if (op == -1) { int inv = qpow(n, mod - 2); for (int i = 0; i < n; ++i) A[i] = 1ll * A[i] * inv % mod; } } int NTT_init(int n) { int lim = 1, l = 0; for (; lim < n; lim <<= 1, ++l); for (int i = 0; i < lim; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); return lim; } void PolyInv(int *F, int *G, int n) { static int A[N], B[N]; if (n == 1) { G[0] = qpow(F[0], mod - 2); return; } PolyInv(F, G, n >> 1); for (int i = 0; i < n; ++i) A[i] = F[i], B[i] = G[i]; int lim = NTT_init(n << 1); NTT(A, lim, 1), NTT(B, lim, 1); for (int i = 0; i < lim; ++i) A[i] = 1ll * A[i] * B[i] % mod * B[i] % mod; NTT(A, lim, -1); for (int i = 0; i < n; ++i) G[i] = (2ll * G[i] - A[i] + mod) % mod; for (int i = 0; i < lim; ++i) A[i] = B[i] = 0; } } void Newton(int *F, int n) { static int F0[N], A[N], B[N], X[N], Y[N], iY[N]; F[0] = 1; for (int lim = 2; lim <= n; lim <<= 1) { for (int i = 0; i < lim; ++i) F0[i] = F[i]; for (int i = 0; i < lim; i += 2) A[i] = F[i / 2]; for (int i = 0; i < lim; i += 3) B[i] = F[i / 3]; Poly::NTT_init(lim << 1); Poly::NTT(F0, lim << 1, 1), Poly::NTT(A, lim << 1, 1); for (int i = 0; i < lim << 1; ++i) X[i] = (1ll * F0[i] * F0[i] % mod * F0[i] + 3ll * A[i] * F0[i]) % mod; Poly::NTT(X, lim << 1, -1), Poly::NTT(F0, lim << 1, -1), Poly::NTT(A, lim << 1, -1); for (int i = 0; i < lim; ++i) X[i] = (X[i] + 2ll * B[i]) % mod; for (int i = lim - 1; i; --i) X[i] = X[i - 1]; X[0] = 6; for (int i = 0; i < lim; ++i) X[i] = (X[i] + 1ll * (mod - 6) * F0[i]) % mod; Poly::NTT(F0, lim << 1, 1); for (int i = 0; i < lim << 1; ++i) Y[i] = 3ll * F0[i] * F0[i] % mod; Poly::NTT(Y, lim << 1, -1), Poly::NTT(F0, lim << 1, -1); for (int i = 0; i < lim; ++i) Y[i] = (Y[i] + 3ll * A[i]) % mod; for (int i = lim - 1; i ; --i) Y[i] = Y[i - 1]; Y[0] = mod - 6; for (int i = lim; i < lim << 1; ++i) X[i] = Y[i] = 0; Poly::PolyInv(Y, iY, lim); Poly::NTT_init(lim << 1); Poly::NTT(X, lim << 1, 1), Poly::NTT(iY, lim << 1, 1); for (int i = 0; i < lim << 1; ++i) X[i] = 1ll * X[i] * iY[i] % mod; Poly::NTT(X, lim << 1, -1); for (int i = 0; i < lim; ++i) F[i] = (F0[i] - X[i] + mod) % mod; for (int i = 0; i < lim << 1; ++i) F0[i] = A[i] = B[i] = X[i] = Y[i] = iY[i] = 0; } } int F[N], G[N], H[N]; int main() { int n = read(), lim = 1; for (; lim <= n; lim <<= 1); Newton(F, lim); Poly::NTT_init(lim << 1); Poly::NTT(F, lim << 1, 1); for (int i = 0; i < lim << 1; ++i) G[i] = 1ll * F[i] * F[i] % mod; Poly::NTT(G, lim << 1, -1), Poly::NTT(F, lim << 1, -1); for (int i = 0; i < lim; i += 2) G[i] = (G[i] + F[i / 2]) % mod; for (int i = 0; i < lim; ++i) G[i] = 1ll * G[i] * inv2 % mod; for (int i = lim - 1; i; --i) G[i] = G[i - 1]; G[0] = 0; for (int i = lim; i < lim << 1; ++i) G[i] = 0; Poly::NTT(G, lim << 1, 1); for (int i = 0; i < lim << 1; ++i) H[i] = 1ll * G[i] * G[i] % mod; Poly::NTT(H, lim << 1, -1), Poly::NTT(G, lim << 1, -1); for (int i = 0; i < lim; i += 2) H[i] = (H[i] + G[i / 2]) % mod; for (int i = 0; i < lim; ++i) H[i] = 1ll * H[i] * inv2 % mod; for (int i = 2; i <= n; ++i) printf("%d\n", H[i]); return 0; }
LOJ 分析 要求的就是每个节点的儿子数都不超过 $3$ 的无标号有根树数。 设其 OGF 为 $F(x)$，考虑列出一个关于 $F(x)$ 的方程。 不妨规定 $[x^0]F(x) = 1$，那么可以看成每个节点的儿子数都恰好为 $3$。 直接得出 $F(x) = 1 + F(x) ^ 3$ 显然是错误的，因为我们同构的情况被重复计算了。 因此可以考虑 Burnside 引理。对于 $(1, 2, 3)$，所有情况都是不动点；对于 $(1, 3, 2)$、$(3, 2, 1)$、$(2, 1, 3)$，在置换中成环的两棵子树同构的情况是不动点；对于 $(3, 1, 2)$、$(2, 3, 1)$，三棵子树都同构的情况是不动点。 这样子我们就可以列出方程 $$ F(x) = 1 + x \frac{F(x) ^ 3 + 3 F(x ^ 2) F(x) + 2 F(x ^ 3)}{6} $$ 考虑使用牛顿迭代解出 $F(x)$。设 $G(F(x)) = 1 + x \frac{F(x) ^ 3 + 3 F(x ^ 2) F(x) + 2 F(x ^ 3)}{6} - F(x) = 0$。 假设当前已知 $F_0(x) = F(x) \bmod x ^ {\frac{n}{2}}$，要求 $F(x) \bmod x ^ n$。那么事实上 $F(x ^ 2) \bmod{x ^ n}$ 和 $F(x ^ 3) \bmod x ^ n$ 已知，因此可以将其看作常量 $A(x)$ 和 $B(x)$。这样子有 $$ G'(F(x)) = x \frac{F(x) ^ 2 + 3 A(x)}{6} - 1 $$ 套牛顿迭代的式子 $$ F(x) = F_0(x) - \frac{G(F_0(x))}{G'(F_0(x))} $$ 这样子就可以求出 $F(x)$ 了。时间复杂度 $\mathcal{O}(n\log n)$。 代码 // ==================================== // author: M_sea // website: https://m-sea-blog.com/ // ==================================== #include <bits/stdc++.h> #define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout) #define debug(...) fprintf(stderr, __VA_ARGS__) using namespace std; typedef long long ll; int read() { int X = 0, w = 1; char c = getchar(); while (c < '0' || c > '9') { if (c == '-') w = -1; c = getchar(); } while (c >= '0' && c <= '9') X = X * 10 + c - '0', c = getchar(); return X * w; } const int N = 262144 + 10; const int mod = 998244353, inv2 = 499122177; int qpow(int a, int b) { int c = 1; for (; b; b >>= 1, a = 1ll * a * a % mod) if (b & 1) c = 1ll * c * a % mod; return c; } namespace Poly { int r[N]; void NTT(int *A, int n, int op) { for (int i = 0; i < n; ++i) if (i < r[i]) swap(A[i], A[r[i]]); for (int i = 1; i < n; i <<= 1) { int rot = qpow(op == 1 ? 3 : 332748118, (mod - 1) / (i << 1)); for (int j = 0; j < n; j += i << 1) for (int k = 0, w = 1; k < i; ++k, w = 1ll * w * rot % mod) { int x = A[j + k], y = 1ll * w * A[j + k + i] % mod; A[j + k] = (x + y) % mod, A[j + k + i] = (x - y + mod) % mod; } } if (op == -1) { int inv = qpow(n, mod - 2); for (int i = 0; i < n; ++i) A[i] = 1ll * A[i] * inv % mod; } } int NTT_init(int n) { int lim = 1, l = 0; for (; lim < n; lim <<= 1, ++l); for (int i = 0; i < lim; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); return lim; } void PolyInv(int *F, int *G, int n) { static int A[N], B[N]; if (n == 1) { G[0] = qpow(F[0], mod - 2); return; } PolyInv(F, G, n >> 1); for (int i = 0; i < n; ++i) A[i] = F[i], B[i] = G[i]; int lim = NTT_init(n << 1); NTT(A, lim, 1), NTT(B, lim, 1); for (int i = 0; i < lim; ++i) A[i] = 1ll * A[i] * B[i] % mod * B[i] % mod; NTT(A, lim, -1); for (int i = 0; i < n; ++i) G[i] = (2ll * G[i] - A[i] + mod) % mod; for (int i = 0; i < lim; ++i) A[i] = B[i] = 0; } } void Newton(int *F, int n) { static int F0[N], A[N], B[N], X[N], Y[N], iY[N]; F[0] = 1; for (int lim = 2; lim <= n; lim <<= 1) { for (int i = 0; i < lim; ++i) F0[i] = F[i]; for (int i = 0; i < lim; i += 2) A[i] = F[i / 2]; for (int i = 0; i < lim; i += 3) B[i] = F[i / 3]; Poly::NTT_init(lim << 1); Poly::NTT(F0, lim << 1, 1), Poly::NTT(A, lim << 1, 1); for (int i = 0; i < lim << 1; ++i) X[i] = (1ll * F0[i] * F0[i] % mod * F0[i] + 3ll * A[i] * F0[i]) % mod; Poly::NTT(X, lim << 1, -1), Poly::NTT(F0, lim << 1, -1), Poly::NTT(A, lim << 1, -1); for (int i = 0; i < lim; ++i) X[i] = (X[i] + 2ll * B[i]) % mod; for (int i = lim - 1; i; --i) X[i] = X[i - 1]; X[0] = 6; for (int i = 0; i < lim; ++i) X[i] = (X[i] + 1ll * (mod - 6) * F0[i]) % mod; Poly::NTT(F0, lim << 1, 1); for (int i = 0; i < lim << 1; ++i) Y[i] = 3ll * F0[i] * F0[i] % mod; Poly::NTT(Y, lim << 1, -1), Poly::NTT(F0, lim << 1, -1); for (int i = 0; i < lim; ++i) Y[i] = (Y[i] + 3ll * A[i]) % mod; for (int i = lim - 1; i ; --i) Y[i] = Y[i - 1]; Y[0] = mod - 6; for (int i = lim; i < lim << 1; ++i) X[i] = Y[i] = 0; Poly::PolyInv(Y, iY, lim); Poly::NTT_init(lim << 1); Poly::NTT(X, lim << 1, 1), Poly::NTT(iY, lim << 1, 1); for (int i = 0; i < lim << 1; ++i) X[i] = 1ll * X[i] * iY[i] % mod; Poly::NTT(X, lim << 1, -1); for (int i = 0; i < lim; ++i) F[i] = (F0[i] - X[i] + mod) % mod; for (int i = 0; i < lim << 1; ++i) F0[i] = A[i] = B[i] = X[i] = Y[i] = iY[i] = 0; } } int F[N], G[N], H[N]; int main() { int n = read(), lim = 1; for (; lim <= n; lim <<= 1); Newton(F, lim); Poly::NTT_init(lim << 1); Poly::NTT(F, lim << 1, 1); for (int i = 0; i < lim << 1; ++i) G[i] = 1ll * F[i] * F[i] % mod; Poly::NTT(G, lim << 1, -1), Poly::NTT(F, lim << 1, -1); for (int i = 0; i < lim; i += 2) G[i] = (G[i] + F[i / 2]) % mod; for (int i = 0; i < lim; ++i) G[i] = 1ll * G[i] * inv2 % mod; for (int i = lim - 1; i; --i) G[i] = G[i - 1]; G[0] = 0; for (int i = lim; i < lim << 1; ++i) G[i] = 0; Poly::NTT(G, lim << 1, 1); for (int i = 0; i < lim << 1; ++i) H[i] = 1ll * G[i] * G[i] % mod; Poly::NTT(H, lim << 1, -1), Poly::NTT(G, lim << 1, -1); for (int i = 0; i < lim; i += 2) H[i] = (H[i] + G[i / 2]) % mod; for (int i = 0; i < lim; ++i) H[i] = 1ll * H[i] * inv2 % mod; for (int i = 2; i <= n; ++i) printf("%d\n", H[i]); return 0; }
darkbzoj 分析 设答案的 OGF 为 $F(x)$，则有 $$ F(x) = x + \sum_{i \in D} F(x) ^i $$ 即 $$ x = F(x) - \sum_{i \in D} F(x)^i $$ 令 $G(x) = F^{-1}(x)$，可知 $$ G(x) = x - \sum_{i \in D} x^i $$ 由拉格朗日反演可得 $$ [x^n]F(x) = \frac{1}{n} [x^{n - 1}] \left(\frac{x}{G(x)}\right) ^ n $$ 注意到 $[x^0] G(x) = 0$，所以我们把上下同时除以 $x$，上式变为 $$ [x^n]F(x) = \frac{1}{n} [x^{n - 1}] \left(\frac{1}{\frac{G(x)}{x}}\right) ^ n $$ 多项式求逆 + 多项式快速幂计算即可。 代码 // ==================================== // author: M_sea // website: https://m-sea-blog.com/ // ==================================== #include <bits/stdc++.h> #define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout) #define debug(...) fprintf(stderr, __VA_ARGS__) using namespace std; typedef long long ll; int read() { int X = 0, w = 1; char c = getchar(); while (c < '0' || c > '9') { if (c == '-') w = -1; c = getchar(); } while (c >= '0' && c <= '9') X = X * 10 + c - '0', c = getchar(); return X * w; } const int N = 262144 + 10; const int mod = 950009857; int qpow(int a, int b) { int c = 1; for (; b; b >>= 1, a = 1ll * a * a % mod) if (b & 1) c = 1ll * c * a % mod; return c; } int r[N]; void NTT(int *A, int n, int op) { for (int i = 0; i < n; ++i) if (i < r[i]) swap(A[i], A[r[i]]); for (int i = 1; i < n; i <<= 1) { int rot = qpow(op == 1 ? 7 : 135715694, (mod - 1) / (i << 1)); for (int j = 0; j < n; j += i << 1) for (int k = 0, w = 1; k < i; ++k, w = 1ll * w * rot % mod) { int x = A[j + k], y = 1ll * w * A[j + k + i] % mod; A[j + k] = (x + y) % mod, A[j + k + i] = (x - y + mod) % mod; } } if (op == -1) { int inv = qpow(n, mod - 2); for (int i = 0; i < n; ++i) A[i] = 1ll * A[i] * inv % mod; } } int NTT_init(int n) { int lim = 1, l = 0; for (; lim < n; lim <<= 1, ++l); for (int i = 0; i < lim; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); return lim; } void PolyInv(int *F, int *G, int n) { static int A[N], B[N]; if (n == 1) { G[0] = qpow(F[0], mod - 2); return; } PolyInv(F, G, n >> 1); for (int i = 0; i < n; ++i) A[i] = F[i], B[i] = G[i]; int lim = NTT_init(n << 1); NTT(A, lim, 1), NTT(B, lim, 1); for (int i = 0; i < lim; ++i) A[i] = 1ll * A[i] * B[i] % mod * B[i] % mod; NTT(A, lim, -1); for (int i = 0; i < n; ++i) G[i] = (2ll * G[i] - A[i] + mod) % mod; for (int i = 0; i < lim; ++i) A[i] = B[i] = 0; } void PolyDeri(int *F, int *G, int n) { for (int i = 1; i < n; ++i) G[i - 1] = 1ll * F[i] * i % mod; G[n - 1] = 0; } void PolyInte(int *F, int *G, int n) { for (int i = 1; i < n; ++i) G[i] = 1ll * F[i - 1] * qpow(i, mod - 2) % mod; G[0] = 0; } void PolyLn(int *F, int *G, int n) { static int A[N], B[N]; PolyDeri(F, A, n), PolyInv(F, B, n); int lim = NTT_init(n << 1); NTT(A, lim, 1), NTT(B, lim, 1); for (int i = 0; i < lim; ++i) A[i] = 1ll * A[i] * B[i] % mod; NTT(A, lim, -1); PolyInte(A, G, n); for (int i = 0; i < lim; ++i) A[i] = B[i] = 0; } void PolyExp(int *F, int *G, int n) { static int A[N], B[N]; if (n == 1) { G[0] = 1; return; } PolyExp(F, G, n >> 1); for (int i = 0; i < n; ++i) A[i] = G[i]; PolyLn(G, B, n); for (int i = 0; i < n; ++i) B[i] = (mod - B[i] + F[i]) % mod; B[0] = (B[0] + 1) % mod; int lim = NTT_init(n << 1); NTT(A, lim, 1), NTT(B, lim, 1); for (int i = 0; i < lim; ++i) A[i] = 1ll * A[i] * B[i] % mod; NTT(A, lim, -1); for (int i = 0; i < n; ++i) G[i] = A[i]; for (int i = 0; i < lim; ++i) A[i] = B[i] = 0; } void PolyPow(int *F, int *G, int n, int k) { static int A[N]; PolyLn(F, A, n); for (int i = 0; i < n; ++i) A[i] = 1ll * A[i] * k % mod; PolyExp(A, G, n); for (int i = 0; i < n; ++i) A[i] = 0; } int n, m, G[N], A[N], B[N]; int main() { n = read(), m = read(); G[1] = 1; for (int i = 1; i <= m; ++i) G[read()] = mod - 1; for (int i = 0; i <= n; ++i) G[i] = G[i + 1]; int lim = 1; for (; lim < n; lim <<= 1); PolyInv(G, A, lim); PolyPow(A, B, lim, n); printf("%d\n", 1ll * B[n - 1] * qpow(n, mod - 2) % mod); return 0; }
球与盒子问题大杂烩 下面我们对这 $12$ 个问题简略分析。 I 球之间互不相同，盒子之间互不相同。 每个球都可以任意选择一个盒子，答案为 $m^n$。 II 球之间互不相同，盒子之间互不相同，每个盒子至多装一个球。 每个球选择一个盒子，已经被选择过的盒子不能再被选择，答案为 $m^{\underline{n}}$。 III 球之间互不相同，盒子之间互不相同，每个盒子至少装一个球。 考虑容斥，枚举空盒子的个数，不难得到答案为 $$ \sum _ {i = 0} ^ m (-1) ^ i {m \choose i} (m - i) ^ n $$ IV 球之间互不相同，盒子全部相同。 枚举非空盒子的个数，不难得到答案为 $$ \sum_{i = 0} ^ m \begin{Bmatrix} n \\ i \end{Bmatrix} $$ 使用 第二类斯特林数·行 的方法计算即可。具体可以看这里。 V 球之间互不相同，盒子全部相同，每个盒子至多装一个球。 答案为 $[n \leq m]$。 VI 球之间互不相同，盒子全部相同，每个盒子至少装一个球。 答案为 $\begin{Bmatrix} n \\ m \end{Bmatrix}$。 VII 球全部相同，盒子之间互不相同。 使用插板法计算，答案为 ${n + m - 1 \choose m - 1}$。 VIII 球全部相同，盒子之间互不相同，每个盒子至多装一个球。 相当于选 $n$ 个盒子放球，答案为 ${m \choose n}$。 IX 球全部相同，盒子之间互不相同，每个盒子至少装一个球。 使用插板法计算，答案为 ${n - 1 \choose m - 1}$。 X 球全部相同，盒子全部相同。 定义划分数 $p_{n, m}$ 为 $n$ 划分为 $m$ 个自然数的方案数。有一个很经典的转移 $$ p_{n, m} = p_{n, m - 1} + p_{n - m, m} $$ 意思是：加一个 $0$，或者把所有数 $+1$。 设生成函数 $F_m(x) = \sum_{i \geq 0} p_{i, m} x ^ i$，那么有 $$ F_m(x) = F_{m-1}(x) (1 + x ^ m + x ^ {2m} + \cdots) = \frac{F_{m - 1}(x)}{1 - x^m} $$ 从而 $$ F_m(x) = \prod_{i = 1} ^ m \frac{1}{1 - x^i} $$ 考虑取对数。通过一些推导可以知道 $$ \ln \frac{1}{1 - x^i} = \sum_{j\geq 1} \frac{x^{ij}}{j} $$ 这样子就可以 $\mathcal{O}(m\ln m)$ 求出 $\ln F_m(x)$，再多项式 exp 回去即可。 XI 球全部相同，盒子全部相同，每个盒子至多装一个球。 答案为 $[n \leq m]$。 XII 球全部相同，盒子全部相同，每个盒子至少装一个球。 相当于先往每个盒子里装一个球，答案为 $p_{n - m, m}$。 代码 // ==================================== // author: M_sea // website: https://m-sea-blog.com/ // ==================================== #include <bits/stdc++.h> #define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout) #define debug(...) fprintf(stderr, __VA_ARGS__) using namespace std; typedef long long ll; int read() { int X = 0, w = 1; char c = getchar(); while (c < '0' || c > '9') { if (c == '-') w = -1; c = getchar(); } while (c >= '0' && c <= '9') X = X * 10 + c - '0', c = getchar(); return X * w; } const int N = 524288 + 10; const int mod = 998244353; int qpow(int a, int b) { int c = 1; for (; b; b >>= 1, a = 1ll * a * a % mod) if (b & 1) c = 1ll * c * a % mod; return c; } int n, m; int fac[N], ifac[N], S[N], P[N]; int r[N]; void NTT(int *A, int n, int op) { for (int i = 0; i < n; ++i) if (i < r[i]) swap(A[i], A[r[i]]); for (int i = 1; i < n; i <<= 1) { int rot = qpow(op == 1 ? 3 : 332748118, (mod - 1) / (i << 1)); for (int j = 0; j < n; j += i << 1) for (int k = 0, w = 1; k < i; ++k, w = 1ll * w * rot % mod) { int x = A[j + k], y = 1ll * A[j + k + i] * w % mod; A[j + k] = (x + y) % mod, A[j + k + i] = (x - y + mod) % mod; } } if (op == -1) { int inv = qpow(n, mod - 2); for (int i = 0; i < n; ++i) A[i] = 1ll * A[i] * inv % mod; } } int NTT_init(int n) { int lim = 1, l = 0; for (; lim < n; lim <<= 1, ++l); for (int i = 0; i < lim; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); return lim; } void PolyInv(int *F, int *G, int n) { static int A[N], B[N]; if (n == 1) { G[0] = qpow(F[0], mod - 2); return; } PolyInv(F, G, n >> 1); for (int i = 0; i < n; ++i) A[i] = F[i], B[i] = G[i]; int lim = NTT_init(n << 1); NTT(A, lim, 1), NTT(B, lim, 1); for (int i = 0; i < lim; ++i) A[i] = 1ll * A[i] * B[i] % mod * B[i] % mod; NTT(A, lim, -1); for (int i = 0; i < n; ++i) G[i] = (2ll * G[i] - A[i] + mod) % mod; for (int i = 0; i < lim; ++i) A[i] = B[i] = 0; } void PolyDeri(int *F, int *G, int n) { for (int i = 1; i < n; ++i) G[i - 1] = 1ll * F[i] * i %mod; G[n - 1] = 0; } void PolyInte(int *F, int *G, int n) { for (int i = 1; i < n; ++i) G[i] = 1ll * F[i-1] * qpow(i, mod - 2) % mod; G[0] = 0; } void PolyLn(int *F, int *G, int n) { static int A[N], B[N]; PolyDeri(F, A, n), PolyInv(F, B, n); int lim = NTT_init(n << 1); NTT(A, lim, 1), NTT(B, lim, 1); for (int i = 0; i < lim; ++i) A[i] = 1ll * A[i] * B[i] % mod; NTT(A, lim, -1); PolyInte(A, G, n); for (int i = 0; i < lim; ++i) A[i] = B[i] = 0; } void PolyExp(int *F, int *G, int n) { static int A[N], B[N]; if (n == 1) { G[0] = 1; return; } PolyExp(F, G, n >> 1); for (int i = 0; i < n; ++i) A[i] = G[i]; PolyLn(G, B, n); for (int i = 0; i < n; ++i) B[i] = (mod - B[i] + F[i]) % mod; B[0] = (B[0] + 1) % mod; int lim = NTT_init(n << 1); NTT(A, lim, 1), NTT(B, lim, 1); for (int i = 0; i < lim; ++i) A[i] = 1ll * A[i] * B[i] % mod; NTT(A, lim, -1); for (int i = 0; i < n; ++i) G[i] = A[i]; for (int i = 0; i < lim; ++i) A[i] = B[i] = 0; } void init() { static int F[N], G[N]; int lim; fac[0] = 1; for (int i = 1; i <= n + m; ++i) fac[i] = 1ll * fac[i - 1] * i % mod; ifac[n + m] = qpow(fac[n + m], mod - 2); for (int i = n + m; i; --i) ifac[i - 1] = 1ll * ifac[i] * i % mod; for (int i = 0; i <= n; ++i) { F[i] = i & 1 ? mod - ifac[i] : ifac[i]; G[i] = 1ll * qpow(i, n) * ifac[i] % mod; } lim = NTT_init(n << 1 | 1); NTT(F, lim, 1), NTT(G, lim, 1); for (int i = 0; i < lim; ++i) F[i] = 1ll * F[i] * G[i] % mod; NTT(F, lim, -1); for (int i = 0; i <= n; ++i) S[i] = F[i]; for (int i = 0; i < lim; ++i) F[i] = G[i] = 0; lim = 1; for (; lim <= n; lim <<= 1); for (int i = 1; i <= m; ++i) for (int j = i; j <= n; j += i) F[j] = (F[j] + qpow(j / i, mod - 2)) % mod; PolyExp(F, G, lim); for (int i = 0; i <= n; ++i) P[i] = G[i]; } int C(int n, int m) { if (n < m) return 0; return 1ll * fac[n] * ifac[m] % mod * ifac[n - m] % mod; } int I() { return qpow(m, n); } int II() { return n > m ? 0 : 1ll * fac[m] * ifac[m - n] % mod; } int III() { int res = 0; for (int i = 0; i <= m; ++i) { int w = 1ll * C(m, i) * qpow(m - i, n) % mod; if (i & 1) res = (res - w + mod) % mod; else res = (res + w) % mod; } return res; } int IV() { int res = 0; for (int i = 1; i <= m; ++i) res = (res + S[i]) % mod; return res; } int V() { return n <= m; } int VI() { return S[m]; } int VII() { return C(n + m - 1, m - 1); } int VIII() { return C(m, n); } int IX() { return C(n - 1, m - 1); } int X() { return P[n]; } int XI() { return n <= m; } int XII() { return n < m ? 0 : P[n - m]; } int main() { n = read(), m = read(); init(); printf("%d\n%d\n%d\n", I(), II(), III()); printf("%d\n%d\n%d\n", IV(), V(), VI()); printf("%d\n%d\n%d\n", VII(), VIII(), IX()); printf("%d\n%d\n%d\n", X(), XI(), XII()); return 0; }
UOJ 分析 设 $f_{i ,j, w}$ 为从 $i$ 走到 $j$、权值和为 $w$ 的方案数，$g_{i, j, w}$ 表示连接 $i,j$ 的权值为 $k$ 的边数。 那么有 $$ f_{i, j, w} = \sum_{k = 1}^n \sum_{v} f_{i, k, v} \times g_{k, j, w - v} + [i = j \land w = 0] $$ 后面那个中括号表示不走的情况。 设 $F_{i,j}(x) = \sum_{k\geq 0} f_{i,j,k} x^k$，$G_{i,j}(x) = \sum_{k\geq 0} g_{i, j, k} x^k$，并设 $F,G$ 表示这两个生成函数构成的矩阵，则有 $$ F = F\times G + I $$ 即 $$ F = \frac{I}{I - G} $$ 可以证明 $I-G$ 一定是可逆的。于是我们只需要求一个每个元素都是一个多项式的矩阵的逆，和普通矩阵求逆类似地做即可。 然而这题还卡常，需要精细地实现来减少 NTT 的次数，最好把预处理原根、取模优化、寻址优化什么的都加上去。 代码 // ==================================== // author: M_sea // website: https://m-sea-blog.com/ // ==================================== #include <bits/stdc++.h> #define file(x) freopen(#x".in", "r", stdin), freopen(#x".out", "w", stdout) #define debug(...) fprintf(stderr, __VA_ARGS__) using namespace std; typedef long long ll; int read() { int X = 0, w = 1; char c = getchar(); while (c < '0' || c > '9') { if (c == '-') w = -1; c = getchar(); } while (c >= '0' && c <= '9') X = X * 10 + c - '0', c = getchar(); return X * w; } const int N = 8 + 5, W = 131072 + 10; const int mod = 998244353; int upd(const int &x) { return x + (x >> 31 & mod); } int qpow(int a, int b) { int c = 1; for (; b; b >>= 1, a = 1ll * a * a % mod) if (b & 1) c = 1ll * c * a % mod; return c; } int r[W], w[W], ow[2][W], Lim; void NTT(int *A, int n, int op) { for (int i = 0; i < n; ++i) if (i < r[i]) swap(A[i], A[r[i]]); for (int i = 1; i < n; i <<= 1) { for (int j = 0, k = 0; j < i; ++j, k += Lim / i) w[j] = ow[op < 0][k]; for (int j = 0; j < n; j += i << 1) for (int k = 0; k < i; ++k) { int x = A[j + k], y = 1ll * w[k] * A[j + k + i] % mod; A[j + k] = upd(x + y - mod), A[j + k + i] = upd(x - y); } } if (op == -1) { int inv = qpow(n, mod - 2); for (int i = 0; i < n; ++i) A[i] = 1ll * A[i] * inv % mod; } } int NTT_init(int n) { int lim = 1, l = 0; for (; lim < n; lim <<= 1, ++l); for (int i = 0; i < lim; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); return lim; } void PolyInv(int *F, int *G, int n) { static int A[W], B[W]; if (n == 1) { G[0] = qpow(F[0], mod - 2); return; } PolyInv(F, G, n >> 1); for (int i = 0; i < n; ++i) A[i] = F[i], B[i] = G[i]; int lim = NTT_init(n << 1); NTT(A, lim, 1), NTT(B, lim, 1); for (int i = 0; i < lim; ++i) A[i] = 1ll * A[i] * B[i] % mod * B[i] % mod; NTT(A, lim, -1); for (int i = 0; i < n; ++i) G[i] = (2ll * G[i] - A[i] + mod) % mod; for (int i = 0; i < lim; ++i) A[i] = B[i] = 0; } int n, m, Q, L; int A[N][N << 1][W]; void MatrixInv() { static int F[W], G[W], H[N << 1][W]; Lim = 1; int l = 0; for (; Lim <= L << 1; Lim <<= 1, ++l); ow[0][0] = 1, ow[0][1] = qpow(3, (mod - 1) / (Lim << 1)); ow[1][0] = 1, ow[1][1] = qpow(ow[0][1], mod - 2); for (int i = 2; i < Lim; ++i) { ow[0][i] = 1ll * ow[0][i - 1] * ow[0][1] % mod; ow[1][i] = 1ll * ow[1][i - 1] * ow[1][1] % mod; } for (int i = 1; i <= n; ++i) A[i][i][0] = A[i][i + n][0] = mod - 1; for (int i = 1; i <= n; ++i) { memset(F, 0, Lim << 2); PolyInv(A[i][i], F, Lim >> 1); for (int j = 0; j < Lim; ++j) r[j] = (r[j >> 1] >> 1) | ((j & 1) << (l - 1)); NTT(F, Lim, 1); for (int j = i; j <= i + n; ++j) { NTT(A[i][j], Lim, 1); for (int k = 0; k < Lim; ++k) A[i][j][k] = 1ll * A[i][j][k] * F[k] % mod; NTT(A[i][j], Lim, -1); for (int k = L + 1; k < Lim; ++k) A[i][j][k] = 0; memcpy(H[j], A[i][j], Lim << 2); NTT(H[j], Lim, 1); } for (int j = 1; j <= n; ++j) { if (i == j) continue; memcpy(F, A[j][i], Lim << 2); NTT(F, Lim, 1); for (int k = i + 1; k <= i + n; ++k) { for (int l = 0; l < Lim; ++l) G[l] = 1ll * F[l] * H[k][l] % mod; NTT(G, Lim, -1); for (int l = 0; l <= L; ++l) A[j][k][l] = upd(A[j][k][l] - G[l]); } } } } int main() { n = read(), m = read(), Q = read(), L = read(); for (int i = 1; i <= m; ++i) { int u = read(), v = read(), w = read(); ++A[u][v][w]; } MatrixInv(); while (Q--) { int u = read(), v = read(), w = read(); printf("%d\n", A[u][v + n][w]); } return 0; }
Luogu 分析 根据二项式定理有 $$ (a_1+a_2+\cdots+a_n)^k=k![x^k]\prod e^{a_ix} $$ 所以我们只要把边权设成 $e^{a_ix}$，然后直接矩阵树定理即可。 因为数据范围很小，所以多项式乘法和求逆可以直接暴力，NTT 甚至会慢。 代码 // ==================================== // author: M_sea // website: https://m-sea-blog.com/ // ==================================== #include <bits/stdc++.h> #define file(x) freopen(#x".in","r",stdin); freopen(#x".out","w",stdout) #define debug(...) fprintf(stderr,__VA_ARGS__) using namespace std; typedef long long ll; int read() { int X=0,w=1; char c=getchar(); while (c<'0'||c>'9') { if (c=='-') w=-1; c=getchar(); } while (c>='0'&&c<='9') X=X*10+c-'0',c=getchar(); return X*w; } const int N=30+5; const int mod=998244353; int qpow(int a,int b) { int c=1; for (;b;b>>=1,a=1ll*a*a%mod) if (b&1) c=1ll*c*a%mod; return c; } int n,k,fac[N],ifac[N]; void init(int n) { fac[0]=1; for (int i=1;i<=n;++i) fac[i]=1ll*fac[i-1]*i%mod; ifac[n]=qpow(fac[n],mod-2); for (int i=n;i;--i) ifac[i-1]=1ll*ifac[i]*i%mod; } struct poly { int a[N]; poly() { memset(a,0,sizeof(a)); } poly(int w) { for (int i=0,pw=1;i<=k;++i,pw=1ll*pw*w%mod) a[i]=1ll*pw*ifac[i]%mod; } int& operator [](int i) { return a[i]; } }; poly operator +(poly a,poly b) { poly c; for (int i=0;i<=k;++i) c[i]=(a[i]+b[i])%mod; return c; } poly operator -(poly a,poly b) { poly c; for (int i=0;i<=k;++i) c[i]=(a[i]-b[i]+mod)%mod; return c; } poly operator *(poly a,poly b) { poly c; for (int i=0;i<=k;++i) for (int j=0;j<=i;++j) c[i]=(c[i]+1ll*a[j]*b[i-j])%mod; return c; } poly polyinv(poly a) { poly b; b[0]=qpow(a[0],mod-2); for (int i=1;i<=k;++i) { int s=0; for (int j=0;j<i;++j) s=(s+1ll*b[j]*a[i-j])%mod; b[i]=1ll*b[0]*(mod-s)%mod; } return b; } poly C[N][N]; poly det(int n) { poly ans(0); for (int i=1;i<=n;++i) { ans=ans*C[i][i]; poly inv=polyinv(C[i][i]); for (int j=i+1;j<=n;++j) { poly t=C[j][i]*inv; for (int k=i;k<=n;++k) C[j][k]=C[j][k]-t*C[i][k]; } } return ans; } int main() { n=read(),k=read(); init(k); for (int i=1;i<=n;++i) for (int j=1;j<=n;++j) { int w=read(); if (i>=j) continue; auto t=poly(w); C[i][i]=C[i][i]+t,C[j][j]=C[j][j]+t; C[i][j]=C[i][j]-t,C[j][i]=C[j][i]-t; } printf("%d\n",1ll*fac[k]*det(n-1)[k]%mod); return 0; }
Luogu 分析 答案显然是 $$ \frac{\sum_{i=1}^n\sum_{j=1}^m(a_i+b_j)^k}{nm} $$ 我们不管分母，只看分子。用二项式定理拆开得到 $$ \sum_{i=1}^n\sum_{j=1}^m\sum_{r=0}^k{k\choose r}a_i\!^rb_j\!^{k-r} $$ 稍微换一下形式得到 $$ \sum_{r=0}^k{k\choose r}\left(\sum_{i=1}^na_i\!^r\right)\left(\sum_{i=1}^mb_i\!^{k-r}\right) $$ 如果我们能够对每个 $r$ 快速求出 $\sum_{i=1}^n a_i\!^r$，那么就可以直接卷了。 这是一个经典问题（鸡贼说叫“等幂和”）。我们相当于要求生成函数 $$ A(x)=\sum_{r\geq 0}x^r\sum_{i=1}^n a_i\!^r $$ 的前若干项系数。 我们把 $A(x)$ 换一个形式写出来 $$ A(x)=\sum_{i=1}^n\frac{1}{1-a_ix} $$ 根据这类问题的方法，设 $$ G(x)=\sum_{i=1}^n\frac{-a_i}{1-a_ix} $$ 那么有 $$ A(x)=n-xG(x) $$ 考虑怎么算 $G(x)$。注意到 $$ \ln(1-a_ix)'=\frac{-a_i}{1-a_ix} $$ 于是可以想到把 $G(x)$ 化成 $$ G(x)=\sum_{i=1}^n\ln(1-a_ix)' $$ 把 $\sum$ 丢到 $\ln$ 里去得到 $$ G(x)=\ln\left(\prod_{i=1}^n(1-a_ix)\right)' $$ 直接分治把所有 $1-a_ix$ 卷起来，然后多项式 $\ln$ + 多项式求导即可求出 $G(x)$，再推回 $A(x)$ 后卷积即可。 代码 不想写 vector，于是从鱼那里蒯了一个分治 /kel // ==================================== // author: M_sea // website: https://m-sea-blog.com/ // ==================================== #include <bits/stdc++.h> #define file(x) freopen(#x".in","r",stdin); freopen(#x".out","w",stdout) #define debug(...) fprintf(stderr,__VA_ARGS__) using namespace std; typedef long long ll; int read() { int X=0,w=1; char c=getchar(); while (c<'0'||c>'9') { if (c=='-') w=-1; c=getchar(); } while (c>='0'&&c<='9') X=X*10+c-'0',c=getchar(); return X*w; } const int N=400000+10; const int mod=998244353; int qpow(int a,int b) { int c=1; for (;b;b>>=1,a=1ll*a*a%mod) if (b&1) c=1ll*c*a%mod; return c; } int r[N]; void NTT(int *A,int n,int op) { for (int i=0;i<n;++i) if (i<r[i]) swap(A[i],A[r[i]]); for (int i=1;i<n;i<<=1) { int rot=qpow(op==1?3:332748118,(mod-1)/(i<<1)); for (int j=0;j<n;j+=i<<1) for (int k=0,w=1;k<i;++k,w=1ll*w*rot%mod) { int x=A[j+k],y=1ll*w*A[j+k+i]%mod; A[j+k]=(x+y)%mod,A[j+k+i]=(x-y+mod)%mod; } } if (op==-1) { int inv=qpow(n,mod-2); for (int i=0;i<n;++i) A[i]=1ll*A[i]*inv%mod; } } int NTT_init(int n) { int lim=1,l=-1; for (;lim<n;lim<<=1,++l); for (int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<l); return lim; } void polyinv(int *F,int *G,int n) { static int A[N],B[N]; if (n==1) { G[0]=qpow(F[0],mod-2); return; } polyinv(F,G,n>>1); for (int i=0;i<n;++i) A[i]=F[i],B[i]=G[i]; int lim=NTT_init(n<<1); NTT(A,lim,1),NTT(B,lim,1); for (int i=0;i<lim;++i) A[i]=1ll*A[i]*B[i]%mod*B[i]%mod; NTT(A,lim,-1); for (int i=0;i<n;++i) G[i]=(2ll*G[i]-A[i]+mod)%mod; for (int i=0;i<lim;++i) A[i]=B[i]=0; } void polyderi(int *F,int *G,int n) { for (int i=1;i<n;++i) G[i-1]=1ll*F[i]*i%mod; G[n-1]=0; } void polyinte(int *F,int *G,int n) { for (int i=1;i<n;++i) G[i]=1ll*F[i-1]*qpow(i,mod-2)%mod; G[0]=0; } void polyln(int *F,int *G,int n) { static int A[N],B[N]; polyderi(F,A,n),polyinv(F,B,n); int lim=NTT_init(n<<1); NTT(A,lim,1),NTT(B,lim,1); for (int i=0;i<lim;++i) A[i]=1ll*A[i]*B[i]%mod; NTT(A,lim,-1); polyinte(A,G,n); for (int i=0;i<lim;++i) A[i]=B[i]=0; } #define ls (o<<1) #define rs (o<<1|1) int len[N]; void divide(int o,int l,int r,int *F,int *G) { len[o]=r-l+1; if (l==r) { G[0]=1,G[1]=mod-F[l]; return; } int L[N],R[N]; int mid=(l+r)>>1; divide(ls,l,mid,F,L),divide(rs,mid+1,r,F,R); int lim=NTT_init(len[o]+1); for (int i=len[ls]+1;i<lim;++i) L[i]=0; for (int i=len[rs]+1;i<lim;++i) R[i]=0; NTT(L,lim,1),NTT(R,lim,1); for (int i=0;i<lim;++i) G[i]=1ll*L[i]*R[i]%mod; NTT(G,lim,-1); for (int i=len[o]+1;i<lim;++i) G[i]=0; } #undef ls #undef rs int n,m,k,a[N],b[N],fac[N],ifac[N],F[N],G[N]; int main() { n=read(),m=read(); for (int i=1;i<=n;++i) a[i]=read(); for (int i=1;i<=m;++i) b[i]=read(); divide(1,1,n,a,F),divide(1,1,m,b,G); k=read(); int l=max({n,m,k})+1; int lim=1; for (;lim<=l;lim<<=1); polyln(F,F,lim),polyderi(F,F,lim); polyln(G,G,lim),polyderi(G,G,lim); for (int i=lim-1;i;--i) F[i]=mod-F[i-1]; F[0]=n; for (int i=lim-1;i;--i) G[i]=mod-G[i-1]; G[0]=m; fac[0]=1; for(int i=1;i<lim;++i) fac[i]=1ll*fac[i-1]*i%mod; ifac[lim-1]=qpow(fac[lim-1],mod-2); for (int i=lim-1;i;--i) ifac[i-1]=1ll*ifac[i]*i%mod; for (int i=0;i<lim;++i) F[i]=1ll*F[i]*ifac[i]%mod; for (int i=0;i<lim;++i) G[i]=1ll*G[i]*ifac[i]%mod; lim=NTT_init(lim<<1); NTT(F,lim,1),NTT(G,lim,1); for (int i=0;i<lim;++i) F[i]=1ll*F[i]*G[i]%mod; NTT(F,lim,-1); int inv=qpow(1ll*n*m%mod,mod-2); for (int i=1;i<=k;++i) printf("%d\n",1ll*F[i]*fac[i]%mod*inv%mod); return 0; }
为了方便，下文中的多项式长度和 $n$ 均为 $2$ 的幂（多项式除法、三角函数和反三角函数除外），模数均为 $998244353$。 下面所有函数中（除了 NTT_init）传进去的 $n,m$ 都是多项式项数而不是次数。 多项式乘法 这是前置知识，如果不会建议「点开洛谷题面单击右边蓝色按钮题解」。 int r[N]; void NTT(int* A,int n,int op) { for (int i=0;i<n;++i) if (i<r[i]) swap(A[i],A[r[i]]); for (int i=1;i<n;i<<=1) { int rot=qpow(op==1?3:332748118,(mod-1)/(i<<1)); for (int j=0;j<n;j+=i<<1) for (int k=0,w=1;k<i;++k,w=1ll*w*rot%mod) { int x=A[j+k],y=1ll*w*A[j+k+i]%mod; A[j+k]=(x+y)%mod,A[j+k+i]=(x-y+mod)%mod; } } if (op==-1) { int inv=qpow(n,mod-2); for (int i=0;i<n;++i) A[i]=1ll*A[i]*inv%mod; } } int NTT_init(int n) { int lim=1,l=-1; for (;lim<n;lim<<=1,++l); for (int i=0;i<lim;++i) r[i]=(r[i>>1]>>1)|((i&1)<<l); return lim; } 多项式求导 $$ (x^a)'=ax^{a-1} $$ void getderi(int* F,int* G,int n) { for (int i=1;i<n;++i) G[i-1]=1ll*F[i]*i%mod; G[n-1]=0; } 多项式不定积分 $$ \int x^a\mathrm{d}x=\frac{x^{a+1}}{a+1} $$ void getinte(int* F,int* G,int n) { for (int i=1;i<n;++i) G[i]=1ll*F[i-1]*qpow(i,mod-2)%mod; G[0]=0; } 多项式牛顿迭代 已知 $G(x)$，求 $F(x)$ 满足 $G(F(x))\equiv 0\pmod{x^n}$。 假设我们已经求出了 $F_o(x)$ 满足 $G(F_0(x))\equiv 0\pmod{x^{\frac{n}{2}}}$，将 $G(x)$ 在 $F_0(x)$ 处泰勒展开得到 $$ G(F(x))=\sum_{i=0}^{+\infty}\frac{G^{(i)}(F_0(x))}{i!}(F(x)-F_0(x))^i $$ 可以发现 $i\geq 2$ 时在模 $x^n$ 意义下都为 $0$，所以 $$ G(F(x))\equiv G(F_0(x))+G'(F_0(x))(F(x)-F_0(x))\pmod{x^n} $$ 又 $G(F(x))=0$，化简得到 $$ F(x)=F_0(x)-\frac{G(F_0(x))}{G'(F_0(x))}\pmod{x^n} $$ 多项式求逆 已知 $F(x)$，求 $G(x)$ 满足 $F(x)G(x)\equiv 1\pmod{x^n}$。 Luogu 假设我们已经求出了 $G_0(x)$ 满足 $F(x)G_0(x)\equiv 1\pmod{x^{\frac{n}{2}}}$，那么 $$ F(x)(G(x)-G_0(x))\equiv 0\pmod{x^\frac{n}{2}} $$ 所以 $$ G(x)-G_0(x)\equiv 0\pmod{x^{\frac{n}{2}}} $$ 平方得到 $$ G(x)^2-2G(x)G_0(x)+G_0(x)^2\equiv 0\pmod{x^n} $$ 乘上 $F(x)$ 即可得到 $$ G(x)\equiv 2G_0(x)-F(x)G_0(x)^2\pmod{x^n} $$ 这样就可以由 $G_0(x)$求出 $G(x)$ 了。 void getinv(int* F,int* G,int n) { static int A[N],B[N]; if (n==1) { G[0]=qpow(F[0],mod-2); return; } getinv(F,G,n>>1); for (int i=0;i<n;++i) A[i]=F[i],B[i]=G[i]; int lim=NTT_init(n<<1); NTT(A,lim,1),NTT(B,lim,1); for (int i=0;i<lim;++i) A[i]=1ll*A[i]*B[i]%mod*B[i]%mod; NTT(A,lim,-1); for (int i=0;i<n;++i) G[i]=(2ll*G[i]-A[i]+mod)%mod; for (int i=0;i<lim;++i) A[i]=B[i]=0; } 多项式对数函数 已知 $F(x)$，求 $G(x)$ 满足 $G(x)\equiv\ln F(x)\pmod{x^n}$。保证 $F(x)$ 常数项为 $1$。 Luogu 两边求导得到 $$ G'(x)\equiv\ln'(F(x))F'(x)\pmod{x^n} $$ 即 $$ G'(x)\equiv\frac{F'(x)}{F(x)}\pmod{x^n} $$ 求导+求逆后乘起来再不定积分回去即可。 void getln(int* F,int* G,int n) { static int A[N],B[N]; getderi(F,A,n),getinv(F,B,n); int lim=NTT_init(n<<1); NTT(A,lim,1),NTT(B,lim,1); for (int i=0;i<lim;++i) A[i]=1ll*A[i]*B[i]%mod; NTT(A,lim,-1); getinte(A,G,n); for (int i=0;i<lim;++i) A[i]=B[i]=0; } 多项式指数函数 已知 $F(x)$，求 $G(x)$ 满足 $G(x)\equiv e^{F(x)}\pmod{x^n}$。保证 $F(x)$ 常数项为 $0$。 Luogu 两边取对数得到 $$ \ln G(x)-F(x)\equiv 0\pmod{x^n} $$ 不妨设 $H(G(x))\equiv \ln G(x)-F(x)\pmod{x^n}$，套牛顿迭代得到 $$ \begin{aligned} G(x)&\equiv G_0(x)-\frac{H(G_0(x))}{H'(G_0(x))}\pmod{x^n}\\ &\equiv G_0(x)-\frac{\ln G_0(x)-F(x)}{\frac{1}{G_0(x)}}\pmod{x^n}\\ &\equiv G_0(x)(1-\ln G_0(x)+F(x))\pmod{x^n} \end{aligned} $$ void getexp(int* F,int* G,int n) { static int A[N],B[N]; if (n==1) { G[0]=1; return; } getexp(F,G,n>>1); for (int i=0;i<n;++i) A[i]=G[i]; getln(G,B,n); for (int i=0;i<n;++i) B[i]=(mod-B[i]+F[i])%mod; B[0]=(B[0]+1)%mod; int lim=NTT_init(n<<1); NTT(A,lim,1),NTT(B,lim,1); for (int i=0;i<lim;++i) A[i]=1ll*A[i]*B[i]%mod; NTT(A,lim,-1); for (int i=0;i<n;++i) G[i]=A[i]; for (int i=0;i<lim;++i) A[i]=B[i]=0; } 多项式幂函数 已知 $F(x)$ 和 $k$，求 $G(x)$ 满足 $G(x)\equiv F(x)^k\pmod{x^n}$。保证 $F(x)$ 常数项为 $1$。 Luogu 求 $\ln$ 后将系数乘上 $k$ 再 $\exp$ 回去即可。 void getpow(int* F,int* G,int n,int k) { static int A[N]; getln(F,A,n); for (int i=0;i<n;++i) A[i]=1ll*A[i]*k%mod; getexp(A,G,n); for (int i=0;i<n;++i) A[i]=0; } 如果不保证常数项为 $1$，可以先提一个 $a\times x^b$ 出来转化成常数项为 $1$ 的情况，最后再乘回去。 多项式开根 已知 $F(x)$，求 $G(x)$ 满足 $G(x)^2\equiv F(x)\pmod{x^n}$，多解保留常数项较小的那一个。 Luogu Luogu * 2 移项得到 $$ G(x)^2-F(x)\equiv 0\pmod{x^n} $$ 不妨设 $H(G(x))\equiv G(x)^2-F(x)\pmod{x^n}$，套牛顿迭代得到 $$ \begin{aligned} G(x)&\equiv G_0(x)-\frac{H(G_0(x))}{H'(G_0(x))}\pmod{x^n}\\ &\equiv G_0(x)-\frac{G_0(x)^2-F(x)}{2G_0(x)}\pmod{x^n}\\ &\equiv \frac{G_0(x)}{2}+\frac{F(x)}{2G_0(x)} \end{aligned} $$ void getsqrt(int* F,int* G,int n) { static int A[N],B[N]; if (n==1) { G[0]=Cipolla::calc(F[0]); return; } getsqrt(F,G,n>>1); for (int i=0;i<n;++i) A[i]=F[i]; getinv(G,B,n); int lim=NTT_init(n<<1); NTT(A,lim,1),NTT(B,lim,1); for (int i=0;i<lim;++i) A[i]=1ll*A[i]*B[i]%mod; NTT(A,lim,-1); for (int i=0;i<n;++i) G[i]=1ll*(A[i]+G[i])*inv2%mod; for (int i=0;i<lim;++i) A[i]=B[i]=0; } 那个 Cipolla::calc 是求二次剩余。 多项式除法 Luogu 根据题意有 $$ F(x)=G(x)Q(x)+R(x) $$ 即 $$ F(\frac{1}{x})=G(\frac{1}{x})Q(\frac{1}{x})+R(\frac{1}{x}) $$ 两边同乘 $x^n$ $$ x^nF(\frac{1}{x})=x^mG(\frac{1}{x})x^{n-m}Q(\frac{1}{x})+x^nR(\frac{1}{x}) $$ 可以发现 $x^nF(\frac{1}{x})$ 就是 $F(x)$ 系数翻转后得到的多项式，于是有 $$ F_R(x)=G_R(x)Q_R(x)+x^{n-m+1}R_R(x) $$ 因为 $Q_R(x)$ 是 $n-m$ 次式，所以我们模上 $x^{n-m+1}$ $$ F_R(x)\equiv G_R(x)Q_R(x)\pmod{x^{n-m+1}} $$ 即 $$ Q_R(x)\equiv\frac{F_R(x)}{G_R(x)}\pmod{x^{n-m+1}} $$ 多项式求逆即可。 求 $R(x)$ 就很简单了 $$ R(x)=F(x)-G(x)Q(x) $$ void getdiv(int *F,int *G,int n,int m,int *Q,int *R) { static int B[N],AR[N],BR[N],iBR[N]; for (int i=0;i<n;++i) AR[n-i-1]=F[i]; for (int i=0;i<m;++i) B[i]=G[i],BR[m-i-1]=G[i]; int lim1=1; for (;lim1<n-m+1;lim1<<=1); getinv(BR,iBR,lim1); int lim=NTT_init(n+lim1); NTT(iBR,lim,1),NTT(AR,lim,1); for (int i=0;i<lim;++i) Q[i]=1ll*iBR[i]*AR[i]%mod; NTT(Q,lim,-1); for (int i=n-m+1;i<lim;++i) Q[i]=0; reverse(Q,Q+n-m+1); lim=NTT_init(n); NTT(B,lim,1),NTT(Q,lim,1); for (int i=0;i<lim;++i) R[i]=1ll*B[i]*Q[i]%mod; NTT(Q,lim,-1),NTT(R,lim,-1); for (int i=0;i<lim;++i) R[i]=(F[i]-R[i]+mod)%mod; for (int i=m-1;i<lim;++i) R[i]=0; for (int i=0;i<lim;++i) B[i]=AR[i]=BR[i]=iBR[i]=0; } 多项式三角函数 Luogu 根据欧拉公式有 $$ e^{ix}=\cos x+i\sin x $$ 同理有 $$ e^{-ix}=\cos x-i\sin x $$ 所以 $$ \begin{aligned} \sin x&=\frac{e^{ix}-e^{-ix}}{2i}\\ \cos x&=\frac{e^{ix}+e^{-ix}}{2} \end{aligned} $$ 多项式 exp 即可。 有一个问题是这个 $i$ 是啥，其实就是 $-1$ 在模 $998244353$ 意义下的二次剩余。 void getsin(int *F,int *G,int n) { static int A[N],B[N],C[N]; for (int i=0;i<n;++i) A[i]=1ll*F[i]*img%mod; getexp(A,B,n),getinv(B,C,n); int inv=qpow(img<<1,mod-2); for (int i=0;i<n;++i) G[i]=1ll*(B[i]-C[i]+mod)*inv%mod; for (int i=0;i<n;++i) A[i]=B[i]=C[i]=0; } void getcos(int *F,int *G,int n) { static int A[N],B[N],C[N]; for (int i=0;i<n;++i) A[i]=1ll*F[i]*img%mod; getexp(A,B,n),getinv(B,C,n); int inv=qpow(2,mod-2); for (int i=0;i<n;++i) G[i]=1ll*(B[i]+C[i])*inv%mod; for (int i=0;i<n;++i) A[i]=B[i]=C[i]=0; } 多项式反三角函数 Luogu 先考虑 $\arcsin$。 $$ G(x)\equiv\arcsin F(x)\pmod{x^n} $$ 两边求导得到 $$ G'(x)\equiv\frac{F'(x)}{\sqrt{1-F(x)^2}}\pmod{x^n} $$ 再不定积分 $$ G(x)\equiv\int\frac{F'(x)}{\sqrt{1-F(x)^2}}\mathrm{d}x\pmod{x^n} $$ 多项式求导+开根+求逆+不定积分即可。 接着考虑 $\arccos$，同理可得 $$ G(x)\equiv-\int\frac{F'(x)}{\sqrt{1-F(x)^2}}\mathrm{d}x\pmod{x^n} $$ $\arctan$ 也一样 $$ G(x)\equiv\int\frac{F'(x)}{1+F(x)^2}\mathrm{d}x\pmod{x^n} $$ void getarcsin(int *F,int *G,int n) { static int A[N],B[N],C[N],D[N]; for (int i=0;i<n;++i) A[i]=F[i]; int lim=NTT_init(n<<1); NTT(A,lim,1); for (int i=0;i<lim;++i) A[i]=1ll*A[i]*A[i]%mod; NTT(A,lim,-1); for (int i=n;i<lim;++i) A[i]=0; for (int i=0;i<n;++i) A[i]=(mod-A[i])%mod; ++A[0]; getsqrt(A,B,n); getinv(B,C,n); getderi(F,D,n); lim=NTT_init(n<<1); NTT(C,lim,1),NTT(D,lim,1); for (int i=0;i<lim;++i) D[i]=1ll*C[i]*D[i]%mod; NTT(D,lim,-1); getinte(D,G,n); for (int i=0;i<lim;++i) A[i]=B[i]=C[i]=D[i]=0; } void getarccos(int *F,int *G,int n) { getarcsin(F,G,n); for (int i=0;i<n;++i) G[i]=(mod-G[i])%mod; } void getarctan(int *F,int *G,int n) { static int A[N],B[N],C[N]; for (int i=0;i<n;++i) A[i]=F[i]; int lim=NTT_init(n<<1); NTT(A,lim,1); for (int i=0;i<lim;++i) A[i]=1ll*A[i]*A[i]%mod; NTT(A,lim,-1); for (int i=n;i<lim;++i) A[i]=0; ++A[0]; getinv(A,B,n); getderi(F,C,n); lim=NTT_init(n<<1); NTT(B,lim,1),NTT(C,lim,1); for (int i=0;i<lim;++i) C[i]=1ll*B[i]*C[i]%mod; NTT(C,lim,-1); getinte(C,G,n); for (int i=0;i<lim;++i) A[i]=B[i]=C[i]=0; } 多项式多点求值 鸽了。 多项式快速插值 鸽了。