支持一下 Qiuly /tyt

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

首先作一步相当重要的转化：
对所有满足恰有 $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$ 记录大于号。
$$
F(t) = 1 + u \int_0^t F(\tau) \,\mathrm d \tau + v \int_t^1 F(\tau) \,\mathrm d \tau
$$

因此
$$
F’ = (u – v)F
$$

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

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

接下来用 $x$ 元计量边数，用 $y$ 元计量大于号的数量。
对于 $p_k$ 的左侧，我们令 $u=x,v=xy$；对于右侧我们令 $u=xy,v=x$。那么
$$
\begin{aligned}
L &= \frac{ (1-y)\mathrm e^{x(1-y)t} }{ 1-y\mathrm e^{x(1-y)} } \\
R &= \frac{ (1-y)\mathrm e^{x(1-y)(1-t)} }{ 1-y\mathrm 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$ 之和应当等于 $\left\langle\begin{matrix}n\\m\end{matrix}\right\rangle$。

令 $\mathrm D$ 是求导算子，我们知道 $f(t+c) = \mathrm e^{c \mathrm 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(\mathrm e^{-\mathrm D}) (t+i)^{n_2}-(t+i)^{n_2} F(\mathrm e^{\mathrm D}) (t+i)^{n_1}\right)
$$

我们知道我们作一次位移是 $\mathrm e^{\mathrm D}$，正向差分是 $1-\mathrm e^{\mathrm D}$，它的逆运算就是 $\frac1{ 1-\mathrm e^{\mathrm D} }$，然后为了能够求逆我们恰好就凑出来了一个伯努利数的形式
$$
\int B(\mathrm D) \left((t+i)^{n_1} F(\mathrm e^{-\mathrm D}) (t+i)^{n_2}-(t+i)^{n_2} F(\mathrm e^{\mathrm D}) (t+i)^{n_1}\right), \quad B(t) = \frac{t}{1-\mathrm 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(\mathrm e^x) = \mathrm e^x f(\mathrm e^x) = \mathrm e^x g(x)$，那么
$$
\begin{aligned}
(n-m+1)f(\mathrm e^x) + \mathrm e^x f'(\mathrm e^x) &= -m\mathrm e^xf(\mathrm e^x) + \mathrm e^{2x}f'(\mathrm e^x) + C \\
(n-m+1)g + g’ &= -m \mathrm e^x g + \mathrm e^x g’ + C \\
(n-m+1+m\mathrm e^x)g + (1-\mathrm e^x)g’ &= C \\
((n-m)\mathrm e^{-x}+m+1)G + (\mathrm 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{ \mathrm e^{p(1-y)t + q(1-y)(1-t)} }{(1-y\mathrm e^{p(1-y)})(1-y\mathrm e^{q(1-y)})}
$$

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

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

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

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

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

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

分别考虑 $\mathrm e^{(p-q)t} B(p-q) G(p),\mathrm e^{(p-q)t} B(p-q) G(q)$。
以前者为例，我们考虑换元：令 $pr = q$，则变为
$$
\mathrm 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}] \mathrm e^{p(1-r)t} B(p(1-r)) G(p) \\
& = \sum\limits_{j \ge 0} ([p^j]B(p) \mathrm 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) \mathrm 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)$。

代码：

#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]);
}