题解 P5408 【第一类斯特林数·行】

· · 题解

题解 P5408 【第一类斯特林数·行】

\mathfrak{View\space it\space on\space my\space Blog}

题意

求\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)

考虑用数学归纳法证明：

  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;
}