PROBLEM LINK:
Practice
Contest: Division 1
Contest: Division 2
Contest: Division 3
Contest: Division 4
Author: gunpoint_88
Tester: mridulahi
Editorialist: iceknight1093
DIFFICULTY:
2585
PREREQUISITES:
Combinatorics
PROBLEM:
Magnus and Hikaru play infinitely many chess games.
Initially, Magnus’ rating is R_m and Hikaru’s is R_h. A player’s winning chance is directly proportional to their rating.
Winning a game increases a player’s rating by 1, losing doesn’t change it.
Find the probability that Magnus has won exactly P games and Hikaru has won exactly Q.
EXPLANATION:
The fact that winning probability is directly proportional to rating means that, if two players with ratings x and y play against each other, their chances of winning are \frac{x}{x+y} and \frac{y}{x+y}, respectively.
For Magnus to win exactly P games and Hikaru to win exactly Q, clearly exactly P+Q games must be played.
So, let’s find the probability that, after P+Q games, Magnus has won P.
Suppose we fix exactly which of the P games were won by Magnus: say, games x_1, x_2, \ldots, x_P,
The probability of this happening is as follows:
- Hikaru wins everything before the x_1-th game.
- So, he wins games 1, 2, 3 ,\ldots, x_1.
- For this to happen has a chance of
\displaystyle\frac{R_h}{R_m + R_h} \cdot \frac{R_h + 1}{R_m + R_h + 1} \cdot \frac{R_h + 2}{R_m + R_h + 2} \cdot \ldots \cdot \frac{R_h + x_i-2}{R_m + R_h + x_i - 2}
- Then, Magnus wins game x_i, with probability \frac{R_m}{R_m + R_h + x_i - 1}
- Then, Hikaru wins everything between x_i+1 and x_2 - 1; once again this can be written out into a product; and so on.
At first glance, this appears to depend on the x_i values chosen - but it actually doesn’t!
Let’s look at the numerators and denominators of the fractions we’re multiplying separately.
- For the i-th game played, the total rating increase of the players before it is exactly i-1.
So, the denominators will always be the sequence (R_m+R_h), (R_m+R_h+1), (R_m+R_h+2), \ldots, (R_m+R_h+P+Q-1).
This is independent of the x_i values.
If factorials are precomputed, their product can be found in \mathcal{O}(1) time since it’s just
\displaystyle \frac{(R_m + R_h + P + Q - 1)!}{(R_m + R_h - 1)!} - As for the numerators, they also follow a pattern.
- Let’s look at the games won by Magnus.
For the first one, his rating is R_m. For the second, it’s (R_m + 1). The third, (R_m+2), and so on — because only when winning does it change.
So, across all games won by Magnus, the product of ratings in the numerator is simply
R_m \cdot (R_m+1)\cdot (R_m+2)\cdot\ldots(R_m+P-1).
Once again, this can be computed quickly with precomputed factorials. - The games won by Hikaru follow a similar pattern, with their product being Q contiguous numbers starting from R_h.
- Let’s look at the games won by Magnus.
So, we’re able to quickly compute the probability that Magnus wins a fixed subset of P matches.
This number is independent of the actual subset, so to find the probability that Magnus wins some subset of P matches, we simply multiply this probability by the number of ways of picking a subset of P matches to win: which is just \binom{P+Q}{P}.
So, the final answer is just
where the last term is inverted because it belongs to the denominator.
If all factorials and inverse factorials are precomputed, this can be found in \mathcal{O}(1) time.
However, \mathcal{O}(\log{MOD}) per testcase is also fast enough (i.e, only precomputing factorials).
Note that the largest value we see is (R_m + R_h + P + Q - 1) which can be almost 4\cdot 10^6, so precomputation needs to be till there.
TIME COMPLEXITY
\mathcal{O}(1) per testcase after precomputation of 4\cdot 10^6 factorials and inverse factorials.
CODE:
Author's code (Python)
import sys
mod=int(1e9+7)
def fastexp(base,exp):
res=1
while exp:
res,base,exp=[res,(res*base)%mod][exp&1],(base*base)%mod,exp>>1
return res
def modinv(num):
return fastexp(num,mod-2)
factorial,inv_factorial=[],[]
def compute_factorial(n=1e6):
global factorial, inv_factorial
factorial,inv_factorial=[1 for i in range(n+1)],[1 for i in range(n+1)]
for i in range(2,n+1):
factorial[i]=factorial[i-1]*i%mod
inv_factorial[n]=modinv(factorial[n])
for i in range(n-1,0,-1):
inv_factorial[i]=inv_factorial[i+1]*(i+1)%mod
def mul(L,R):
return factorial[R-1]*inv_factorial[L-1]%mod
def solve(p,q,rm,rh):
num=mul(rm,rm+p)*mul(rh,rh+q)%mod
den=mul(rm+rh,rm+rh+p+q);
prob=num*modinv(den)%mod;
ways=(mul(1,p+q+1)*modinv(mul(1,p+1))%mod)*modinv(mul(1,q+1))%mod
return prob*ways%mod;
def read(T=int):
return [T(i) for i in sys.stdin.readline().split()]
def main():
def work():
p,q,rm,rh=read()
print(solve(p,q,rm,rh))
compute_factorial(int(6e6))
for i in range(read(int)[0]):
work()
main()
Tester's code (C++)
//#pragma GCC optimize("O3")
//#pragma GCC optimize("Ofast")
//#pragma GCC optimize("unroll-loops")
//#pragma GCC target("avx2,bmi,bmi2,lzcnt,popcnt")
#include<bits/stdc++.h>
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
using namespace std;
using namespace __gnu_pbds;
struct custom_hash {
static uint64_t splitmix64(uint64_t x) {
// http://xorshift.di.unimi.it/splitmix64.c
x += 0x9e3779b97f4a7c15;
x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9;
x = (x ^ (x >> 27)) * 0x94d049bb133111eb;
return x ^ (x >> 31);
}
size_t operator()(uint64_t x) const {
static const uint64_t FIXED_RANDOM = chrono::steady_clock::now().time_since_epoch().count();
return splitmix64(x + FIXED_RANDOM);
}
};
mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
template<typename T>
T rand(T a, T b){
return uniform_int_distribution<T>(a, b)(rng);
}
typedef tree<int, null_type, less<int>, rb_tree_tag, tree_order_statistics_node_update> ordered_set;
typedef tree<pair<int,int>, null_type, less<pair<int,int>>, rb_tree_tag, tree_order_statistics_node_update> ordered_multiset;
typedef long long ll;
typedef vector<ll> vl;
typedef vector<int> vi;
#define rep(i, a, b) for(int i = a; i < b; i++)
#define all(x) begin(x), end(x)
#define sz(x) static_cast<int>((x).size())
#define int long long
const ll mod = 1e9 + 7;
const ll INF = 1e18;
/* ----------------------------------------------------- GO DOWN ---------------------------------------------------------------------- */
int norm (int x) {
if (x < 0) {
x += mod;
}
if (x >= mod) {
x -= mod;
}
return x;
}
template<class T>
T power(T a, int b) {
T res = 1;
for (; b; b /= 2, a *= a) {
if (b % 2) {
res *= a;
}
}
return res;
}
struct Z {
int x;
Z(int x = 0) : x(norm(x)) {}
int val() const {
return x;
}
Z operator-() const {
return Z(norm(mod - x));
}
Z inv() const {
assert(x != 0);
return power(*this, mod - 2);
}
Z &operator*=(const Z &rhs) {
x = x * rhs.x % mod;
return *this;
}
Z &operator+=(const Z &rhs) {
x = norm(x + rhs.x);
return *this;
}
Z &operator-=(const Z &rhs) {
x = norm(x - rhs.x);
return *this;
}
Z &operator/=(const Z &rhs) {
return *this *= rhs.inv();
}
friend Z operator*(const Z &lhs, const Z &rhs) {
Z res = lhs;
res *= rhs;
return res;
}
friend Z operator+(const Z &lhs, const Z &rhs) {
Z res = lhs;
res += rhs;
return res;
}
friend Z operator-(const Z &lhs, const Z &rhs) {
Z res = lhs;
res -= rhs;
return res;
}
friend Z operator/(const Z &lhs, const Z &rhs) {
Z res = lhs;
res /= rhs;
return res;
}
friend std::istream &operator>>(std::istream &is, Z &a) {
int v;
is >> v;
a = Z(v);
return is;
}
friend std::ostream &operator<<(std::ostream &os, const Z &a) {
return os << a.val();
}
};
const int maxn = 4e6;
Z fact[maxn];
Z ifact[maxn];
Z C (int n, int r) {
return fact[n] * ifact[r] * ifact[n - r];
}
signed main() {
ios::sync_with_stdio(0);
cin.tie(0);
cout.tie(0);
fact[0] = 1;
for (int i = 1; i < maxn; i++) fact[i] = fact[i - 1] * i;
ifact[maxn - 1] = Z(1) / fact[maxn - 1];
for (int i = maxn - 2; i >= 0; i--) ifact[i] = ifact[i + 1] * (i + 1);
int t;
cin >> t;
while (t--) {
int p, q, a, b;
cin >> p >> q >> a >> b;
Z ans = 1;
ans *= C(p + q, p);
ans *= fact[a + p - 1];
ans /= fact[a - 1];
ans *= fact[b + q - 1];
ans /= fact[b - 1];
ans /= fact[p + q + a + b - 1];
ans *= fact[a + b - 1];
cout << ans << "\n";
}
}
Editorialist's code (Python)
mod = 10**9 + 7
maxN = 4* 10**6 + 100
fac = [1]*maxN
for i in range(2, maxN): fac[i] = fac[i-1] * i % mod
def C(n, r):
if n < r or r < 0: return 0
return fac[n] * pow(fac[r] * fac[n-r], mod-2, mod) % mod
for _ in range(int(input())):
p, q, m, h = map(int, input().split())
ans = fac[p + m - 1] * fac[q + h - 1] % mod
ans = ans * pow(fac[m-1] * fac[h-1], mod-2, mod) % mod
ans = ans * C(p + q, p) % mod
ans = ans * pow(fac[p + q + m + h - 1], mod-2, mod) % mod
ans = ans * fac[m + h - 1] % mod
print(ans)