pick K random integers without repetition

Question

Suppose we've to pick $K (\le N)$ random integers in the range $[0, N - 1]$ for very large $N$ such that there is no repetition, while also deterministically minimizing the number of calls made to the $rand()$ function. How would we do it?

Assume $rand(X)$ produces a random integer in the range $[0,X - 1]$; calling it with any argument counts as one call.

Picking at random till we get a unique value may not be efficient for smaller N, and $O(N)$ algorithms, i.e. ones which require looping through the entire range of integers is also not feasible, as $N$ could be up to $10^{18}$.

Is there any way to accomplish this?

Answer 1

Ṃųỻịgǻňạcểơửṩ mentioned that this is the reservoir problem. The reservoir sampling problem, though, is explicitly a very strict requirement for a streaming algorithm, i.e. it works in just one pass over a stream of values whose eventual size we don't know.

In this case, we want to do it in one pass, but we do know both N and K, so the problem is a little different and it permits a simpler solution.

Floyd's algorithm is designed for this task. It has some resemblance to both the Fisher-Yates shuffle and to reservoir sampling.

A couple of minor differences from the framing in your question:

• K < N
• we use randint(i, j) which produces a random integer between i and j inclusive
• the samples will be generated in the range 1 to N inclusive (but of course it is easy to adjust this at the end)

Here's the Python code:

import random
n = 200
k = 5
s = set()
for j in range(n-k+1, n+1):
    t = random.randint(1, j)
    if t not in s:
        s.add(t)
    else:
        s.add(j)
print(s)

A good source for the origin of and mathematical justification for Floyd's algorithm is this article by Jon Bentley from 1987. It's worth noting that the algorithm only calls randint() K times.

A very brief description of how it works is that, on iteration t, it draws a number randomly from [1..n-k+t]. If the number has not been drawn before, it is added to the sample. Otherwise, it adds the number n-k+t to the sample. This works to give the proper probabilities of inclusion in the final set.

In Bentley's article there is a description of the recursive formulation of the algorithm:

We can appreciate the correctness of Floyd’s algorithm anecdotally. When M = 5 and N = 10, the algorithm first recursively computes in S a 4-element random sample from 1..9. Next it assigns to T a random integer in the range 1..10. Of the 10 values that T can assume, exactly 5 result in inserting 10 into S: the four values already in S, and the value 10 itself. Thus element 10 is inserted into the set with the correct probability of 5/10.

I can't do better than that explanation, right now!

The algorithm has been mentioned on StackOverflow and proofs are here on Math SE.

Answer 2

What quality of randomness do you need? You can construct a Linear Congruential Generator as a PRNG with period of N (or the next prime above N) and discard out-of-range samples.

Use rand() only to seed your LCG PRNG.

Of course you get the same sequence every time, just from different start points, unless you also try to randomize your LCG's multiplier and adder parameters. Or as John Bollinger puts it, for any given K it can produce only about N distinct K-element samples.

Re-randomizing the multiplier (a) and adder (c) constants while still satisfying the conditions for period m (the modulus) can make multiples sets of sequences possible. And/or raise m so there are more numbers between N and m you'd have to discard if generated. But there's limited freedom so this only goes so far. (And depending on LCG parameters there can be significant correlation between numbers.)

This requires O(1) storage and negligible time per result, just a 64-bit multiply, add, and modulo.
If you need bigints for this, I guess O(log2(N)) storage and time per query with a small constant factor. But 10^18 < 2^64. Initial setup time requires finding a prime above N.

---

Quality of randomness is low: even the best LCGs are not great pseudo-random generators by modern standards, although some of their downsides come from choices like a power-of-2 modulus which makes the low-order bits highly correlated. (Other downsides include a short period, which we're taking advantage of here.) But automated selection of LCG parameters to fit a given period can result in generators much worse than a "good" LCG. For some small N, pairs of consecutive numbers are common in the output sequence. Maybe different LCG parameters with the same or similar period could avoid that.

I used this in practice for a partly-finished subtree-pruning-regrafting library that I never did much with, GPLed code on Github with comments explaining the algorithm, citing Knuth TAOCP for the conditions that produce period m, and Numerical Recipes for some PRNG quality folklore and best practices. The LCG-selection part does work correctly, but basically gives up and makes the multiplier 1 and adder 0 in some corner cases. Also for tiny maxval = N, 6 or less, since for this library's purposes, 6 was few enough to just brute-force try all possibilities and see which tree gave the maximum likelihood.

struct lcg {
    unsigned int state;
    unsigned int a, c, m;
    unsigned int startstate;
};
/******************** Linear Congruential Generator setup ***********/
/ to generate all possible SPRs in a pseudo-random order, we generate

all the numbers between 0 and the number of possible SPRs once each
without repetition using an LCG of the form: x_n+1 = x_n*a + c mod m.

/
/

Knuth: TAOCP 3.2.1: ex 2: if a and m are relatively prime,
the number X_0 will always appear in the period.  (will return to start?)

3.2.1.2: Theorem A: An LCG will have period m iff:

c is relatively prime to m



b = a-1 is a multiple of p, for every prime p dividing m



b is a multiple of 4, if m is a multiple of 4



3.2.1.3: ex 4:  m = 2**e >= 8  ->  maximum potency when a mod 8 = 5.
small multipliers are to be avoided.

Numerical recipies suggests c = a prime close to (1/2 - sqrt(3)/6)*m

*/
/* make up some parameters for an LCG that will have the maximum period

equal to the range, so every value is generated once.
When maxval doesn't have any repeated prime factors, a = m+1,
which is the same as a=1.  It's not exactly random, but it does still
mix up which SPRs are done.

TODO: take advantage of the fact that maxval = floor(sqrt(maxval))*ceil(sqrt(maxval))
could do that, but then the code would be less general-purpose
successfully brute-force tested for maxval=1..1000.

/
void findlcg(struct lcg lcg_params, int maxval)
{
    unsigned int a, b, c, m = maxval;
    int i;
primesetup (maxval+maxval/2);

if (m<=6){ // will be either 6 or 2.  Just loop in order
    b=0;
    c=1;
}else{
    int divlimit = m;
    b=1; // b must be a multiple of all of m's prime factors
    if (!(m%2)){
        b=2;
        while (divlimit%2 == 0) divlimit /= 2;
    }
    for (i=3 ; i <= divlimit ; i+=2){
        if (is_prime(i) && m%i == 0){
            b *= i;
            while (divlimit%i == 0) divlimit /= i;
        }
    }

    if (!(m%4)){    // if m is a mult of 4, b must be.
        while (b%4) b *= 2;
    }

    /* make sure a isn't too small */
    while (b<sqrtf(m)) b*=7;

    if (b == m) b=0;  // just give up and avoid overflow


// Numerical Recipies says there is "lore" behind this... :)
// TAOCP says it's useless unless the multiplier sucks (section 3.3.3, eq. 40)
// That would be us.
        c = next_prime(max(5, (0.5 - sqrtf(3)/6.0)m - 2));
        while (m%c == 0) c = next_prime(c+1);
        // Luckily we don't have to test for c>m, because it doesn't
        // happen with any m<100, and there are enough primes later...
/ I've observed that when a == m, (e.g. a=13, c=13, m=72) you often get

two consecutive numbers...  Do something to avoid that if it's a problem */
 }
a = b+1;
unsigned long long l = (unsigned long long)a * m;
 if (l > ULONG_MAX){
     fprintf(stderr, "spr: chosen Linear Congruential Generator is bogus\n"
     "   x_n+1 = x_n%u+%u mod %u\n"
     "   am > ULONG_MAX, so it would overflow :(\n", a, c, m);
 }
lcg_params->a = a;
 lcg_params->c = c;
 lcg_params->m = m;
 lcg_params->startstate = UINT_MAX;
 lcg_params->state = rand() % m;


}

My use-case was smallish trees so for prime finding I just used a straightforward Sieve to find true primes, not just relatively-prime which would be sufficient. Quality of pseudo-randomness was not a priority at the time, before moving on to other work. There might be room to spend more time choosing LCG parameters better than the algorithm shown here.

Answer 3

If you care only about having uniform probability for each of the $N$ elements to be included in the $K$-element sample, then @PeterCordes's suggestion to use just one call to $rand()$ to seed a suitably tuned LCG is very attractive.

If it is additionally important to have equal probability of each of the $\tbinom{N}{K}$ distinct possible $K$-element samples, without regard to order, then Floyd's algorithm for sampling without replacement, as suggested by @TheoH, is hard to beat.

But if the order of the selections is important, such that for every draw, each as-yet unselected element should have an equal probability of being selected, then you need something more. A partial Fischer-Yates shuffle could do this job in $K$ calls to $rand()$, but that requires an impractical amount of storage for large $N$. Or does it?

You can perform a Fischer-Yates-like selection of $K$ elements with $K$ calls to $rand()$ and with only $O(K)$ overhead by maintaining a map from previously selected values to replacement values. Such a map can be constructed with $O(K)$ memory for which storage and retrieval operations each have $O(1)$ amortized asymptotic complexity, so this does not asymptotically increase either computational complexity or the storage required ($O(K)$ storage already being required for the result). Having first initialized an empty map, $M$, one proceeds as follows on the $i$^th draw ($i$ starting from 0):

• compute $x = \mathrm{rand}(N - i)$
• compute the selected element $y$ as $y = \begin{cases} M[x] & \text{if it exists} \\ x & \text{otherwise} \end{cases}$
• set $M[x] = \begin{cases} M[N - i - 1] & \text{if it exists} \\ N - i - 1 & \text{otherwise} \end{cases}$

Answer 4

You can do it in one call to the rand() function, with a perfect uniform distribution, if you're willing to use large integers.

Here's the basic idea:

1. List out all of the possible ways to choose K distinct integers from the range [0, N-1].
2. Throw a dart at the list, choosing one of those ways

Once you've done this, you'll realize that you don't actually need that list of possibilities: You can use a mathematical formula to compute the selection choice directly from the id, in O(1) time, without the list.

For example, if you're trying to select k = 2 distinct items from n = 3 options, then you could use the list:

L(3, 2) = ["011", "101", "110"]

...where "110" means that you would select the first two items, and ignore the third. So, if rand() generates "id = 2", then you'd look up the binary sequence at position 2 ("110"), and select those first two items.

In general, there are M = (N K) = N!/(K!(N-K)!) possible ways to select K distinct integers from the set [0, N-1]. So a single call to the random number generator, rand(M), will randomly choose one of those possibilities.

All you need to get started is a canonical list of the possibilities--that is, a canonical ordering of the list.

Really, any ordering will work. But here's one that has a nice recursive formula:

L(n, 0) = ["0...0"] (where "0...0" has length n), and
L(n, n) = ["1...1"] (where "1...1" has length n), and
L(n, k) = [], for n < k, and
L(n, k) = ["0" + L(n-1, k), "1" + L(n-1, k-1)]
...where "1" + "001" = "1001" is string concatenation.

This lets you build up a canonical list of any size.

For example, if you'd like to generate the list from earlier, L(n, k) where n = 3 and k = 2, then you could do it like this:

L(1, 1) = ["1"]
L(2, 2) = ["11"]
L(1, 0) = ["0"]
L(2, 0) = ["00"]
L(3, 0) = ["000"]
L(2, 1)
 = ["0" + L(1, 1), "1" + L(1, 0)]
 = ["0" + ["1"], "1" + ["0"]]
 = ["01", "10"]
L(3, 1)
= ["0" + L(2, 1), "1" + L(2, 0)]
= ["0" + ["01", "10"], "1" + ["00"]]
= ["001", "010", "100"]
L(3, 2)
= ["0" + L(2,2), "1" + L(2,1)]
= ["0" + ["11"], "1" + ["01", "10"]]
= ["011", "101", "110"]

This allows you to compute the canonical list for any problem size. This already enough to get a solution using dynamic programming.

But you can do better. The recursive formula also yields a closed-form solution, and its inverse, g(id), directly computes the binary sequence.

This is a high-effort approach, but it does yield a fantastic O(1) solution that involves only one call to the random number generator.

Answer 5

Theo H's suggestion of using Floyd's algorithm is a good one, but requires knowing the number of elements we wish to sample in advance. However, there's a fairly similar algorithm that allows the sampled elements to be generated one at a time "on demand".

The basic idea is to run a partial Fisher–Yates-Durstenfeld–Knuth shuffle on a sparse array, storing only the elements that have changed from their initial values.

To understand this algorithm, it may be useful to start with a quick recap of the basic FYDK shuffle (zero-based ascending-index version with initially sorted array):

• Let a be an n element array initialized with a[k] := k for all 0 ≤ k < n.
• For each i from 0 to n-1:
  • Let j be a uniformly chosen random integer with i ≤ j < n.
  • If i ≠ j, swap a[i] and a[j].
• Now a is a random permutation of the integers from 0 to n-1.

Looking at the FYDK shuffle algorithm above, we can see a few notable features and invariants:

• In each loop iteration, the new value of a[i] is chosen uniformly at random from among the values not yet assigned to any a[k] for some k < i.
• The value assigned to a[i] in each loop iteration is never overwritten by later iterations.
• By induction on the above, at the end of each loop iteration, a[0] … a[i] is a uniformly chosen (and uniformly shuffled) random sample of i+1 distinct integers from the range [0, n).

In particular, this leads to the following observations:

• If we stop the FYDK shuffle (as given above) after k ≤ n iterations, the first k elements of the array a form a uniformly chosen (and uniformly shuffled) random sample of k distinct integers from the range [0, n).
• During the first k iterations, only at most 2k elements of a are changed from their initial values. Thus, if k is much smaller than n, we can save time and storage space by using a sparse map to only store the changed elements of the array a.

Also, since the first k-th elements of a are never changed after the first k iterations of the FYDK shuffle loop, we can easily transform the algorithm into a coroutine that yields the k-th element to the caller immediately after the k-th iteration:

• Let a be a sparse map with default values a[k] = k for all 0 ≤ k < n.
• For each i from 0 to n-1 (or until caller stops iteration):
  • Let j be a uniformly chosen random integer with i ≤ j < n.
  • Yield a[j].
  • Let a[j] := a[i].
  • Unset a[i].

…or, equivalently (using a map without implicit default values):

• Let a be a sparse map.
• For each i from 0 to n-1 (or until caller stops iteration):
  • Let j be a uniformly chosen random integer with i ≤ j < n.
  • If j is a key in a, yield a[j], else yield j.
  • If i is a key in a, let a[j] := a[i], else let a[j] := i.
  • Unset a[i].

(The "unset a[i]" step at the end of the loop is not strictly necessary for correctness, but since no later iteration will access a[i], there is no point in storing it any longer.)

Anyway, here's an example implementation in Python:

import random
def sample_range(start, stop=None, step=1):
    """
    Yields a stream of random unique integers from range(start, stop[, step]) using only
    O(k) storage for k outputs. This is done using the Fisher-Yates-Durstenfeld-Knuth
    shuffle on a sparse array (implemented using a dict) that only stores elements that
    have been changed and not yet yielded to the caller. The first k outputs form a
    uniformly chosen (and shuffled) random sample of k distinct elements from the range.
    """
    if stop is None:
        # Python convention: range(n) == range(0, n)
        start, stop = 0, start
# Conceptually this dictionary is initialized so that a[k] = k for all k in the range,
# but we only explicitly store values for entries that are changed from this default.
a = dict()

for i in range(start, stop, step):
    j = random.randrange(i, stop, step)
    # Conceptually we swap a[i] and a[j] and then yield a[i] to the caller. However,
    # as a[i] will never again be accessed after this, we can save some time and
    # storage space by leaving it unset.
    yield a.get(j, j)   # a.get(key, default) returns a[key] if set, else default
    a[j] = a.pop(i, i)  # a.pop(key, default) does the same, but also unsets a[key]

(And here's an online demo with a basic test of correctness and uniformity.)

FWIW, I do not know if this particular algorithm / coroutine has a name. I've used and written about it before, but for all I know I could be the first person to invent it. Given its simplicity, however, I very much doubt that.

In particular, I rather suspect that, had someone described this algorithm to Knuth 50 years ago, he would've considered it a natural and obvious extension of the FYDK shuffle. Which of course doesn't mean that it's necessarily quite as obvious to us lesser mortals.

(Edit: In fact, except for trivial indexing differences, I believe this is essentially the same algorithm as the one sketched out in John Bollinger's answer.)