実装が地味に面倒だった問題。
https://atcoder.jp/contests/arc113/tasks/arc113_f
問題
数直線上に原点およびN個の点の座標X[i]が与えられる。
N個の区間[X[i],X[i+1]]内それぞれにランダムな実数座標上に1個ずつ点を置く。
点同士の距離の最小値の期待値を求めよ。
解法
f(x) := 距離の最小値がx以上である確率
とすると、f(x)を積分すればよいことになる。
そこでf(x)を頑張ってxの多項式で表現しよう。
各区間を[X[i]-x*i,X[i+1]-x*i]のように区間ごとにxだけ左にずらしていこう。
この各区間に点を置いた時、xが昇順になっていれば良いことになる。
そこで登場する座標を座標圧縮し、手前の点から順に、圧縮後のどの区間に入るかを数え上げて行こう。
その際、いくつかの連続する点が座標圧縮後の同じ区間内に入るケースがあるが、その場合n個の点が同じ区間にある場合にそれらが昇順に並ぶ確率は1/(n!)と考えてよい。
xは実数の値を取りうるが、考えるべきはずらした後の複数の区間[X[i]-x*i,X[i+1]-x*i]の大小関係が入れ替わるケースだけである。
それ以外の区間では、f(x)は同じ式で計算することができる。
そこで、区間の大小関係が入れ替わるxでのみ、上記数え上げを行なおう。
int N; int X[22]; ll A[55],B[55][2]; ll P[201][51]; const ll mo=998244353; int id[2][21]; ll pre[30]; ll modpow(ll a, ll n = mo-2) { ll r=1;a%=mo; while(n) r=r*((n%2)?a:1)%mo,a=a*a%mo,n>>=1; return r; } int cmp(pair<ll,ll> L,pair<ll,ll> R) { return L.first*R.second<L.second*R.first; } ll PP,QQ; int cmp2(pair<ll,ll>& a,pair<ll,ll>& b) { ll ap=a.first*PP+a.second*QQ; ll bp=b.first*PP+b.second*QQ; if(ap!=bp) return ap<bp; return a.first>b.first; } vector<ll> hoge(ll p, ll q) { PP=p,QQ=q; vector<pair<ll,ll>> C; int i,j,x; FOR(i,N) { C.push_back({-i,X[i]}); C.push_back({-i,X[i+1]}); } sort(ALL(C),cmp2); FOR(i,N) { FOR(j,C.size()) { if(C[j].first==-i&&C[j].second==X[i]) id[i][0]=j; if(C[j].first==-i&&C[j].second==X[i+1]) id[i][1]=j; } } ll from[44][22][22]={}; from[0][0][0]=1; FOR(i,N) { ll to[44][22][22]={}; ll len=X[i+1]-X[i]; for(int range=id[i][0];range<id[i][1];range++) { ll X=(C[range+1].first-C[range].first+mo)*modpow(len)%mo; ll Y=(C[range+1].second-C[range].second+mo)*modpow(len)%mo; //別のrangeから来るケース for(int order=0;order<=i+1;order++) { ll t=0; for(int num=0;num<=i+1;num++) { for(int sr=0;sr<range;sr++) { t+=from[sr][order][num]; } } (to[range][order][1]+=t%mo*Y)%=mo; (to[range][order+1][1]+=t%mo*X)%=mo; } for(int order=0;order<=i+1;order++) { for(int num=0;num<=i+1;num++) if(from[range][order][num]) { (to[range][order][num+1]+=from[range][order][num]*Y%mo*pre[num+1])%=mo; (to[range][order+1][num+1]+=from[range][order][num]*X%mo*pre[num+1])%=mo; } } } swap(from,to); } vector<ll> R(22); FOR(i,21) { FOR(j,21) { FOR(x,44) { R[i]+=from[x][i][j]; } } R[i]%=mo; } return R; } void solve() { int i,j,k,l,r,x,y; string s; for(i=1;i<=21;i++) pre[i]=modpow(i); cin>>N; FOR(i,N+1) cin>>X[i]; FOR(i,N) { A[i]=-i; B[i][0]=X[i]; B[i][1]=X[i+1]; } vector<pair<ll,ll>> cand; FOR(y,N) FOR(x,y) { cand.push_back({B[y][0]-B[x][0],A[x]-A[y]}); cand.push_back({B[y][0]-B[x][1],A[x]-A[y]}); cand.push_back({B[y][1]-B[x][0],A[x]-A[y]}); cand.push_back({B[y][1]-B[x][1],A[x]-A[y]}); } FORR2(a,b,cand) { ll g=__gcd(a,b); a/=g; b/=g; } sort(ALL(cand),cmp); cand.erase(unique(ALL(cand)),cand.end()); ll ret=0; vector<ll> A; FOR(i,cand.size()) if(i) { ll a=cand[i-1].first*modpow(cand[i-1].second)%mo; ll b=cand[i].first*modpow(cand[i].second)%mo; ll x1=1,x2=1; int order=0; A=hoge(cand[i].first,cand[i].second); FORR(p,A) { x1=x1*a%mo; x2=x2*b%mo; order++; (ret+=p*(x2+mo-x1)%mo*modpow(order))%=mo; } } cout<<ret<<endl; }
まとめ
有理数で多項式の積分をやらせるとかいう非常に面倒な問題。
本番、数値積分に持ち込むことは思いついたものの実装が全く間に合わず。