北大集训2019 n扇门问题（有趣的概率题）

在2020 CCPC-Wannafly Winter Camp Day3发现了这题，应该是一个出题人。

https://ac.nowcoder.com/acm/contest/4114/I

n扇门问题：

有一个猜奖者和一个主持人，一共有n扇门，只有一扇门后面有奖，主持人事先知道哪扇门后有奖，而猜奖者不知道。

每一轮，猜奖者选择它认为的有奖概率最大（如果有多个最大，随机选一个）的一扇门，主持人从剩下的且门后没有奖的门中随机打开一扇。

直到剩两扇门时，猜奖者做出的选择就是他最后的选择。

现在由你来安排主持人每次打开哪一扇门，猜奖者不知道有内幕，他还认为主持人是从可以打开的门中随机一扇打开。

你要使猜奖者获奖概率最低，求这个概率。

step1：

如何实时计算：从猜奖者的视角，每一扇门有奖的概率？

假设现在还有(n)扇门，第(i)扇有奖的概率是(p[i])。

猜奖者选了第(x)扇门，主持人打开了第(y)扇门：

1.有(p[x])的概率第(x)扇门就是有奖的，经过这次操作，(p[x])显然不会变。

2.有(1-p[x])的概率奖不在第(x)扇门，现在又多排除了第(y)扇门，

所以对(z≠x且z≠y,p[z]*={1-p[x]over 1-p[x]-p[y]})，意义为剩下的门均分这个这个多出来的概率。

step2：

每次选择的门，这次操作后，都会成为概率最小且独一无二的门。

这个可以归纳证明。

考虑对于任意时候的(n)扇门，设(p[x])为最大概率，(p[y])为最小概率。

上面的命题相当于证明任意时候有：
(p[x]<p[y]*{1-p[x]over 1 - p[x] - p[y]})

(px(1-px)<py)

第一轮显然满足，对于下一轮：
(px'<=px*{1-p[x]over 1-p[x]-p[y]})

(py'=px)

代进去应该就可以证明(px')和(py')也满足。

step3：

问题变成了：

有一个队列，一开始队列头一个位置上有(n)个球。

每次从队列头随机一个球(x)，然后从剩下的且不是答案的球去掉一个，把球x单独加到队尾。

最后剩两个球就清晰了，当一开始的(n>2)时，队列头是答案球就赢了。

那么设(f[i][j][k])表示一共还剩(i)个球，队列头有(j)个球，答案球在队列第k个位置上。

得到一个转移(O(1))，总复杂度(O(n^3))的dp。

step4：

不难想到(n)比较大时候，主持人可以通过控制奇偶，使得剩两个球的时候，答案球一定在队尾。

实际上n>10猜奖者赢的概率就是0了，其实直接搜索也能搜出10以内的答案。

step5：

原题强行加了个扩展中国剩余定理。

Code：

#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, _b = y; i <= _b; i ++)
#define ff(i, x, y) for(int i = x, _b = y; i <  _b; i ++)
#define fd(i, x, y) for(int i = x, _b = y; i >= _b; i --)
#define ll long long
#define pp printf
#define hh pp("
")
using namespace std;
 
#define db double
 
const int N = 205;
 
db f[N][N][N];
 
void dp() {
    f[2][0][1] = 1; f[2][0][2] = 0;
    f[2][1][1] = 1; f[2][1][2] = 0;
    f[2][2][1] = 0.5; f[2][2][2] = 0.5;
    fo(i, 3, 200) {
        fo(j, 0, i) {
            fo(k, 1, i) {
                if(j > 0) {
                    if(k == 1) {
                        db u = 1;
                        if(j > 1) u = min(u, f[i - 1][j - 2][i - 1]);
                        if(j < i) u = min(u, f[i - 1][j - 1][i - 1]);
                         
                        db v = 1;
                        if(j > 2) v = min(v, f[i - 1][j - 2][1]);
                        if(j < i) v = min(v, f[i - 1][j - 1][1]);
                         
                        f[i][j][k] = u / j + v / j * (j - 1);
                    } else {
                        db v = 1;
                        if(k - 1 > j) v = min(v, f[i - 1][j - 1][k - 2]);
                        if(k < i) v = min(v, f[i - 1][j - 1][k - 1]);
                        if(j > 1) v = min(v, f[i - 1][j - 2][k - 2]);
                        f[i][j][k] = v;
                    }
                } else {
                    if(k == 1) f[i][j][k] = f[i - 1][j][i - 1]; else
                    if(k == i) f[i][j][k] = f[i - 1][j][i - 2]; else {
                        f[i][j][k] = f[i - 1][j][k - 1];
                        if(k > 2) f[i][j][k] = min(f[i][j][k], f[i - 1][j][k - 2]);
                    }
                }
            }
        }
    }
}
 
int T;
ll m1, c1, m2, c2;
 
ll mul(ll x, ll y, ll mo) {
    x %= mo, y %= mo;
    ll z = (long double) x * y / mo;
    z = x * y - z * mo;
    if(z < 0) z += mo; else if(z >= mo) z -= mo;
    return z;
}
 
ll gcd(ll x, ll y) {
    return !y ? x : gcd(y, x % y);
}
void exgcd(ll a, ll b, ll &x, ll &y) {
    if(!b) { x = 1, y = 0; return;}
    exgcd(b, a % b, y, x); y -= (a / b) * x;
}
ll inv(ll a, ll b) {
    ll x, y;
    exgcd(a, b, x, y);
    x = (x % b + b) % b;
    return x;
}
 
void mer() {
    ll d = gcd(m1, m2);
    if((c2 - c1) % d != 0) {
        pp("error
");
        exit(0);
    }
    m1 /= d, m2 /= d;
    ll c = (c2 - c1) / d, M = m1 * m2 * d;
    c1 = (mul(inv(m1, m2), (c2 - c1) / d, m2) * (m1 * d) % M + c1) % M;
    c1 = (c1 % M + M) % M; m1 = M;
}
 
int main() {
    dp();
    m1 = 1, c1 = 0;
    for(scanf("%d", &T); T; T --) {
        scanf("%lld %lld", &c2, &m2);
        mer();
    }
    if(c1 < 2) {
        pp("error
");
        exit(0);
    }
    if(c1 <= 10) pp("%.6lf
", f[c1][c1][1]); else
        pp("0.000000
");
}