CF947E Perpetual Subtraction

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^8\]

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^j\sum_{j=0}^n\frac{f_j}{j+1}\\ &=\sum_{j=0}^n\frac{f_j}{j+1}\sum_{i=0}^j x^j\\ &=\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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
#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;
}