Project Euler 430

题目

Range flips

N disks are placed in a row, indexed \(1\) to \(N\) from left to right.

Each disk has a black side and white side. Initially all disks show their white side.

At each turn, two, not necessarily distinct, integers \(A\) and \(B\) between \(1\) and \(N\) (inclusive) are chosen uniformly at random.

All disks with an index from \(A\) to \(B\) (inclusive) are flipped.

The following example shows the case \(N = 8\). At the first turn \(A = 5\) and \(B = 2\), and at the second turn \(A = 4\) and \(B = 6\).

Let \(E(N, M)\) be the expected number of disks that show their white side after \(M\) turns.

We can verify that \(E(3, 1) = 10/9, E(3, 2) = 5/3, E(10, 4) \approx 5.157\) and \(E(100, 10) \approx 51.893\).

Find \(E(10^{10}, 4000)\).

Give your answer rounded to \(2\) decimal places behind the decimal point.

解决方案

我们将独立考虑每一个棋子的情况。

令\(I_{i,m}\)表示示性随机变量：经过\(m\)轮翻转后，棋子仍然白色面朝上。（也就是说，发生了偶数次的翻转）

对于第\(i\)个棋子，它在某一轮被翻转一次的概率为\(p_i=\dfrac{2i(n-i+1)-1}{n^2}\).

那么第\(i\)个棋子在这\(m\)次翻转中，其翻转次数\(X_{i}\)服从二项分布\(B(m,p_i)\)。令\(q_i=1-p_i\)

结合

\(\begin{aligned} &(p+q)^m=\sum_{i=0}^m\dbinom{m}{i}p^iq^{m-i}\\ &(-p+q)^m=\sum_{i=0}^m\dbinom{m}{i}(-p)^iq^{m-i}\\ \end{aligned}\)

这两个式子可以得到\(X\)为偶数时的概率，也就是有

\(\begin{aligned} I_{i,m}&=P\{X\text{ is even}\}\\ &=\dfrac{1}{2}((q_i+p_i)^m+(q_i-p_i)^m)\\ &=\dfrac{1}{2}(1+(1-2p_i)^m) \end{aligned}\)

那么

\(\begin{aligned} E(n,m)&=\sum_{i=1}^{n}I_{i,m}\\ &=\sum_{i=1}^{n}\dfrac{1}{2}(1+(1-2p_i)^m)\\ &=\dfrac{n}{2}+\dfrac{1}{2}\sum_{1=1}^n(1-2p_i)^m \end{aligned}\)

这时我们考虑计算\(\sum_{1=1}^n(1-2p_i)^m\)的值。当\(i\approx\dfrac{n}{2}\)时，\(p_i\approx\dfrac{1}{2}\)，此时级数靠中央的所有项几乎可以忽略。因此，直接暴力进行计算即可。

不难证明\(p_i=p_{n+1-i}\)，因此我们只需要计算一侧的结果即可。

代码

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

const ll N=10000000000;
const int M=4000;
const int F=2;
double eps= exp10(-F-20);
double g(ll i,ll n) {
    return (2.0 * i * (n - i + 1) - 1) / n / n;

}
int main(){
    double ans=0.5*N;
    double s=0;
    if(N<=1000000){
        for(int i=1;i<=N;i++){
            s+=pow(1.0-g(i,N)*2,M);
        }
        ans+=s/2;
    }
    else{
        ll i;
        for(i=1;i+i<=N;i++){
            double w=pow(1.0-g(i,N)*2,M);
            if(w<eps) break;
            s+=w;
        }
        ans+=s;
        if(i+i==N+1) ans+=0.5*pow(1.0-g(i,N)*2,M);
    }
    printf("%.*f\n",F,ans);
}