好吧 我承认我只是在这里存一下EXCRT板子QAQ。
EXCRT推导过程:
其中(温爷爷起这种奇葩名字我也很绝望)
注意一下 因为exgcd
算出来的是的一组解,所以如果需要求原方程的解,需要给和都乘上一个。
#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;
}