This post is inspired by this tweet:

Eight years ago I found this way to fairly sample k distinct integers below n: no hash set, no extra memory, O(k log k). Still can't believe it's novel. The innovation is a step that deliberately favours duplicates...https://t.co/x6qmj1yrAI pic.twitter.com/54sCyUGcMn

— Paul Crowley (@ciphergoth) October 12, 2026

It is such an elegant little algorithm. It is somewhat similar to Floyd’s sampling algorithm, but without the need for a set data structure. In a way, this algorithm replaces the set data structure with sorting, which is nice because sorting can use less memory and is more cache-friendly.

The problem we are trying to solve is this: Given integers $n$ and $k$, where $0 \le k \le n$, uniformly sample $k$ distinct integers from the range $[0,n)$. The algorithm works in two phases. First, it uniformly samples a multiset of $k$ integers from the range $[0, n-k]$. Second, it transforms that multiset into $k$ distinct integers from the range $[0, n)$.

Let’s assume we already have a way to sample a multiset of size $k$ from the range $[0, n-k]$. To transform it into $k$ distinct integers in the range $[0, n)$, we can do this:

1def transform_multiset(d):
2    d.sort()
3    k = len(d)
4    for i in range(k):
5        d[i] += i
6    return d

Notice that the above code is similar to the final four lines of the X post. It is intuitive to see that every valid multiset of size $k$ in the range $[0, n-k]$ can be transformed into $k$ distinct integers from the range $[0, n)$, and every set of $k$ distinct integers from the range $[0, n)$ can be transformed back into a multiset of size $k$ in the range $[0, n-k]$, or in other words, transform_multiset is a bijective function. Since it is a bijective function, choosing the multiset uniformly is the same as choosing distinct integers uniformly.

Now, the question becomes: how do we uniformly sample the multiset? To make it more concrete, the problem is reduced to: Given integers $n$ and $k$, where $0 \le k \le n$, uniformly sample a multiset of size $k$ from the range $[0,n-k]$. This is where the Pólya urn model comes in.

To understand how it works, consider this algorithm to sample a multiset of $k$ integers in the range $[0,n)$:

1def sample_multiset(n, k):
2    d = [None] * k
3    for i in range(k):
4        d[i] = random.randrange(n)
5    return d

The algorithm above generates a valid multiset, but not uniformly. As an example, suppose $n = 3\, k = 3$, in which case the algorithm above may return these possible values:

Each of the elements above has a $\frac{1}{27}$ chance to be chosen. The problem is that some of the outcomes above represent the same multiset. For example, $[0, 0, 1]$, $[0, 1, 0]$, and $[1, 0, 0]$ refer to the same multiset $\{0, 0, 1\}$, which means $\{0, 0, 1\}$ has a $\frac{3}{27}$ chance of being chosen. In contrast, there is only one representation of the multiset $\{0, 0, 0\}$, which means it has only a $\frac{1}{27}$ chance of being chosen.

One observation that you might see is that a multiset with more repeated elements, such as $\{0, 0, 0\}$, is less likely to be chosen than a multiset with more distinct elements, such as $\{0, 1, 2\}$. The idea of the Pólya urn model is to compensate for this bias. When we choose $0$ as the first element, we want to make it more likely to choose $0$ again to compensate for the probability imbalance.

Now, imagine you have a bag of $3$ numbers, i.e., $\{0, 1, 2\}$. Next, you will pick one of the numbers from the bag randomly, record it, and put it back. Then, add one more copy of the number you picked to the bag. Now, you have $4$ numbers inside the bag, one for each number that you didn’t pick, and two for the number that you picked. Repeat this process $3$ times to get the multiset. This method, magically, gives us a uniform distribution over all possible multisets.

Let me demonstrate with a few examples. Suppose $n = 3$ and $k = 3$. The multiset $\{0, 0, 0\}$ has a $\frac{1}{10}$ chance of being chosen:

The multiset $\{1, 1, 2\}$ also has a $\frac{1}{10}$ chance of being chosen. It can be produced by three ordered sequences. For $[1,1,2]$, the process is:

This holds for other configurations as well. To implement this scheme, we can do this:

1def sample_multiset(n, k):
2  d = [None] * k
3  for i in range(k):
4      r = random.randrange(n + i)
5      if r < n:
6          d[i] = r
7      else:
8          d[i] = d[r - n]
9  return d

Notice that the above code is similar to the beginning of the X post.

It might not be obvious that the code above actually implements the Pólya urn model. The intuition is quite simple: $r$ is the random choice that represents our action of picking a number from the bag. Notice that it first chooses from $[0, n)$, then $[0, n+1)$, $[0, n+2)$, $...$, $[0, n+k-1)$. It keeps increasing one by one because we add one number to the bag in each step. The r < n check determines whether we choose one of the original numbers that were in the bag from the beginning. Otherwise, r represents one of the additional copies. r - n ranges from 0 to i - 1, which we can use to choose one of the previously added copies.

Combining phase one (sampling the multiset) and phase two (transforming the multiset into distinct integers), we get this:

1def sample(n, k):
2  d = sample_multiset(n - k + 1, k)
3  d.sort()
4  for i in range(k):
5      d[i] += i
6  return d

If we inline the sample_multiset, we get:

 1def sample(n, k):
 2  t = n - k + 1
 3  d = [None] * k
 4  for i in range(k):
 5      r = random.randrange(t + i)
 6      if r < t:
 7          d[i] = r
 8      else:
 9          d[i] = d[r - t]
10  d.sort()
11  for i in range(k):
12      d[i] += i
13  return d

This is the algorithm posted on X.