p498.ans1
======================
0.7142857142
0.5000000000
0.5333333333
0.8000000000

=================
p498.cpp
======================
#include <bits/stdc++.h>

using namespace std;

template<typename T1, typename T2>
ostream& operator<<(ostream& out, const pair<T1, T2>& x) {
    return out << x.first << ' ' << x.second;
}

template<typename T1, typename T2>
istream& operator>>(istream& in, pair<T1, T2>& x) {
    return in >> x.first >> x.second;
}

template<typename T>
istream& operator>>(istream& in, vector<T>& a) {
    for(auto& x: a) {
        in >> x;
    }
    return in;
};

template<typename T>
ostream& operator<<(ostream& out, const vector<T>& a) {
    for(auto x: a) {
        out << x << ' ';
    }
    return out;
};

int n1, m1, n2, m2;

void read() { cin >> n1 >> m1 >> n2 >> m2; }

void solve() {
    // We can solve this problem with a Bayesian approach:
    //
    // The prior is: p, q ~ Uniform(0, 1), independent
    // This means f(p) = f(q) = 1 for p, q in [0, 1].
    //
    // Then the likelihood is:
    //    P(data | p, q) = P(m1 heads in n1 tosses | p) *
    //                     P(m2 heads in n2 tosses | q)
    //                   = C(n1, m1) * p^{m1} * (1-p)^{n1-m1} *
    //                     C(n2, m2) * q^{m2} * (1-q)^{n2-m2}
    //
    // Then by Bayes' theorem the posterior will be (ignoring normalizing
    // constants):
    //    P(p, q | data) ~ p^{m1} * (1-p)^{n1-m1} * q^{m2} * (1-q)^{n2-m2}
    //
    // The normalizing constant is Z = P(data), and can be calculated as:
    //     P(data) = integral_0^1 integral_0^1 p^{m1}(1-p)^{n1-m1}
    //                                         q^{m2}(1-q)^{n2-m2} dp dq
    //             = B(m1+1, n1-m1+1) * B(m2+1, n2-m2+1)
    //
    // where B(a, b) is the Beta function:
    //     B(a, b) = integral_0^1 t^{a-1}(1-t)^{b-1} dt
    //             = (a-1)!(b-1)! / (a+b-1)!
    //               (for positive integers)
    //
    // See the Wikipedia page for the identity:
    // https://en.wikipedia.org/wiki/Beta_function
    //
    // Then the integral we want is:
    //   P(p < q | data) = (1/Z) * integral_{p<q} p^{m1}(1-p)^{n1-m1}
    //                                            q^{m2}(1-q)^{n2-m2} dp dq
    //
    // The inner integral is the CDF of Beta(m1+1, n1-m1+1) at q (times Z_p).
    // For integer parameters, the Beta CDF has a closed form:
    //   F_p(q) = sum_{j=m1+1}^{n1+1} C(n1+1, j) * q^j * (1-q)^{n1+1-j}
    //
    // So the inner integral (unnormalized) is:
    //   integral_0^q p^{m1}(1-p)^{n1-m1} dp = B(m1+1, n1-m1+1) * F_p(q)
    //
    // Substituting back and using B(m1+1, n1-m1+1) * B(m2+1, n2-m2+1) = Z:
    //   P(p < q | data) = (1/B(m2+1, n2-m2+1)) *
    //                      sum_{j=m1+1}^{n1+1} C(n1+1, j) *
    //                                          integral_0^1 q^{j+m2}
    //                                                       (1-q)^{n1+n2-m2-j+1}
    //                                                       dq
    //
    //                   = (1/B(m2+1, n2-m2+1)) *
    //                     sum_{j=m1+1}^{n1+1} C(n1+1, j) *
    //                                         B(j+m2+1, n1+n2-m2-j+2)
    //
    // This gives us a O(N) solution per test case.

    long double term = 1.0;
    // C(n1+1, m1+1) = prod_{i=0}^{m1} (n1+1-i)/(i+1)
    for(int i = 0; i <= m1; i++) {
        term *= (long double)(n1 + 1 - i) / (i + 1);
    }

    // B(m1+m2+2, n1+n2-m1-m2+1) / B(m2+1, n2-m2+1)
    // = [(m1+m2+1)!/m2!] * [(n1+n2-m1-m2)!/(n2-m2)!] * [(n2+1)!/(n1+n2+2)!]
    for(int i = 1; i <= m1 + 1; i++) {
        term *= (long double)(m2 + i);
    }
    for(int i = 1; i <= n1 - m1; i++) {
        term *= (long double)(n2 - m2 + i);
    }
    for(int i = 1; i <= n1 + 1; i++) {
        term /= (long double)(n2 + 1 + i);
    }

    // T(j+1)/T(j) = (n1+1-j)*(j+m2+1) / [(j+1)*(n1+n2-m2-j+1)]
    long double ans = term;
    for(int j = m1 + 1; j <= n1; j++) {
        term *= (long double)(n1 + 1 - j) * (j + m2 + 1);
        term /= (long double)(j + 1) * (n1 + n2 - m2 - j + 1);
        ans += term;
    }

    cout << fixed << setprecision(10) << ans << '\n';
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    int T = 1;
    cin >> T;
    for(int test = 1; test <= T; test++) {
        read();
        // cout << "Case #" << test << ": ";
        solve();
    }

    return 0;
}

=================
p498.in1
======================
4
2 1 4 3
8 4 16 8
2 0 6 1
2 0 2 1

=================
statement.txt
======================
498. Coins
Time limit per test: 0.75 second(s)
Memory limit: 262144 kilobytes
input: standard
output: standard

Andrew has just made a breakthrough in gambling: he has invented a way to make coins which have probability p of landing heads-up and 1 - p of landing tails-up when tossed, where p is some number between 0 and 1. He can't, however, control the number p itself, and creates coins with any p from 0 to 1 with equal probability. That is, the value of p is a random value uniformly distributed on . Andrew has generated two coins independently. One with probability p and another with probability q of landing heads-up. Random values p and q are both uniformly distributed on  and are independent. Of course, neither Andrew nor we know the numbers p and q, we can only try to guess them using our observations. The observations are the following: the first coin was tossed n1 times, and m1 of them landed heads-up. The second coin was tossed n2 times, and m2 of them landed heads-up. Your task is to compute the probability that p < q.
Input
The first line of the input file contains one integer T () — the number of test cases to solve. Each of the following T lines contains 4 integers each: n1, m1, n2, m2. 1 ≤ n1, n2 ≤ 1000, 0 ≤ m1, m2 ≤ 50, 0 ≤ m1 ≤ n1, 0 ≤ m2 ≤ n2.
Output
For each test case output one line with a floating-point number, the probability of p < q. Your answer will be considered correct if it is within 10-4 of the right answer.
Example(s)
sample input
sample output
4
2 1 4 3
8 4 16 8
2 0 6 1
2 0 2 1
0.7142857142
0.5000000000
0.5333333333
0.8000000000

=================
