Solution 「CF 923E」Perpetual Subtraction

\(\mathcal{Description}\)

  Link.

  有一个整数 \(x\in[0,n]\)，初始时以 \(p_i\) 的概率取值 \(i\)。进行 \(m\) 轮变换，每次均匀随机取整数 \(r\in[0,x]\)，令 \(x\leftarrow r\)。求变换完成后 \(x=i~(i=0..n)\) 的概率。答案模 \(998244353\)。

\(\mathcal{Solution}\)

  令向量 \(\boldsymbol p\) 为此时 \(x\) 的取值概率，显然一次变换是对 \(\boldsymbol p\) 的线性变换 \(A\)，有

\[A=\begin{bmatrix} 1&\frac{1}{2}&\cdots&\frac{1}{n}&\frac{1}{n+1}\\ &\frac{1}{2}&\cdots&\frac{1}{n}&\frac{1}{n+1}\\ &&\ddots&\vdots&\vdots\\ &&&\frac{1}{n}&\frac{1}{n+1}\\ &&&&\frac{1}{n+1} \end{bmatrix}. \]

  我们的目标即为求出 \(A^m\boldsymbol p\)。注意到 \(A\) 的特征值很明显——\(\lambda_i=\frac{1}{i+1},i\in[0,n]\)，所以可以考虑将其对角化来加速矩阵幂的计算。手算一下 \(\lambda_i\) 所对应的特征向量 \(\boldsymbol v_i\)，发现一组特解

\[\boldsymbol v_i=\begin{bmatrix} \binom{i}{0}\\ -\binom{i}{1}\\ \binom{i}{2}\\ \vdots\\ (-1)^n\binom{i}{n} \end{bmatrix}, \]

那么 \(V=\begin{bmatrix}\boldsymbol v_0&\boldsymbol v_1&\cdots&\boldsymbol v_n\end{bmatrix}\) 有

\[V_{ij}=(-1)^i\binom{j}{i}. \]

  尝试对 \(V\) 求逆，由于

\[\begin{aligned} V_{ij}^2 &= \sum_{k=0}^{n-1}(-1)^i\binom{k}{i}\cdot(-1)^k\binom{j}{k}\\ &= \sum_{k=0}^{n-1}(-1)^{i+k}\binom{j}{i}\binom{j-i}{k-i}\\ &= (-1)^i\binom{j}{i}\sum_{k=0}^{n=1}(-1)^k\binom{j-i}{k-i}\\ &= (-1)^i\binom{j}{i}(1-1)^{j-i}\\ &= [i=j], \end{aligned} \]

所以 \(V=V^{-1}\)。因此答案为

\[V\begin{bmatrix} \lambda_0^m\\ &\lambda_1^m\\ &&\ddots\\ &&&\lambda_n^m \end{bmatrix}V\boldsymbol p. \]

其中 \(V\) 的变换效果是一个差卷积，NTT 实现即可。复杂度 \(\mathcal O(n\log n)\)。

\(\mathcal{Code}\)

/*+Rainybunny+*/

#include <bits/stdc++.h>

#define rep(i, l, r) for (int i = l, rep##i = r; i <= rep##i; ++i)
#define per(i, r, l) for (int i = r, per##i = l; i >= per##i; --i)

typedef long long LL;

const int MAXL = 1 << 18, MOD = 998244353, G = 3;
int n, p[MAXL + 5], fac[MAXL + 5], ifac[MAXL + 5];
LL m;

inline int sgn(const int u) { return u & 1 ? MOD - 1 : 1; }
inline int mul(const int u, const int v) { return 1ll * u * v % MOD; }
inline int sub(int u, const int v) { return (u -= v) < 0 ? u + MOD : u; }
inline int add(int u, const int v) { return (u += v) < MOD ? u : u - MOD; }
inline int mpow(int u, int v) {
    int ret = 1;
    for (; v; u = mul(u, u), v >>= 1) ret = mul(ret, v & 1 ? u : 1);
    return ret;
}

inline void init() {
    fac[0] = 1;
    rep (i, 1, n) fac[i] = mul(i, fac[i - 1]);
    ifac[n] = mpow(fac[n], MOD - 2);
    per (i, n - 1, 0) ifac[i] = mul(i + 1, ifac[i + 1]);
}

inline void ntt(const int n, int* u, const int tp) {
    static int rev[MAXL + 5]; int lgn = 31 - __builtin_clz(n);
    rep (i, 1, n - 1) rev[i] = rev[i >> 1] >> 1 | (i & 1) << lgn >> 1;
    rep (i, 0, n - 1) if (i < rev[i]) std::swap(u[i], u[rev[i]]);
    for (int stp = 1; stp < n; stp <<= 1) {
        int wi = mpow(G, (MOD - 1) / (stp << 1));
        for (int j = 0; j < n; j += stp <<1 ) {
            for (int wk = 1, k = j; k < j + stp; ++k, wk = mul(wk, wi)) {
                int ev = u[k], ov = mul(wk, u[k + stp]);
                u[k] = add(ev, ov), u[k + stp] = sub(ev, ov);
            }
        }
    }
    if (!~tp) {
        std::reverse(u + 1, u + n);
        int inv = mpow(n, MOD - 2);
        rep (i, 0, n - 1) u[i] = mul(u[i], inv);
    }
}

inline void transV() {
    static int f[MAXL + 5], g[MAXL + 5];
    int len = 1 << 32 - __builtin_clz((n << 1) - 1);
    rep (i, 0, len - 1) f[i] = g[i] = 0;
    rep (i, 0, n - 1) f[i] = mul(fac[i], p[i]), g[i] = ifac[n - 1 - i];
    ntt(len, f, 1), ntt(len, g, 1);
    rep (i, 0, len - 1) f[i] = mul(f[i], g[i]);
    ntt(len, f, -1);
    rep (i, 0, n - 1) p[i] = mul(mul(sgn(i), ifac[i]), f[n - 1 + i]);
}

int main() {
    scanf("%d %lld", &n, &m);
    ++n, m %= MOD - 1, init();
    rep (i, 0, n - 1) scanf("%d", &p[i]);

    transV();
    rep (i, 0, n - 1) p[i] = mul(p[i], mpow(i + 1, MOD - 1 - m));
    transV();
    rep (i, 0, n - 1) printf("%d%c", p[i], i < repi ? ' ' : '\n');
    return 0;
}