# CHESSELO - Editorial

Author: gunpoint_88
Tester: mridulahi
Editorialist: iceknight1093

2585

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.

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

\binom{P+Q}{Q}\cdot \frac{(R_m + P - 1)!}{(R_m-1)!} \cdot \frac{(R_h + Q - 1)!}{(R_h-1)!} \cdot \left(\frac{(R_m + R_h + P + Q - 1)!}{(R_m + R_h - 1)!}\right)^{-1}

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;

return [T(i) for i in sys.stdin.readline().split()]

def main():
def work():
print(solve(p,q,rm,rh))

compute_factorial(int(6e6))
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);
}
};
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 = 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 = *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)

1 Like

Can someone provide me blog/ material to read about taking X * Y^(-1) mod (M), would really appreciate that. I wanna undersand the Z struct is it same as mint?

There are many articles, the term you’re looking for is “modular inverse”.
For example, cp-algorithms has a page on it.

Assuming you’re talking about the tester’s code, yes.

1 Like