LUCMAT - Editorial


Contest: Division 1
Contest: Division 2
Contest: Division 3
Contest: Division 4

Author: indreshsingh
Tester: airths
Editorialist: iceknight1093




Probability and expected value, binomial coefficients


There’s an M\times N boolean matrix B.
Each entry of B is 1 with probability p = \frac{a}{b}, and 0 otherwise.

A boolean array A of length N is said to be good if \sum_{j=1}^N B_{ij}A_j \equiv 0 \pmod 2 for every i from 1 to M.
Let E_s be the expected number of good arrays whose sum is s.
Given k, find \sum_{s=1}^k E_s.


Both A and B hold boolean values. So, B_{ij}A_j = 1 if and only if B_{ij} and A_j are both 1 themselves; otherwise the product is 0 mod 2.

That means the sum \sum_{j=1}^N B_{ij}A_j denotes the number of positions where the arrays B_i and A both have ones.
For this to be 0\pmod 2 means that the number of such positions must be even.

Now, let’s consider some boolean array A that has s ones.
Rather than the entire matrix B, let’s look only at its first row, B_1.
In this first row,

  • There are N-s positions of A that contain zeros. The corresponding values in B_1 can be anything at all - they won’t affect anything.
  • Among the s positions of A that are ones, B_1 must contain ones at an even number of them.
    Let’s fix which 2r positions of B_1 contain these ones (where 0 \leq 2r \leq s).
    • There are \binom{s}{2r} ways to choose these positions 2r positions.
    • The probability that each of them ends up being 1 is p^{2r}.
    • The probability that the other s-2r positions are all zeros is (1-p)^{s-2r}.

So, the probability that the first row satisfies the required condition is

\sum_{r=0}^{\left\lfloor \frac{s}{2} \right\rfloor} \binom{s}{2r} p^{2r} (1-p)^{s-2r}

Using the binomial theorem and a bit of algebra, this summation can be simplified into a closed form.


Recall the binomial theorem:

(x+y)^N = \sum_{k=0}^N \binom{N}{k} x^k y^{N-k}

Choosing x = p and y = 1-p along with N = s gives us

(p + (1-p))^s = \sum_{k=0}^s \binom{s}{k} p^k (1-p)^{s-k}

Note that the left side evaluates to just 1, while the right side looks rather similar to the quantity we want to compute: just that we only care about the even-indexed values of that summation, rather than all of them.

To get rid of odd-indexed values, we use another standard trick to bring in subtraction via powers of -1: consider a = -p and b = 1-p to get

(-p + (1-p))^s = \sum_{k=0}^s \binom{s}{k} (-p)^k (1-p)^{s-k} \\ = \sum_{k=0}^s \binom{s}{k} (-1)^k p^k (1-p)^{s-k}

Note that on adding the two summations we have, the even-indexed terms are the same, while the odd-indexed terms cancel out!
So, we obtain

1 + (1 - 2p)^s = \sum_{\substack{k = 0 \\ k \text{ even}}}^{s} 2\cdot \binom{s}{k} p^k (1-p)^{s-k}

The right-hand side is now (twice) the summation we want, while the left hand side is easily computed, being a simple expression.

So, the value we’re looking for is exactly

\frac{1}{2} \cdot \left(1 + (1 - 2p)^s \right)

Now, observe that while this computation was done for the first row, it can really be repeated for any row - they’re all independent of each other, since the only condition we care about is row-specific.
So, the probability that A is good for the entire matrix is simply the probability for a single row, raised to the power M - that is,

\left(\frac{1}{2} \cdot \left(1 + (1 - 2p)^s \right) \right)^M

Next, to compute E_s,

  • There are \binom{N}{s} ways to choose a set of s positions of A that should contain ones.
  • For each of them, the probability that the matrix B makes it good is the quantity above.

So, we simply obtain

E_s = \binom{N}{s} \left(\frac{1}{2} \cdot \left(1 + (1 - 2p)^s \right) \right)^M

which is easily computed in \mathcal{O}(\log s + \log M) time using binary exponentiation.

Finding the sum of E_s for all s \leq k thus multiplies the above complexity by k, and the sum of k across all tests is bounded so it’s fast enough.
The \log s factor can be further removed since we only require it to find (1-2p)^s, but since s is increasing iteratively this can instead be computed in constant time by multiplying (1-2p) to (1-2p)^{s-1} each time.
The overall complexity is thus \mathcal{O}(k\log M).
Depending on implementation and how you compute modular inverses/binomial coefficients, the \log M factor might be heightened to \log{MOD} instead.


\mathcal{O}(k\log M) or \mathcal{O}(k\log{MOD}) per testcase.


Author's code (C++)
#include <bits/stdc++.h>
using namespace std;
#define int               long long
#define pb                push_back
#define ppb               pop_back
#define pf                push_front
#define ppf               pop_front
#define all(x)            (x).begin(),(x).end()
#define fr                first
#define sc                second
#define rep(i,a,b)        for(int i=a;i<b;i++)
#define ppc               __builtin_popcount
#define ppcll             __builtin_popcountll
#define debug(x)  cout<<(x)<<'\n';
#define vi                vector<int>

const long long INF=1e18;
const int32_t M=1e9+7;
const int32_t MM=998244353;

const int N=1e5+5;

int  fac[N + 1],p=M;
void fact()
    fac[0] = 1;
    for (int i = 1; i <= N; i++)
        fac[i] = (fac[i - 1] * i) % p;

int binpow(int a, int b,int m) {
    a %= m;
    int res = 1;
    while (b > 0) {
        if (b & 1)
            res = res * a % m;
        a = a * a % m;
        b >>= 1;
    return res;

int modInverse(int n, int p)
    return binpow(n, p - 2, p);

int calc(int n,int r)

if (n < r)
        return 0;
    // Base case
if (r == 0)
return 1;
return (fac[n] * modInverse(fac[r], p) % p
            * modInverse(fac[n - r], p) % p)% p;


void solve(){
int n,m,k,a,b;

int ans=0;
for(int s=1;s<=k;s++)
  int num=(b-2*a+M)%M;
  int sum1=(1+binpow(num,s,M)*binpow(b,(M-2)*s,M))%M;
  int sum=sum1*binpow(2,M-2,M)%M;




signed main(){

    int t=1;
    while(t--) solve();
    return 0;

Tester's code (C++)
 * 	^v^
#include <iostream>
#include <numeric>
#include <set>
#include <cctype>
#include <iomanip>
#include <chrono>
#include <queue>
#include <string>
#include <vector>
#include <functional>
#include <tuple>
#include <map>
#include <bitset>
#include <algorithm>
#include <array>
#include <random>
#include <cassert>

using namespace std;

using ll = long long int;
using ld = long double;

#define iamtefu ios_base::sync_with_stdio(false); cin.tie(0);

const ll mod = 1'000'000'007LL, N = 100'005;

mt19937 rng(chrono::high_resolution_clock::now().time_since_epoch().count());

void ext(ll a, ll b, ll &x, ll &y){
	if (b==0){
		x = 1, y = 0;
	} else {
		ext(b, a%b, x, y);
		ll tem = x;
		x = y;
		y = tem - (a/b)*y;

ll moi(ll a, ll m = mod){
	ll x = 0, y = 0;
	ext(a, m, x, y);
	return (x+m)%m;

ll pw(ll a, ll b, ll m){
	ll res=1;
	while (b){
		if (b&1){
	return res;

vector <ll> fac(N+1, 1ll), ifac(N+1, 1ll);

ll com(ll n, ll r){
	if (n<r){
		return 0ll;
	return ((fac[n]*ifac[r])%mod*ifac[n-r])%mod;

void scn(){
	// not necessarily distinct
	// right down ytdm
	ll n, m, k, a, b; cin>>n>>m>>k>>a>>b;
	ll even = pw(2, m, mod);
	ll p = (a*moi(b))%mod, q = (1-p+mod)%mod;
	vector <pair<ll,ll>> arr(n+1, {0,0});
	arr[1] = {q, p};
	for (int i=2; i<=n; i++){
		arr[i].first = ((arr[i-1].second*p)%mod+(arr[i-1].first*q)%mod)%mod;
		arr[i].second = ((arr[i-1].first*p)%mod+(arr[i-1].second*q)%mod)%mod;
	ll ans = 0;
	for (int i=1; i<=k; i++){
		(ans+=(com(n, i)*(pw(arr[i].first, m, mod)))%mod)%=mod;

int main(){
#if defined(airths)
	auto t1=chrono::high_resolution_clock::now();
	freopen("input.txt", "r", stdin);
	freopen("output.txt", "w", stdout);
	for (int i = 1; i<=N; i++){
		fac[i] = (fac[i-1]*i)%mod;
		ifac[i] = moi(fac[i], mod);
	int _; for(cin>>_; _--;)
#if defined(airths)
	auto t2=chrono::high_resolution_clock::now();
	ld ti=chrono::duration_cast<chrono::nanoseconds>(t2-t1).count();
	cerr<<"Time: "<<setprecision(12)<<ti;
	return 0;
Editorialist's code (Python)
mod = 10**9 + 7
for _ in range(int(input())):
    n, m, k, a, b = map(int, input().split())
    p = a * pow(b, mod-2, mod) % mod
    ans, comb = 0, 1
    for s in range(1, k+1):
        comb = (comb * (n+1-s) * pow(s, mod-2, mod)) % mod
        ans += comb * pow(pow(2, mod-2, mod) * (1 + pow(1 - 2*p, s, mod)), m, mod) % mod
    print(ans % mod)
1 Like