Investigating a conjecture from Stack Overflow 1/2

math
stack-overflow
exercise

A user on Stack Overflow proposes an interesting problem featuring modulos (and thus arithmetic) and sampling without replacement. I attempt to understand it better by identifying easy cases and performing some simulations.

Published

December 10, 2024

1 Problem statement

Initial question on SO.

There are \(n\) cards which have numbers \(1\)~\(n\) on each. You pick \(m\) cards from it, and you don’t put it back once you pick from it. Is the probability that their sum is divisible by \(k\) always \(\dfrac 1 k\), while \(k|n\)? If not, how do we generalize the probability?

The key features of this problem are:

  • picking without replacement \(X_i \in [1, n]\).
  • \(k\) possible values for the output of \(\sum X_i \bmod k\).

We would like to show that a single value has probability \(1/k\). This hints at the distribution being uniform over the \(k\) possibilities, but maybe something more complex is going on.

To investigate such questions, unless I happen to glance exactly the right analysis straight away, my process is to:

  • look for easy cases
  • perform simulations and try to identify patterns.

2 Simulation code

The situation is straightforward to implement in python using numpy.

import numpy as np
from scipy.stats import binomtest

repetitions = 10000

def single_simulation(n, m, k):
    """Compute the sum modulo k of m elements in [1, n]."""
    # Check input
    assert n % k == 0, f"Expected n % k == 0 but got {n % k} with {n=} and {k=}"
    assert m <= n, f"Excpected m <= n but got {m=} and {n=}"

1    picks = np.random.choice(1 + np.arange(n), size=(m,), replace=False)
    result = np.sum(picks) % k

    return result

def compute_success_probability(n, m, k, repetitions=10000, confidence_level=0.95):
    """Compute a Monte-Carlo approximation of the probability of success."""
    num_success = 0
    for idx in range(repetitions):
        num_success += single_simulation(n, m, k) == 0
    
    probability = num_success / repetitions

2    result = binomtest(k=num_success, n=repetitions).proportion_ci(confidence_level=confidence_level)
    low = result.low
    high = result.high

    return probability, low, high
1
Sample without replacement: cf numpy docs.
2
Exact confidence interval: cf scipy docs.

3 A simple counter-example

There are \(n \choose m\) possibilities for the selection of the cards. These possibilities are equally likely. Thus, the final probability must be some integer multiple of \(\left[ n \choose m \right]^{-1}\). Thus, if \(n \choose m\) is not divisible by \(k\), then the conjecture is false.

For example, for the triplet \(n=6, m=3, k=3\), we have \({n \choose m} = 20\) and the conjecture must be false.

This is confirmed by a simulation:

from scipy.special import comb

n = 6
m = 3
k = 3
probability, low, high = compute_success_probability(n=n, m=m, k=k)

print(f"Number of combinations: {comb(n, m) = }")
print(f"Number of combinations makes conjecture impossible: {comb(n, m) % k = }")
print(f"Estimated probability: {np.round(probability, 3)} (0.95 interval: {np.round(low,3), np.round(high,3)})")
print(f"Conjectured probability: {np.round(1/k, 3)}")
Number of combinations: comb(n, m) = np.float64(20.0)
Number of combinations makes conjecture impossible: comb(n, m) % k = np.float64(2.0)
Estimated probability: 0.404 (0.95 interval: (np.float64(0.395), np.float64(0.414)))
Conjectured probability: 0.333

4 The simple cases: \(m=1\), \(m=n\), \(m=n-1\)

Some cases are immediate.

If \(m=1\), then the result is a uniform variable with \(k\) possibilities. The conjecture is thus immediately true. [^If the sampling was with replacement instead of without, then the result would be immediately true due to this argument.]

If \(m=n\), there is no randomness: we just select all cards. The final result depends on the parity of \(k\).

  • If \(k\) is odd, then the terms cancel out in the sum by pairing them: \(1\) with \(k-1\), etc. The final sum is thus always 0.
  • If \(k\) is even, then we again pair the terms, but the term \(k/2\) is left standing alone. The final sum is thus:

\[ ((n / k) \bmod 2) * k / 2 \]

If \(m=n-1\), we can represent this situation as selecting a single card to exclude from the total sum. Thus, by the same argument as with \(m=1\), we know that the probability is \(1/k\).

5 An exhaustive exploration

We can also run a simple exhaustive experiment:

  • test all values of \(k\) in \([2,7]\),
  • test \(n\) in \(k * [1, 5]\),
  • test \(m\) in \([2, n-2]\) (excluding the easy cases discussed above),
  • if \({n \choose m} \bmod k = 0\) but the probability is different (which we measure using an exact 0.99 confidence interval) from \(1/k\), print out the details of the results.

Note that we expect there to be roughly 1 percent of false positives.

from IPython.display import Markdown
from tabulate import tabulate

table = []
for k in range(2, 7 + 1):
    for n in range(k, k * 5 + 1, k):
        for m in range(0, n+1):
            if m in [0, 1, n-1, n]:
                continue
            p, l, h = compute_success_probability(n=n,m=m,k=k,confidence_level=0.99)
            c = int(comb(n, m))
            conjecture = 1/k
            conjecture_validated = l <= conjecture <= h
            if c % k == 0 and not conjecture_validated:
                line = [
                    (n, m, k),
                    np.round(conjecture,4),
                    (np.round(l,4), np.round(h,4)),
                ]
                table.append(line)

Markdown(
    tabulate(
        table, 
        headers=[
            "(n, m, k)",
            "Conjectured proba $1/k$",
            "Probability 0.99 confidence interval"
        ],
    )
)
(n, m, k) Conjectured proba \(1/k\) Probability 0.99 confidence interval
(4, 2, 2) 0.5 (np.float64(0.3284), np.float64(0.3529))
(8, 2, 2) 0.5 (np.float64(0.4101), np.float64(0.4356))
(8, 4, 2) 0.5 (np.float64(0.5255), np.float64(0.5513))
(8, 6, 2) 0.5 (np.float64(0.4155), np.float64(0.4411))
(10, 4, 2) 0.5 (np.float64(0.5171), np.float64(0.5429))
(10, 6, 2) 0.5 (np.float64(0.4662), np.float64(0.492))
(9, 3, 3) 0.3333 (np.float64(0.3424), np.float64(0.3671))
(9, 6, 3) 0.3333 (np.float64(0.3465), np.float64(0.3713))
(8, 2, 4) 0.25 (np.float64(0.1993), np.float64(0.2204))
(8, 6, 4) 0.25 (np.float64(0.2047), np.float64(0.226))
(16, 2, 4) 0.25 (np.float64(0.2254), np.float64(0.2474))
(16, 14, 4) 0.25 (np.float64(0.2178), np.float64(0.2395))
(12, 2, 6) 0.1667 (np.float64(0.145), np.float64(0.1637))
(12, 10, 6) 0.1667 (np.float64(0.1363), np.float64(0.1546))
(18, 12, 6) 0.1667 (np.float64(0.1682), np.float64(0.188))
(24, 22, 6) 0.1667 (np.float64(0.1468), np.float64(0.1656))
(35, 3, 7) 0.1429 (np.float64(0.1434), np.float64(0.162))

We find many examples where the \({n \choose m} \bmod k = 0\) conjecture is not true. For some of these, it is straightforward to prove that the probability is not \(1/k\).

For example, consider the triplet \(n=4, m=2, k=2\).

  • The winning pairs are (2, 4) and (1, 3).
  • There are 6 total pairs.
  • The true probability is thus \(1/3\).

Other cases with \(m=2\) are similarly easy to analyze.

6 Mathew Spam improvement

In of the answers, Matthew Spam proposes the following observation:

If \(\operatorname{gcd}(m,n)=1\) then the probability is equal to \(\frac{1}{k}\).

and offers a proof: there is a symmetry that enables us to show that the distribution of \(\sum X_i \bmod k\) is uniform over \([1,k]\).

This is straightforward to extend to the case \(\operatorname{gcd}(m, k) = 1\).

In view of this, let’s rerun the analysis focusing only on cases where \(\operatorname{gcd}(m,k) \neq 1\).

table = []
for k in range(2, 7 + 1):
    for n in range(k, k * 5 + 1, k):
        for m in range(1, (n-1)+1):
            divisor = np.gcd(m, k)
            if divisor == 1:
                continue

            p, l, h = compute_success_probability(n=n,m=m,k=k,confidence_level=0.99)
            c = int(comb(n, m))
            conjecture = 1/k
            conjecture_validated = l <= conjecture <= h
            if not conjecture_validated:
                line = [
                    (n, m, k),
                    divisor,
                    np.round(conjecture,4),
                    (np.round(l,4), np.round(h,4)),
                ]
                table.append(line)

Markdown(
    tabulate(
        table, 
        headers=[
            "(n, m, k)",
            "gcd(m, k)",
            "Conjectured proba $1/k$",
            "Probability 0.99 confidence interval"
        ],
    )
)
(n, m, k) gcd(m, k) Conjectured proba \(1/k\) Probability 0.99 confidence interval
(4, 2, 2) 2 0.5 (np.float64(0.3153), np.float64(0.3396))
(6, 2, 2) 2 0.5 (np.float64(0.3885), np.float64(0.4138))
(6, 4, 2) 2 0.5 (np.float64(0.58), np.float64(0.6054))
(8, 2, 2) 2 0.5 (np.float64(0.4199), np.float64(0.4455))
(8, 4, 2) 2 0.5 (np.float64(0.5306), np.float64(0.5564))
(8, 6, 2) 2 0.5 (np.float64(0.4162), np.float64(0.4418))
(10, 2, 2) 2 0.5 (np.float64(0.4292), np.float64(0.4549))
(10, 4, 2) 2 0.5 (np.float64(0.5156), np.float64(0.5414))
(10, 6, 2) 2 0.5 (np.float64(0.4724), np.float64(0.4982))
(10, 8, 2) 2 0.5 (np.float64(0.5486), np.float64(0.5743))
(6, 3, 3) 3 0.3333 (np.float64(0.3908), np.float64(0.4161))
(9, 3, 3) 3 0.3333 (np.float64(0.346), np.float64(0.3708))
(9, 6, 3) 3 0.3333 (np.float64(0.3469), np.float64(0.3717))
(12, 9, 3) 3 0.3333 (np.float64(0.3346), np.float64(0.3592))
(15, 3, 3) 3 0.3333 (np.float64(0.3339), np.float64(0.3585))
(4, 2, 4) 2 0.25 (np.float64(0.1565), np.float64(0.1758))
(8, 2, 4) 2 0.25 (np.float64(0.2102), np.float64(0.2317))
(8, 6, 4) 2 0.25 (np.float64(0.2011), np.float64(0.2222))
(12, 2, 4) 2 0.25 (np.float64(0.2201), np.float64(0.2419))
(12, 8, 4) 4 0.25 (np.float64(0.2532), np.float64(0.276))
(12, 10, 4) 2 0.25 (np.float64(0.2113), np.float64(0.2328))
(16, 14, 4) 2 0.25 (np.float64(0.221), np.float64(0.2428))
(20, 2, 4) 2 0.25 (np.float64(0.2262), np.float64(0.2482))
(15, 10, 5) 5 0.2 (np.float64(0.2005), np.float64(0.2216))
(6, 2, 6) 2 0.1667 (np.float64(0.1217), np.float64(0.1391))
(6, 3, 6) 3 0.1667 (np.float64(0.1901), np.float64(0.2108))
(6, 4, 6) 2 0.1667 (np.float64(0.1821), np.float64(0.2025))
(12, 2, 6) 2 0.1667 (np.float64(0.142), np.float64(0.1605))
(12, 10, 6) 2 0.1667 (np.float64(0.1449), np.float64(0.1636))
(24, 2, 6) 2 0.1667 (np.float64(0.1466), np.float64(0.1654))
(24, 22, 6) 2 0.1667 (np.float64(0.1434), np.float64(0.162))
(30, 2, 6) 2 0.1667 (np.float64(0.1465), np.float64(0.1653))

In these results, we observe that, when \(\operatorname{gcd}(m, k) \neq 1\):

  • We often have that the probability does not match \(1/k\).
  • The pattern is less apparent when \(n\) grows, but it could be due to the missmatch being smaller.

7 Conclusions

  1. The conjecture is false, as is: we identified several counter-examples.

  2. If \({n \choose m} \bmod k != 0\) then the conjecture is clearly wrong.

  3. There are cases where \({n \choose m} \bmod k = 0\) but the conjecture is wrong. For example: \(n=4, m=2, k=2\) where the probability is provably \(1/3\) instead of \(1/2\).

  4. If \(\operatorname{gcd}(m, k)=1\) the conjecture is provably true.

  5. From a simulation study, it seems possible that when \(\operatorname{gcd}(m, k) \neq 1\), the probability is always different from \(1/k\).

These final observations has put me on the right track, and I am now able to move on to rigorous mathematical analysis. I will do so in a second post to keep things organized.

Back to top