[IAMOI R2] 污水博弈 题解

[IAMOI R2] 污水博弈

有 \(n\) 个水槽，初始互不连通，内部水位高度为 \(a_i\)，有 \(n-1\) 个闸门连接相邻的水槽，现在随机开启一些闸门，会形成一些连通块，一个连通块内的水位高度等于每个位置水位的平均值，然后又会在连通块中随机抽取一个，询问抽取到的连通块水位的期望高度，\(n \le 10^6\)。

思路：

求随机抽取到的连通块的水位高度的期望，我们不能一一求出每一种方案的贡献，那就考虑每一个连通块对答案的贡献，最后除以总方案数即可。

钦定我们会选中 \([l,r]\) 的位置，且这些位置构成一整个连通块，那这个连通块内部的阀门的开关情况是固定了的。考虑一下此时贡献长什么样子，大概就是形如 \(\frac{ht(l,r)}{x} \times f(x)\)，其中 \(ht(l,r)\) 表示这一块高度的平均数，\(x\) 表示一共有几个连通块，\(f(x)\) 则表示此时形成 \(x\) 个连通块的方案数。

由于高度的平均值比较好算，我们考虑一下 \(\frac{f(x)}{x}\) 怎么表示，令 \(m\) 表示当前还有几个随意开关的阀门，先假设他们是全部打开的，如果我们关闭了其中 \(j\) 个，就会在原本的基础上新增 \(j\) 个连通块，这里原本的基础就是 \([l,r]\) 把整个序列分成了几块，也可以说是当前最少有几个连通块，假设它为 \(k\) 吧，待会我们来分讨一下这个 \(k\) 的取值情况。

那我们要求的式子除了高度外的部分就变成了 \(\sum_{j=0}^m\binom{m}{j}\frac{1}{j+k}\)。这个 \(\frac{1}{j+k}\) 是比较讨厌的，考虑把它转化掉。由积分的知识可知，函数 \(x^{n-1}\) 在 \([0,1]\) 内与 \(x\) 轴围成的面积大小为 \(\frac{1}{n}\)，所以我们可以把 \(\frac{1}{j+k}\) 变成 \(\int_0^1 x^{j+k-1} \mathrm{d}x\)。则设 \(g(k,m)\) 等于：

\[\int_0^1 \sum_{j=0}^m\binom{m}{j}x^{j+k-1} \mathrm{d}x \]

\[=\int_0^1x^{k-1}\sum_{j=0}^m\binom{m}{j}x^j \mathrm{d}x \]

\[=\int_0^1x^{k-1}(x+1)^m \mathrm{d}x \]

又因为牛莱公式，我们知道 \(\int_a^b f'(x) \mathrm{d}x = f(b)-f(a)\)，所以我们可以 \(O(1)\) 算出这个式子，前提是要对 \(k\) 进行一些分讨。

分类讨论如下：

1. 当 \(l=1,r=n\) 时，此时 \(m=0\)，\(k=1\)，权重为 \(g(1,m)=\int_0^1(1+x)^m \mathrm{d}x\)。
2. 当 \(l=1,r<n\) 或者 \(l>1,r=n\) 时，此时 \(m=n-r-1\) 或 \(l-2\)，而 \(k=2\)，因为此时最少被分成两个固定的连通块。权重为 \(g(2,m) = \int_0^1 x(1+x)^m \mathrm{d}x\)。
3. 当 \(l>1,r<n\) 时，此时令 \(d = r-l+1\)，则 \(m=n-d-2\)，\(k=3\)，相当于把序列从中取了一段，那就是至少三个连通块。权重为 \(g(3,m) = \int_0^1 x^2(1+x)^m \mathrm{d}x\)。

对每一种情况我们考虑怎么计算，令 \(s_i\) 为前缀和序列。第一种很好算，\(g(1,m)\) 为其权重，用牛莱公式转化一下得 \(\int_0^1(x+1)^m \mathrm{d}x=\frac{2^{m+1}-1}{m+1}\)，记为 \(h(m)\)，因为 \(m=0\)，贡献就是 \(1\times \frac{s_n}{n}\)。

第二种情况，固定一段，枚举另一端，此时可以得到平均高度，权重为 \(g(2,m)=\int_0^1x(1+x)^m \mathrm{d}x\)，乍一看好像看不出那个函数的导数是 \(x(x+1)^m\)，但我们把 \(x\) 表示为 \((1+x)-1\)，则原式等于 \(\int_0^1 (1+x)^{m+1}-(1+x)^m \mathrm{d}x\)，就是 \(h(m+1)-h(m)\)。

第三种情况比较难搞，我们毕竟不可能枚举每一个区间。但是我们发现对于同一个长度 \(d\) 的区间其权值 \(g(3,n-d-2)\) 都是一样的，那我们就考虑对于一个长度算出所有该长度的区间的高度的平均数之和，即为 \(\sum_{l=2}^{n-d} \dfrac{s_{l+d-1}-s_{l-1}}{d}\)，展开之后发现分子等于 \(\sum_{i=d+1}^{n-1}s_i-\sum_{i=0}^{n-d-1}s_i\)，这可以求前缀和的前缀和快速求出。考虑一下 \(g(3,m)\) 怎么算，其实和上面差不多，\(\int_0^1x^2(1+x)^m \mathrm{d}x\)，把 \(x^2\) 化成 \((1+x-1)^2\)，则原式等于 \(h(m+2)+h(m)-2\times h(m+1)\)。

至此就做完了，语言描述可能难懂，可以对照代码理解。时间复杂度 \(O(n)\)。

代码：

#include<bits/stdc++.h>
using namespace std;
#define int long long
#define _int __int128
#define ull unsigned long long
#define pii pair<int,int>
#define fst first
#define scd second
#define pq priority_queue
#define mkp make_pair
#define popcount(x) __builtin_popcount(x)
#define endl '\n'
int n,ans;
const int N = 1e6+10,mod=998244353;
int a[N],s[N],ss[N],pw[N],inv[N];
int qpow(int a,int b){
    int res=1;
    while(b){
        if(b&1)res=(res*a)%mod;
        a=(a*a)%mod;
        b>>=1;
    }
    return res;
}
void add(int &a,int b){
    a%=mod,b%=mod;
    a+=b;
    if(a>=mod)a-=mod;
}
int h(int m){
    int res=pw[m+1]-1;
    res=(res%mod*inv[m+1]%mod)%mod;
    return res;
}
int g1(int m){
    return h(m);
}
int g2(int m){
    return (h(m+1)-h(m)+mod)%mod;
}
int g3(int m){
    return ((h(m+2)+h(m))%mod-2*h(m+1)%mod+mod)%mod;
}
signed main(){
    ios::sync_with_stdio(0);
    cin.tie(0),cout.tie(0);
    cin>>n;
    pw[0]=1;
    for(int i=1;i<=n;i++){
        pw[i]=(pw[i-1]*2)%mod;
    }
    inv[1]=1;
    for(int i=2;i<=n;i++){
        inv[i]=(mod-(mod/i)*inv[mod%i]%mod)%mod;
    }
    for(int i=1;i<=n;i++){
        cin>>a[i];
        s[i]=(s[i-1]+a[i])%mod;
    }
    for(int i=1;i<=n;i++){
        ss[i]=(ss[i-1]+s[i])%mod;
    }
    add(ans,s[n]%mod*inv[n]%mod);
    for(int r=1;r<n;r++){
        int ht=(s[r]%mod*inv[r]%mod)%mod;
        int res=(ht%mod*g2(n-r-1)%mod)%mod;
        add(ans,res);
    }
    for(int l=2;l<=n;l++){
        int ht=(s[n]-s[l-1]);
        ht=(ht%mod*inv[n-l+1]%mod)%mod;
        int res=(ht%mod*g2(l-2)%mod)%mod;
        add(ans,res);
    }
    for(int d=1;d<=n-2;d++){
        int ht=(ss[n-1]-ss[d]+mod)%mod;
        ht=(ht-ss[n-d-1]+mod)%mod;
        ht=(ht%mod*inv[d]%mod)%mod;
        int res=(ht%mod*g3(n-d-2)%mod)%mod;
        add(ans,res);
    }
    ans=(ans%mod*qpow(pw[n-1],mod-2)%mod)%mod;
    cout<<ans;
    return 0;
}