<|instruction|>
Solve the below problem. The solution should start with an abridged problem statement. Then key observations. Then full solution based on the observations. Then C++ and Python implementations with comments.

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

<|response|>
## 1) Abridged problem statement (concise)

You observe results of tossing two different coins:

- Coin 1: unknown head probability \(p\), tossed \(n_1\) times, got \(m_1\) heads
- Coin 2: unknown head probability \(q\), tossed \(n_2\) times, got \(m_2\) heads

Before observing data, \(p\) and \(q\) are independent and uniformly random on \([0,1]\).

For each test case, compute the posterior probability:
\[
\Pr(p < q \mid \text{observations})
\]
Output as a floating-point number (absolute error \(\le 10^{-4}\)).

---

## 2) Key observations

1. **Uniform prior + binomial observations ⇒ Beta posterior**
   Uniform\([0,1]\) is \(\text{Beta}(1,1)\). After observing \(m\) heads in \(n\) tosses:
   \[
   p \mid \text{data} \sim \text{Beta}(m+1,\, n-m+1)
   \]

2. **Independence remains after conditioning**
   Since the two experiments are independent and priors are independent:
   \[
   (p \mid \text{coin1 data}) \perp (q \mid \text{coin2 data})
   \]

3. **We need \(\Pr(p<q)\) for two independent Beta variables**
   This can be expressed as:
   \[
   \Pr(p<q)=\int_0^1 F_p(x)\, f_q(x)\, dx
   \]
   where \(F_p\) is the Beta CDF of coin 1, and \(f_q\) is the Beta PDF of coin 2.

4. **For integer parameters, Beta CDF becomes a binomial-tail polynomial**
   For \(p\sim \text{Beta}(a,b)\) with integers \(a,b\):
   \[
   F_p(x)=I_x(a,b)=\sum_{j=a}^{a+b-1} \binom{a+b-1}{j} x^j(1-x)^{a+b-1-j}
   \]
   Here \(a=m_1+1\), \(b=n_1-m_1+1\), so \(a+b-1=n_1+1\).

5. **After swapping sum and integral, each integral is a Beta function**
   Result becomes a single sum with Beta functions, computable in \(O(n_1)\).

6. **Avoid factorial overflow using a multiplicative recurrence**
   Compute the first term carefully by products of ratios, then iterate:
   \[
   \frac{T(j+1)}{T(j)}=
   \frac{(n_1+1-j)(j+m_2+1)}{(j+1)(n_1+n_2-m_2-j+1)}
   \]
   This is stable with `long double`.

---

## 3) Full solution approach

### Step A: Posterior distributions
With uniform priors:
- Coin 1:
  \[
  p \mid data \sim \text{Beta}(a_1,b_1),\quad a_1=m_1+1,\; b_1=n_1-m_1+1
  \]
- Coin 2:
  \[
  q \mid data \sim \text{Beta}(a_2,b_2),\quad a_2=m_2+1,\; b_2=n_2-m_2+1
  \]

We need \(\Pr(p<q)\) where these two Beta variables are independent.

### Step B: Turn into one-dimensional integral
\[
\Pr(p<q)=\int_0^1 \left(\int_0^q f_p(p)\,dp\right) f_q(q)\,dq
= \int_0^1 F_p(q)\, f_q(q)\, dq
\]

### Step C: Use integer-parameter Beta CDF expansion
Since \(a_1,b_1\) are integers,
\[
F_p(q)=\sum_{j=m_1+1}^{n_1+1} \binom{n_1+1}{j} q^j(1-q)^{n_1+1-j}
\]

Plug into integral and swap sum/integral:
\[
\Pr(p<q)=\sum_{j=m_1+1}^{n_1+1} \binom{n_1+1}{j}
\int_0^1 q^j(1-q)^{n_1+1-j}\, f_q(q)\,dq
\]

### Step D: Each integral is a Beta function
Beta PDF of \(q\):
\[
f_q(q)=\frac{1}{B(m_2+1,n_2-m_2+1)} q^{m_2}(1-q)^{n_2-m_2}
\]

Combine exponents:
\[
\int_0^1 q^{j+m_2}(1-q)^{(n_1+1-j)+(n_2-m_2)}\,dq
= B(j+m_2+1,\; n_1+n_2-m_2-j+2)
\]

Therefore:
\[
\Pr(p<q)=\frac{1}{B(m_2+1,n_2-m_2+1)}
\sum_{j=m_1+1}^{n_1+1}\binom{n_1+1}{j}\,
B(j+m_2+1,\; n_1+n_2-m_2-j+2)
\]

### Step E: Compute in \(O(n_1)\) via recurrence
Compute \(T(m_1+1)\) by multiplicative products (no factorials), then:
\[
T(j+1)=T(j)\cdot
\frac{(n_1+1-j)(j+m_2+1)}{(j+1)(n_1+n_2-m_2-j+1)}
\]
Sum all terms.

**Complexity:** \(O(n_1)\) per test case, with \(n_1\le 1000\).
**Precision:** use `long double` and print ~10 decimals.

---

## 4) C++ Solution

```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;
}
```

---

## 5) Python implementation (detailed comments)

```python
import sys

def solve_case(n1, m1, n2, m2) -> float:
    """
    Computes P(p < q | data) where:
      p|data ~ Beta(m1+1, n1-m1+1)
      q|data ~ Beta(m2+1, n2-m2+1)

    Uses the same O(n1) summation and recurrence as the C++ code,
    keeping everything as floating-point products (no factorials).
    """
    j0 = m1 + 1

    term = 1.0

    # Build C(n1+1, j0) = C(n1+1, m1+1) multiplicatively
    for i in range(m1 + 1):  # i=0..m1
        term *= (n1 + 1 - i) / (i + 1)

    # Multiply by product expansion of:
    # B(j0+m2+1, n1+n2-m2-j0+2) / B(m2+1, n2-m2+1)

    # Multiply by (m2+1)(m2+2)...(m2+m1+1)
    for i in range(1, m1 + 2):
        term *= (m2 + i)

    # Multiply by (n2-m2+1)...(n2-m2 + (n1-m1))
    for i in range(1, n1 - m1 + 1):
        term *= (n2 - m2 + i)

    # Divide by (n2+2)(n2+3)...(n2+n1+2)
    for i in range(1, n1 + 2):
        term /= (n2 + 1 + i)

    ans = term

    # Recurrence:
    # T(j+1)/T(j) = (n1+1-j)(j+m2+1) / ((j+1)(n1+n2-m2-j+1))
    for j in range(j0, n1 + 1):  # j=j0..n1
        term *= (n1 + 1 - j) * (j + m2 + 1)
        term /= (j + 1) * (n1 + n2 - m2 - j + 1)
        ans += term

    return ans


def main() -> None:
    data = sys.stdin.read().strip().split()
    if not data:
        return
    t = int(data[0])
    idx = 1
    out = []
    for _ in range(t):
        n1 = int(data[idx]); m1 = int(data[idx+1])
        n2 = int(data[idx+2]); m2 = int(data[idx+3])
        idx += 4
        out.append(f"{solve_case(n1, m1, n2, m2):.10f}")
    sys.stdout.write("\n".join(out))


if __name__ == "__main__":
    main()
```
