「洛谷P5408」第一类斯特林数·行 题解

题目

题意

求$\left[\begin{matrix}n\\i\end{matrix}\right],i=0,1,\dots,n$

其中$\left[\begin{matrix}n\\i\end{matrix}\right]$为第一类斯特林数,组合意义为将$n$个不同元素放入$i$个相同圆排列(非空)的方案数,有递推式$\left[\begin{matrix}n\\m\end{matrix}\right]=\left[\begin{matrix}n-1\\m-1\end{matrix}\right]+\left[\begin{matrix}n-1\\m\end{matrix}\right]\times(n-1)$

推导

一个结论:

$$\boxed{\sum\limits^n_{i=0}\left[\begin{matrix}n\\i\end{matrix}\right]x^i=x^{\overline{n}}}$$

$x^{\overline{n}}$为$x$的$n$次上升幂,等于$\prod\limits^{n-1}_{i=0}(x+i)$

考虑用数学归纳法证明:

  1. $n=1$时该结论显然成立

  2. 证明为$n$时成立则为$n+1$时也成立:

    $$\begin{aligned}x^{\overline{n+1}}&=x^{\overline{n}}\times(x+n)\\&=(x+n)\times\sum\limits^n_{i=0}\left[\begin{matrix}n\\i\end{matrix}\right]x^i\\&=\sum\limits^{n+1}_{i=1}\left[\begin{matrix}n\\i-1\end{matrix}\right]x^i+\sum\limits^n_{i=0}\left(n\times\left[\begin{matrix}n\\i\end{matrix}\right]x^i\right)\\&=\sum\limits^{n+1}_{i=0}\left(\left[\begin{matrix}n\\i-1\end{matrix}\right]+n\times\left[\begin{matrix}n\\i\end{matrix}\right]\right)x^i\\&=\sum\limits^{n+1}_{i=0}\left[\begin{matrix}n+1\\i\end{matrix}\right]x^i\end{aligned}$$

于是原命题成立

直接求上升幂是$O(n^2)$,考虑如何优化:

假设$n$是偶数:

$$x^{\overline{n}}=x^{\overline{\frac{n}{2}}}\times(x+\frac{n}{2})^{\overline{\frac{n}{2}}}$$

如果我们已经得到了前半部分对应的多项式$f(x)=x^{\overline{\frac{n}{2}}}=\sum\limits^{\frac{n}{2}}_{i=0}a_ix^i$

考虑如何通过前半部分的结果求出后半部分:

$$\begin{aligned}(x+\frac{n}{2})^{\overline{\frac{n}{2}}}&=f(x+\frac{n}{2})\\&=\sum\limits^{\frac{n}{2}}_{i=0}a_i(x+\frac{n}{2})^i\\&=\sum\limits^{\frac{n}{2}}_{i=0}\left(a_i\times\sum\limits^{i}_{j=0}\left(\binom{i}{j}\times x^j\times\left(\frac{n}{2}\right)^{i-j}\right)\right)\\&=\sum\limits^{\frac{n}{2}}_{j=0}\frac{1}{j!}\times\left(\sum\limits^{\frac{n}{2}}_{i=j}\left(i!\cdot a_i\right)\left(\frac{\left(\frac{n}{2}\right)^{i-j}}{\left(i-j\right)!}\right)\right)\times x^j\end{aligned}$$

到这里可能还不能发现它是卷积的形式,但是注意到$i-(i-j)=j$所以将后面一部分的$i$次项改为$-i$次项即可,但是负次数多项式不好实现,可以先将后半部分对应的多项式反过来(即将负次数多项式乘上$x^{\frac{n}{2}}$),即将$\sum\limits^{\frac{n}{2}}_{i=0}a_ii!x^i$和$\sum\limits^{\frac{n}{2}}_{i=0}\frac{\left(\frac{n}{2}\right)^i}{i!}x^{\frac{n}{2}-i}$乘起来的多项式第$i$项的系数为要求的多项式的第$i-\frac{n}{2}$项的系数

于是就可以由$x^{\overline{\frac{n}{2}}}$在$O(n\log_2n)$时间内求出$(x+\frac{n}{2})^{\overline{\frac{n}{2}}}$

而$x^{\overline{\frac{n}{2}}}$可以递归求解

若$n$是奇数,则$x^{\overline{n}}=x^{\overline{n-1}}\times(x+n-1)$

时间复杂度$T(n)=T(\frac{n}{2})+\Theta(nlog_2n)=\Theta(nlog_2n)$

代码

本题卡常,NTT的常数需要写得较小

const int N = 530000;
const int Mod = 167772161, g = 3;

class Polynomial {
public:
  int deg, *a;
  Polynomial() {
    a = new int[N];
    memset(a, 0, sizeof(int) * N);
  }
  void NTT(int, bool);
};
int r[N];
Polynomial solve(int);

int fac[N], ifac[N];
int omega[26];
int siz;

int main () {
  int n;
  read(n);
  int lim = 1;
  while (lim <= (n << 1)) {
    lim <<= 1;
  }
  siz = __builtin_ctz(lim);
  omega[25] = pow(g, (Mod - 1) >> 25, Mod);
  for (int i = 24; i >= 0; --i) {
    omega[i] = omega[i + 1] * 1ll * omega[i + 1] % Mod;
  }
  for (int i = 0; i < lim; ++i) {
    r[i] = (r[i >> 1] >> 1) | ((i & 1) * (lim >> 1));
  }
  fac[0] = 1;
  for (int i = 1; i <= lim; ++i) {
    fac[i] = fac[i - 1] * 1ll * i % Mod;
  }
  ifac[lim] = pow(fac[lim], Mod - 2, Mod);
  for (int i = lim - 1; i >= 0; --i) {
    ifac[i] = ifac[i + 1] * 1ll * (i + 1) % Mod;
  }
  Polynomial ans = solve(n);
  for (int i = 0; i <= n; ++i) {
    write(ans.a[i]), SP;
  }
  EL;
  return 0;
}

void Polynomial::NTT(int lim, bool opt) {
  if (opt) {
    std::reverse(a + 1, a + lim);
  }
  int tp = siz - __builtin_ctz(lim);
  for (int i = 0; i < lim; ++i) {
    if (i < (r[i] >> tp)) {
      std::swap(a[i], a[r[i] >> tp]);
    }
  }
  for (int l = 1, t = 1; l < lim; l <<= 1, ++t) {
    int gw = omega[t];
    for (int i = 0; i < lim; i += (l << 1)) {
      for (int j = 0, gn = 1; j < l; ++j, gn = gn * 1ll * gw % Mod) {
        int x = a[i + j], y = a[i + j + l] * 1ll * gn % Mod;
        a[i + j] = (x + y) % Mod;
        a[i + j + l] = (x - y + Mod) % Mod;
      }
    }
  }
  if (opt) {
    int iv = pow(lim, Mod - 2, Mod);
    for (int i = 0; i < lim; ++i) {
      a[i] = a[i] * 1ll * iv % Mod;
    }
  }
}
void calc(const Polynomial &f, int n, int lim, Polynomial &g) {
  Polynomial h;
  memset(g.a + n + 1, 0, sizeof(int) * (lim - n));
  g.a[0] = 1, h.a[0] = f.a[0];
  for (int i = 1; i <= n; ++i) {
    g.a[i] = g.a[i - 1] * 1ll * n % Mod;
    h.a[i] = f.a[i] * 1ll * fac[i] % Mod;
  }
  for (int i = 1; i <= n; ++i) {
    g.a[i] = g.a[i] * 1ll * ifac[i] % Mod;
  }
  std::reverse(g.a, g.a + n + 1);
  g.NTT(lim, false), h.NTT(lim, false);
  for (int i = 0; i < lim; ++i) {
    g.a[i] = g.a[i] * 1ll * h.a[i] % Mod;
  }
  g.NTT(lim, true);
  for (int i = 0; i <= n; ++i) {
    g.a[i] = g.a[i + n] * 1ll * ifac[i] % Mod;
  }
  memset(g.a + n + 1, 0, sizeof(int) * (lim - n));
}
Polynomial solve(int n) {
  std::stack<bool> S;
  while (n) {
    S.push(n & 1);
    n >>= 1;
  }
  int lim = 1;
  Polynomial A, B;
  n = A.a[1] = S.top();
  S.pop();
  while (!S.empty()) {
    while (lim <= (n << 1)) {
      lim <<= 1;
    }
    calc(A, n, lim, B);
    A.NTT(lim, false), B.NTT(lim, false);
    for (int i = 0; i < lim; ++i) {
      A.a[i] = A.a[i] * 1ll * B.a[i] % Mod;
    }
    A.NTT(lim, true);
    n <<= 1;
    if (!S.top()) {
      S.pop();
      continue;
    }
    S.pop();
    A.a[n + 1] = A.a[n];
    for (int i = n; i > 0; --i) {
      A.a[i] = (A.a[i - 1] + A.a[i] * 1ll * n % Mod) % Mod;
    }
    A.a[0] = A.a[0] * 1ll * n % Mod;
    ++n;
  }
  A.deg = n;
  return A;
}

本博客采用知识共享署名-相同方式共享 4.0 国际许可协议进行许可
本文链接:https://blog.seniorious.cc/2020/luogu-5408/