[LibreOJ 6268]分拆数

danihao123 posted @ 2018年9月01日 17:25 in 题解 with tags loj NTT FFT 多项式牛顿迭代 多项式求逆 多项式ln 多项式exp , 1013 阅读
转载请注明出处:http://danihao123.is-programmer.com/

定义分拆数\(f(x)\)表示将\(x\)拆为若干正整数的本质不同方案数。对于\(i = 1\ldots n\),输出\(f(i)\)。

\(1\leq n\leq 10^5\)。

啊啊啊多项式牛迭又一题。

考虑定义分拆数数列的生成函数\(F(x) = \sum_{i = 1}^{+\infty} f(i)x^i\),那么显然有:

\[F(x) = \prod_{i = 1}^{+\infty} \frac{1}{1 - x^i}\]

两边取对数,得:

\[
\begin{aligned}
\ln F(x) &= \sum_{i = 1}^{+\infty}\ln(\frac{1}{1 - x^i})\\
&= -\sum_{i = 1}^{+\infty}\ln(1 - x^i)\\
&= \sum_{i = 1}^{+\infty} \sum_{j = 1}^{+\infty} \frac{x^{ij}}{j}
\end{aligned}
\]

其中一步应该用了\(\ln(1 - x)\)的麦克劳林展开。

不过注意到我们只需要求\(F(x)\bmod{x^{n + 1}}\),换言之我们也只需要求\(\ln F(x)\bmod{x^{n + 1}}\),因此我们只需要用上面的那个式子就能在\(O(n\ln n)\)的时间内搞出来\(\ln F(x)\bmod{x^{n + 1}}\),然后多项式exp一下就能得到\(F(x)\bmod{x^{n + 1}}\)了。

然后下面考虑怎么去弄多项式exp……套用多项式牛迭那一套理论,可以得知设\(G(t) = e^{P(x)} - t\)(\(P(x)\)为要求exp的多项式),则要解\(G(F(x))\equiv 0\pmod{x^{p}}\)这个方程。我们将两边取一下对数再乘下\(-1\),得到:

\[\ln F(x) - P(x)\equiv 0\pmod{x^{p}}\]

因此可以设\(G(t) = \ln(t) - P(x)\),这样的话可以得到倍增的式子是:

\[F(x) = F_0(x)(1 - \ln F_0(x) + P(x))\]

不过这就牵扯出来一个多项式ln的问题了……不过这个相当好弄,观察到\((\ln(A(x)))' = \frac{A'(x)}{A(x)}\),所以说\(\ln A(x) = \int\frac{A'(x)}{A(x)}\mathrm{d}x\),这个只要弄上多项式求导、积分和求逆就行了,都是传统艺能。

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <cassert>
#include <algorithm>
#include <functional>
#include <utility>
using ll = long long;
const ll ha = 998244353LL;
ll pow_mod(ll a, ll b) {
  ll ans = 1, res = a;
  while(b) {
    if(1LL & b) ans = ans * res % ha;
    res = res * res % ha; b >>= 1;
  }
  return ans;
}
ll inv(ll x) {
  return pow_mod(x, ha - 2LL);
}

const int maxn = 400050;
ll invp[maxn];
void process() {
  invp[1] = 1;
  for(int i = 2; i < maxn; i ++) {
    invp[i] = (ha - (ha / (ll)i));
    invp[i] = invp[i] * invp[ha % (ll)i] % ha;
  }
}

int flip(int bi, int x) {
  int ans = 0;
  for(int j = 0; j < bi; j ++) {
    if((1 << j) & x) {
      ans |= (1 << (bi - j - 1));
    }
  }
  return ans;
}
void NTT(ll *A, int bi, bool flag = false) {
  int n = 1 << bi;
  for(int i = 0; i < n; i ++) {
    int v = flip(bi, i);
    if(v < i) std::swap(A[v], A[i]);
  }
  for(int L = 1; L < n; L <<= 1) {
    ll xi_n = pow_mod(3LL, (ha - 1LL) / (ll(L << 1)));
    if(flag) xi_n = inv(xi_n);
    for(int i = 0; i < n; i += (L << 1)) {
      ll xi = 1LL;
      for(int j = i; j < i + L; j ++) {
        ll x = A[j], y = A[j + L];
        ll tmp = y * xi % ha;
        A[j] = (x + tmp) % ha;
        A[j + L] = (x - tmp + ha) % ha;
        xi = xi * xi_n % ha;
      }
    }
  }
  if(flag) {
    ll inv_n = inv(n);
    for(int i = 0; i < n; i ++) {
      A[i] = A[i] * inv_n % ha;
    }
  }
}

void poly_inv(int mod, ll *B, ll *BB) {
  if(mod == 1) {
    BB[0] = inv(B[0]);
  } else {
    poly_inv((mod + 1) >> 1, B, BB);
    int bi = 0, sz = 1;
    while(sz <= ((mod * 2) + 1)) {
      bi ++; sz <<= 1;
    }
    static ll tmp[maxn];
    std::copy(B, B + mod, tmp);
    std::fill(tmp + mod, tmp + sz, 0LL);
    NTT(tmp, bi); NTT(BB, bi);
    for(int i = 0; i < sz; i ++) {
      tmp[i] = (tmp[i] * BB[i]) % ha;
      tmp[i] = (tmp[i] * (ha - 1LL)) % ha;
      tmp[i] = (tmp[i] + 2LL) % ha;
      tmp[i] = (tmp[i] * BB[i]) % ha;
    }
    NTT(tmp, bi, true);
    std::copy(tmp, tmp + mod, BB);
    std::fill(BB + mod, BB + sz, 0LL);
  }
}
void poly_ln(int mod, ll *B, ll *BB) {
  ll *tmp = (ll*)calloc(maxn, sizeof(ll));
  poly_inv(mod, B, BB);
  for(int i = 0; i < mod - 1; i ++) tmp[i] = B[i + 1] * (ll(i + 1)) % ha;
  
  int bi = 0, sz = 1;
  while(sz <= ((mod * 2) + 1)) {
    bi ++; sz <<= 1;
  }
  NTT(tmp, bi); NTT(BB, bi);
  for(int i = 0; i < sz; i ++) BB[i] = BB[i] * tmp[i] % ha;
  NTT(BB, bi, true);
  for(int i = mod - 1; i > 0; i --) BB[i] = BB[i - 1] * invp[i] % ha;
  BB[0] = 0; std::fill(BB + mod, BB + sz, 0LL);
  free(tmp);
}
void poly_exp(int mod, ll *B, ll *BB) {
  if(mod == 1) {
    BB[0] = 1;
  } else {
    poly_exp((mod + 1) >> 1, B, BB);
    int bi = 0, sz = 1;
    while(sz <= ((mod * 2) + 1)) {
      bi ++; sz <<= 1;
    }
    
    std::fill(BB + mod, BB + sz, 0LL);
    ll *tmp = (ll*)calloc(maxn, sizeof(ll));
    poly_ln(mod, BB, tmp);
    for(int i = 0; i < mod; i ++) {
      tmp[i] = (B[i] - tmp[i] + ha) % ha;
    }
    tmp[0] = (tmp[0] + 1LL) % ha;
    NTT(tmp, bi); NTT(BB, bi);
    for(int i = 0; i < sz; i ++) BB[i] = tmp[i] * BB[i] % ha;
    NTT(BB, bi, true);
    
    std::fill(BB + mod, BB + sz, 0LL);
    free(tmp);
  }
}

ll A[maxn], B[maxn];
int main() {
  process();
  int n; scanf("%d", &n);
  for(int i = 1; i <= n; i ++) {
    for(int j = 1; j * i <= n; j ++) {
      A[j * i] = (A[j * i] + invp[j]) % ha;
    }
  }
  poly_exp(n + 1, A, B);
  for(int i = 1; i <= n; i ++) {
    printf("%lld\n", B[i]);
  }
  return 0;
}

登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter