UOJ 593 新年的军队

发表于 | 分类于 |

新年的军队题解人话重置版！

首先作一步相当重要的转化：
对所有满足恰有 \(m\) 个位置 \(i\) 满足 \(\alpha_i > \alpha_{i+1}\) 的实数序列 \(\alpha_1,\alpha_2,\dots,\alpha_n\)，计算 \(\alpha_k\) 的排名的分布。其中 \(\alpha_i\) 在 \([0,1]\) 上均匀分布。

我们知道若干随机实数相等的概率是 \(0\)，所以这里其实区间的开闭也不影响结果。那么若干个实数就可以排出排名，并且显然每种排列出现的概率是相等的。
所以这个问题的结果和原问题是等价的。

然后稍微插入一下对概率密度函数的普及。
（我为啥要开个引用块？为了让会的神仙方便区分我的瞎扯和正经内容）
假装我们有一个在 \([0,4]\) 上均匀分布的实数变量，考虑一个非常憨批的问题：其落在区间 \([1,2]\) 内的概率。
我们发现这就是长度之比的事。
不过这个时候我们假装有一个在函数图像上的长方形，横坐标的范围是 \([0,4]\)，纵坐标都是 \(\frac14\)，然后我们在 \([1,2]\) 上作定积分就得到了它的概率。
这个 \(\frac14\) 其实就是个类似概率密度函数的东西。它用来描述连续型随机变量的输出值。我们想要知道它落在某个区间上的概率，只需要在这个区间上对概率密度函数作积分即可。

我们考虑用概率密度来考察特定 \(\alpha_k\) 的分布。若 \(\alpha_k\) 的排名为 \(j\)，则其贡献的概率密度即为 \[ \frac{ x^{j-1}(1-x)^{n-j} }{ (j-1)!(n-j)! } \]

这是因为当 \(\alpha_k = x\) 时，有 \(j-1\) 个数小于其的概率各是 \(x\)，然后它们有 \((j-1)!\) 种排列方法，大于其的情况类似。
这个地方我们需要考虑除以 \((j-1)!\) 的原因在于这些情况是本质相同的。

这样的话，刻画两个数之间的大小关系就变得容易了。因为我们可以直接在 \([0,x]\) 或 \([x,1]\) 上积分。

接下来我们建立所求的东西和概率密度函数的关系。
考虑令 \(a_i\) 为 \(p_i=i+1\) 的方案数，设 \(f(x)\) 为 \(\alpha_k\) 最终的概率密度函数，那么我们有 \[ f(x) = \sum\limits_{i=0}^{n-1} a_i \frac{ x^i (1-x)^{n-1-i} }{i!(n-1-i)!} \]

容易反演 \[ \sum\limits_{i=0}^{n-1} a_i \frac{x^i}{i!(n-1-i)!} = \sum\limits_{i=0}^{n-1} f_i x^i (1+x)^{n-1-i} \]

具体地 \[ \begin{aligned} \sum\limits_{i=0}^{n-1} f_i x^i (1+x)^{n-1-i} &= \sum\limits_{i=0}^{n-1} (1+x)^i f_{n-1-i} x^{n-1-i} \\ &= \sum\limits_{i=0}^{n-1} f_{n-1-i} \sum\limits_{j=0}^i \binom ij x^{n-1-(i-j)} \\ &= \sum\limits_{i=0}^{n-1} f_{n-1-i} \sum\limits_{j=0}^i \binom ij x^{n-1-j} \\ &= \sum\limits_{j=0}^{n-1} x^{n-1-j} \sum\limits_{i=j}^{n-1} \binom ij f_{n-1-i} \end{aligned} \]

接下来正式开始处理原问题。
考虑在概率密度函数中加入两个元区分大于号和小于号：设 \(F(u,v,t)\) 表示排列末尾的分布，其中 \(u\) 记录小于号，\(v\) 记录大于号。 \[ \def\d{ {\rm d} } \def\e{ {\rm e} } \def\D{ {\rm D} } F(t) = 1 + u \int_0^t F(\tau) \d \tau + v \int_t^1 F(\tau) \d \tau \]

因此 \[ F' = (u - v)F \]

显然可以待定系数得到其解 \[ F = C\e^{(u-v)t} \]

代入 \(t=0,1\) 可以确定 \[ F = \frac{ (u - v)\e^{(u - v)t} }{ u - v\e^{u-v} } \]

接下来用 \(x\) 元计量边数，用 \(y\) 元计量大于号的数量。
对于 \(p_k\) 的左侧，我们令 \(u=x,v=xy\)；对于右侧我们令 \(u=xy,v=x\)。那么 \[ \begin{aligned} L &= \frac{ (1-y)\e^{x(1-y)t} }{ 1-y\e^{x(1-y)} } \\ R &= \frac{ (1-y)\e^{x(1-y)(1-t)} }{ 1-y\e^{x(1-y)} } \end{aligned} \]

接下来有两条路。

离散微积分

我们直接考虑提取其系数。
我们知道 \(([x^{n_1}] L) ([x^{n_2} R])\) 就等于 \[ \frac1{n_1!n_2!} (1-y)^{n_1+n_2+2}\left(\sum\limits_{j\ge 0} y^j (t+j)^{n_1}\right) \left(\sum\limits_{j\ge 0} y^j ((j+1)-t)^{n_2}\right) \]

事实上在这个地方我们容易求出结果（关于 \(t\) 的多项式）在 \(0,1,\dots,n\) 处的点值，能否通过取决于插值的常数。
不过这个做法不太有趣，此处略过。

去掉常数 \(\frac{(-1)^{n_2}}{n_1!n_2!}\)，我们把它变成 \[ \begin{aligned} & \quad\; [y^m] (1-y)^{n+1} \sum\limits_{i\ge 0} \sum\limits_{j\ge 0} y^{i+j} (t+i)^{n_1} (t-(j+1))^{n_2} \\ &= \sum\limits_{s=0}^m [y^m](1-y)^{n+1} y^s \sum\limits_{i=0}^s (t+i)^{n_1} (t+i-(s+1))^{n_2} \end{aligned} \]

我们将其拆成正方向的无穷求和加上一个差分 \(\Delta_{s+1} f(t) = f(t) - f(t+s+1)\)。
设 \(f_s = [y^m](1-y)^{n+1} y^s\)。
\[ \sum\limits_{s=0}^m f_s \Delta_{s+1} \left(\sum\limits_{i\ge 0} (t+i)^{n_1} (t+i-(s+1))^{n_2}\right) \]

这个地方我们考虑把差分换到无穷求和里面，看起来我们求和的东西是不变的，但是实际上无穷求和各项的关系并不这么单纯，正像你作不定积分后会多出一个常数。
所以实际上变成了 \[ \sum\limits_{s=0}^m f_s \left(\sum\limits_{i\ge 0} \Delta_{s+1} (t+i)^{n_1} (t+i-(s+1))^{n_2}\right) + C \]

不管常数项，即 \[ \sum\limits_{i\ge 0} \left( \sum\limits_{s=0}^m f_s [(t+i)^{n_1} (t+i-(s+1))^{n_2}-(t+i+(s+1))^{n_1} (t+i)^{n_2}]\right) \]

但是这个 \(C\) 一看就不像是人工分析可以搞出来的东西……
不过回到原先的概率密度函数，我们会发现令 \(f_0\) 加上 \(C\) 会导致所有 \(a_i\) 之和加上 \(C\) 倍的某个常数（而这个常数是容易算出的），而我们知道 \(a_i\) 之和应当等于 \(\genfrac{\langle}{\rangle}{0pt}{}nm\)。

令 \(\D\) 是求导算子，我们知道 \(f(t+c) = \e^{c \D} f(t)\)，因而不妨设 \(F(x) = \sum\limits_{s=0}^m f_s x^{s+1}\)，就有 \[ \sum\limits_{i\ge 0} \left((t+i)^{n_1} F(\e^{-\D}) (t+i)^{n_2}-(t+i)^{n_2} F(\e^{\D}) (t+i)^{n_1}\right) \]

我们知道我们作一次位移是 \(\e^{\D}\)，正向差分是 \(1-\e^{\D}\)，它的逆运算就是 \(\frac1{ 1-\e^{\D} }\)，然后为了能够求逆我们恰好就凑出来了一个伯努利数的形式 \[ \int B(\D) \left((t+i)^{n_1} F(\e^{-\D}) (t+i)^{n_2}-(t+i)^{n_2} F(\e^{\D}) (t+i)^{n_1}\right), \quad B(t) = \frac{t}{1-\e^t} \]

接下来问题落在计算 \(G(x)\)。我们知道 \[ F(x) = \sum\limits_{i=0}^m \binom{n+1}{m-i} (-1)^{m-i} x^{i+1} \]

不妨设 \(F = xf\)，令 \(C = (n-m+1)f_0 = (n-m+1) \binom{n+1}m (-1)^m\)，那么我们知道 \[ \begin{aligned} (n-m+1+i)f_i &= -(m-i+1)f_{i-1} + [i=0]C \\ (n-m+1)f + xf' &= -mxf + x^2f' + C \end{aligned} \]

令 \(G(x) = F(\e^x) = \e^x f(\e^x) = \e^x g(x)\)，那么 \[ \begin{aligned} (n-m+1)f(\e^x) + \e^x f'(\e^x) &= -m\e^xf(\e^x) + \e^{2x}f'(\e^x) + C \\ (n-m+1)g + g' &= -m \e^x g + \e^x g' + C \\ (n-m+1+m\e^x)g + (1-\e^x)g' &= C \\ ((n-m)\e^{-x}+m+1)G + (\e^{-x}-1)G' &= C \end{aligned} \]

提取 \([x^i]\) 得 \[ \begin{aligned} (n-m)\sum\limits_{j=1}^i \frac{(-1)^j}{j!} G_{i-j} &+ (n+1) G_i \\ +\sum\limits_{j=1}^i \frac{(-1)^j}{j!} (i-j+1) G_{i-j+1} &= [i=0]C \end{aligned} \]

整理得 \[ \begin{aligned} (n-i+1)G_i &= [i=0]C \\ -(n-m)\sum\limits_{j=0}^{i-1} G_j \frac{ (-1)^{i-j} }{(i-j)!} &-\sum\limits_{j=1}^{i-1} j G_j \frac{ (-1)^{i-j+1} }{(i-j+1)!} \end{aligned} \]

四元生成函数

我们先不提取系数，而是先综合考虑整体的生成函数。

左侧和右侧的 \(x\) 元分别用 \(p,q\) 代替，令 \(n_1 = k-1,n_2 = n-k\)，那么问题变为提取 \[ [y^m p^{n_1} q^{n_2}] (1-y)^2 \frac{ \e^{p(1-y)t + q(1-y)(1-t)} }{(1-y\e^{p(1-y)})(1-y\e^{q(1-y)})} \]

对 \(p,q\) 换元，提出 \((1-y)\)，则 \[ (1-y)^{n+1} \frac{ \e^{pt + q(1-t)} }{(1-y\e^p)(1-y\e^q)} \]

作部分分式，则 \[ (1-y)^{n+1} \e^{pt+q(1-t)} \frac 1{\e^p - \e^q}\left(\frac{\e^p}{1-y\e^p} - \frac{\e^q}{1-y\e^q}\right) \]

这看起来不太能计算，但是某些奇思妙想让我们试图对 \(t\) 求导， \[ \begin{aligned} & \quad\; (1-y)^{n+1} \e^{pt+q(1-t)} \frac{p-q}{\e^p-\e^q}\left(\frac{\e^p}{1-y\e^p} - \frac{\e^q}{1-y\e^q}\right) \\ & = (1-y)^{n+1} \e^{pt+q(1-t)} \frac{ (p-q)\e^{-q} }{\e^{p-q}-1}\left(\frac{\e^p}{1-y\e^p} - \frac{\e^q}{1-y\e^q}\right) \\ & = (1-y)^{n+1} \e^{(p-q)t} B(p-q) \left(\frac{\e^p}{1-y\e^p} - \frac{\e^q}{1-y\e^q}\right) \end{aligned} \]

其中 \(B(x) = \frac x{\e^x-1}\)。

损失的常数可以使用上面提到的方法类似处理。

设 \(F(x) = \sum\limits_{i\ge 0} [y^{m-i}] (1-y)^{n+1} x^{i+1},G(x) = F(\e^x)\)，提取 \(y^m\) 的系数 \[ \e^{(p-q)t} B(p-q) (G(p)-G(q)) \]

分别考虑 \(\e^{(p-q)t} B(p-q) G(p),\e^{(p-q)t} B(p-q) G(q)\)。
以前者为例，我们考虑换元：令 \(pr = q\)，则变为 \[ \e^{p(1-r)t} B(p(1-r)) G(p) \]

对其提取 \(p^{n-1} r^{n_2}\) 的系数，有 \[ \begin{aligned} & \quad\; [p^{n-1} r^{n_2}] \e^{p(1-r)t} B(p(1-r)) G(p) \\ & = \sum\limits_{j \ge 0} ([p^j]B(p) \e^{pt})([r^{n_2}](1-r)^j)([p^{n-1-j}]G(p)) \end{aligned} \]

设 \(P(p) = \sum\limits_{j\ge 0} p^{n-1-j} ([r^{n_2}](1-r)^j)([p^{n-1-j}]G(p))\)，那么 \[ [p^{n-1}] P(p) B(p) \e^{pt} = \sum\limits_{j \ge 0} \frac{t^j}{j!} [p^{n-1-j}] P(p) B(p) \]

类似地，设 \(Q(q) = \sum\limits_{j\ge 0} q^{n-1-j} ([r^{n_1}](r-1)^j)([q^{n-1-j}] G(q))\)，就有 \[ \sum\limits_{j\ge 0} \frac{t^j}{j!} [q^{n-1-j}] Q(q) B(q) \]

\(G\) 的求算已经在上文提到过。

时间复杂度 \(O(n \log^2 n)\) 或 \(O\left(\frac{n \log^2 n}{\log \log n}\right)\)。

代码： 

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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
#include <cstdio>
#include <vector>
#include <cstring>
#include <algorithm>
#define add(a,b) (a + b >= mod ? a + b - mod : a + b)
#define dec(a,b) (a < b ? a - b + mod : a - b)
using namespace std;
const int N = 5e5;
const int mod = 998244353;
int n,m,k;
int n1,n2;
inline int fpow(int a,int b)
{
    int ret = 1;
    for(;b;b >>= 1)
        (b & 1) && (ret = (long long)ret * a % mod),a = (long long)a * a % mod;
    return ret;
}
namespace Poly
{
    const int LG = 19;
    const int N = 1 << LG + 1;
    const int G = 3;
    int lg2[N + 5];
    int rev[N + 5],fac[N + 5],ifac[N + 5],inv[N + 5];
    int rt[N + 5];
    inline void init()
    {
        for(register int i = 2;i <= N;++i)
            lg2[i] = lg2[i >> 1] + 1;
        rt[0] = 1,rt[1 << LG] = fpow(G,(mod - 1) >> LG + 2);
        for(register int i = LG;i;--i)
            rt[1 << i - 1] = (long long)rt[1 << i] * rt[1 << i] % mod;
        for(register int i = 1;i < N;++i)
            rt[i] = (long long)rt[i & i - 1] * rt[i & -i] % mod;
        fac[0] = 1;
        for(register int i = 1;i <= N;++i)
            fac[i] = (long long)fac[i - 1] * i % mod;
        ifac[N] = fpow(fac[N],mod - 2);
        for(register int i = N;i;--i)
            ifac[i - 1] = (long long)ifac[i] * i % mod;
        for(register int i = 1;i <= N;++i)
            inv[i] = (long long)ifac[i] * fac[i - 1] % mod;
    }
    struct poly
    {
        vector<int> a;
        inline poly(int x = 0)
        {
            x && (a.push_back(x),1);
        }
        inline poly(const vector<int> &o)
        {
            a = o,shrink();
        }
        inline poly(const poly &o)
        {
            a = o.a,shrink();
        }
        inline void shrink()
        {
            for(;!a.empty() && !a.back();a.pop_back());
        }
        inline int size() const
        {
            return a.size();
        }
        inline void resize(int x)
        {
            a.resize(x);
        }
        inline int operator[](int x) const
        {
            if(x < 0 || x >= size())
                return 0;
            return a[x];
        }
        inline void clear()
        {
            vector<int>().swap(a);
        }
        inline poly rever() const
        {
            return poly(vector<int>(a.rbegin(),a.rend()));
        }
        inline void dif()
		{
            int n = size();
            for(register int i = 0,len = n >> 1;len;++i,len >>= 1)
                for(register int j = 0,*w = rt;j < n;j += len << 1,++w)
                    for(register int k = j,R;k < j + len;++k)
                        R = (long long)*w * a[k + len] % mod,
                        a[k + len] = dec(a[k],R),
                        a[k] = add(a[k],R);
        }
        inline void dit()
        {
            int n = size();
            for(register int i = 0,len = 1;len < n;++i,len <<= 1)
                for(register int j = 0,*w = rt;j < n;j += len << 1,++w)
                    for(register int k = j,R;k < j + len;++k)
                        R = add(a[k],a[k + len]),
                        a[k + len] = (long long)(a[k] - a[k + len] + mod) * *w % mod,
                        a[k] = R;
            reverse(a.begin() + 1,a.end());
            for(register int i = 0;i < n;++i)
                a[i] = (long long)a[i] * inv[n] % mod;
        }
        inline void ntt(int type = 1)
        {
            type == 1 ? dif() : dit();
        }
        friend inline poly operator+(const poly &a,const poly &b)
        {
            vector<int> ret(max(a.size(),b.size()));
            for(register int i = 0;i < ret.size();++i)
                ret[i] = add(a[i],b[i]);
            return poly(ret);
        }
        friend inline poly operator-(const poly &a,const poly &b)
        {
            vector<int> ret(max(a.size(),b.size()));
            for(register int i = 0;i < ret.size();++i)
                ret[i] = dec(a[i],b[i]);
            return poly(ret);
        }
        friend inline poly operator*(poly a,poly b)
        {
            if(a.a.empty() || b.a.empty())
                return poly();
            if(a.size() < 40 || b.size() < 40)
            {
                if(a.size() > b.size())
                    swap(a,b);
                poly ret;
                ret.resize(a.size() + b.size() - 1);
                for(register int i = 0;i < ret.size();++i)
                    for(register int j = 0;j <= i && j < a.size();++j)
                        ret.a[i] = (ret[i] + (long long)a[j] * b[i - j]) % mod;
                ret.shrink();
                return ret;
            }
            int lim = 1,tot = a.size() + b.size() - 1;
            for(;lim < tot;lim <<= 1);
            a.resize(lim),b.resize(lim);
            a.ntt(),b.ntt();
            for(register int i = 0;i < lim;++i)
                a.a[i] = (long long)a[i] * b[i] % mod;
            a.ntt(-1),a.shrink();
            return a;
        }
        poly &operator+=(const poly &o)
        {
            resize(max(size(),o.size()));
            for(register int i = 0;i < o.size();++i)
                a[i] = add(a[i],o[i]);
            return *this;
        }
        poly &operator-=(const poly &o)
        {
            resize(max(size(),o.size()));
            for(register int i = 0;i < o.size();++i)
                a[i] = dec(a[i],o[i]);
            return *this;
        }
        poly &operator*=(poly o)
        {
            return (*this) = (*this) * o;
        }
        poly deriv() const
        {
            if(a.empty())
                return poly();
            vector<int> ret(size() - 1);
            for(register int i = 0;i < size() - 1;++i)
                ret[i] = (long long)(i + 1) * a[i + 1] % mod;
            return poly(ret);
        }
        poly integ() const
        {
            if(a.empty())
                return poly();
            vector<int> ret(size() + 1);
            for(register int i = 0;i < size();++i)
                ret[i + 1] = (long long)a[i] * inv[i + 1] % mod;
            return poly(ret);
        }
        inline poly modxn(int n) const
        {
            if(a.empty())
                return poly();
            n = min(n,size());
            return poly(vector<int>(a.begin(),a.begin() + n));
        }
        inline poly inver(int m) const
        {
            poly ret(fpow(a[0],mod - 2)),f,g;
            for(register int k = 1;k < m;)
            {
                k <<= 1,f.resize(k),g.resize(k);
                for(register int i = 0;i < k;++i)
                    f.a[i] = (*this)[i],g.a[i] = ret[i];
                f.ntt(),g.ntt();
                for(register int i = 0;i < k;++i)
                    f.a[i] = (long long)f[i] * g[i] % mod;
                f.ntt(-1);
                for(register int i = 0;i < (k >> 1);++i)
                    f.a[i] = 0;
                f.ntt();
                for(register int i = 0;i < k;++i)
                    f.a[i] = (long long)f[i] * g[i] % mod;
                f.ntt(-1);
                ret.resize(k);
                for(register int i = (k >> 1);i < k;++i)
                    ret.a[i] = dec(0,f[i]);
            }
            return ret.modxn(m);
        }
        inline pair<poly,poly> div(poly o) const
        {
            if(size() < o.size())
                return make_pair(poly(),*this);
            poly f,g;
            f = (rever().modxn(size() - o.size() + 1) * o.rever().inver(size() - o.size() + 1)).modxn(size() - o.size() + 1).rever();
            g = (modxn(o.size() - 1) - o.modxn(o.size() - 1) * f.modxn(o.size() - 1)).modxn(o.size() - 1);
            return make_pair(f,g);
        }
        inline poly log(int m) const
        {
            return (deriv() * inver(m)).integ().modxn(m);
        }
        inline poly exp(int m) const
        {
            poly ret(1),iv,it,d = deriv(),itd,itd0,t1;
            if(m < 70)
            {
                ret.resize(m);
                for(register int i = 1;i < m;++i)
                {
                    for(register int j = 1;j <= i;++j)
                        ret.a[i] = (ret[i] + (long long)j * operator[](j) % mod * ret[i - j]) % mod;
                    ret.a[i] = (long long)ret[i] * inv[i] % mod;
                }
                return ret;
            }
            for(register int k = 1;k < m;)
            {
                k <<= 1;
                it.resize(k >> 1);
                for(register int i = 0;i < (k >> 1);++i)
                    it.a[i] = ret[i];
                itd = it.deriv(),itd.resize(k >> 1);
                iv = ret.inver(k >> 1),iv.resize(k >> 1);
                it.ntt(),itd.ntt(),iv.ntt();
                for(register int i = 0;i < (k >> 1);++i)
                    it.a[i] = (long long)it[i] * iv[i] % mod,
                    itd.a[i] = (long long)itd[i] * iv[i] % mod;
                it.ntt(-1),itd.ntt(-1),it.a[0] = dec(it[0],1);
                for(register int i = 0;i < k - 1;++i)
                    itd.a[i % (k >> 1)] = dec(itd[i % (k >> 1)],d[i]);
                itd0.resize((k >> 1) - 1);
                for(register int i = 0;i < (k >> 1) - 1;++i)
                    itd0.a[i] = d[i];
                itd0 = (itd0 * it).modxn((k >> 1) - 1);
                t1.resize(k - 1);
                for(register int i = (k >> 1) - 1;i < k - 1;++i)
                    t1.a[i] = itd[(i + (k >> 1)) % (k >> 1)];
                for(register int i = k >> 1;i < k - 1;++i)
                    t1.a[i] = dec(t1[i],itd0[i - (k >> 1)]);
                t1 = t1.integ();
                for(register int i = 0;i < (k >> 1);++i)
                    t1.a[i] = t1[i + (k >> 1)];
                for(register int i = (k >> 1);i < k;++i)
                    t1.a[i] = 0;
                t1.resize(k >> 1),t1 = (t1 * ret).modxn(k >> 1),t1.resize(k);
                for(register int i = (k >> 1);i < k;++i)
                    t1.a[i] = t1[i - (k >> 1)];
                for(register int i = 0;i < (k >> 1);++i)
                    t1.a[i] = 0;
                ret -= t1;
            }
            return ret.modxn(m);
        }
        inline poly sqrt(int m) const
        {
            poly ret(1),f,g;
            for(register int k = 1;k < m;)
            {
                k <<= 1;
                f = ret,f.resize(k >> 1);
                f.ntt();
                for(register int i = 0;i < (k >> 1);++i)
                    f.a[i] = (long long)f[i] * f[i] % mod;
                f.ntt(-1);
                for(register int i = 0;i < k;++i)
                    f.a[i % (k >> 1)] = dec(f[i % (k >> 1)],(*this)[i]);
                g = (2 * ret).inver(k >> 1),f = (f * g).modxn(k >> 1),f.resize(k);
                for(register int i = (k >> 1);i < k;++i)
                    f.a[i] = f[i - (k >> 1)];
                for(register int i = 0;i < (k >> 1);++i)
                    f.a[i] = 0;
                ret -= f;
            }
            return ret.modxn(m);
        }
        inline poly pow(int m,int k1,int k2 = -1) const
        {
            if(a.empty())
                return poly();
            if(k2 == -1)
                k2 = k1;
            int t = 0;
            for(;t < size() && !a[t];++t);
            if((long long)t * k1 >= m)
                return poly();
            poly ret;
            ret.resize(m);
            int u = fpow(a[t],mod - 2),v = fpow(a[t],k2);
            for(register int i = 0;i < m - t * k1;++i)
                ret.a[i] = (long long)operator[](i + t) * u % mod;
            ret = ret.log(m - t * k1);
            for(register int i = 0;i < ret.size();++i)
                ret.a[i] = (long long)ret[i] * k1 % mod;
            ret = ret.exp(m - t * k1),t *= k1,ret.resize(m);
            for(register int i = m - 1;i >= t;--i)
                ret.a[i] = (long long)ret[i - t] * v % mod;
            for(register int i = 0;i < t;++i)
                ret.a[i] = 0;
            return ret;
        }
    };
}
using Poly::init;
using Poly::poly;
inline int C(int n,int m)
{
    return n < m ? 0 : (long long)Poly::fac[n] * Poly::ifac[m] % mod * Poly::ifac[n - m] % mod;
}
poly f,b,g,p,q;
int rest,c;
void solve(int l,int r)
{
    if(l == r)
    {
        g.a[l] = (long long)Poly::inv[n - l + 1] * g[l] % mod;
        return ;
    }
    int mid = l + r >> 1;
    solve(l,mid);
    poly t1,t2,t3,t4,t;
    int lim = 1,tot = mid - l + r - l + 2;
    for(;lim < tot;lim <<= 1);
    t1.resize(lim),t2.resize(lim),t3.resize(lim),t4.resize(lim),t.resize(lim);
    for(register int i = 0;i < lim;++i)
        t1.a[i] = t2.a[i] = t3.a[i] = t4.a[i] = t.a[i] = 0;
    for(register int i = 0;i <= mid - l;++i)
        t1.a[i] = (long long)(n - m) * g[i + l] % mod,
        t3.a[i] = (long long)(i + l) * g[i + l] % mod;
    for(register int i = 0;i <= r - l;++i)
        t2.a[i] = i & 1 ? dec(0,Poly::ifac[i]) : Poly::ifac[i],
        t4.a[i] = i & 1 ? Poly::ifac[i + 1] : dec(0,Poly::ifac[i + 1]);
    if(mid - l <= 60 || r - l <= 100)
    {
        for(register int i = mid - l + 1;i <= r - l;++i)
            for(register int j = 0;j <= i && j <= mid - l;++j)
                t.a[i] = (t[i] + (long long)t1[j] * t2[i - j] + (long long)t3[j] * t4[i - j]) % mod;
    }
    else
    {
        t1.ntt(),t2.ntt(),t3.ntt(),t4.ntt();
        for(register int i = 0;i < lim;++i)
            t.a[i] = ((long long)t1[i] * t2[i] + (long long)t3[i] * t4[i]) % mod;
        t.ntt(-1);
    }
    for(register int i = mid + 1;i <= r;++i)
        g.a[i] = dec(g[i],t[i - l]);
    solve(mid + 1,r);
}
poly t1,t2,a;
int main()
{
    init();
    scanf("%d%d%d",&n,&m,&k),n1 = k - 1,n2 = n - k;
    for(register int i = 0;i < n - m;++i)
        rest = (rest + (long long)(i & 1 ? mod - 1 : 1) * Poly::ifac[i] % mod * Poly::ifac[n - i] % mod * (fpow(n - m - i,n) - fpow(n - m - i - 1,n) + mod)) % mod;
    rest = (long long)rest * Poly::fac[n] % mod;
    g.resize(n),
    g.a[0] = (long long)(m & 1 ? (mod - C(n + 1,m)) % mod : C(n + 1,m)) * (n - m + 1) % mod,
    solve(0,n - 1);
    b.resize(n);
    for(register int i = 0;i < n;++i)
        b.a[i] = Poly::ifac[i + 1];
    b = b.inver(n);
    p.resize(n),q.resize(n);
    for(register int i = n2;i < n;++i)
        p.a[n - 1 - i] = (long long)C(i,n2) * (n2 & 1 ? mod - 1 : 1) % mod * g[n - 1 - i] % mod;
    for(register int i = n1;i < n;++i)
        q.a[n - 1 - i] = (long long)C(i,n1) * (i - n1 & 1 ? mod - 1 : 1) % mod * g[n - 1 - i] % mod;
    p = (p * b).modxn(n),q = (q * b).modxn(n),f.resize(n);
    for(register int i = 0;i < n;++i)
        f.a[i] = (long long)Poly::ifac[i] * (p[n - 1 - i] - q[n - 1 - i] + mod) % mod;
    f = f.integ().modxn(n);
    t1.resize(n),t2.resize(n);
    for(register int i = 0;i < n;++i)
        t1.a[n - 1 - i] = (long long)Poly::fac[i] * f[n - 1 - i] % mod,
        t2.a[i] = Poly::ifac[i];
    t1 *= t2,a.resize(n),c = 0;
    for(register int i = 0;i < n;++i)
        a.a[i] = (long long)Poly::ifac[n - 1 - i] * t1[i] % mod * Poly::fac[i] % mod * Poly::fac[n - 1 - i] % mod,
        rest = dec(rest,a[i]),
        c = (c + (long long)C(n - 1,i) * Poly::fac[i] % mod * Poly::fac[n - 1 - i]) % mod;
    rest = (long long)rest * fpow(c,mod - 2) % mod;
    for(register int i = 0;i < n;++i)
        a.a[i] = (a[i] + (long long)C(n - 1,i) * Poly::fac[i] % mod * Poly::fac[n - 1 - i] % mod * rest) % mod;
    for(register int i = 0;i < n;++i)
        printf("%d%c",a[i]," \n"[i == n - 1]);
}