CF947E Perpetual Subtraction

2018 年 04 月 28 日发布.

Description

有一个非负整数 $x$ 。你要执行 $m$ 次操作，每次操作是 x = randint(0, x) （ = 表示赋值， randint(0, x) 均匀随机地在 $[0, x]$ 中取一个整数）。

现在已知初始的 $x$ 会随机在 $[0, n]$ 取值，且取 $i$ 的概率是 $p_i$ ，求最后取到 $[0, n]$ 每个数的概率。

$$n\leqslant10^5, m\leqslant10^{18}$$

Solution

首先我们把概率都去掉：

如果我当前取到 $i\in[0, n]$ 的概率是 $f_i$ ，变换之后概率是 $f^*_i$ ，那么易知

$$f^*_i=\sum_{j=i}^n\frac{f_j}{j+1}$$

搞一个生成函数 $F(x)=\sum_{i=0}^n f_ix^i$ ，则

$$\begin{aligned} F^*(x)&=\sum_{i=0}^nx^i\sum_{j=0}^n\frac{f_j}{j+1}\\ &=\sum_{j=0}^n\frac{f_j}{j+1}\sum_{i=0}^j x^i\\ &=\sum_{j=0}^n\frac{f_j}{j+1}\frac{x^{j+1}-1}{x-1}\\ &=\frac1{x-1}\sum_{j=0}^n f_j\int_1^x t^j{\rm d}t\\ &=\frac{\int_1^x F(t){\rm d}t}{x-1} \end{aligned}$$

这个东西很不好：首先它积分下限是 $1$ 而不是 $0$ ；其次它除以 $x-1$ 而不是 $x$ 。

那么我们可以想到令 $G(x)=F(x+1)$ ，则

$$G^(x)=F^(x+1)=\frac{\int_1^{x+1}F(t){\rm d}t}{x+1-1}=\frac{\int_0^x G(t){\rm d}t}x$$

这个就很好：换成数列形式就是 $g^*_i=\frac{g_i}{i+1}$ 。

所以进行 $m$ 次操作就是 $g^*_i=\frac{g_i}{(i+1)^m}$ 。

于是我们只需要考虑如何由 $f$ 变换成 $g$ 和如何逆变换即可。

由于

$$\sum_{i=0}^ng_ix^i=\sum_{j=0}^nf_j(x+1)^j=\sum_{j=0}^nf_j\sum_{i=0}^j{j\choose i}x^i$$

所以

$$g_i=\sum_{j=i}^n{j\choose i}f_j$$

即

$$i!g_i=\sum_{j=i}^n\frac{j!f_j}{(j-i)!}$$

反过来卷积一下即可。同样也可得到

$$i!f_i=\sum_{j=1}^n\frac{j!f_j}{(-1)^{j-i}(j-i)!}$$

总时间复杂度 $O(n\log n)$ （快速幂常数小忽略）

Code

#include <algorithm>
#include <cstdio>
#include <cstring>

const int N = 100050;
const int mod = 998244353;
const int g = 3;

typedef long long LL;

LL pow_mod(LL a, LL p) {
  LL ans = 1;
  ((p %= (mod - 1)) += mod - 1) %= mod - 1;
  for (; p > 0; p >>= 1, (a *= a) %= mod)
    if (p & 1) (ans *= a) %= mod;
  return ans;
}

void NTT(LL *A, int len, int opt) {
  for (int i = 1, j = 0; i < len; ++i) {
    for (int k = len; ~j & k; j ^= (k >>= 1));
    if (i < j) std::swap(A[i], A[j]);
  }
  for (int h = 2; h <= len; h <<= 1) {
    LL wn = pow_mod(g, (mod - 1) / h * opt);
    for (int j = 0; j < len; j += h) {
      LL w = 1;
      for (int i = 0; i < h >> 1; ++i) {
        LL _t1 = A[i + j], _t2 = A[i + j + (h >> 1)] * w;
        A[i + j] = (_t1 + _t2) % mod;
        A[i + j + (h >> 1)] = (_t1 - _t2) % mod;
        (w *= wn) %= mod;
      }
    }
  }
  if (opt == -1)
    for (int i = 0, inv = -(mod - 1) / len; i < len; ++i)
      (A[i] *= inv) %= mod;
}

int n;
LL m, fac[N], ifac[N], F[N << 2], G[N << 2];

int main() {
  scanf("%d%lld", &n, &m);
  fac[0] = fac[1] = ifac[0] = ifac[1] = 1;
  for (int i = 2; i <= n; ++i)
    ifac[i] = -(mod / i) * ifac[mod % i] % mod;
  for (int i = 2; i <= n; ++i) {
    fac[i] = i * fac[i - 1] % mod;
    (ifac[i] *= ifac[i - 1]) %= mod;
  }
  for (int i = 0; i <= n; ++i) {
    scanf("%lld", &F[n - i]);
    (F[n - i] *= fac[i]) %= mod;
    G[i] = ifac[i];
  }
  int len = 1;
  while (len < (n + 1) * 2) len <<= 1;
  NTT(F, len, 1);
  NTT(G, len, 1);
  for (int i = 0; i < len; ++i)
    (G[i] *= F[i]) %= mod;
  NTT(G, len, -1);
  for (int i = 0; i <= n; ++i)
    (G[n - i] *= pow_mod(pow_mod(i + 1, m), mod - 2)) %= mod;
  for (int i = n + 1; i < len; ++i) G[i] = 0;
  NTT(G, len, 1);
  memset(F, 0, sizeof F);
  for (int i = 0; i <= n; ++i)
    F[i] = ifac[i] * (i & 1 ? -1 : 1);
  NTT(F, len, 1);
  for (int i = 0; i < len; ++i)
    (F[i] *= G[i]) %= mod;
  NTT(F, len, -1);
  for (int i = 0; i <= n; ++i)
    printf("%lld ", (F[n - i] * ifac[i] % mod + mod) % mod);
  return 0;
}