[北大集训2019 D4T2] N门问题

好吧 我承认我只是在这里存一下EXCRT板子QAQ。
EXCRT推导过程：


其中$orz=k1$（温爷爷起这种奇葩名字我也很绝望）
注意一下 因为exgcd算出来的是$ax+by=gcd(a,b)$的一组解，所以如果需要求原方程的解，需要给$x$和$y$都乘上一个$\frac{c}{gcd(a,b)}$。

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <vector>
#include <queue>
#include <set>
#define vi vector<int>
#define pb push_back
#define mk make_pair
#define pii pair<int, int>
#define rep(i, a, b) for (int i = (a), i##end = (b); i <= i##end; i++)
#define fi first
#define se second
#define ldb long double
typedef long long ll;
using namespace std;
const ldb eps = 1e-8;
void exgcd(ll a, ll b, ll &d, ll &x, ll &y) {
    if (!b) {
        x = 1, y = 0, d = a;
        return;
    }
    exgcd(b, a % b, d, y, x);
    y -= x * (a / b);
}
void magic(vector<ldb> &p) {
    ldb tmp = 0;
    for (int i = 0; i < p.size(); i++) tmp += p[i];
    for (int i = 0; i < p.size(); i++) p[i] /= tmp;
}
ldb dfs(vector<ldb> &p, int remain) {
    //	for(int i=0;i<p.size();i++)cout<<p[i]<<' ';
    //	cout<<endl;
    //	cout<<p[0]<<endl;
    if (remain == 2)
        if (fabs(p[0] - p[1]) < eps)
            return 0.5;
        else
            return (p[0] < p[1]) ? 0 : 1;

    ldb mx = *max_element(p.begin(), p.end());
    // cout<<">>"<<p.size()<<' '<<remain<<' '<<mx<<endl;
    int ct = 0, pos;
    rep(i, 0, remain - 1) if (fabs(p[i] - mx) < eps) ct++, pos = i;
    ldb mn1 = 1, mn2 = 1, res = 0;
    if (fabs(mx - p[0]) < eps) {  //火山哥选的是第一个
        rep(i, 1, remain - 1) {   //主持人开哪个
            vector<ldb> tmp;
            tmp.pb(p[0] / (remain - 1));  //除了第一个节点 剩下的都可能开
            rep(j, 1, remain - 1) if (i != j) {
                tmp.pb(p[j] / (remain - 2));  //？
            }
            magic(tmp);
            mn1 = min(mn1, dfs(tmp, remain - 1));
        }
        res += mn1 / ct;
    }
    if (pos) {
        rep(i, 1, remain - 1) if (i != pos) {
            vector<ldb> tmp;
            rep(j, 0, remain - 1) if (i != j) {
                if (j == pos)
                    tmp.pb(p[j] / (remain - 1));
                else
                    tmp.pb(p[j] / (remain - 2));
            }
            magic(tmp);
            mn2 = min(mn2, dfs(tmp, remain - 1));
        }
        if (fabs(mx - p[0]) < eps)
            res += mn2 * (ct - 1) / ct;
        else
            res += mn2;
    }
    return res;
}
struct gg {
    ll num, mod;
} equ[50000 + 100];
ll mul(ll x, ll y, ll mod) {
    ll tmp = (ldb)x * y / mod;
    return ((x * y - tmp * mod) % mod + mod) % mod;
}
bool excrt(int n) {
    rep(i, 2, n) {
        gg st = equ[1];
        ll x, y, d;
        exgcd(st.mod, equ[i].mod, d, x, y);
        y = -y;
        if ((equ[i].num - st.num) % d)
            return 1;
        equ[1].mod = st.mod / d * equ[i].mod;
        equ[1].num =
            (st.num + mul(mul((equ[i].num - st.num) / d, x, equ[1].mod), st.mod, equ[1].mod)) % equ[1].mod;
    }
    return 0;
}
int main() {
    // freopen("in","r",stdin);
    ll n;
    scanf("%lld", &n);
    rep(i, 1, n) scanf("%lld%lld", &equ[i].num, &equ[i].mod);
    if (excrt(n)) {
        printf("error");
        return 0;
    }
    n = equ[1].num;
    if (n < 2) {
        printf("error");
        return 0;
    }
    if (n > 10) {
        printf("%lf", 0);
        return 0;
    }
    vector<ldb> tmp;
    rep(i, 1, n) tmp.pb(1.0 / n);
    //	cout<<tmp[0]<<' '<<n<<endl;
    printf("%.6lf", (double)dfs(tmp, n));
    return 0;
}