Atcoder Regular Contest 113F. Social Distance

给定一个长度为 \(n+1\) 的非负整数序列 \(X_0,X_1,\ldots,X_{n}\)，满足 \(0=X_0<X_1<\ldots<X_{n}\)。
在每个区间 \([X_0,X_1],[X_1,X_2],\ldots,[X_{n-1},X_{n}]\) 中均匀随机一个点，求相邻两个点之间的最短距离的期望。
答案对 \(998244353\) 取模。
\(2\le n\le 20,0\le X_i\le 10^6\)。

---

首先，此题中的连续随机变量的期望可以用下式计算

\[\text{E}(X)=\int_{0}^{+\infty}\text{P}(t<X)\text{d}t \]

考虑一般的离散型连续变量 \(X\)，它有 \(n+1\) 个取值 \(0,\Delta x,2\Delta x,\ldots,n\Delta x\)，\(\Delta x\) 满足 \(\Delta x=\frac{V}{n}\)，期望为

\[\begin{aligned}\text{E}(X)&=\sum_{i=0}^{n}i\Delta x\cdot \text{P}(X=i\Delta x)\\ &=\sum_{i=0}^{n}\text{P}(X=i\Delta x)\sum_{j=0}^{i-1}\Delta x\\&=\sum_{j=0}^{n}\Delta x\sum_{i=j+1}^{n}\text{P}(X=j\Delta x)\\ &=\sum_{j=0}^{n}\Delta x\cdot\text{P}(j\Delta x<X)\end{aligned} \]

令 \(n\rightarrow \infty\)，可得

\[\text{E}(X)=\int_{0}^{V}\text{P}(t<X)\text{d}t \]

考虑固定 \(t\) 后计算 \(\text{P}(t<X)\)。设第 \(i\) 个人点的位置为 \(p_i\)，那么有

\[\forall i\in[0,n-2],p_{i+1}-p_{i}\ge t \]

令 \(p'_{i}=p_{i}-it\)，则有

\[\forall i\in[0,n-2],p'_{i+1}\ge p'_{i} \]

根据 \(X\) 序列可以求出 \(p'\) 的取值范围，进而问题转化成

有 \(n\) 个连续随机变量，第 \(i\) 个连续随机变量 \(x_i\) 的取值范围为 \([l_i,r_i]\)，求 \(x_1\le x_2\le \ldots\le x_{n}\) 的概率。

先根据所有取值范围的左右边界放到坐标轴上，此时形成 \(2n-1\) 个区间，设第 \(i\) 个区间为 \([L_i,R_i]\)。

设 \(dp_{i,j}\) 为第 \(i\) 个连续随机变量在第 \(j\) 个区间且 \(x_1\le x_2\le\ldots\le x_{i}\) 的概率，则有转移

\[dp_{i,j}=\sum_{k=1}^{i}\left(\left(\sum_{l=1}^{j-1}dp_{k,l}\right)\times \frac{1}{(i-k)!}\left(\prod_{p=k+1}^i[[L_j,R_j)\sube [l_p,r_p]]\frac{R_j-L_j}{r_p-l_p}\right)\right) \]

即枚举最后一个在第 \(j-1\) 个区间的随机变量 \(x_k\)，于是 \(x_{k+1},x_{k+2},\ldots,x_i\) 都在第 \(j\) 个区间，满足 \(x_{k+1},x_{k+2},\ldots,x_i\) 单调不减的概率显然是 \(\frac{1}{(i-k)!}\)。

维护 \(dp_{i,j}\) 关于 \(j\) 的前缀和以及上式的乘积式就可以做到 \(O(n^3)\) 求出答案。

但是 \(t\) 是变量，而不是固定的。考虑到 \(p'_{i}=-i\cdot t+p_i\) 是关于 \(t\) 的一次函数，\(\text{dp}\) 的转移方式只与 \(L_i,R_i\) 的相对顺序有关，而在若干个段之间的 \(t\) 形成的 \(L_i,R_i\) 的相对顺序是相同的，可以发现不同的相对顺序是 \(O(n^2)\) 级别的。维护关于 \(t\) 的多项式进行转移，最后再对结果多项式进行积分即可。

总时间复杂度 \(O(n^6)\)。

#include<bits/stdc++.h>
#define fi first
#define se second
using namespace std;
typedef long long ll;
const double eps=1e-6;
const int N=45,M=5005,mod=998244353;
struct Poly{
	int n;ll coe[N];
	Poly(){n=0;memset(coe,0,sizeof coe);}
	Poly operator+(const Poly&p)const{
		Poly ret;ret.n=max(n,p.n);
		for(int k=0;k<=ret.n;++k)ret.coe[k]=(coe[k]+p.coe[k])%mod;
		while(ret.n&&!ret.coe[ret.n])--ret.n;
		return ret;
	}
	Poly operator*(const Poly&p)const{
		Poly ret;ret.n=n+p.n;
		for(int j=0;j<=n;++j)for(int k=0;k<=p.n;++k)
			(ret.coe[j+k]+=coe[j]*p.coe[k]%mod)%=mod;
		while(ret.n&&!ret.coe[ret.n])--ret.n;
		return ret;
	}
	Poly operator*(const int c)const{
		Poly ret;ret.n=n;
		for(int k=0;k<=n;++k)ret.coe[k]=coe[k]*c%mod;
		return ret;
	}
}Len[M],dp[N][N];
struct Frac{
	ll x,y;
	bool operator<(const Frac&a){return x*a.y<a.x*y;}
}seq[M];
int n,X[N],cnt,idL[N],idR[N];ll ans,fac[M],inv[M],Inv[M];
pair<double,int>p[N];
inline ll qpow(ll a,ll b){ll ans=1;for(;b;b>>=1,a=a*a%mod)if(b&1)ans=ans*a%mod;return ans;}
inline void precalc(int n){
	fac[0]=inv[0]=Inv[0]=fac[1]=inv[1]=Inv[1]=1;
	for(int i=2;i<=n;++i)
		fac[i]=fac[i-1]*i%mod,
		Inv[i]=(mod-mod/i)*Inv[mod%i]%mod,
		inv[i]=inv[i-1]*Inv[i]%mod;
}
inline ll Integral(Poly f,ll x){
	ll ans=0,pw=x;
	for(int i=0;i<=f.n;++i,pw=pw*x%mod)(ans+=pw*Inv[i+1]%mod*f.coe[i]%mod)%=mod;
	return ans;
}
inline ll calc(Frac sL,Frac sR){
	int tot=0;ll FR=sL.x*qpow(sL.y,mod-2)%mod,RF=sR.x*qpow(sR.y,mod-2)%mod;
	double tL=1.0*sL.x/sL.y+eps;
	for(int i=0;i<n;++i)
		p[++tot]=make_pair(X[i]-i*tL,i+1),p[++tot]=make_pair(X[i+1]-i*tL,-i-1);
	sort(p+1,p+1+tot);
	for(int i=1;i<tot;++i)
		Len[i].coe[0]=(X[p[i+1].se>0?p[i+1].se-1:-p[i+1].se]-X[p[i].se>0?p[i].se-1:-p[i].se]+mod)%mod,
		Len[i].coe[1]=(abs(p[i].se)-abs(p[i+1].se)+mod)%mod,Len[i].n=1;
	for(int i=1;i<=tot;++i){
		if(p[i].se>0)idL[p[i].se-1]=i;
		else idR[-p[i].se-1]=i-1;
	}
	for(int i=0;i<n;++i)for(int j=0;j<tot;++j)dp[i][j]=Poly();
	for(int i=idL[0];i<=idR[0];++i)dp[0][i]=Len[i]*qpow(X[1],mod-2);
	for(int i=1;i<tot;++i)dp[0][i]=dp[0][i]+dp[0][i-1];
	for(int i=1;i<n;++i){
		for(int j=idL[i];j<=idR[i];++j){
			Poly sufp;sufp.coe[0]=1;
			for(int k=i;~k&&idL[k]<=j&&j<=idR[k];--k){
				Poly tmp;
				if(!k)tmp.coe[0]=1;else tmp=dp[k-1][j-1];
				sufp=sufp*Len[j]*qpow(X[k+1]-X[k],mod-2);
				dp[i][j]=dp[i][j]+sufp*tmp*inv[i-k+1];
			}
		}
		for(int j=1;j<tot;++j)dp[i][j]=dp[i][j]+dp[i][j-1];
	}
	ll Ans=(Integral(dp[n-1][tot-1],RF)-Integral(dp[n-1][tot-1],FR)+mod)%mod;
	return Ans;
}
int main(){
	precalc(5000);
	scanf("%d",&n);
	for(int i=0;i<=n;++i)scanf("%d",X+i);
	for(int i=0;i<n;++i)
		for(int j=i+1;j<n;++j)
			seq[++cnt]=(Frac){X[j]-X[i],j-i},seq[++cnt]=(Frac){X[j+1]-X[i+1],j-i},
			seq[++cnt]=(Frac){X[j]-X[i+1],j-i},seq[++cnt]=(Frac){X[j+1]-X[i],j-i};
	sort(seq+1,seq+1+cnt);
	for(int i=1;i<=cnt;++i)(ans+=calc(seq[i],seq[i+1]))%=mod;
	return printf("%lld\n",ans),0;
}