# LUCMAT - Editorial

Author: indreshsingh
Tester: airths
Editorialist: iceknight1093

TBD

# PREREQUISITES:

Probability and expected value, binomial coefficients

# PROBLEM:

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.

# EXPLANATION:

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.

How?

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.

# TIME COMPLEXITY:

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

# CODE:

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;
cin>>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;

ans=(ans+calc(n,s)*binpow(sum,m,M)%M)%M;

}

cout<<ans<<'\n';

}
signed main(){
ios_base::sync_with_stdio(false);
cin.tie(0);cout.tie(0);

int t=1;
fact();
cin>>t;
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;
a%=m;
while (b){
if (b&1){
res=(res*a)%m;
}
a=(a*a)%m;
b/=2;
}
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;
}
cout<<(ans)%mod<<'\n';
}

int main(){
iamtefu;
#if defined(airths)
auto t1=chrono::high_resolution_clock::now();
freopen("input.txt", "r", stdin);
freopen("output.txt", "w", stdout);
#endif
for (int i = 1; i<=N; i++){
fac[i] = (fac[i-1]*i)%mod;
ifac[i] = moi(fac[i], mod);
}
int _; for(cin>>_; _--;)
{
scn();
}
#if defined(airths)
auto t2=chrono::high_resolution_clock::now();
ld ti=chrono::duration_cast<chrono::nanoseconds>(t2-t1).count();
ti*=1e-6;
cerr<<"Time: "<<setprecision(12)<<ti;
cerr<<"ms\n";
#endif
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