「VK Cup 2018 - Round 1」E. Perpetual Subtraction

Posted on | Edited on | In , ,

Problem

Description

有一个初值以 p_{i} 概率取值为 i\left(0 \le i \le n\right) 的离散随机变量 x ，定义对其的一次操作为将其等概率随机赋值为 \left[0, x\right] 中的一个整数。
现给出 n, m, p_{i} ，要求对于 \left[0, n\right] 中的每个 i 求出对 x 进行 m 次操作后其取值为 i 的概率。

答案对 998244353 取模。

Constraints

1 \le n \le 10^{5}, 0 \le m \le 10^{18}, \sum_{i = 0} ^ {n} p_{i} \equiv 1 \pmod{998244353}

Solution

Analysis

对于单次变换，设变换前 x 的 PGF 为 f(x) ，变换后 x 的 PGF 为 f^{*}(x) ，稍加推导可得：

\begin{aligned} [x^{i}] f^{*}(x) & = \sum_{j = i}^{n} \frac{[x^{j}] f(x)}{j + 1} \\ f^{*}(x) &= \sum_{i = 0}^{n} x^{i} \sum_{j = i}^{n} \frac{[x^{j}] f(x)}{j + 1} \\ &= \sum_{i = 0}^{n} \frac{[x^{i}] f(x)}{i + 1} \sum_{j = 0}^{i} x^{j} \\ &= \sum_{i = 0}^{n} \frac{[x^{i}] f(x)}{i + 1} \frac{x^{i + 1} - 1}{x - 1} \\ &= \frac{1}{x - 1} \sum_{i = 0}^{n} [x^{i}] f(x) \int_{1}^{x} t^{i} \mathrm{d} t \\ &= \frac{1}{x - 1} \int_{1}^{x} f(t) \mathrm{d} t \end{aligned}

但是我们发现这个东西的积分下限是 1 而不是 0 ，分母是 x - 1 而不是 x ，这就很难办了。
于是考虑令 g(x) = f\left(x + 1\right), g^{*}(x) = f^{*}\left(x + 1\right) ，有：

g^{*}(x) = f^{*}\left(x + 1\right) = \frac{1}{x + 1 - 1} \int_{1}^{x + 1} f(t) \mathrm{d} t = \frac{1}{x} \int_{0}^{x} g(t) \mathrm{d} t

比较系数可得：

[x^{i}] g^{*}(x) = \frac{[x^{i}] g(x)}{i + 1}

所以 m 次变换后就有：

[x^{i}] g^{*}(x) = \frac{[x^{i}] g(x)}{\left(i + 1\right)^{m}}

于是现在只需要考虑从 f 到 g 的变换及逆变换即可。

\begin{aligned} g(x) &= f\left(x + 1\right) \\ \sum_{i = 0}^{n} [x^{i}] g(x) x^{i} &= \sum_{i = 0}^{n} [x^{i}] f(x) \left(x + 1\right)^{i} \\ \sum_{i = 0}^{n} [x^{i}] g(x) x^{i} &= \sum_{i = 0}^{n} [x^{i}] f(x) \sum_{k = 0}^{i} \binom{i}{k} x^{k} \\ [x^{i}] g(x) &= \sum_{k = 0}^{n} \binom{k}{i} [x^{k}] f(x) \\ i![x^{i}] g(x) &= \sum_{k = 0}^{n} \frac{k! [x^{k}] f(x)}{\left(k - i\right)!} \end{aligned}

可以卷积解决。逆变换同理：

\begin{aligned} f(x) &= g\left(x - 1\right) \\ \sum_{i = 0}^{n} [x^{i}] f(x) x^{i} &= \sum_{i = 0}^{n} [x^{i}] g(x) \left(x - 1\right)^{i} \\ \sum_{i = 0}^{n} [x^{i}] f(x) x^{i} &= \sum_{i = 0}^{n} [x^{i}] g(x) \sum_{k = 0}^{i} (-1)^{i - k} \binom{i}{k} x^{k} \\ [x^{i}] f(x) &= \sum_{k = 0}^{n} (-1)^{k - i} \binom{k}{i} [x^{k}] g(x) \\ i![x^{i}] f(x) &= \sum_{k = 0}^{n} \frac{(-1)^{k - i} k! [x^{k}] g(x)}{\left(k - i\right)!} \end{aligned}

时间复杂度 \displaystyle{O\left(n\log{n} + \frac{n}{\ln{n}}\log{m}\right)} 。

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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
#include <cstdio>
#include <bitset>
#include <algorithm>

using i64 = long long;

constexpr int maxn = 100000;
constexpr int mod = 998244353;

template <class x_tp, class y_tp>
inline void inc(x_tp &x, const y_tp &y) {
	x += y; (mod <= x) && (x -= mod);
}
template <class x_tp, class y_tp>
inline int sub(const x_tp &x, const y_tp &y) {
	return x < y ? x - y + mod : (x - y);
}
template <typename _Tp>
inline int fpow(_Tp v, int n) {
	_Tp pw = 1;
	for (; n; n >>= 1, v = (i64)v * v % mod)
		if (n & 1) pw = (i64)pw * v % mod;
	return pw;
}
template <typename _Tp>
inline void swap(_Tp &x, _Tp &y) {
	_Tp z = x; x = y; y = z;
}

namespace IOManager {
	constexpr int FILESZ = 131072;
	char buf[FILESZ];
	const char *ibuf = buf, *tbuf = buf;

	struct IOManager {
		inline char gc() {
			return (ibuf == tbuf) && (tbuf = (ibuf = buf) + fread(buf, 1, FILESZ, stdin), ibuf == tbuf) ? EOF : *ibuf++;
		}

		template <typename _Tp>
			inline operator _Tp() {
				_Tp s = 0u; char c = gc();
				for (; c < 48; c = gc());
				for (; c > 47; c = gc())
					s = (_Tp)(s * 10u + c - 48u);
				return s;
			}
	};
} IOManager::IOManager io;

namespace poly {

constexpr int maxn = 262144;
constexpr int g = 3;

using poly_t = int[maxn];
using poly = int *const;

inline int calcpw2(const int &n) {
	for (int t = 1; ; t <<= 1)
		if (n < t) return t;
}

poly_t wn, iwn, cw;

void polyinit() {
	const int wnv = fpow(g, (mod - 1) / maxn);
	wn[0] = iwn[0] = 1;
	for (int i = 1; i != maxn; ++i)
		wn[i] = (i64)wn[i - 1] * wnv % mod;
	std::reverse_copy(wn + 1, wn + maxn, iwn + 1);
}

void DFT(poly &a, const int &n) {
	for (int i = 0, j = 0; i != n; ++i) {
		if (i < j) swap(a[i], a[j]);
		for (int k = n >> 1; (j ^= k) < k; k >>= 1);
	}
	poly endpos = a + n, wendpos = wn + maxn; cw[0] = 1;
	for (int l = 1, tp = maxn >> 1; l != n; l <<= 1, tp >>= 1) {
		for (int *w = wn + tp, *i = cw; w != wendpos; w += tp)
			*++i = *w;
		for (int *i = a, z; i != endpos; i += l + l)
			for (int j = 0; j != l; ++j)
				z = (i64)i[j + l] * cw[j] % mod,
				i[j + l] = sub(i[j], z),
				inc(i[j], z);
	}
}

void IDFT(poly &a, const int &n) {
	for (int i = 0, j = 0; i != n; ++i) {
		if (i < j) swap(a[i], a[j]);
		for (int k = n >> 1; (j ^= k) < k; k >>= 1);
	}
	poly endpos = a + n, wendpos = iwn + maxn; cw[0] = 1;
	for (int l = 1, tp = maxn >> 1; l != n; l <<= 1, tp >>= 1) {
		for (int *w = iwn + tp, *i = cw; w != wendpos; w += tp)
			*++i = *w;
		for (int *i = a, z; i != endpos; i += l + l)
			for (int j = 0; j != l; ++j)
				z = (i64)i[j + l] * cw[j] % mod,
				i[j + l] = sub(i[j], z),
				inc(i[j], z);
	}
	const int invn = fpow(n, mod - 2);
	for (int *i = a; i != endpos; ++i)
		*i = (i64)*i * invn % mod;
}

}

int fac[maxn + 1],
	ifac[maxn + 1];
std::bitset<maxn + 1> notp;
int pri[maxn / 7],
	ipw[maxn + 1];
poly::poly_t F, G;

int main() {
	poly::polyinit();
	const int n = io, n1 = n + 1, deg = poly::calcpw2(n1 + n1); const i64 m = io;
	const int t = (m % (mod - 1)) * (mod - 2ll) % (mod - 1);

	fac[0] = fac[1] = ifac[0] = ifac[1] = 1;
	for (int i = 2; i != n1; ++i)
		fac[i] = (i64)fac[i - 1] * i % mod;
	ifac[n] = fpow(fac[n], mod - 2);
	for (int i = n; 2 < i; --i)
		ifac[i - 1] = (i64)ifac[i] * i % mod;

	/* pre-calculate x^{-m} by linear sieve */
	int cnt = 0;
	ipw[1] = 1;
	for (int i = 2; i <= n1; ++i) {
		if (!notp[i])
			pri[++cnt] = i, ipw[i] = fpow(i, t);
		for (int j = 1, v; j <= cnt && (v = i * pri[j]) <= n1; ++j) {
			notp[v] = true;
			ipw[v] = (i64)ipw[i] * ipw[pri[j]] % mod;
			if (i % pri[j] == 0)
				break;
		}
	}

	for (int i = 0; i != n1; ++i) {
		F[n - i] = (i64)fac[i] * (int)io % mod;
		G[i] = ifac[i];
	}

	poly::DFT(F, deg); poly::DFT(G, deg);
	for (int i = 0; i != deg; ++i)
		F[i] = (i64)F[i] * G[i] % mod;
	poly::IDFT(F, deg);

	for (int i = 0; i != n1; ++i)
		F[n - i] = (i64)F[n - i] * ipw[i + 1] % mod;
	std::fill(F + n1, F + deg, 0);
	for (int i = 0; i <= n; i += 2)
		G[i] = ifac[i],
		G[i + 1] = mod - ifac[i + 1];
	std::fill(G + n1, G + deg, 0);

	poly::DFT(F, deg); poly::DFT(G, deg);
	for (int i = 0; i != deg; ++i)
		F[i] = (i64)F[i] * G[i] % mod;
	poly::IDFT(F, deg);

	for (int i = 0; i != n1; ++i)
		printf("%lld ", (i64)F[n - i] * ifac[i] % mod);
	return 0;
}