Random Permutation, in depth: Lehmer-code ranking and unranking, keyed permutations you never store, exact bucket shuffles and the statistics that test them

What uniform means, and the four jobs

Start by counting. There are n! orderings, so choosing one uniformly needs about log2(n!) random bits. That is about n log2 n - 1.44 n bits. For a deck of 52 cards it is 225.6 bits. For a billion training examples it is roughly 2.8 times 10^10 bits. Two consequences follow. First, a generator seeded with 64 bits can reach at most 2^64 of those orderings. For anything larger than about 20 items, most permutations are unreachable. That is fine for training-data order and wrong for a card game that has to be fair. Second, you never need all those bits at once. Every method below consumes randomness in pieces.

There are four jobs, and each wants a different tool:

Ranking and unranking with the Lehmer code

The factorial number system gives every permutation of n items a unique integer in [0, n!). Write the permutation as a sequence of choices. The first element is chosen from n remaining items, the second from n-1, and so on. Record each choice as its index among the items still unused. That sequence of indices is the Lehmer code. Its digits have place values (n-1)!, (n-2)!, ..., 0!.

Worked example. Unrank r = 13 with n = 4. 13 divided by 3! = 6 gives 2 remainder 1, so take item 2 from [0,1,2,3] and leave [0,1,3]. 1 divided by 2! gives 0 remainder 1, so take item 0 and leave [1,3]. 1 divided by 1! gives 1 remainder 0, so take 3. The last item is 1. The permutation is [2,0,3,1] with Lehmer code (2,0,1,0). Running the steps backwards gives 2*6 + 0*2 + 1*1 = 13.

import math

def unrank(r, n):
    """Map an integer 0 <= r < n! to a permutation of range(n)."""
    items, out = list(range(n)), []
    for i in range(n, 0, -1):
        d, r = divmod(r, math.factorial(i - 1))
        out.append(items.pop(d))
    return out

def rank(perm):
    n, r, items = len(perm), 0, sorted(perm)
    for i, x in enumerate(perm):
        d = items.index(x)
        r += d * math.factorial(n - 1 - i)
        items.pop(d)
    return r

assert [rank(unrank(r, 5)) for r in range(120)] == list(range(120))
assert unrank(13, 4) == [2, 0, 3, 1]

To get a uniform permutation this way, draw r uniformly from [0, n!) with an exact bounded draw such as Python's secrets.randbelow(math.factorial(n)), then unrank it. Taking a 64-bit integer modulo n! is already biased for n of 3 or more, and for n of 21 or more it cannot even reach most permutations, because 21! is larger than 2^64. The list pop makes this O(n^2). A Fenwick tree over used and unused flags finds the d-th unused item with binary lifting in O(log n), so the whole conversion takes O(n log n).

Ranking is useful well beyond shuffling. One integer can serve as a compact key for a permutation in a cache or database. A test can name a failing order as "permutation 4,117" and rebuild it later. And for small n you can split the n! range into chunks and enumerate them across workers without coordination.

Keyed permutations you never store

A training data loader over 10^10 samples cannot afford an 80 GB int64 index array for each epoch. What it needs is a function perm(i) that is a bijection on range(n), costs O(1) memory, changes with the epoch and is cheap to compute. A block cipher is exactly a keyed bijection on fixed-width integers. The Feistel network is the easiest one to build. Split x into halves L and R and repeat L, R = R, L ^ F(key, round, R). Each round is invertible whatever F is, so the result is a bijection on [0, 2^b). F only has to mix well.

When n is not a power of two, use cycle-walking. Encrypt i. If the result is n or more, encrypt that result again, and repeat until it lands below n. Because the cipher is a permutation of the larger domain, following its cycle from a value below n must return below n. Restricted to range(n), the walk is therefore still a bijection.

import hashlib

class FeistelPermutation:
    """Keyed bijection on range(n), computed one index at a time."""
    def __init__(self, n, key: bytes, rounds=8):
        self.n, self.key, self.rounds = n, key, rounds
        bits = max(2, (n - 1).bit_length())
        bits += bits % 2                       # two equal halves
        self.half = bits // 2
        self.mask = (1 << self.half) - 1

    def _f(self, rnd, x):
        h = hashlib.blake2b(x.to_bytes(8, "little") + bytes([rnd]),
                            key=self.key, digest_size=8)
        return int.from_bytes(h.digest(), "little") & self.mask

    def _encrypt(self, x):
        left, right = x >> self.half, x & self.mask
        for rnd in range(self.rounds):
            left, right = right, left ^ self._f(rnd, right)
        return (left << self.half) | right

    def __call__(self, i):
        if not 0 <= i < self.n:
            raise IndexError(i)
        x = self._encrypt(i)
        while x >= self.n:                     # cycle-walk back into range
            x = self._encrypt(x)
        return x

for n in [1, 2, 3, 10, 1000, 4097]:
    P = FeistelPermutation(n, b"epoch-7")
    assert sorted(P(i) for i in range(n)) == list(range(n))

With n = 10 and key epoch-7, this gives [4, 3, 0, 1, 6, 9, 8, 7, 2, 5]. With epoch-8 it gives [8, 5, 2, 4, 9, 7, 6, 3, 0, 1]. To use it in a loader, derive the key as a hash of (run seed, epoch). Worker r of W reads positions r, r+W, r+2W and so on. To resume at step s, recompute the key and start from the right position; there is no state to save beyond the step counter. It gives the stateless, replayable order that the training data pipeline article requires, without storing an index per epoch.

Be clear about what you get. The family has at most as many members as there are keys, so it is pseudorandom: it is not a uniform draw over all n! orderings. For data order that does not matter. For anything adversarial, use a vetted format-preserving cipher. Rounding the bit width up to an even number can make the domain almost 4n (n = 4097 gets a domain of 16,384), so the expected walk can approach four encryptions per index. If that cost matters, use unbalanced halves. Cheaper bijections exist, such as an affine map modulo a prime, or the ZMap scanner's walk through a cyclic group modulo a prime just above 2^32. They are fast but visibly structured: consecutive outputs differ by a constant factor.

Shuffling data that does not fit

For data larger than one machine there are two standard designs.

Random keys plus sort. Give each item an independent random 64-bit key and sort by key. This is how a MapReduce or Spark job shuffles, because the framework already distributes sorting. The order is uniform only if the keys are distinct. By the birthday bound, a tie among n keys of 64 bits has probability about n^2 / 2^65, which is about 2.7 percent at n = 10^9. Breaking ties by input position then leaks the input order into tied pairs. Use 128-bit keys, or draw fresh keys for tied runs.

Rao-Sandelius buckets. Send each item to one of k buckets uniformly at random, shuffle each bucket independently (recursively, or with Fisher-Yates once it fits in memory) and concatenate. The result is exactly uniform. For a fixed final ordering, the bucket sizes s_1..s_k determine which bucket each item had to enter. Summing over size vectors gives k^-n times the sum of n!/(s_1!...s_k!) times 1/n!, and that equals 1/n!. Out of core, this becomes a two-pass shuffle: stream the input once and scatter it into k files on disk, then load each file, shuffle it in RAM and append it to the output.

import random

def rao_sandelius(items, k=2, rng=random):
    if len(items) <= 1:
        return list(items)
    buckets = [[] for _ in range(k)]
    for x in items:
        buckets[rng.randrange(k)].append(x)
    out = []
    for b in buckets:
        out.extend(rao_sandelius(b, k, rng))
    return out

Measured: in 240,000 shuffles of four items, all 24 orderings appeared, with counts between 9,770 and 10,281 against an expected 10,000. That spread is ordinary sampling noise.

Statistics of a random permutation, and a test harness

A uniform random permutation has well-known statistics. They are your test oracle. Use linearity of expectation with indicator variables to derive the mean of each one.

These values also produce algorithms. A uniform random derangement comes from shuffling and rejecting any result with a fixed point. Acceptance tends to 1/e, so you need about e = 2.72 shuffles on average, and the output is exactly uniform over derangements. Sattolo's algorithm is not a substitute: it produces only single n-cycles, a much smaller set.

import collections, math

def test_generator(shuffle, n=1000, trials=20000):
    fixed, cycles = collections.Counter(), 0
    for _ in range(trials):
        p = shuffle(list(range(n)))
        fixed[sum(p[i] == i for i in range(n))] += 1
        seen = [False] * n
        for i in range(n):
            if not seen[i]:
                cycles += 1
                while not seen[i]:
                    seen[i], i = True, p[i]
    h_n = sum(1 / k for k in range(1, n + 1))
    print("P(derangement)", fixed[0] / trials, "expect", 1 / math.e)
    print("mean cycles", cycles / trials, "expect", h_n)

Add a position test for small n: count how often each item lands in each slot and run a chi-square test on the n by n table. For n up to about 8, enumerate every ordering and compare its frequency with 1/n!. Passing these tests does not prove uniformity. Failing one proves bias.

Failure modes

• Modulo bias in bounded draws. rand64() % m over-weights small values when m does not divide 2^64. Use rejection or a library's exact bounded draw. With unranking the bias is huge, because m = n! is enormous.
• Sorting with a random comparator. The result depends on the sort algorithm and is biased. Generate random keys first, then sort by key.
• Key ties in distributed sorts. At billions of items, 64-bit keys collide. Measure the tie count in the job and fail if it is not zero.
• Same key every epoch. A Feistel loader that forgets to fold the epoch into its key repeats the same order every epoch. The loss curve looks fine and generalization suffers. Log a hash of the first 100 positions per epoch and assert that they differ.
• Too few Feistel rounds or a weak F. With two or three rounds and a linear F, outputs keep visible structure. Use at least six to eight rounds and a real hash as F, then run the statistics above on the output.
• Seed space smaller than the job. A 32-bit seed for a reproducible shuffle gives only four billion distinct orders. That is plenty for training. It is not enough when users can benefit from guessing the order.

Trade-offs

What to do next

1. Run unrank and rank above, then replace the list pops with a Fenwick tree and time n = 10^6.
2. Build the Feistel permutation for your dataset size. Assert that it is a bijection on a small n, and log the mean cycle-walk length for your real n.
3. Point test_generator at your production shuffle, including the distributed one, and keep it in CI.
4. Check whether your loader folds the epoch into its key, and read dataset ordering for how much mixing training actually needs.
5. If you shuffle with random keys at scale, add a tie counter to the job.