北大集训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("\n")
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\n");
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\n");
exit(0);
}
if(c1 <= 10) pp("%.6lf\n", f[c1][c1][1]); else
pp("0.000000\n");
}
