多项式与生成函数学习笔记

作者是个民科 OIer，本文并不能作为严谨、完备的多项式学习资料，也不保证不存在理论错误。本文主要是对多项式在 OI 层面的运用做了简单的整理，是一份偏实用主义的学习笔记，因此并不侧重于介绍一些重要但非 OI 考察重点的数学原理，而是侧重于如何处理生成函数、如何解题上。

FFT / NTT

以 \(\Theta(n \log n)\) 的复杂度，计算多项式在 \(n\) 个单位根处的点值，以及通过 \(n\) 个单位根处的点值还原多项式的算法。常用于 \(\Theta(\mathsf{M}(n)) = \Theta(n \log n)\) 地计算多项式乘法。

由于这个算法的原理在 OI 中是相当板的存在，就不在这里列出了。OI 中计算多项式乘法基本只需用到 NTT，即使模数非 NTT 模数也可以结合 NTT 与中国剩余定理求解。下面给出作者计算多项式乘法的模板：

Code

inline void NTT(int type) {
  int n = size(); rev_init(n);
  for(int i = 1; i < n; i++) {
    if(i < rev[i]) swap(a[i], a[rev[i]]);
  }
  static int gr[N]; gr[0] = 1;
  for(int d = 1; d < n; d <<= 1) {
    int gw = qpow(type == 1 ? 3 : invg, (mod - 1) / (d << 1));
    for(int i = 1; i < d; i++) {
      gr[i] = 1ll * gr[i - 1] * gw % mod;
    }
    for(int i = 0; i < n; i += d << 1) {
      for(int j = 0; j < d; j++) {
        int x = a[i + j], y = 1ll * gr[j] * a[i + j + d] % mod;
        a[i + j] = add(x, y);
        a[i + j + d] = dec(x, y);
      }
    }
  }
  if(type == -1) {
    for(int i = 0, invn = inv(n); i < n; i++) {
      a[i] = 1ll * a[i] * invn % mod;
    }
  }
}
inline friend Poly operator * (Poly A, Poly B) {
  int n = A.size() + B.size() - 1, m = 1;
  while(m < n) m <<= 1;
  A.resize(m), B.resize(m);
  A.NTT(1), B.NTT(1);
  for(int i = 0; i < m; i++) {
    A[i] = 1ll * A[i] * B[i] % mod;
  }
  A.NTT(-1), A.resize(n);
  return A;
}

二维 NTT：将二元 GF \(F(x, y)\) 的 \(x^n y^m\) 项映射到一个一元 GF \(G(z)\) 的 \(z ^ {n \times M + m}\) 上，\(M\) 为卷积后 \(y\) 的最高次数加一。
分治 NTT：计算出 \(f_{[l, mid]}\) 后，处理其对 \(f_{(mid, r]}\) 的贡献。多用于处理 \(f_n = a_n + c_n \sum \limits_{i = 1} ^ n f_{n - i} g_i\) 这一类半在线卷积。

Newton 迭代

可以用它实现大部分多项式基本计算（如 Inv, Sqrt, Exp 等），当然更广泛的用途是求解多项式方程。Newton 迭代的原理（即泰勒展开等）此处略过，只记录一下其通用的一种使用方式。

记 \(g(f(x)) = 0\) 为一关于 \(f(x)\) 的方程。则我们可以通过以下方式倍增地求出 \(f(x)\)。

假设我们已经在模 \(x ^ n\) 意义下得到了解为 \(f_0(x)\)，那么模 \(x ^ {2n}\) 意义下的解即为 \(f(x) = f_0(x) - \dfrac{g(f(x))}{g'(f(x))}\)，其中 \(g'(f(x))\) 表示将 \(g\) 对 \(f(x)\) 求偏导。由此，只需先求出解的常数项，就可以倍增得到需要的解了。

以下举几个多项式基本操作的例子：

多项式求逆：记原多项式为 \(h(x)\)，有 \(g(f(x)) = \dfrac{1}{f(x)} - h(x)\)，\(g'(f(x)) = -\dfrac{1}{f ^ 2(x)}\)，故每次令 \(f(x) \gets 2f(x) - f^2(x)h(x)\) 即可。
多项式 exp：记原多项式为 \(h(x)\)，有 \(g(f(x)) = \ln f(x) - h(x)\)，\(g'(f(x)) = \dfrac{1}{f(x)}\)，故每次令 \(f(x) \gets f(x) - f(x) \ln f(x) + h(x)f(x)\) 即可。
多项式开根：记原多项式为 \(h(x)\)，有 \(g(f(x)) = f^2(x) - h(x)\)，\(g'(f(x)) = 2f(x)\)，故每次令 \(f(x) \gets \dfrac{1}{2}(f(x) + \dfrac{h(x)}{f(x)})\) 即可。

一些数学基础及常用封闭形式

纯属想到多少写多少。太水的就不出示了。

广义二项式系数：定义 \(r ^ {\underline{k}} = \prod \limits_{i = 0} ^ {k - 1} (r - i), r \in \mathbf{C}\)。那么可以定义广义二项式系数 \(\dbinom{r}{k} = \dfrac{r ^ {\underline k}}{k!}, r \in \mathbf C, k \in \mathbf N\)。
广义二项式定理（牛顿二项式定理）：对于 \(\alpha \in \mathbf C\)，有 \((1 + x) ^ \alpha = \sum \limits_{n \geq 0} \dbinom{\alpha}{n} x ^ n\)。
\(\dfrac{1}{(1 - x) ^ {c + 1}} = \sum \limits_{n \geq 0} \dbinom{n + c}{c} x ^ n\)，可自行使用广义二项式定理或组合意义推导。
\(\dfrac 1 {\sqrt{1 - 4x}} = \sum \limits_{n \geq 0} \dbinom{- \frac 1 2}{n} (-4x) ^ n = 1 + \sum \limits_{n \geq 1} \dfrac{(2n - 1)!! 2 ^ n}{n!} x ^ n = \sum \limits_{n \geq 0} \dbinom{2n}{n} x ^ n\)，你可能会在含 \(\dbinom{2n}{n}\) 的 GF 中见到它。
令 \(G(x)\) 为一常数项为 \(0\) 的多项式，对于一类令 \(F(x) \gets \dfrac{F(x)}{1 - G(x)}\) 的操作，可以被等价于从 \(0\) 到 \(n\) 令 \(f_n \gets f_n + \sum \limits_{i = 1} ^ n f_{n - i}g_i\)。当 \(G(x)\) 只含 \(\mathcal{O}(1)\) 项时会采取这种方式。例如乘 \(\dfrac{1}{1 + x}\) 等价于 \(f_n \gets f_n - f_{n - 1}\)。
\(\sum \limits_{n \geq 0} [x ^ n]F(x) = F(1)\)。
在一些 EGF 有关问题（如期望停时问题）中可能会遇到 \(\sum \limits_{n \geq 0} n![x ^ n]ax ^ b e ^ {cx}\) 这种形式的求和。此时有 \(\sum \limits_{n \geq 0} n![x ^ n]a x^b e^{cx} = a \sum\limits_{n \geq b} \dfrac{n!c ^ {n - b}}{(n - b)!} = \dfrac{ab!}{(1 - c) ^ {b + 1}}\)。
\(-\ln(1 - x) = \sum\limits_{n \geq 1} \dfrac{x^n}{n}\)。在使用 ln/exp 技巧时常用，如 Euler 变换。
记 \(F(x) = \sum \limits_{n \geq 1} \dfrac{f_n}{n!} x ^ n\)，则 \(\exp(F(x))\) 具有的组合意义为：将 \(n\) 个元素无序地划分为若干集合，每个集合内部的方案数为 \(f_n\)，则方案数的 EGF 即为 \(\exp(F(x))\)。可以用递推的方式证明，即枚举 \(1\) 号元素所在集合的大小，可知结果与代数意义上的式子吻合。也有更直白的理解方式：\(\exp(F(x)) = \sum \limits_{i \geq 0} \dfrac{F ^ i(x)}{i!}\)，意义显然。

Lagrange 反演

解决一定条件下多项式复合中单项系数提取的问题。

首先引入 Lagrange 反演前需要了解形式 Laurent 级数，这里就不多介绍了。

称两个多项式 \(F(x), G(x)\) 互为复合逆，当且仅当 \(F(G(x)) = x\)（可以证明这与 \(G(F(x)) = x\) 等价）。默认下文用到复合逆的多项式均满足常数项为 \(0\)，一次项非 \(0\)。

对于两互为复合逆的多项式 \(F(x), G(x)\)，有 Lagrange 反演如下：

\[n[x ^ n]F^k(x) = k[x ^ {n - k}](\dfrac{x}{G(x)})^{n} \]

将此式推及其线性组合，可以得到扩展 Lagrange 反演（\(H(x)\) 可以认为是任意多项式）：

\[[x ^ n]H(F(x)) = \frac{1}{n} [x ^ {n - 1}] H'(x) (\frac{x}{G(x)}) ^ n \]

容易发现，此式是最通用的 Lagrange 反演公式。由此式可以进一步推导得到，若 \(A(B(x)) = C(x)\)，那么（用 \(B ^ {<-1>}(x)\) 代换 \(x\)，然后使用上面的式子）

\[[x ^ n]A(x) = \dfrac{1}{n} [x ^ {n - 1}] C'(x) (\dfrac{x}{B(x)}) ^ n \]

做题记录

标 * 的是经典或相对优质的题目。基本按照作者完成题目的时间顺序编排。前面的题目多来自
Aleph1022 的题单生成函数之 OGF & EGF - 从入门到入土，后面的题目可能会混乱邪恶一些（

P5110 块速递推

用 GF 解决数列递推通项的问题。算是一道开胃前菜。

设数列 \(a\) 的 OGF 为 \(F(x)\)，立得 \(F(x) = 233xF(x) + 666x^2F(x) + x\)，整理得 \(F(x) = \dfrac{x}{1 - 233x - 666x^2}\)。这个形式不便于提取系数，考虑我们一定可以将分母分解为 \((1 - px)(1 - qx)\) 的形式，这样可以将 \(F(x)\) 转写为 \(\dfrac{\frac{1}{1 - px} - \frac{1}{1 - qx}}{p - q}\)。此式的系数是容易提取的，有 \([x ^ n]F(x) = \dfrac{p ^ n - q ^ n}{p - q}\)。容易解出 \(p, q\)，暴力求二次剩余后光速幂即可。时间复杂度 \(\Theta(\sqrt V + T)\)。

Code

#include <bits/stdc++.h>
using namespace std;
using i64 = unsigned long long;

namespace Data {
  i64 SA, SB, SC;
  inline void init() {cin >> SA >> SB >> SC;}
  inline i64 rand() {
    SA ^= SA << 32, SA ^= SA >> 13, SA ^= SA << 1;
    i64 t = SA;
    SA = SB, SB = SC, SC ^= t ^ SA;
    return SC;
  }
}

constexpr int mod = 1e9 + 7;
namespace basic {
  inline int add(int x, int y) {return (x + y >= mod ? x + y - mod : x + y);}
  inline int dec(int x, int y) {return (x - y < 0 ? x - y + mod : x - y);}
  inline void ad(int &x, int y) {x = add(x, y);}
  inline void de(int &x, int y) {x = dec(x, y);}

  inline int qpow(int a, int b) {
    int r = 1;
    while(b) {
      if(b & 1) r = 1ll * r * a % mod;
      a = 1ll * a * a % mod; b >>= 1;
    }
    return r;
  }
  inline int inv(int x) {return qpow(x, mod - 2);}
}
using namespace basic;

constexpr int G = 188305837, invG = 233230706, inv2 = 500000004;
constexpr int A = 94153035, B = 905847205, A_ = 168644639, B_ = 560842502;
int A0[1 << 15], A1[1 << 15];
int B0[1 << 15], B1[1 << 15];

constexpr int U = (1 << 15) - 1;

inline int powA(int n) {
  return 1ll * A1[n >> 15] * A0[n & U] % mod;
}
inline int powB(int n) {
  return 1ll * B1[n >> 15] * B0[n & U] % mod;
}

int main() {
  ios::sync_with_stdio(false);
  cin.tie(nullptr), cout.tie(nullptr);
  
  A0[0] = A1[0] = B0[0] = B1[0] = 1;
  for(int i = 1; i < (1 << 15); i++) {
    A0[i] = 1ll * A0[i - 1] * A % mod;
    B0[i] = 1ll * B0[i - 1] * B % mod;
  }
  for(int i = 1; i < (1 << 15); i++) {
    A1[i] = 1ll * A1[i - 1] * A_ % mod;
    B1[i] = 1ll * B1[i - 1] * B_ % mod;
  }
  int T; cin >> T; Data::init();
  int ret = 0;
  while(T--) {
    i64 n = Data::rand(); n %= (mod - 1);
    ret ^= 1ll * dec(powA(n), powB(n)) * invG % mod;
  }
  cout << ret << "\n";

  // find sqrt(56953):
  // for(int i = 1; i < mod; i++) {
  //   if(1ll * i * i % mod == 56953) {
  //     cout << i << "\n";
  //     return 0;
  //   }
  // }
}

P5320 [BJOI2019] 勘破神机

用 GF 解决数列递推通项，斯特林数实现下降幂转普通幂，扩域技巧。

两行网格下枚举最后是摆了一个 \(2 \times 1\) 的 piece 还是两个 \(1 \times 2\) 的 piece 可得 \(f_{n} = f_{n - 1} + f_{n - 2}, f_0 = 1, f_1 = 1\)。

三行网格下 \(n\) 为奇数时无解，故记录 \(g_n\) 表示 \(2n\) 的答案。最后的形态可能是摆了三个 \(1 \times 2\) 的 piece，否则一定是两个 \(2 \times 1\) 的 piece 夹了若干 \(1 \times 2\) 的 piece，其上方 / 下方铺满 \(1 \times 2\) 的 piece。可得 \(g_{n} = g_{n - 1} + 2\sum \limits_{i \geq 1} g_{n - i}\)，可化简为 \(g_{n} = 4g_{n - 1} - g_{n - 2}, g_0 = 1, g_1 = 3\)。

两者的做法基本一致，这里就只讲述三行网格的情况了。问题所求即为 \(\sum \limits_{i = 0} ^ n \dbinom{g_i}{k} = \dfrac{1}{k!} \sum \limits_{i = 0} ^ n g_i ^ {\underline{k}}\)。使用第一类斯特林数将下降幂转为普通幂，即 \(\alpha ^ {\underline{k}} = \sum \limits_{i = 0} ^ k (-1) ^ {k - i} {k \brack i} \alpha ^ i\)，问题所求可转写为 \(\sum \limits_{j = 0} ^ k (-1) ^ {k - j} {k \brack j} \sum \limits_{i = 0} ^ n g_i ^ j\)。

至此容易得到 \(\Theta(k ^ 3 \log n)\) 的矩阵做法，但本题的数据范围并不允许。考虑在这种类型的递推下，\(g_n\) 一定可以表示为 \(A \alpha ^ n + B \beta ^ n\) 的形式，这是一个值得尝试的突破口。使用上一题的方法可以解出 \(g_n = \dfrac{(p ^ {n + 1} - q ^ {n + 1}) - (p ^ n - q ^ n)}{p - q}, p = 2 + \sqrt 3, q = 2 - \sqrt 3\)，整理得 \(g_n = \dfrac{(3 + \sqrt 3) (2 + \sqrt 3) ^ n}{6} + \dfrac{(3 - \sqrt 3) (2 - \sqrt 3) ^ n}{6}\)，即 \(A = \dfrac{3 + \sqrt 3}{6}, \alpha = 2 + \sqrt 3, B = \dfrac{3 - \sqrt 3}{6}, \beta = 2 - \sqrt 3\)。

这样就可以进一步化简答案了：

\[\begin{aligned} \sum \limits_{i = 0} ^ n g_i ^ k &= \sum \limits_{i = 0} ^ n (A\alpha^i + B\beta^i) ^ k \\ &= \sum \limits_{i = 0} ^ n \sum \limits_{j = 0} ^ k \binom{k} j (A\alpha^i) ^ j (B\beta^i) ^ {k - j} \\ &= \sum \limits_{j = 0} ^ k \binom k j A ^ j B ^ {k - j} \sum \limits_{i = 0} ^ n \alpha ^ {ij} \beta ^ {ik - ij} \end{aligned} \]

后面的式子是等比数列求和的形式，可以 \(\Theta(\log n)\) 完成。那么此式就可以 \(\Theta(k \log n)\) 地计算，总时间复杂度即为 \(\Theta(k ^ 2 \log n)\)。

有一个问题是，\(\sqrt 3\) 在模 \(998244353\) 意义下是无解的。考虑将运算过程中的数扩域为 \(a + b \sqrt 3\) 的形式即可解决（在此基础上定义模 \(998244353\) 意义下的四则运算，具体可参考代码），显然答案一定会是整数。

Code

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

constexpr int mod = 998244353, N = 500 + 10;
namespace basic {
  inline int add(int x, int y) {return (x + y >= mod ? x + y - mod : x + y);}
  inline int dec(int x, int y) {return (x - y < 0 ? x - y + mod : x - y);}
  inline void ad(int &x, int y) {x = add(x, y);}
  inline void de(int &x, int y) {x = dec(x, y);}

  inline int qpow(int a, int b) {
    int r = 1;
    while(b) {
      if(b & 1) r = 1ll * r * a % mod;
      a = 1ll * a * a % mod; b >>= 1;
    }
    return r;
  }
  inline int inv(int x) {return qpow(x, mod - 2);}

  int fac[N], ifac[N];
  inline void fac_init(int n = N - 1) {
    fac[0] = 1;
    for(int i = 1; i <= n; i++)
      fac[i] = 1ll * fac[i - 1] * i % mod;
    ifac[n] = inv(fac[n]);
    for(int i = n - 1; i >= 0; i--)
      ifac[i] = 1ll * ifac[i + 1] * (i + 1) % mod;
  }
  int invx[N];
  inline void inv_init(int n = N - 1) {
    invx[1] = 1;
    for(int i = 2; i <= n; i++)
      invx[i] = 1ll * (mod - mod / i) * invx[mod % i] % mod;
  }
  inline int binom(int n, int m) {
    if(n < m || m < 0) return 0;
    return 1ll * fac[n] * ifac[m] % mod * ifac[n - m] % mod;
  }
}
using namespace basic;

template <int P>
struct Sq {
  int a, b;

  inline friend Sq operator + (Sq x, Sq y) {
    return {add(x.a, y.a), add(x.b, y.b)};
  }
  inline friend Sq operator - (Sq x, Sq y) {
    return {dec(x.a, y.a), dec(x.b, y.b)};
  }
  inline friend Sq operator * (Sq x, Sq y) {
    return {add(1ll * x.a * y.a % mod, 1ll * x.b * y.b % mod * P % mod), 
            add(1ll * x.a * y.b % mod, 1ll * x.b * y.a % mod)};
  }
  inline Sq Inv() {
    int coef = inv(dec(1ll * a * a % mod, 1ll * b * b % mod * P % mod));
    return {1ll * a * coef % mod, 1ll * dec(0, b) * coef % mod};
  }
  inline friend Sq operator / (Sq x, Sq y) {
    return x * y.Inv();
  }

  inline Sq power(i64 n) {
    Sq ret = {1, 0}, t = *this;
    while(n) {
      if(n & 1) ret = ret * t;
      n >>= 1, t = t * t;
    }
    return ret;
  }
  inline Sq Getsum(i64 n) {
    if(a == 1 && b == 0) return {(n + 1) % mod, 0};
    Sq x = (*this).power(n + 1); de(x.a, 1);
    Sq y = *this; de(y.a, 1);
    return x / y;
  }
};

int type;
int S[N][N];

inline int calc(i64 n, int k) {
  int ans = 0;
  if(type == 2) {
    Sq<5> A = {0, inv(5)}, B = {0, dec(0, inv(5))};
    Sq<5> x = {inv(2), inv(2)}, y = {inv(2), dec(0, inv(2))};
    for(int j = 0; j <= k; j++) {
      Sq<5> g = {0, 0};
      for(int t = 0; t <= j; t++) {
        Sq<5> c = x.power(t) * y.power(j - t); c = c.Getsum(n + 1), de(c.a, 1);
        Sq<5> res = A.power(t) * B.power(j - t) * c;
        res.a = 1ll * res.a * binom(j, t) % mod;
        res.b = 1ll * res.b * binom(j, t) % mod;
        g = g + res;
      }
      assert(g.b == 0);
      int res = g.a;
      int coef = 1ll * ((k - j) % 2 == 0 ? 1 : mod - 1) * S[k][j] % mod;
      ad(ans, 1ll * res * coef % mod);
    }
  } else {
    Sq<3> A = {inv(2), inv(6)}, B = {inv(2), dec(0, inv(6))};
    Sq<3> x = {2, 1}, y = {2, mod - 1};
    for(int j = 0; j <= k; j++) {
      Sq<3> g = {0, 0};
      for(int t = 0; t <= j; t++) {
        Sq<3> c = x.power(t) * y.power(j - t);
        Sq<3> res = A.power(t) * B.power(j - t) * c.Getsum(n);
        res.a = 1ll * res.a * binom(j, t) % mod;
        res.b = 1ll * res.b * binom(j, t) % mod;
        g = g + res;
      }
      assert(g.b == 0);
      int res = g.a;
      int coef = 1ll * ((k - j) % 2 == 0 ? 1 : mod - 1) * S[k][j] % mod;
      ad(ans, 1ll * res * coef % mod);
    }
  }
  return ans;
}

int k; i64 l, r;

void Main() {
  cin >> l >> r >> k;
  int coef = inv((r - l + 1) % mod);
  if(type == 3) r = r / 2, l = (l + 1) / 2;
  int res = dec(calc(r, k), calc(l - 1, k));
  res = 1ll * res * ifac[k] % mod;
  res = 1ll * res * coef % mod;
  cout << res << "\n";
}

int main() {
  ios::sync_with_stdio(false);
  cin.tie(nullptr), cout.tie(nullptr);

  fac_init();
  S[0][0] = 1;
  for(int i = 1; i < N; i++) {
    for(int j = 1; j <= i; j++) {
      S[i][j] = add(S[i - 1][j - 1], 1ll * (i - 1) * S[i - 1][j] % mod);
    }
  }

  int T; cin >> T >> type;
  while(T--) {
    Main();
  }
}

P3978 [TJOI2015] 概率论

纯粹的推式子。

\(n\) 个点的二叉树数量当然是人尽皆知的，不过这里我们需要它的 OGF。令它的 OGF 为 \(F(x)\)，有 \(F(x) = 1 + xF^2(x)\)。通过解二次方程，可得 \(F(x) = \dfrac{1 \pm \sqrt{1 - 4x}}{2x}\)，这里需要讨论取哪一个根。考虑将其变形为 \(\dfrac{2}{1 \mp \sqrt{1 - 4x}}\)，由于常数项收敛，代入 \(x = 0\) 可知应取 \(F(x) = \dfrac{1 - \sqrt{1 - 4x}}{2x}\)。

令所有 \(n\) 个点的二叉树的叶子数总和为 \(g_n\)，建立 \(g\) 的 OGF \(G(x)\)。考虑拆开两棵子树的贡献，可得 \(G(x) = 2xF(x)G(x) + x\)。

\[\begin{aligned} G(x) &= \frac{x}{1 - 2xF(x)} \\ &= \frac{x}{\sqrt{1 - 4x}} \\ &= \sum \limits_{n \geq 1} \binom{2n - 2}{n - 1} x ^ n \end{aligned} \]

因此答案为 \(\dfrac{\binom{2n - 2}{n - 1}}{\frac{\binom{2n}{n}}{n + 1}} = \dfrac{n(n + 1)}{2(2n - 1)}\)。

Code

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

int main() {
  ios::sync_with_stdio(false);
  cin.tie(nullptr), cout.tie(nullptr);
  
  int n; cin >> n;
  double res = 1. * n * (n + 1) / 2 / (2 * n - 1);
  cout << setprecision(9) << fixed << res << "\n";
}

CF438E The Child and Binary Tree

记 \(W(x) = \sum \limits_{i = 1} ^ {n} x ^ {c_i}\)，答案的 OGF 为 \(F(x)\)。考虑根的权值立得 \(F(x) = W(x)F ^ 2(x) + 1\)，即 \(F(x) = \dfrac{1 \pm \sqrt{1 - 4W(x)}}{2W(x)}\)。用类似上一题的方式讨论一下常数项是否收敛，最终可得 \(F(x) = \dfrac{1 - \sqrt{1 - 4W(x)}}{2W(x)}\)。注意这里分母的常数项为 \(0\)，不能直接求逆，一种处理方式是将分子分母同时除以 \(x ^ {n_0}\)，另一种是直接转写为 \(\dfrac{2}{1 + \sqrt{1 - 4W(x)}}\)。

Code

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

constexpr int mod = 998244353, N = 1 << 20, invg = (mod + 1) / 3;
namespace basic {
  inline int add(int x, int y) {return (x + y >= mod ? x + y - mod : x + y);}
  inline int dec(int x, int y) {return (x - y < 0 ? x - y + mod : x - y);}
  inline void ad(int &x, int y) {x = add(x, y);}
  inline void de(int &x, int y) {x = dec(x, y);}

  inline int qpow(int a, int b) {
    int r = 1;
    while(b) {
      if(b & 1) r = 1ll * r * a % mod;
      a = 1ll * a * a % mod; b >>= 1;
    }
    return r;
  }
  inline int inv(int x) {return qpow(x, mod - 2);}

  int fac[N], ifac[N];
  inline void fac_init(int n = N - 1) {
    fac[0] = 1;
    for(int i = 1; i <= n; i++)
      fac[i] = 1ll * fac[i - 1] * i % mod;
    ifac[n] = inv(fac[n]);
    for(int i = n - 1; i >= 0; i--)
      ifac[i] = 1ll * ifac[i + 1] * (i + 1) % mod;
  }
  int invx[N];
  inline void inv_init(int n = N - 1) {
    invx[1] = 1;
    for(int i = 2; i <= n; i++)
      invx[i] = 1ll * (mod - mod / i) * invx[mod % i] % mod;
  }
  inline int binom(int n, int m) {
    if(n < m || m < 0) return 0;
    return 1ll * fac[n] * ifac[m] % mod * ifac[n - m] % mod;
  }

  int rev[N];
  inline void rev_init(int n) {
    for(int i = 1; i < n; i++) {
      rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? n >> 1 : 0);
    }
  }
}
using namespace basic;

struct Poly {
  vector<int> a;
  inline int & operator [] (int x) {return a[x];}
  inline int size() {return a.size();}
  inline void resize(int n) {a.resize(n, 0);}
  Poly() {}
  explicit Poly(int n) {resize(n);}
  Poly(vector<int> b) {a = b;}
  inline void clear() {a.clear();}

  inline void NTT(int type) {
    int n = size();
    rev_init(n);
    for(int i = 1; i < n; i++) {
      if(i < rev[i]) {
        swap(a[i], a[rev[i]]);
      }
    }
    static int gr[N]; gr[0] = 1;
    for(int d = 1; d < n; d <<= 1) {
      int gw = qpow(type == 1 ? 3 : invg, (mod - 1) / (d << 1));
      for(int i = 1; i < d; i++) {
        gr[i] = 1ll * gr[i - 1] * gw % mod;
      }
      for(int i = 0; i < n; i += d << 1) {
        for(int j = 0; j < d; j++) {
          int x = a[i + j], y = 1ll * a[i + j + d] * gr[j] % mod;
          a[i + j] = add(x, y);
          a[i + j + d] = dec(x, y);
        }
      }
    }
    if(type == -1) {
      for(int i = 0, invn = inv(n); i < n; i++) {
        a[i] = 1ll * a[i] * invn % mod;
      }
    }
  }
  inline friend Poly operator * (Poly x, Poly y) {
    int m = x.size() + y.size() - 1;
    int n = 1;
    while(n < m) {
      n <<= 1;
    }
    x.resize(n), y.resize(n);
    x.NTT(1), y.NTT(1);
    for(int i = 0; i < n; i++) {
      x[i] = 1ll * x[i] * y[i] % mod;
    }
    x.NTT(-1);
    x.resize(m);
    return x;
  }

  inline friend Poly operator + (Poly x, Poly y) {
    Poly ret(max(x.size(), y.size()));
    for(int i = 0; i < ret.size(); i++) {
      if(i < x.size()) ad(ret[i], x[i]);
      if(i < y.size()) ad(ret[i], y[i]);
    }
    return ret;
  }
  inline friend Poly operator - (Poly x, Poly y) {
    Poly ret(max(x.size(), y.size()));
    for(int i = 0; i < ret.size(); i++) {
      if(i < x.size()) ad(ret[i], x[i]);
      if(i < y.size()) de(ret[i], y[i]);
    }
    return ret;
  }

  inline Poly operator += (Poly x) {
    return (*this) = (*this) + x;
  }
  inline Poly operator -= (Poly x) {
    return (*this) = (*this) - x;
  }
  inline Poly operator *= (Poly x) {
    return (*this) = (*this) * x;
  }

  inline Poly operator * (int x) {
    for(int i = 0; i < size(); i++) {
      a[i] = 1ll * a[i] * x % mod;
    }
    return (*this);
  }

  inline Poly Inv() {
    assert(a[0] != 0);
    Poly f0 = vector<int>{inv(a[0])};
    while(f0.size() < size()) {
      int n = f0.size() << 1;
      Poly h = (*this); h.resize(n);
      f0 = f0 * (vector<int>{2} - f0 * h); f0.resize(n);
    }
    f0.resize(size());
    return f0;
  }
  inline Poly Sqrt() {
    assert(a[0] != 0);
    Poly f0 = vector<int>{1};
    while(f0.size() < size()) {
      int n = f0.size() << 1; f0.resize(n);
      Poly h = (*this); h.resize(n);
      f0 = (f0 * f0 + h) * f0.Inv(); f0.resize(n);
      int inv2 = (mod + 1) / 2;
      for(int i = 0; i < n; i++) {
        f0[i] = 1ll * f0[i] * inv2 % mod;
      }
    }
    f0.resize(size());
    return f0;
  }
};

int main() {
  ios::sync_with_stdio(false);
  cin.tie(nullptr), cout.tie(nullptr);
  
  int m, n; cin >> m >> n;
  Poly W(n + 1);
  for(int i = 1; i <= m; i++) {
    int x; cin >> x;
    if(x <= n) {
      W[x] = 1;
    }
  }
  Poly F = vector<int>{1} - W * 4;
  F = vector<int>{1} + F.Sqrt();
  F = F.Inv() * 2;
  for(int i = 1; i <= n; i++) {
    cout << F[i] << "\n";
  }
}

P2767 树的数量

【模板】Lagrange 反演。

直接记有标号 \(m\) 叉树的 EGF 为 \(F(x)\)，那么 \(F(x) = x (1 + F(x)) ^ m\)。记 \(G(x)\) 为 \(F(x)\) 的复合逆，\(x \to G(x)\) 得 \(G(x) = \dfrac{x}{(1 + x) ^ m}\)。由 Lagrange 反演，\([x ^ n]F(x) = \dfrac{1}{n} [x ^ {n - 1}] (\dfrac{x}{G(x)}) ^ n = \dfrac{1}{n} [x ^ {n - 1}] (1 + x) ^ {nm} = \dfrac{\binom{nm}{n - 1}}{n}\)。

Code

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

constexpr int mod = 10007, N = mod;
namespace basic {
  inline int add(int x, int y) {return (x + y >= mod ? x + y - mod : x + y);}
  inline int dec(int x, int y) {return (x - y < 0 ? x - y + mod : x - y);}
  inline void ad(int &x, int y) {x = add(x, y);}
  inline void de(int &x, int y) {x = dec(x, y);}

  inline int qpow(int a, int b) {
    int r = 1;
    while(b) {
      if(b & 1) r = 1ll * r * a % mod;
      a = 1ll * a * a % mod; b >>= 1;
    }
    return r;
  }
  inline int inv(int x) {return qpow(x, mod - 2);}

  int fac[N], ifac[N];
  inline void fac_init(int n = N - 1) {
    fac[0] = 1;
    for(int i = 1; i <= n; i++)
      fac[i] = 1ll * fac[i - 1] * i % mod;
    ifac[n] = inv(fac[n]);
    for(int i = n - 1; i >= 0; i--)
      ifac[i] = 1ll * ifac[i + 1] * (i + 1) % mod;
  }
  int invx[N];
  inline void inv_init(int n = N - 1) {
    invx[1] = 1;
    for(int i = 2; i <= n; i++)
      invx[i] = 1ll * (mod - mod / i) * invx[mod % i] % mod;
  }
  inline int binom(int n, int m) {
    if(n < m || m < 0) return 0;
    return 1ll * fac[n] * ifac[m] % mod * ifac[n - m] % mod;
  }
}
using namespace basic;

inline int lucas(int n, int m) {
  if(n < mod) {
    return binom(n, m);
  }
  return 1ll * lucas(n / mod, m / mod) * binom(n % mod, m % mod) % mod;
}

int main() {
  ios::sync_with_stdio(false);
  cin.tie(nullptr), cout.tie(nullptr);
  
  fac_init();

  int n, m; cin >> n >> m;
  cout << 1ll * lucas(n * m, n - 1) * inv(n) % mod << "\n";
}

P4841 [集训队作业2013] 城市规划

考虑 EGF exp 的组合意义是将元素作无序划分，而一张无向图中的若干个连通块恰巧可以被视为将 \(n\) 个点作无序划分后，每个子集内部连通。令 \(F(x)\) 为有标号简单无向图的 EGF，\(G(x)\) 为连通图的 EGF，则 \(\exp(G(x)) = F(x)\)，\(G(x) = \ln F(x)\)。

Code

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

constexpr int mod = 1004535809, N = 1 << 20, invg = (mod + 1) / 3;
namespace basic {
  inline int add(int x, int y) {return (x + y >= mod ? x + y - mod : x + y);}
  inline int dec(int x, int y) {return (x - y < 0 ? x - y + mod : x - y);}
  inline void ad(int &x, int y) {x = add(x, y);}
  inline void de(int &x, int y) {x = dec(x, y);}

  inline int qpow(int a, int b) {
    int r = 1;
    while(b) {
      if(b & 1) r = 1ll * r * a % mod;
      a = 1ll * a * a % mod; b >>= 1;
    }
    return r;
  }
  inline int inv(int x) {return qpow(x, mod - 2);}

  int fac[N], ifac[N];
  inline void fac_init(int n = N - 1) {
    fac[0] = 1;
    for(int i = 1; i <= n; i++)
      fac[i] = 1ll * fac[i - 1] * i % mod;
    ifac[n] = inv(fac[n]);
    for(int i = n - 1; i >= 0; i--)
      ifac[i] = 1ll * ifac[i + 1] * (i + 1) % mod;
  }
  int invx[N];
  inline void inv_init(int n = N - 1) {
    invx[1] = 1;
    for(int i = 2; i <= n; i++)
      invx[i] = 1ll * (mod - mod / i) * invx[mod % i] % mod;
  }
  inline int binom(int n, int m) {
    if(n < m || m < 0) return 0;
    return 1ll * fac[n] * ifac[m] % mod * ifac[n - m] % mod;
  }

  int rev[N];
  inline void rev_init(int n) {
    for(int i = 1; i < n; i++) {
      rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? n >> 1 : 0);
    }
  }
}
using namespace basic;

struct Poly {
  vector<int> a;
  inline int & operator [] (int x) {return a[x];}
  inline int size() {return a.size();}
  inline void resize(int n) {a.resize(n, 0);}
  Poly() {}
  explicit Poly(int n) {resize(n);}
  Poly(vector<int> b) {a = b;}
  inline void clear() {a.clear();}

  inline void NTT(int type) {
    int n = size();
    rev_init(n);
    for(int i = 1; i < n; i++) {
      if(i < rev[i]) {
        swap(a[i], a[rev[i]]);
      }
    }
    static int gr[N]; gr[0] = 1;
    for(int d = 1; d < n; d <<= 1) {
      int gw = qpow(type == 1 ? 3 : invg, (mod - 1) / (d << 1));
      for(int i = 1; i < d; i++) {
        gr[i] = 1ll * gr[i - 1] * gw % mod;
      }
      for(int i = 0; i < n; i += d << 1) {
        for(int j = 0; j < d; j++) {
          int x = a[i + j], y = 1ll * a[i + j + d] * gr[j] % mod;
          a[i + j] = add(x, y);
          a[i + j + d] = dec(x, y);
        }
      }
    }
    if(type == -1) {
      for(int i = 0, invn = inv(n); i < n; i++) {
        a[i] = 1ll * a[i] * invn % mod;
      }
    }
  }
  inline friend Poly operator * (Poly x, Poly y) {
    int m = x.size() + y.size() - 1;
    int n = 1;
    while(n < m) {
      n <<= 1;
    }
    x.resize(n), y.resize(n);
    x.NTT(1), y.NTT(1);
    for(int i = 0; i < n; i++) {
      x[i] = 1ll * x[i] * y[i] % mod;
    }
    x.NTT(-1);
    x.resize(m);
    return x;
  }

  inline friend Poly operator + (Poly x, Poly y) {
    Poly ret(max(x.size(), y.size()));
    for(int i = 0; i < ret.size(); i++) {
      if(i < x.size()) ad(ret[i], x[i]);
      if(i < y.size()) ad(ret[i], y[i]);
    }
    return ret;
  }
  inline friend Poly operator - (Poly x, Poly y) {
    Poly ret(max(x.size(), y.size()));
    for(int i = 0; i < ret.size(); i++) {
      if(i < x.size()) ad(ret[i], x[i]);
      if(i < y.size()) de(ret[i], y[i]);
    }
    return ret;
  }

  inline Poly operator += (Poly x) {
    return (*this) = (*this) + x;
  }
  inline Poly operator -= (Poly x) {
    return (*this) = (*this) - x;
  }
  inline Poly operator *= (Poly x) {
    return (*this) = (*this) * x;
  }

  inline Poly operator * (int x) {
    for(int i = 0; i < size(); i++) {
      a[i] = 1ll * a[i] * x % mod;
    }
    return (*this);
  }

  inline Poly Inv() {
    assert(a[0] != 0);
    Poly f0 = vector<int>{inv(a[0])};
    while(f0.size() < size()) {
      int n = f0.size() << 1;
      Poly h = (*this); h.resize(n);
      f0 = f0 * (vector<int>{2} - f0 * h); f0.resize(n);
    }
    f0.resize(size());
    return f0;
  }
  inline Poly Sqrt() {
    assert(a[0] != 0);
    Poly f0 = vector<int>{1};
    while(f0.size() < size()) {
      int n = f0.size() << 1; f0.resize(n);
      Poly h = (*this); h.resize(n);
      f0 = (f0 * f0 + h) * f0.Inv(); f0.resize(n);
      int inv2 = (mod + 1) / 2;
      for(int i = 0; i < n; i++) {
        f0[i] = 1ll * f0[i] * inv2 % mod;
      }
    }
    f0.resize(size());
    return f0;
  }

  inline Poly Deriv() {
    Poly ret(size() - 1);
    for(int i = 1; i < size(); i++) {
      ret[i - 1] = 1ll * i * a[i] % mod;
    }
    return ret;
  }
  inline Poly Integ() {
    Poly ret(size() + 1);
    for(int i = 0; i < size(); i++) {
      ret[i + 1] = 1ll * inv(i + 1) * a[i] % mod;
    }
    return ret;
  }

  inline Poly Ln() {
    assert(a[0] == 1);
    Poly ret = ((*this).Deriv() * (*this).Inv()).Integ();
    ret.resize(size());
    return ret;
  }
};

int main() {
  ios::sync_with_stdio(false);
  cin.tie(nullptr), cout.tie(nullptr);
  
  fac_init();

  int n; cin >> n;
  Poly F(n + 1);
  for(int i = 0; i <= n; i++) {
    F[i] = qpow(2, 1ll * i * (i - 1) / 2 % (mod - 1));
    F[i] = 1ll * ifac[i] * F[i] % mod;
  }
  Poly G = F.Ln();
  cout << 1ll * fac[n] * G[n] % mod << "\n";
}

P5219 无聊的水题 I

最大值为 \(m\) 可以容斥为 \(\leq m\) 的方案数减去 \(\leq m - 1\) 的方案数。然后直接使用 prufer 序列就好了：\([\dfrac{x ^ {n - 2}}{(n - 2)!}](\sum \limits_{i = 0} ^ {m - 1} \dfrac{x ^ i}{i!}) ^ n\)。

Code

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

constexpr int mod = 998244353, N = 1 << 20, invg = (mod + 1) / 3;
namespace basic {
  inline int add(int x, int y) {return (x + y >= mod ? x + y - mod : x + y);}
  inline int dec(int x, int y) {return (x - y < 0 ? x - y + mod : x - y);}
  inline void ad(int &x, int y) {x = add(x, y);}
  inline void de(int &x, int y) {x = dec(x, y);}

  inline int qpow(int a, int b) {
    int r = 1;
    while(b) {
      if(b & 1) r = 1ll * r * a % mod;
      a = 1ll * a * a % mod; b >>= 1;
    }
    return r;
  }
  inline int inv(int x) {return qpow(x, mod - 2);}

  int fac[N], ifac[N];
  inline void fac_init(int n = N - 1) {
    fac[0] = 1;
    for(int i = 1; i <= n; i++)
      fac[i] = 1ll * fac[i - 1] * i % mod;
    ifac[n] = inv(fac[n]);
    for(int i = n - 1; i >= 0; i--)
      ifac[i] = 1ll * ifac[i + 1] * (i + 1) % mod;
  }
  int invx[N];
  inline void inv_init(int n = N - 1) {
    invx[1] = 1;
    for(int i = 2; i <= n; i++)
      invx[i] = 1ll * (mod - mod / i) * invx[mod % i] % mod;
  }
  inline int binom(int n, int m) {
    if(n < m || m < 0) return 0;
    return 1ll * fac[n] * ifac[m] % mod * ifac[n - m] % mod;
  }

  int rev[N];
  inline void rev_init(int n) {
    for(int i = 1; i < n; i++) {
      rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? n >> 1 : 0);
    }
  }
}
using namespace basic;

struct Poly {
  vector<int> a;
  inline int & operator [] (int x) {return a[x];}
  inline int size() {return a.size();}
  inline void resize(int n) {a.resize(n, 0);}

  Poly() {}
  explicit Poly(int n) {resize(n);}
  Poly(vector<int> b) {a = b;}

  inline void NTT(int type) {
    int n = size(); rev_init(n);
    for(int i = 1; i < n; i++) {
      if(i < rev[i]) {
        swap(a[i], a[rev[i]]);
      }
    }
    static int gr[N]; gr[0] = 1;
    for(int d = 1; d < n; d <<= 1) {
      int gw = qpow(type == 1 ? 3 : invg, (mod - 1) / (d << 1));
      for(int i = 1; i < d; i++) {
        gr[i] = 1ll * gr[i - 1] * gw % mod;
      }
      for(int i = 0; i < n; i += d << 1) {
        for(int j = 0; j < d; j++) {
          int x = a[i + j], y = 1ll * a[i + j + d] * gr[j] % mod;
          a[i + j] = add(x, y);
          a[i + j + d] = dec(x, y);
        }
      }
    }
    if(type == -1) {
      for(int i = 0, invn = inv(n); i < n; i++) {
        a[i] = 1ll * a[i] * invn % mod;
      }
    }
  }
  inline friend Poly operator * (Poly x, Poly y) {
    int m = x.size() + y.size() - 1;
    int n = 1;
    while(n < m) {
      n <<= 1;
    }
    x.resize(n), y.resize(n);
    x.NTT(1), y.NTT(1);
    for(int i = 0; i < n; i++) {
      x[i] = 1ll * x[i] * y[i] % mod;
    }
    x.NTT(-1);
    x.resize(m);
    return x;
  }

  inline friend Poly operator + (Poly x, Poly y) {
    Poly ret(max(x.size(), y.size()));
    for(int i = 0; i < ret.size(); i++) {
      if(i < x.size()) ad(ret[i], x[i]);
      if(i < y.size()) ad(ret[i], y[i]);
    }
    return ret;
  }
  inline friend Poly operator - (Poly x, Poly y) {
    Poly ret(max(x.size(), y.size()));
    for(int i = 0; i < ret.size(); i++) {
      if(i < x.size()) ad(ret[i], x[i]);
      if(i < y.size()) de(ret[i], y[i]);
    }
    return ret;
  }

  inline Poly operator += (Poly x) {
    return (*this) = (*this) + x;
  }
  inline Poly operator -= (Poly x) {
    return (*this) = (*this) - x;
  }
  inline Poly operator *= (Poly x) {
    return (*this) = (*this) * x;
  }

  inline friend Poly operator * (Poly x, int y) {
    for(int i = 0; i < x.size(); i++) {
      x[i] = 1ll * x[i] * y % mod;
    }
    return x;
  }

  inline Poly Deriv() {
    Poly ret(size() - 1);
    for(int i = 1; i < size(); i++) {
      ret[i - 1] = 1ll * i * a[i] % mod;
    }
    return ret;
  }
  inline Poly Integ() {
    Poly ret(size() + 1);
    for(int i = 0; i < size(); i++) {
      ret[i + 1] = 1ll * inv(i + 1) * a[i] % mod;
    }
    return ret;
  }

  inline Poly Inv() {
    assert(a[0] != 0);
    Poly f0 = vector<int>{inv(a[0])};
    while(f0.size() < size()) {
      int n = f0.size() << 1; f0.resize(n);
      Poly h = (*this); h.resize(n);
      f0 = f0 * (vector<int>{2} - f0 * h); f0.resize(n);
    }
    f0.resize(size());
    return f0;
  }

  inline Poly Ln() {
    assert(a[0] == 1);
    Poly ret = ((*this).Deriv() * (*this).Inv()).Integ();
    ret.resize(size());
    return ret;
  }

  inline Poly Exp() {
    assert(a[0] == 0);
    Poly f0 = vector<int>{1};
    while(f0.size() < size()) {
      int n = f0.size() << 1; f0.resize(n);
      Poly h = (*this); h.resize(n);
      f0 = f0 * (vector<int>{1} - f0.Ln() + h); f0.resize(n);
    }
    f0.resize(size());
    return f0;
  }
};

int n;
inline int calc(int m) {
  Poly F(n - 1);
  for(int i = 0; i < m; i++) {
    F[i] = ifac[i];
  }
  F = (F.Ln() * n).Exp();
  return 1ll * fac[n - 2] * F[n - 2] % mod;
}

int main() {
  ios::sync_with_stdio(false);
  cin.tie(nullptr), cout.tie(nullptr);
  
  fac_init();

  int m; cin >> n >> m;
  if(m >= n) {
    cout << 0 << "\n";
    return 0;
  }
  cout << dec(calc(m), calc(m - 1)) << "\n";
}

不过在题解区瞎翻的时候，发现了一种代数推导，与 prufer 序列得到的结果相似，在此记录一下，就当供大家理性愉悦了。

如果不用 prufer 序列，还能怎么做这类题呢？

设序列 \([d_1, d_2, \dots, d_m]\) 为一两两不同的序列，求 \(n\) 个点的有标号无根树数量，使得 \(\deg_u \in \{d\}\)。

记序列 \(d' = [d_1 - 1, d_2 - 1, \dots, d_m - 1]\)。记 \(F(x)\) 表示有标号有根树的 EGF，满足每个结点的儿子数均属于 \(d'\)，有 \(F(x) = x \sum \limits_{i = 1} ^ m \dfrac{F ^ {d_i - 1}(x)}{(d_i - 1)!}\)。这么定义的动机源于非根结点均会在父亲结点处再获得一个度数。而我们真正需要的有标号无根树可以通过 \(\dfrac{F ^ 2(x)}{2}\) 得到，即枚举无根树上的一条边并去序。不难发现，此时一棵树在每条边处都被统计了一次，因此 \([x ^ n]\dfrac{F ^ 2(x)}{2}\) 是 \(n\) 个点的有标号无根树数量的 \(n - 1\) 倍。

回到如何求 \(F(x)\) 的系数上。设 \(G(x)\) 为 \(F(x)\) 的复合逆，\(x \to G(x)\) 可得 \(G(x) = \dfrac{x}{\sum \limits_{i = 1} ^ {m} \dfrac{x ^ {d_i - 1}}{(d_i - 1)!}}\)。由 Lagrange 反演，

\[\begin{aligned} [x ^ n]F^2(x) &= \dfrac{2}{n} [x ^ {n - 2}] (\dfrac{x}{G(x)}) ^ n \\ &= \dfrac{2}{n} [x ^ {n - 2}] (\sum \limits_{i = 1} ^ m \dfrac{x ^ {d_i - 1}}{(d_i - 1)!}) ^ n \end{aligned} \]

即 \([\dfrac{x ^ n}{n!}]\dfrac{F^2(x)}{2} = (n - 1)[\dfrac{x ^ {n - 2}}{(n - 2)!}] (\sum \limits_{i = 1} ^ m \dfrac{x ^ {d_i - 1}}{(d_i - 1)!}) ^ n\)。

可以发现，这一结果与 prufer 序列的一些结论是吻合的。当然这只是一些 useless 的小想法罢了，遇到这种题肯定是直接使用轻便又具有双射意义的 prufer 序列啦。

P4233 射命丸文的笔记

接下来有一火车的图论计数问题。

先解决分子，即所有 \(n\) 个点的竞赛图的哈密顿回路数量的总和。考虑拆贡献计算之：枚举回路，求有多少个竞赛图满足存在这样的回路。这等价于钦定了图中 \(n\) 条边的方向，因此方案数为 \(2 ^ {\binom n 2 - n}\)。每条回路对应的方案数显然相同，而回路的数量有 \((n - 1)!\)。因此分子为 \((n - 1)! 2 ^ {\frac{n(n - 3)}{2}}\)。

考虑分母，即 \(n\) 个点的强连通竞赛图数量。考虑竞赛图强连通分量缩点后会形成一条链的结构，因此将 \(n\) 个点作划分，要求子集内部强连通，并钦定这若干个强连通分量在链上的顺序，即为 \(n\) 个点任意竞赛图的数量。将上述内容形式化，设任意竞赛图的 EGF 为 \(F(x)\)，强连通竞赛图的 EGF 为 \(G(x)\)，那么 \(\sum \limits_{i \geq 0} G^i(x) = F(x)\)。容易推得 \(G(x) = 1 - \dfrac{1}{F(x)}\)。直接计算即可，时间复杂度 \(\Theta(n \log n)\)。注意特判 \(n \leq 2\) 的情况。

Code

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

constexpr int mod = 998244353, N = 1 << 20, invg = (mod + 1) / 3;
namespace basic {
  inline int add(int x, int y) {return (x + y >= mod ? x + y - mod : x + y);}
  inline int dec(int x, int y) {return (x - y < 0 ? x - y + mod : x - y);}
  inline void ad(int &x, int y) {x = add(x, y);}
  inline void de(int &x, int y) {x = dec(x, y);}

  inline int qpow(int a, int b) {
    int r = 1;
    while(b) {
      if(b & 1) r = 1ll * r * a % mod;
      a = 1ll * a * a % mod; b >>= 1;
    }
    return r;
  }
  inline int inv(int x) {return qpow(x, mod - 2);}

  int fac[N], ifac[N];
  inline void fac_init(int n = N - 1) {
    fac[0] = 1;
    for(int i = 1; i <= n; i++)
      fac[i] = 1ll * fac[i - 1] * i % mod;
    ifac[n] = inv(fac[n]);
    for(int i = n - 1; i >= 0; i--)
      ifac[i] = 1ll * ifac[i + 1] * (i + 1) % mod;
  }
  int invx[N];
  inline void inv_init(int n = N - 1) {
    invx[1] = 1;
    for(int i = 2; i <= n; i++)
      invx[i] = 1ll * (mod - mod / i) * invx[mod % i] % mod;
  }
  inline int binom(int n, int m) {
    if(n < m || m < 0) return 0;
    return 1ll * fac[n] * ifac[m] % mod * ifac[n - m] % mod;
  }

  int rev[N];
  inline void rev_init(int n) {
    for(int i = 1; i < n; i++) {
      rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? n >> 1 : 0);
    }
  }
}
using namespace basic;

struct Poly {
  vector<int> a;
  inline int & operator [] (int x) {return a[x];}
  Poly() {}
  inline int size() {return a.size();}
  inline void resize(int n) {a.resize(n, 0);}
  explicit Poly(int n) {resize(n);}
  Poly(vector<int> b) {a = b;}
  inline void clear() {a.clear();}

  inline friend Poly operator + (Poly x, Poly y) {
    Poly ret(max(x.size(), y.size()));
    for(int i = 0; i < ret.size(); i++) {
      if(i < x.size()) {ad(ret[i], x[i]);}
      if(i < y.size()) {ad(ret[i], y[i]);}
    }
    return ret;
  }
  inline friend Poly operator - (Poly x, Poly y) {
    Poly ret(max(x.size(), y.size()));
    for(int i = 0; i < ret.size(); i++) {
      if(i < x.size()) {ad(ret[i], x[i]);}
      if(i < y.size()) {de(ret[i], y[i]);}
    }
    return ret;
  }

  inline void NTT(int type) {
    int n = size();
    rev_init(n);
    for(int i = 1; i < n; i++) {
      if(i < rev[i]) {
        swap(a[i], a[rev[i]]);
      }
    }
    static int gr[N]; gr[0] = 1;
    for(int d = 1; d < n; d <<= 1) {
      int gw = qpow(type == 1 ? 3 : invg, (mod - 1) / (d << 1));
      for(int i = 1; i < d; i++) {
        gr[i] = 1ll * gr[i - 1] * gw % mod;
      }
      for(int i = 0; i < n; i += d << 1) {
        for(int j = 0; j < d; j++) {
          int x = a[i + j], y = 1ll * gr[j] * a[i + j + d] % mod;
          a[i + j] = add(x, y);
          a[i + j + d] = dec(x, y);
        }
      }
    }
    if(type == -1) {
      for(int i = 0, invn = inv(n); i < n; i++) {
        a[i] = 1ll * a[i] * invn % mod;
      }
    }
  }
  inline friend Poly operator * (Poly x, Poly y) {
    int m = x.size() + y.size() - 1;
    int n = 1;
    while(n < m) {
      n <<= 1;
    }
    x.resize(n), y.resize(n);
    x.NTT(1), y.NTT(1);
    for(int i = 0; i < n; i++) {
      x[i] = 1ll * x[i] * y[i] % mod;
    }
    x.NTT(-1);
    x.resize(m);
    return x;
  }

  inline Poly operator += (Poly x) {
    return (*this) = (*this) + x;
  }
  inline Poly operator -= (Poly x) {
    return (*this) = (*this) - x;
  }
  inline Poly operator *= (Poly x) {
    return (*this) = (*this) * x;
  }

  inline Poly operator * (int x) {
    for(int i = 0; i < size(); i++) {
      a[i] = 1ll * a[i] * x % mod;
    }
    return (*this);
  }

  inline Poly Deriv() {
    Poly ret(size() - 1);
    for(int i = 1; i < size(); i++) {
      ret[i - 1] = 1ll * a[i] * i % mod;
    }
    return ret;
  }
  inline Poly Integ() {
    Poly ret(size() + 1); 
    for(int i = 0; i < size(); i++) {
      ret[i + 1] = 1ll * a[i] * inv(i + 1) % mod;
    }
    return ret;
  }

  // Implement followings by Newton's Method.
  // f(x) = f0(x) - g(f0(x)) / g'(f0(x))
  inline Poly Inv() {
    assert(a[0] != 0);
    Poly f0 = vector<int>{inv(a[0])};
    while(f0.size() < size()) {
      int n = f0.size() << 1;
      Poly h = (*this); h.resize(n), h.resize(n << 1); h.NTT(1);
      f0.resize(n << 1); f0.NTT(1);
      for(int i = 0; i < (n << 1); i++) {
        f0[i] = 1ll * f0[i] * dec(2, 1ll * f0[i] * h[i] % mod) % mod;
      }
      f0.NTT(-1), f0.resize(n);
    }
    f0.resize(size());
    return f0;
  }
  
  inline Poly Ln() {
    assert(a[0] == 1);
    Poly ret = ((*this).Deriv() * (*this).Inv()).Integ();
    ret.resize(size());
    return ret;
  }

  inline Poly Exp() {
    assert(a[0] == 0);
    Poly f0 = vector<int>{1};
    while(f0.size() < size()) {
      int n = f0.size() << 1; f0.resize(n);
      Poly t = f0.Ln();
      Poly h = (*this); h.resize(n);
      for(int i = 0; i < n; i++) {
        t[i] = dec(h[i], t[i]);
      }
      ad(t[0], 1);
      f0 = f0 * t; f0.resize(n);
    }
    f0.resize(size());
    return f0;
  }

  inline Poly Sqrt() {
    Poly f0 = vector<int>{1};
    while(f0.size() < size()) {
      int n = f0.size() << 1; f0.resize(n);
      Poly h = (*this); h.resize(n);
      f0 = f0 + h * f0.Inv(); f0.resize(n);
      int inv2 = inv(2);
      for(int i = 0; i < n; i++) {
        f0[i] = 1ll * f0[i] * inv2 % mod;
      }
    }
    return f0;
  }
};

int n, H[N];

int main() {
  ios::sync_with_stdio(false);
  cin.tie(nullptr), cout.tie(nullptr);
  
  fac_init();

  cin >> n;
  for(int i = 3; i <= n; i++) {
    H[i] = 1ll * fac[i - 1] * qpow(2, 1ll * i * (i - 3) / 2 % (mod - 1)) % mod;
  }

  Poly F(n + 1);
  for(int i = 0; i <= n; i++) {
    F[i] = qpow(2, 1ll * i * (i - 1) / 2 % (mod - 1));
    F[i] = 1ll * F[i] * ifac[i] % mod;
  }
  Poly G = vector<int>{1} - F.Inv();
  // for(int i = 1; i <= n; i++) {
  //   cout << 1ll * G[i] * fac[i] % mod << "\n";
  // }
  for(int i = 1; i <= n; i++) {
    if(i == 1) {
      cout << 1 << "\n";
      continue;
    }
    if(i == 2) {
      cout << -1 << "\n";
      continue;
    }
    cout << 1ll * H[i] * inv(1ll * G[i] * fac[i] % mod) % mod << "\n";
  }
}

P5900 无标号无根树计数

先考虑无标号有根树怎么做。对于 OGF \(F(x)\)，有 Euler 变换：\(\mathcal{E}(F(x)) = \prod \limits_{i \geq 1} \dfrac{1}{(1 - x ^ i) ^ {f_i}}\)。考虑当大小为 \(i\) 的组合结构有 \(f_i\) 种时，任意选取这些组合结构组成一个大小为 \(n\) 的组合结构的方案数即 \([x ^ n]\mathcal{E}(F(x))\)。

记无标号有根树的 OGF 为 \(F(x)\)，有 \(F(x) = x \mathcal{E}(F(x))\)。直接代用上述 Euler 变换的形式并不便于求解 \(F(x)\)，考虑转写一下 \(\mathcal{E}(F(x))\)。

\[\begin{aligned} \mathcal{E}(F(x)) &= \prod \limits_{i \geq 1} \dfrac{1}{(1 - x ^ i) ^ {f_i}} \\ &= \exp(\sum \limits_{i \geq 1} -f_i\ln(1 - x ^ i)) \\ &= \exp(\sum \limits_{i \geq 1} f_i \sum \limits_{j \geq 1} \frac{x ^ {ij}}{j}) \\ &= \exp(\sum \limits_{j \geq 1} \frac{1}{j} \sum \limits_{i \geq 1} f_i x ^ {ij}) \\ &= \exp(\sum \limits_{j \geq 1} \frac{F(x ^ j)}{j}) \end{aligned} \]

这样就便于用牛顿迭代解出 \(F(x)\) 了。一个问题是，\(F(x ^ k), k \geq 2\) 应该如何处理。考虑牛顿迭代的过程，在解出模 \(x ^ n\) 意义下的解 \(F_0(x)\) 后，我们希望求出模 \(x ^ {2n}\) 意义下的解。这意味着，\(F_0(x ^ k), k \geq 2\) 在模 \(x ^ {2n}\) 意义下其实是一个完全已知的多项式。令 \(P(x) = x \exp(\sum \limits_{i \geq 2} \dfrac{F(x ^ i)}{i})\)，则有 \(g(F(x)) = P(x)\exp(F(x)) - F(x)\)，\(g'(F(x)) = P(x)\exp(F(x)) - 1\)，直接使用牛顿迭代即可。

得到 \(F(x)\) 后，考虑怎么处理无根的情况。为简化分析，先不考虑树具有两个重心的情况。如果我们能在一棵树的重心处统计到这棵树，则可以保证每棵无标号无根树被统计恰好一次。一个点是一棵树的重心当且仅当其每个子树的大小都不超过 \(\lfloor \dfrac{n}{2} \rfloor\)。考虑怎么计算根非重心的情况数：钦定一棵子树的大小为 \(s > \lfloor \dfrac{n}{2} \rfloor\)，将其接到一棵 \(n - s\) 个点的有根树的根下即可。由于大小超过一半的子树至多只有一棵，故可以保证不重不漏。综上所述，无标号有根树且根为重心的情况数，即无标号无根树的数量，等于 \(f_n - \sum \limits_{s > \lfloor \frac{n}{2} \rfloor} f_sf_{n - s}\)。

当 \(n\) 为偶数且存在一大小为 \(\dfrac n 2\) 的子树时，树会有两个重心。考虑这样的树会被统计两次，额外减掉一次它们的贡献即可，为 \(\binom{f_{n / 2}}{2}\)。（如果两个 \(\dfrac{n}{2}\) 大小的子树形态相同，那么这棵树只会被统计一次，故不需要减掉这种情况的贡献。）

Code

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

constexpr int mod = 998244353, N = 1 << 20, invg = (mod + 1) / 3;
namespace basic {
  inline int add(int x, int y) {return (x + y >= mod ? x + y - mod : x + y);}
  inline int dec(int x, int y) {return (x - y < 0 ? x - y + mod : x - y);}
  inline void ad(int &x, int y) {x = add(x, y);}
  inline void de(int &x, int y) {x = dec(x, y);}

  inline int qpow(int a, int b) {
    int r = 1;
    while(b) {
      if(b & 1) r = 1ll * r * a % mod;
      a = 1ll * a * a % mod; b >>= 1;
    }
    return r;
  }
  inline int inv(int x) {return qpow(x, mod - 2);}

  int fac[N], ifac[N];
  inline void fac_init(int n = N - 1) {
    fac[0] = 1;
    for(int i = 1; i <= n; i++)
      fac[i] = 1ll * fac[i - 1] * i % mod;
    ifac[n] = inv(fac[n]);
    for(int i = n - 1; i >= 0; i--)
      ifac[i] = 1ll * ifac[i + 1] * (i + 1) % mod;
  }
  int invx[N];
  inline void inv_init(int n = N - 1) {
    invx[1] = 1;
    for(int i = 2; i <= n; i++)
      invx[i] = 1ll * (mod - mod / i) * invx[mod % i] % mod;
  }
  inline int binom(int n, int m) {
    if(n < m || m < 0) return 0;
    return 1ll * fac[n] * ifac[m] % mod * ifac[n - m] % mod;
  }

  int rev[N];
  inline void rev_init(int n) {
    for(int i = 1; i < n; i++) {
      rev[i] = (rev[i >> 1] >> 1) | (i & 1 ? n >> 1 : 0);
    }
  }
}
using namespace basic;

struct Poly {
  vector<int> a;
  inline int & operator [] (int x) {return a[x];}
  Poly() {}
  inline int size() {return a.size();}
  inline void resize(int n) {a.resize(n, 0);}
  explicit Poly(int n) {resize(n);}
  Poly(vector<int> b) {a = b;}
  inline void clear() {a.clear();}

  inline friend Poly operator + (Poly x, Poly y) {
    Poly ret(max(x.size(), y.size()));
    for(int i = 0; i < ret.size(); i++) {
      if(i < x.size()) {ad(ret[i], x[i]);}
      if(i < y.size()) {ad(ret[i], y[i]);}
    }
    return ret;
  }
  inline friend Poly operator - (Poly x, Poly y) {
    Poly ret(max(x.size(), y.size()));
    for(int i = 0; i < ret.size(); i++) {
      if(i < x.size()) {ad(ret[i], x[i]);}
      if(i < y.size()) {de(ret[i], y[i]);}
    }
    return ret;
  }

  inline void NTT(int type) {
    int n = size();
    rev_init(n);
    for(int i = 1; i < n; i++) {
      if(i < rev[i]) {
        swap(a[i], a[rev[i]]);
      }
    }
    static int gr[N]; gr[0] = 1;
    for(int d = 1; d < n; d <<= 1) {
      int gw = qpow(type == 1 ? 3 : invg, (mod - 1) / (d << 1));
      for(int i = 1; i < d; i++) {
        gr[i] = 1ll * gr[i - 1] * gw % mod;
      }
      for(int i = 0; i < n; i += d << 1) {
        for(int j = 0; j < d; j++) {
          int x = a[i + j], y = 1ll * gr[j] * a[i + j + d] % mod;
          a[i + j] = add(x, y);
          a[i + j + d] = dec(x, y);
        }
      }
    }
    if(type == -1) {
      for(int i = 0, invn = inv(n); i < n; i++) {
        a[i] = 1ll * a[i] * invn % mod;
      }
    }
  }
  inline friend Poly operator * (Poly x, Poly y) {
    int m = x.size() + y.size() - 1;
    int n = 1;
    while(n < m) {
      n <<= 1;
    }
    x.resize(n), y.resize(n);
    x.NTT(1), y.NTT(1);
    for(int i = 0; i < n; i++) {
      x[i] = 1ll * x[i] * y[i] % mod;
    }
    x.NTT(-1);
    x.resize(m);
    return x;
  }

  inline Poly operator += (Poly x) {
    return (*this) = (*this) + x;
  }
  inline Poly operator -= (Poly x) {
    return (*this) = (*this) - x;
  }
  inline Poly operator *= (Poly x) {
    return (*this) = (*this) * x;
  }

  inline Poly operator * (int x) {
    for(int i = 0; i < size(); i++) {
      a[i] = 1ll * a[i] * x % mod;
    }
    return (*this);
  }

  inline Poly Deriv() {
    Poly ret(size() - 1);
    for(int i = 1; i < size(); i++) {
      ret[i - 1] = 1ll * a[i] * i % mod;
    }
    return ret;
  }
  inline Poly Integ() {
    Poly ret(size() + 1); 
    for(int i = 0; i < size(); i++) {
      ret[i + 1] = 1ll * a[i] * inv(i + 1) % mod;
    }
    return ret;
  }

  // Implement followings by Newton's Method.
  // f(x) = f0(x) - g(f0(x)) / g'(f0(x))
  inline Poly Inv() {
    assert(a[0] != 0);
    Poly f0 = vector<int>{inv(a[0])};
    while(f0.size() < size()) {
      int n = f0.size() << 1;
      Poly h = (*this); h.resize(n), h.resize(n << 1); h.NTT(1);
      f0.resize(n << 1); f0.NTT(1);
      for(int i = 0; i < (n << 1); i++) {
        f0[i] = 1ll * f0[i] * dec(2, 1ll * f0[i] * h[i] % mod) % mod;
      }
      f0.NTT(-1), f0.resize(n);
    }
    f0.resize(size());
    return f0;
  }
  
  inline Poly Ln() {
    assert(a[0] == 1);
    Poly ret = ((*this).Deriv() * (*this).Inv()).Integ();
    ret.resize(size());
    return ret;
  }

  inline Poly Exp() {
    assert(a[0] == 0);
    Poly f0 = vector<int>{1};
    while(f0.size() < size()) {
      int n = f0.size() << 1; f0.resize(n);
      Poly t = f0.Ln();
      Poly h = (*this); h.resize(n);
      for(int i = 0; i < n; i++) {
        t[i] = dec(h[i], t[i]);
      }
      ad(t[0], 1);
      f0 = f0 * t; f0.resize(n);
    }
    f0.resize(size());
    return f0;
  }

  inline Poly Sqrt() {
    Poly f0 = vector<int>{1};
    while(f0.size() < size()) {
      int n = f0.size() << 1; f0.resize(n);
      Poly h = (*this); h.resize(n);
      f0 = f0 + h * f0.Inv(); f0.resize(n);
      int inv2 = inv(2);
      for(int i = 0; i < n; i++) {
        f0[i] = 1ll * f0[i] * inv2 % mod;
      }
    }
    return f0;
  }
};

inline Poly Euler(int m) {
  Poly f0 = vector<int>{0};
  while(f0.size() < m) {
    int n = f0.size() << 1; f0.resize(n);
    Poly P(n);
    for(int i = 2; i < n; i++) {
      int coef = inv(i);
      for(int j = 0; j <= (n - 1) / i; j++) {
        ad(P[i * j], 1ll * f0[j] * coef % mod);
      }
    }
    P = P.Exp();
    Poly C = vector<int>{0, 1} * P * f0.Exp();
    Poly p = f0 - C; p.resize(n);
    Poly q = vector<int>{1} - C; q.resize(n);
    f0 = f0 - p * q.Inv(); f0.resize(n);
  }
  f0.resize(m);
  return f0;
}

int n, H[N];

int main() {
  ios::sync_with_stdio(false);
  cin.tie(nullptr), cout.tie(nullptr);
  
  fac_init();

  cin >> n;
  Poly F = Euler(n + 1);
  int ans = F[n];
  for(int i = n / 2 + 1; i < n; i++) {
    de(ans, 1ll * F[i] * F[n - i] % mod);
  }
  if(n % 2 == 0) {
    de(ans, 1ll * F[n / 2] * dec(F[n / 2], 1) % mod * inv(2) % mod);
  }
  cout << ans << "\n";
}

写这个有点累的，咕咕中。

posted @ 2024-05-15 11:56 ChroneZ 阅读(281) 评论(0) 收藏举报

刷新页面返回顶部

chroneZ

(ง •_•)ง

多项式与生成函数学习笔记

FFT / NTT

Newton 迭代

一些数学基础及常用封闭形式

Lagrange 反演

做题记录

P5110 块速递推

P5320 [BJOI2019] 勘破神机

P3978 [TJOI2015] 概率论

CF438E The Child and Binary Tree

P2767 树的数量

P4841 [集训队作业2013] 城市规划

P5219 无聊的水题 I

P4233 射命丸文的笔记

P5900 无标号无根树计数

公告