「2025CCPC哈尔滨」H. 匹配题解

[题目链接](https://qoj.ac/contest/2575/problem/14821)

面对巨大的状态空间，直接进行期望 DP 显然是不现实的。我们引入鞅与停时定理，为状态 $S$ 构造一个“势能函数” $\Phi(S)$。

在非终止状态下，每一次操作后整个系统的期望势能恰好减少 $1$：

$$
E[\Phi(S')] = \Phi(S)-1
$$

根据停时定理，只要满足上述条件，从初始状态走到终止状态的期望步数 $E[T]$，就等于两者的势能差。若我们令终止状态（所有元素都在一个集合）的势能为 $0$，那么答案就直接等于初始状态的势能：

$$
E[T] = \Phi(\{a_1,a_2,\cdots a_n\})
$$

因为每个集合经过一次操作之后的变化是独立的（只和它自己的元素数量有关）

我们可以试着把 $\Phi(S)$ 分成多个的和，但还需要保证转移式子成立

令 $\Phi(S) = \sum_{i=1}^n f(a_i)$

$E[\Phi(S')] = E\left[ \sum_{i=1}^n f(a_i') \right]= \sum_{i=1}^n E[f(a_i')]$

所以 $\sum_{i=1}^n \left( f(a_i) - E[f(a_i')] \right) = 1$

为了求 $f$，我们需要把右边的 $1$ 拆成和左边可对应的多项，并使得上面的式子永远成立，而一个不变的量是 $\sum_{i=1}^n a_i=2m$

$\sum_{i=1}^n \left( f(a_i) - E[f(a_i')] \right) = \sum_{i=1}^n \frac{a_i}{2m}$

然后就可以得到转移了：

$f(a) - E[f(a')] = \frac{a}{2m}$，对于 $a\in[1,2m-1]$ 成立，边界情况是 $f(0)=f(2m)=0$

要求 $E[f(a')]$，我们可以枚举 $i$ 个配对中元素数量 $+1$，$j$ 个配对中 $-1$，剩下的内部配对不影响，但内部配对元素数量需要是偶数

可以得到

$f(a) - \sum_{i, j} P(a \to a+i-j) \cdot f(a+i-j) = \frac{a}{2m}$

其中

$P(a \to a+i-j) = [\text{令 k=i+j, a-k 为偶数}] \frac{1}{(2m-1)!!} \times \binom{a}{k} \binom{2m-a}{k} k!\binom{k}{i} \left(\frac{1}{2}\right)^k (a-k-1)!! (2m-a-k-1)!!$

就可以得到 $2m+1$ 个方程和未知数

使用高斯消元即可得到每一个 $f$， $\sum_{i=1}^n f(a_i)$ 即为我们要求的答案！

---

CODE

```cpp
#include <bits/stdc++.h>
#define int long long

using namespace std;

const int N=810, P=998244353;
int n, m, a[N][N], f[N];

int fac[N], inv[N], facc[N], invv[N], inv2[N];
int power(int x, int y) {
    int r=1;
    while (y) {
        if (y&1) r=r*x%P;
        x=x*x%P; y>>=1;
    }
    return r;
}

void init(int n) {
    fac[0]=1;
    for (int i=1; i<=n; i++) fac[i]=fac[i-1]*i%P;
    inv[n]=power(fac[n], P-2);
    for (int i=n-1; ~i; i--) inv[i]=inv[i+1]*(i+1)%P;
    facc[1]=1;
    for (int i=3; i<=n; i++) facc[i]=facc[i-2]*i%P;
    for (int i=1; i<=n; i++) invv[i]=power(facc[i], P-2);
    inv2[0]=1;
    int Inv2=(P+1)/2;
    for (int i=1; i<=n; i++) inv2[i]=inv2[i-1]*Inv2%P;
}

int C(int n, int m) {
    return fac[n]*inv[m]%P*inv[n-m]%P;
}

int getP(int a, int i, int j) {
    int k=i+j;
    int res=invv[2*m-1]*C(a, k)%P*C(2*m-a, k)%P*fac[k]%P*C(k, i)%P*inv2[k]%P;
    if (a-k>0) res=res*facc[a-k-1]%P;
    if (2*m-a-k>0) res=res*facc[2*m-a-k-1]%P;
    return res;
}

void gauss() {
    for (int A=0; A<=m*2; A++) {
        fill(a[A], a[A]+m*2+1+1, 0);
        if (A==0 || A==m*2) {
            a[A][A]=1;
            continue;
        }
        a[A][A]=1;
        for (int i=0; i<=A; i++) {
            for (int j=0; j<=A; j++) {
                if (i+j>A || i+j>2*m-A) continue;
                if ((A-i-j)%2==1) continue;
                (a[A][A+i-j]+=P-getP(A, i, j))%=P;
            }
        }
        a[A][m*2+1]=A*power(m*2, P-2)%P;
    }
    for (int i=0; i<=m*2; i++) {
        int pos=-1;
        for (int j=i; j<=m*2; j++)
            if (a[j][i]) {pos=j; break;}
        if (pos==-1) continue;
        if (pos!=i) swap(a[pos], a[i]);
        for (int j=0; j<=m*2; j++) {
            if (i==j) continue;
            int t=a[j][i]*power(a[i][i], P-2)%P;
            for (int k=0; k<=m*2+1; k++)
                (a[j][k]+=P-a[i][k]*t%P)%=P;
        }
    }
    for (int i=0; i<=m*2; i++)
        f[i]=a[i][m*2+1]*power(a[i][i], P-2)%P;
}

void solve() {
    cin>>n>>m;
    gauss();
    int ans=0;
    for (int i=1; i<=n; i++) {
        int x;
        cin>>x;
        ans=(ans+f[x])%P;
    }
    cout<<ans<<'\n';
}

signed main() {
    ios::sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    int T=1;
    cin>>T;
    init(801);
    while (T--) solve();
    return 0;
}
```