AGC007 C-Pushing Balls

ボールは1-indexed、穴は0-indexedで左側から番号を付ける
つまり、 $i$ 番目のボールの左には $i-1$ 番目の穴、右には $i$ 番目の穴があるとする( $1\leq i\leq n$ )

立式

同様に
とおくと明らかに $Pr(i,S)=Pl(i,S)$ でありそれらを $P(i,S)$ とおき直すと答えは

[tex: \sum{i=1}^n \left( \sum{S=1}^{i+S=n+1} P(i,S)(R(i,S)+L(i,S) \right)]

具体的な値を求める

P(i,S)を求める

とすると
1.Xが右に転がる
2.XがXの右S-1個よりも後に転がる
3.Xの右S-1個がYにたどり着かない
の3つを満たす必要があり、それらの確率を考えると
[tex:a_S=\dfrac{1}{2S}\left(1-\sum{i=1}^{S-1}a_i\right)]
が導かれる
$P(i,S)$ はさらにYよりも右にある $n-i-S+1$ 個のボールがYに落ちない必要があり、Yより $k$ 個右のボールがXより先にYにたどり着く確率は適当に考えると
$a_k\dfrac{S}{S+k}$
であるから
[tex:P(i,S)=a_S\left(1-\sum{k=1}^{n-i-S+1}a_k\dfrac{S}{S+k}\right)]

R(i,S)とL(i,S)を求める

tex:R(i,S)は
初項: $d+(2i-1)x$
公差: $x$
項数: $2S-1$
の等差数列の和であるから
$R(i,S)=(d+x(2i+S-2))(2S-1)$
tex:L(i,S)は
初項: $d+(2n-1)x+(2i-1)(-x)$
公差: $-x$
項数: $2S-1$
の等差数列の和であるから
$L(i,S)=(d+x(2n-2i-S+1))(2S-1)$
$∴R(i,S)+L(i,S)=(2d+x(2n-1))(2S-1)$
$i$ に寄らなくなったのでこれを $L(S)$ と置き直して和の順序を入れ替えれば求める答えは
[tex:\sum{S=1}^{n}\left(L(S)\sum{i=1}^{i+S=n+1}P(i,S)\right)]
$L(S)$ と $P(i,S)$ に具体的な値を代入すると
[tex:(2d+x(2n-1))\sum{S=1}^n(2S-1)\left(\sum{i=1}^{i+S=n+1}a_S\left(1-\sum{k=1}^{n-i-S+1}\dfrac{a_kS}{S+k}\right)\right)]
$C=(2x+x(2n-1))$ と置き、総和部分を適切に見ると
[tex:C\sum{S=1}^n(2S-1)(n-S+1)a_S+C\sum{i,S,k=1}^{i+S+k=n+1}\dfrac{S(2S-1)a_Sa_k}{S+k}]
よく考えると二項目の $i$ は消せて
[tex:answer=C\sum{S=1}^n(2S-1)(n-S+1)a_S+C\sum_{S,k=1}^{S+k=n}(n+1-S-k)\dfrac{S(2S-1)a_Sa_k}{S+k}]

tex:a_Sの具体的な値を求める

[tex:a_S=\dfrac{1}{2S}\left(1-\sum{i=1}^{S-1}a_i\right)]
を見ると総和の形をした漸化式なので、
[tex:A_S=\sum{i=1}^Sa_i]
とおいてみると
[tex:A_S-A{S-1}=\dfrac{1}{2S}\left(1-A{S-1}\right)]
特性方程式等を使って変形すると
[tex:A_S-1=\left(1-\dfrac{1}{2S}\right)(A{S-1}-1)]
$∴A_S-1=-\dfrac{(2S-1)!!}{(2S)!!}$
[tex:a_S=A_S-A{S-1}=(A_S-1)-(A{S-1}-1)]に代入して整理すると
$a_S=\dfrac{(2S-3)!!}{(2S)!!}$
これは順番に求めれば $a_1$ から $a_n$ まで $O(n)$ で求めることが出来るのは明らか
よってanswerの一項目はO(n)で求められる
また二項目に関しては
$F_S=S(2S-1)a_S=\dfrac{S(2S-1)!!}{(2S)!!}$
$G_k=a_k=\dfrac{(2k-3)!!}{(2k)!!}$
に対して
[tex:H_T=\sum{S+k=T}F_SG_k]
となるHが求められれば
$C\sum_{T=2}^n\dfrac{n+1-T}{T}H_T$
このようなHはFFTで $O(n\log n)$ で求められる

実装上の注意

・FFTの定数倍が重くlong doubleでやるとTLEが（俺の実力では）取れない
・二重階乗を求めてから分数に直そうとすると値が大きくなりすぎてオーバーフローするので分数のまま扱う

#include <bits/stdc++.h>
using namespace std;

using ld=double;
#define REP(i,s,n) for(long long i=s;i<n;i++)
ld ans;
//[競プロのための高速フーリエ変換](https://www.creativ.xyz/fast-fourier-transform/)
// Cooley–Tukey FFT algorithm
const ld PI=acos(-1);
complex<ld> W,S,T;
inline vector<complex<ld>> fft(vector<complex<ld>> a,bool inverse=false) {
  int n=a.size(),h=0;
  for(int i=0;(1<<i)<n;i++)h++;
  REP(i,0,n){
    int j=0;
    REP(k,0,h)j|=(i>>k&1)<<(h-1-k);
    if(i<j)swap(a[i],a[j]);
  }
  for(int b=1;b<n;b*=2){
    REP(j,0,b){
      W=polar((ld)1.0,(2*PI)/(2*b)*j*(inverse?1:-1));
      for(int k=0;k<n;k+=b*2){
        S=a[j+k];
        T=a[j+k+b]*W;
        a[j+k]=S+T;
        a[j+k+b]=S-T;
      }
    }
  }
  if(inverse)REP(i,0,n)a[i]/=n;
  return a;
}
inline vector<complex<ld>> fft(vector<ld> a,bool inverse=false){
  vector<complex<ld>> a_complex(a.size());
  REP(i,0,a.size())a_complex[i]=complex<ld>(a[i],0);
  return fft(a_complex,inverse);
}
vector<ld> convolve(vector<ld> a,vector<ld> b){
  int s=a.size()+b.size()-1,t=1;
  while(t<s)t*=2;
  a.resize(t);
  b.resize(t);
  vector<complex<ld>> A=fft(a);
  vector<complex<ld>> B=fft(b);
  REP(i,0,t)A[i]*=B[i];
  A=fft(A,true);
  a.resize(s);
  REP(i,0,s)a[i]=A[i].real();
  return a;
}

int main(){
  long long n,d,x;cin>>n>>d>>x;
  ld temp=1;
  REP(i,1,n+1){
    temp*=(ld)(2*i-1)/(2*i);
    ans+=temp*(n-i+1);
  }
  vector<ld> a(n+1,0),b(n+1,0);
  a[1]=b[1]=0.5;
  REP(i,2,n+1)a[i]=a[i-1]*(2*i-1)*i/(ld)((2*i)*(i-1));
  REP(i,2,n+1)b[i]=b[i-1]*(2*i-3)/(ld)(2*i);
  vector<ld> c=convolve(a,b);
  REP(i,2,n+1)ans-=c[i]*(n-i+1)/(ld)i;
  ans*=(2*d+x*(2*n-1));
  cout<<fixed<<setprecision(12)<<ans<<endl;
}