Randomised algorithms: expected time, sampling

A randomised algorithm makes some of its choices at random, so two runs on the same input can take different steps. There are two kinds. A Las Vegas algorithm always returns a correct answer, and only its running time varies: quicksort with random pivots sorts every list correctly and averages about 2n ln n comparisons on every input. A Monte Carlo algorithm runs for a fixed time and may be wrong, with a probability you choose: the Miller–Rabin primality test can mistake a composite number for a prime, but every extra round divides the bound on that chance by 4. Expected time is the average over the algorithm’s own random choices, taken for the worst input, not an average over inputs, so no input is reliably slow. Randomness also builds useful things directly: a fair shuffle and a fair sample of a stream. This lesson measures each of them and makes every run repeatable with a seed.

Las Vegas: always right, the time varies

Quicksort picks one item of a range as the pivot, moves the smaller items in front of it and the larger ones behind it, and sorts the two sides the same way. With the first item as pivot, a list that is already sorted is its worst case: every pivot is the smallest item left, and the work grows to n(n − 1)/2 comparisons. Choosing each pivot uniformly at random from its range removes that weakness, because a bad split now needs bad luck rather than a particular input. The program sorts the same sorted list of 1,000 numbers 500 times, each time with a different seed, and counts the comparisons:

# Randomised quicksort is a Las Vegas algorithm: every run sorts correctly, and only the work
# depends on the random pivots. The same sorted input is sorted 500 times, each run with its own seed.
import random
import statistics


def quicksort(a, rng):
    """Sorts a in place with a uniformly random pivot in every range; returns the comparisons made."""
    comparisons = 0
    ranges = [(0, len(a) - 1)]
    while ranges:
        lo, hi = ranges.pop()
        if lo >= hi:
            continue
        p = rng.randint(lo, hi)
        a[p], a[hi] = a[hi], a[p]
        pivot, i = a[hi], lo
        for j in range(lo, hi):
            comparisons += 1
            if a[j] < pivot:
                a[i], a[j] = a[j], a[i]
                i += 1
        a[i], a[hi] = a[hi], a[i]
        ranges.append((lo, i - 1))
        ranges.append((i + 1, hi))
    return comparisons


n, runs = 1000, 500
sorted_input = list(range(n))
counts = []
correct = 0
for seed in range(runs):
    a = sorted_input[:]
    counts.append(quicksort(a, random.Random(seed)))
    correct += a == sorted_input

expected = 2 * (n + 1) * sum(1 / k for k in range(1, n + 1)) - 4 * n
print(f"runs that sorted correctly: {correct} of {runs}")
print(f"comparisons per run: fewest {min(counts):,}, mean {statistics.mean(counts):,.0f}, most {max(counts):,}")
print(f"expected value 2(n + 1)H(n) - 4n: {expected:,.0f}")
for factor in (1.1, 1.25, 1.5):
    over = sum(c > factor * expected for c in counts)
    print(f"runs above {factor} x expected: {over}")
print(f"first item as pivot on this input: {n * (n - 1) // 2:,} comparisons, every run")

runs that sorted correctly: 500 of 500
comparisons per run: fewest 9,787, mean 10,984, most 13,390
expected value 2(n + 1)H(n) - 4n: 10,986
runs above 1.1 x expected: 36
runs above 1.25 x expected: 0
runs above 1.5 x expected: 0
first item as pivot on this input: 499,500 comparisons, every run

All 500 runs sorted correctly; that part never depends on luck. The work did: between 9,787 and 13,390 comparisons, with a mean of 10,984, against 499,500 on every run for the first-item pivot. Only 36 runs went above 1.1 times the expected value, and none above 1.25 times.

The expected value printed, 2(n + 1)H(n) − 4n with H(n) = 1 + 1/2 + … + 1/n, has a short derivation. Name the items by their places in sorted order. Items i and j, with i below j, are compared exactly when the first pivot drawn from the j − i + 1 items between them, both included, is i or j: a pivot strictly between them sends them to different sides before they ever meet. Each of those items is equally likely to be drawn first, so the chance is 2/(j − i + 1). Adding it over all pairs gives 2(n + 1)H(n) − 4n, which is about 2n ln n, or 1.39 n log₂ n. For n = 1,000 that is 10,986, two comparisons from the measured mean.

Monte Carlo: fixed time, a small chance of error

Trying divisors is far too slow to test a number of hundreds of digits for primality. Fermat’s little theorem gives a quick test: when n is prime, a raised to the power n − 1 leaves remainder 1 when divided by n, for every base a from 1 to n − 1. Most composite numbers fail that for most bases, but some, the Carmichael numbers, pass it for every base that shares no factor with them. The Miller–Rabin test checks a stronger condition, from Gary Miller’s work, which Michael Rabin turned into a randomised test, proving that an odd composite number passes it for at most a quarter of the bases. Each round draws a fresh random base, so k rounds let a composite through with probability at most (1/4)^k. The program counts the lying bases of six composite numbers, then tests every number up to 200,000 against a sieve:

# Miller-Rabin is a Monte Carlo test: a fixed number of rounds, and a small chance of a wrong
# "probably prime". First, for six composite numbers, how many bases would let each one pass;
# then every n up to 200,000, tested with 1, 2 and 5 random bases and checked against a sieve.
import random


def passes_fermat(a, n):
    """n looks prime to base a by Fermat's test: a^(n-1) leaves remainder 1."""
    return pow(a, n - 1, n) == 1


def passes_strong(a, n):
    """n looks prime to base a by the strong test that Miller-Rabin repeats (n odd, above 3)."""
    d, s = n - 1, 0
    while d % 2 == 0:
        d //= 2
        s += 1
    x = pow(a, d, n)
    if x == 1 or x == n - 1:
        return True
    for _ in range(s - 1):
        x = x * x % n
        if x == n - 1:
            return True
    return False


def probably_prime(n, rounds, rng):
    """False is always right (n is composite); True is wrong when every base drawn was a liar."""
    if n < 4:
        return n in (2, 3)
    if n % 2 == 0:
        return False
    return all(passes_strong(rng.randint(2, n - 2), n) for _ in range(rounds))


print("composite   bases   Fermat liars   strong liars")
for n in (221, 561, 1105, 1729, 2047, 8911):
    bases = range(2, n - 1)
    fermat = sum(passes_fermat(a, n) for a in bases)
    strong = sum(passes_strong(a, n) for a in bases)
    print(f"{n:>9,}   {len(bases):>5,}   {fermat:>12,}   {strong:>5,} = {strong / len(bases):.1%}")

limit = 200_000
sieve = bytearray([1]) * (limit + 1)
sieve[0] = sieve[1] = 0
for p in range(2, int(limit**0.5) + 1):
    if sieve[p]:
        sieve[p * p :: p] = bytearray(len(range(p * p, limit + 1, p)))
composites = limit - 1 - sum(sieve)
for rounds in (1, 2, 5):
    rng = random.Random(rounds)
    wrong = sum(probably_prime(n, rounds, rng) and not sieve[n] for n in range(2, limit + 1))
    print(f"{rounds} round(s): {wrong} of {composites:,} composites up to {limit:,} called prime")

composite   bases   Fermat liars   strong liars
      221     218             14       4 = 1.8%
      561     558            318       8 = 1.4%
    1,105   1,102            766      28 = 2.5%
    1,729   1,726          1,294     160 = 9.3%
    2,047   2,044            482     240 = 11.7%
    8,911   8,908          7,126   1,780 = 20.0%
1 round(s): 48 of 182,015 composites up to 200,000 called prime
2 round(s): 0 of 182,015 composites up to 200,000 called prime
5 round(s): 0 of 182,015 composites up to 200,000 called prime

The Carmichael number 561 = 3 × 11 × 17 fooled Fermat’s test for 318 of its 558 bases, but only 8 bases passed the strong test. 8,911 was the hardest of the six, with 1,780 lying bases, 20.0 per cent, close to the bound of a quarter. With one random base per number, 48 of the 182,015 composite numbers up to 200,000 were called prime; with two or five bases, none were. A round costs one modular exponentiation, O(log n) multiplications, and 40 rounds bring the error bound down to 4⁻⁴⁰, about 8 × 10⁻²⁵. A “composite” answer is never wrong, because a base that fails the test is proof.

Karger’s minimum-cut algorithm is another classic Monte Carlo algorithm. It merges the two ends of randomly chosen edges of a graph until only two groups of vertices are left and reports the edges between them, which form a minimum cut with a probability that repeating the run pushes as close to 1 as needed.

Shuffling fairly with Fisher–Yates

A shuffle is fair when each of the n! orders of the items is equally likely. Fisher–Yates, in the in-place form that Richard Durstenfeld published, walks from the last position down to the second. At position i it draws a position j from 0 to i, each equally likely, and swaps the two items. After that step the item at position i stays where it is, and the next step chooses among the items in front of it.

The choices multiply to n × (n − 1) × … × 2 = n! equally likely sequences, and each sequence produces a different order, so every order has probability exactly 1/n!. The shuffle draws n − 1 random numbers and makes n − 1 swaps: O(n) time and no extra memory.

A shortcut that looks just as good swaps every position with a position drawn from the whole list. It is biased, and the reason is arithmetic: it makes nⁿ equally likely sequences of choices, and for n of 3 or more that number is not a multiple of n!, so the orders cannot all get the same share. The program finds the exact chances for three items by trying every sequence of choices, then shuffles 600,000 times each way:

# Three ways to shuffle the items A, B, C: the shortcut "swap each position with any position", Fisher-Yates,
# and Fisher-Yates with an off-by-one bug. First the exact chances, from every equally likely sequence of random
# choices; then 600,000 real shuffles of the shortcut and of Fisher-Yates.
import itertools
import random
from collections import Counter
from fractions import Fraction


def naive_shuffle(a, choose):
    for i in range(len(a)):
        j = choose(len(a))  # any position from 0 to n - 1: the biased shortcut
        a[i], a[j] = a[j], a[i]


def fisher_yates(a, choose):
    for i in range(len(a) - 1, 0, -1):
        j = choose(i + 1)  # a position from 0 to i: Durstenfeld's version
        a[i], a[j] = a[j], a[i]


def off_by_one(a, choose):
    for i in range(len(a) - 1, 0, -1):
        j = choose(i)  # a position from 0 to i - 1: the item at i can never stay
        a[i], a[j] = a[j], a[i]


def exact(shuffle, sizes):
    """Runs the shuffle once for every sequence of choices; each sequence is equally likely."""
    orders = Counter()
    for picks in itertools.product(*(range(s) for s in sizes)):
        it = iter(picks)
        a = list("ABC")
        shuffle(a, lambda bound: next(it))
        orders["".join(a)] += 1
    total = sum(orders.values())
    return {order: Fraction(count, total) for order, count in orders.items()}, total


def measured(shuffle, runs, seed):
    rng = random.Random(seed)
    orders = Counter()
    for _ in range(runs):
        a = list("ABC")
        shuffle(a, rng.randrange)
        orders["".join(a)] += 1
    return orders


naive, naive_sequences = exact(naive_shuffle, [3, 3, 3])
fy, fy_sequences = exact(fisher_yates, [3, 2])
bug, bug_sequences = exact(off_by_one, [2, 1])
print(f"equally likely choice sequences: shortcut {naive_sequences}, Fisher-Yates {fy_sequences}, off-by-one {bug_sequences}")
print()
print("order   shortcut   Fisher-Yates   off-by-one")
for order in sorted(fy):
    print(f"{order}     {str(naive.get(order, 0)):>8}   {str(fy.get(order, 0)):>12}   {str(bug.get(order, 0)):>10}")

runs = 600_000
naive_runs = measured(naive_shuffle, runs, seed=5)
fy_runs = measured(fisher_yates, runs, seed=5)
print()
print(f"in {runs:,} shuffles (a fair shuffle gives each order about {runs // 6:,})")
print("order   shortcut   Fisher-Yates")
for order in sorted(fy):
    print(f"{order}     {naive_runs[order]:>8,}   {fy_runs[order]:>12,}")

equally likely choice sequences: shortcut 27, Fisher-Yates 6, off-by-one 2

order   shortcut   Fisher-Yates   off-by-one
ABC         4/27            1/6            0
ACB         5/27            1/6            0
BAC         5/27            1/6            0
BCA         5/27            1/6          1/2
CAB         4/27            1/6          1/2
CBA         4/27            1/6            0

in 600,000 shuffles (a fair shuffle gives each order about 100,000)
order   shortcut   Fisher-Yates
ABC       88,691        100,074
ACB      111,408        100,018
BAC      111,208        100,206
BCA      111,410         99,857
CAB       88,836        100,102
CBA       88,447         99,743

The shortcut’s 27 sequences fall 4 or 5 to an order, so three orders came up more than 111,000 times and three fewer than 89,000 times, where a fair shuffle gives 100,000. Fisher–Yates stayed within 260 of 100,000 for every order. The last column of exact chances is a common slip: drawing j from 0 to i − 1 instead of 0 to i, so the item at position i can never stay put. Only 2 of the 6 orders can then appear, each half the time. Sandra Sattolo showed that this version produces exactly the orders that form a single cycle, which is useful when a cycle is the goal and wrong for a shuffle. Sorting with a comparator that returns a random sign is no shuffle either: the ECMAScript standard leaves the order of a sort implementation-defined when the comparator is not consistent, so the result depends on the engine.

A shuffle can also produce no more orders than its generator has starting states:

# A generator can only produce as many different shuffles as it has starting states.
import math

deck = math.factorial(52)
print(f"orders of a 52-card deck: {deck:.3e}, about 2^{math.log2(deck):.1f}")
print(f"orders a 32-bit seed can reach: 2^32 = {2**32:,}")
print(f"share of deck orders reachable from a 32-bit seed: {2**32 / deck:.1e}")

period = 2**19937 - 1  # the period of Python's Mersenne Twister
print("2080! below the period:", math.factorial(2080) < period)
print("2081! below the period:", math.factorial(2081) < period)

orders of a 52-card deck: 8.066e+67, about 2^225.6
orders a 32-bit seed can reach: 2^32 = 4,294,967,296
share of deck orders reachable from a 32-bit seed: 5.3e-59
2080! below the period: True
2081! below the period: False

A deck of 52 cards has about 8.07 × 10⁶⁷ orders, about 2²²⁵·⁶. A generator started from a 32-bit seed reaches at most 2³² of them, a share of about 5 × 10⁻⁵⁹, however fair the shuffle. Python’s Mersenne Twister has a period of 2¹⁹⁹³⁷ − 1, and its documentation notes that 2,080 items is the longest list whose orders all fit in that period; the program confirms that 2,080! is below the period and 2,081! is not.

Sampling a stream with reservoir sampling

Choosing k winners uniformly from n entries is easy when all n sit in a list. It is harder when entries arrive one at a time and nobody knows n until the last one: the answers to a live quiz, the lines of a log too large to keep, the readings of a sensor. Reservoir sampling does it in one pass, with memory for only k items:

1. Keep the first k items.
2. When item number m arrives, for m greater than k, draw j uniformly from 0 to m − 1. If j is below k, the new item replaces the one in place j; otherwise the new item is dropped.

Each item ends up in the sample with probability exactly k/n. Item number m gets in with probability k/m. Each later item number t then replaces it with probability (k/t) × (1/k) = 1/t, so it survives all of them with probability m/(m + 1) × (m + 1)/(m + 2) × … × (n − 1)/n = m/n, and k/m × m/n = k/n. Jeffrey Vitter’s paper calls this method Algorithm R and adds Algorithm Z, which jumps over the items it would drop and runs in O(k(1 + log(n/k))) expected time. The program draws five winners from 48,213 entries made by a generator, so the stream is never in memory, and then checks that the draw is fair:

# Five prize winners from a stream of quiz entries whose length nobody knows in advance:
# reservoir sampling keeps exactly five entries in memory however long the stream runs.
import random
from collections import Counter


def reservoir_sample(stream, k, rng):
    """A uniform sample of k items of an iterable of unknown length, in one pass and O(k) memory."""
    sample = []
    for seen, item in enumerate(stream, start=1):
        if seen <= k:
            sample.append(item)
        else:
            j = rng.randrange(seen)  # 0 .. seen - 1, so the new item is kept with probability k / seen
            if j < k:
                sample[j] = item
    return sample


def quiz_entries(count):
    """Entries as they arrive, one at a time: a generator, so the stream is never held in memory."""
    for number in range(1, count + 1):
        yield f"entry-{number:05d}"


print("winners:", reservoir_sample(quiz_entries(48_213), 5, random.Random(3)))

# Is it fair? Draw 5 of 40 entries 40,000 times and count how often each entry wins.
rng = random.Random(8)
wins = Counter()
draws = 40_000
for _ in range(draws):
    wins.update(reservoir_sample(range(1, 41), 5, rng))
low, high = min(wins.values()), max(wins.values())
print(f"each entry should win {draws * 5 // 40:,} times; every entry won between {low:,} and {high:,} times")
first, last = sum(wins[e] for e in range(1, 6)), sum(wins[e] for e in range(36, 41))
print(f"entries 1-5, kept at the start, won {first:,} times; entries 36-40, seen last, won {last:,}")

winners: ['entry-07150', 'entry-45289', 'entry-22755', 'entry-31899', 'entry-39202']
each entry should win 5,000 times; every entry won between 4,827 and 5,120 times
entries 1-5, kept at the start, won 25,015 times; entries 36-40, seen last, won 24,964

Over 40,000 draws of 5 from 40 entries, each entry should win 5,000 times, and every one won between 4,827 and 5,120 times. The first five entries, which start in the reservoir, won 25,015 times between them and the last five 24,964, so arriving early or late gives no advantage. Python’s random.sample() needs the whole population as a sequence, so it cannot read a generator at all; reservoir sampling can.

Seeds: the same random numbers on every run

A pseudo-random generator is a deterministic program: from the same starting state, its seed, it produces the same numbers. That makes randomised code testable, and the lesson on testing algorithms shows how a seed lets a failing random test be repeated. Two details decide whether a seed really repeats a run:

• Python guarantees that random.Random(seed) gives the same random() values in future versions, but not that shuffle(), randint() or sample() keep their algorithms. Code that must reproduce exact results on every version takes its generator as a parameter.
• JavaScript’s Math.random() cannot be seeded at all: the engine chooses both the seed and the algorithm, and the ECMAScript standard asks only for approximately uniform values.

A test that must give the same numbers in both languages, on every version, can bring its own generator. Park and Miller’s minimal standard generator multiplies its state by a constant modulo 2³¹ − 1; with the multiplier 48,271, which Park, Miller and Stockmeyer recommended later, C++ ships it as std::minstd_rand, and the C++ standard requires the 10,000th value from the default seed 1 to be 399,268,537. Every product stays below 2⁵³, so JavaScript’s ordinary numbers compute it exactly:

# A seeded generator small enough to copy into any language: Park and Miller's "minimal standard"
# with the multiplier 48271, which C++ calls std::minstd_rand. Same seed, same numbers, everywhere.
M = 2_147_483_647  # 2**31 - 1, a prime


class MinStd:
    def __init__(self, seed):
        self.state = seed % M or 1  # the state must never be 0

    def next(self):
        """The next value, from 1 to M - 1."""
        self.state = self.state * 48271 % M
        return self.state

    def below(self, n):
        """A whole number from 0 to n - 1, every one equally likely (values in the last, short block are redrawn)."""
        limit = (M - 1) // n * n
        while True:
            x = self.next() - 1
            if x < limit:
                return x % n


g = MinStd(1)
values = [g.next() for _ in range(10_000)]
print("first three values from seed 1:", values[:3])
print("the 10,000th value:", values[-1], "(the C++ standard requires 399268537)")

g = MinStd(42)
deck = list(range(1, 11))
for i in range(len(deck) - 1, 0, -1):  # Fisher-Yates with the seeded generator
    j = g.below(i + 1)
    deck[i], deck[j] = deck[j], deck[i]
print("ten cards shuffled with seed 42:", deck)

first three values from seed 1: [48271, 182605794, 1291394886]
the 10,000th value: 399268537 (the C++ standard requires 399268537)
ten cards shuffled with seed 42: [4, 1, 7, 3, 9, 10, 8, 5, 6, 2]

// The same generator in JavaScript. Math.random() cannot be seeded, so tests bring their own.
// Every product below stays under 2**53, so plain numbers are exact and no BigInt is needed.
const M = 2147483647; // 2**31 - 1, a prime

class MinStd {
  constructor(seed) {
    this.state = seed % M || 1; // the state must never be 0
  }

  next() {
    this.state = (this.state * 48271) % M;
    return this.state;
  }

  below(n) {
    const limit = Math.floor((M - 1) / n) * n;
    for (;;) {
      const x = this.next() - 1;
      if (x < limit) return x % n;
    }
  }
}

let g = new MinStd(1);
const values = Array.from({ length: 10000 }, () => g.next());
console.log(`first three values from seed 1: [${values.slice(0, 3).join(', ')}]`);
console.log(`the 10,000th value: ${values[9999]} (the C++ standard requires 399268537)`);

g = new MinStd(42);
const deck = Array.from({ length: 10 }, (_, k) => k + 1);
for (let i = deck.length - 1; i > 0; i--) {
  const j = g.below(i + 1);
  [deck[i], deck[j]] = [deck[j], deck[i]];
}
console.log(`ten cards shuffled with seed 42: [${deck.join(', ')}]`);

first three values from seed 1: [48271, 182605794, 1291394886]
the 10,000th value: 399268537 (the C++ standard requires 399268537)
ten cards shuffled with seed 42: [4, 1, 7, 3, 9, 10, 8, 5, 6, 2]

Both versions printed 399,268,537 at the 10,000th step, as the standard requires, and both shuffled ten cards into 4, 1, 7, 3, 9, 10, 8, 5, 6, 2 with seed 42. below(n) draws again when a value lands in the last, incomplete block of n values, because x % n on its own favours small results. Here x can take 2,147,483,646 values, so for n = 10⁹ every result from 0 to 147,483,645 would come from three values of x and every other result from two: one and a half times as likely.

Randomness that must not be guessed

A generator that can be repeated can also be predicted. Python’s documentation calls the Mersenne Twister completely unsuitable for cryptographic purposes and points to the secrets module; MDN says the same of Math.random() and points to crypto.getRandomValues(). The program shows how little it takes: it reads 624 consecutive 32-bit outputs of a generator seeded from the operating system, rebuilds the generator’s whole internal state from them, and predicts what comes next:

# Why the random module must never make secrets: from 624 consecutive 32-bit outputs, the whole state
# of Python's Mersenne Twister can be rebuilt, and every later output predicted. The generator below is
# seeded from the operating system, so not even its seed is known.
import random
import secrets

MASK = 0xFFFFFFFF


def undo_right(y, shift):
    """Inverts y ^= y >> shift."""
    x = y
    for _ in range(32 // shift):
        x = y ^ (x >> shift)
    return x & MASK


def undo_left(y, shift, mask):
    """Inverts y ^= (y << shift) & mask."""
    x = y
    for _ in range(32 // shift):
        x = y ^ ((x << shift) & mask)
    return x & MASK


def untemper(y):
    """Turns one output back into the state word it came from (the tempering steps, in reverse)."""
    y = undo_right(y, 18)
    y = undo_left(y, 15, 0xEFC60000)
    y = undo_left(y, 7, 0x9D2C5680)
    return undo_right(y, 11)


service = random.Random()  # seeded from os.urandom
seen = [service.getrandbits(32) for _ in range(624)]  # values an observer collected

copy = random.Random()
copy.setstate((3, tuple(untemper(v) for v in seen) + (624,), None))
predicted = [copy.getrandbits(32) for _ in range(1000)]
actual = [service.getrandbits(32) for _ in range(1000)]
print("the next 1,000 outputs predicted exactly:", predicted == actual)
print("a 'random' 6-digit code predicted too:", copy.randrange(10**6) == service.randrange(10**6))

token = secrets.token_urlsafe(32)  # 32 bytes from the operating system's secure source
print("secrets.token_urlsafe(32) gives", len(token), "characters")

the next 1,000 outputs predicted exactly: True
a 'random' 6-digit code predicted too: True
secrets.token_urlsafe(32) gives 43 characters

The copy predicted the next 1,000 outputs exactly, and a six-digit code drawn after them as well. Nobody knew the seed; the outputs gave the state away, because every step that turns the twister’s state into an output can be undone. Login codes, password-reset links, session identifiers and API keys must therefore come from the operating system’s secure source: secrets.token_urlsafe(32) in Python, which turned 32 random bytes into 43 characters, or crypto.getRandomValues() in browsers and Node.js. Seeded generators stay the right tool for simulations, tests and the algorithms of this lesson, where repeatability helps and nobody gains by guessing.

Practise: a seeded shuffle and a reservoir

Write Fisher–Yates and reservoir sampling against the lesson’s seeded generator. The sample tests check exact results for fixed seeds, then count thousands of draws to check that every item is chosen about equally often.

M = 2_147_483_647  # 2**31 - 1


class MinStd:
    """The seeded generator of the lesson (leave it as it is): below(n) is a whole number from 0 to n - 1."""

    def __init__(self, seed):
        self.state = seed % M or 1

    def next(self):
        self.state = self.state * 48271 % M
        return self.state

    def below(self, n):
        limit = (M - 1) // n * n
        while True:
            x = self.next() - 1
            if x < limit:
                return x % n


def shuffle(items, rng):
    """Shuffles items in place with Fisher-Yates, using rng.below; returns nothing."""
    # Write the shuffle here.


def reservoir(stream, k, rng):
    """Returns k items of the iterable stream, chosen uniformly in one pass, using rng.below."""
    # Write reservoir sampling here.
    return []

from collections import Counter

from sampling import MinStd, reservoir, shuffle


def test_shuffle_exact():
    """with seed 42, ten cards come out as 4, 1, 7, 3, 9, 10, 8, 5, 6, 2"""
    cards = list(range(1, 11))
    assert shuffle(cards, MinStd(42)) is None
    assert cards == [4, 1, 7, 3, 9, 10, 8, 5, 6, 2]


def test_shuffle_in_place_and_complete():
    """shuffles the list it is given and keeps every item, for lists of 0 to 12 items"""
    for size in range(13):
        items = [f"item{i}" for i in range(size)]
        same = items
        shuffle(items, MinStd(size + 1))
        assert same is items
        assert sorted(items) == sorted(f"item{i}" for i in range(size)), items


def test_shuffle_fair():
    """6,000 shuffles of A, B, C: every order comes up between 930 and 1,070 times"""
    rng = MinStd(7)
    orders = Counter()
    for _ in range(6000):
        a = list("ABC")
        shuffle(a, rng)
        orders["".join(a)] += 1
    assert len(orders) == 6, orders
    assert all(930 <= count <= 1070 for count in orders.values()), orders


def test_reservoir_exact():
    """with seed 3, five of the numbers 1 to 100 are 73, 49, 86, 15, 52"""
    assert reservoir(range(1, 101), 5, MinStd(3)) == [73, 49, 86, 15, 52]


def test_reservoir_one_pass():
    """reads a generator once: four squares of 0 to 999 with seed 17"""
    squares = (x * x for x in range(1000))
    assert reservoir(squares, 4, MinStd(17)) == [6241, 197136, 962361, 841]


def test_reservoir_short_stream():
    """a stream shorter than k gives all its items in order; k = 0 gives nothing"""
    assert reservoir("abc", 5, MinStd(1)) == ["a", "b", "c"]
    assert reservoir("abc", 0, MinStd(1)) == []
    assert reservoir([], 3, MinStd(1)) == []


def test_reservoir_fair():
    """20,000 draws of 2 of 10 items: every item is drawn between 3,800 and 4,200 times"""
    rng = MinStd(9)
    drawn = Counter()
    for _ in range(20000):
        drawn.update(reservoir(range(10), 2, rng))
    assert sorted(drawn) == list(range(10)), drawn
    assert all(3800 <= count <= 4200 for count in drawn.values()), drawn

const M = 2147483647; // 2**31 - 1

/** The seeded generator of the lesson (leave it as it is): below(n) is a whole number from 0 to n - 1. */
export class MinStd {
  constructor(seed) {
    this.state = seed % M || 1;
  }

  next() {
    this.state = (this.state * 48271) % M;
    return this.state;
  }

  below(n) {
    const limit = Math.floor((M - 1) / n) * n;
    for (;;) {
      const x = this.next() - 1;
      if (x < limit) return x % n;
    }
  }
}

/** Shuffles items in place with Fisher-Yates, using rng.below; returns nothing. */
export function shuffle(items, rng) {
  // Write the shuffle here.
}

/** Returns k items of the iterable stream, chosen uniformly in one pass, using rng.below. */
export function reservoir(stream, k, rng) {
  // Write reservoir sampling here.
  return [];
}

import { test, assert } from 'mysmartcopilot:test';
import { MinStd, reservoir, shuffle } from './sampling.mjs';

test('with seed 42, ten cards come out as 4, 1, 7, 3, 9, 10, 8, 5, 6, 2', () => {
  const cards = Array.from({ length: 10 }, (_, i) => i + 1);
  assert.equal(shuffle(cards, new MinStd(42)), undefined);
  assert.deepEqual(cards, [4, 1, 7, 3, 9, 10, 8, 5, 6, 2]);
});

test('shuffles the array it is given and keeps every item, for arrays of 0 to 12 items', () => {
  for (let size = 0; size <= 12; size++) {
    const items = Array.from({ length: size }, (_, i) => `item${i}`);
    const same = items;
    shuffle(items, new MinStd(size + 1));
    assert.ok(same === items);
    assert.deepEqual(items.slice().sort(), Array.from({ length: size }, (_, i) => `item${i}`).sort());
  }
});

test('6,000 shuffles of A, B, C: every order comes up between 930 and 1,070 times', () => {
  const rng = new MinStd(7);
  const orders = new Map();
  for (let n = 0; n < 6000; n++) {
    const a = ['A', 'B', 'C'];
    shuffle(a, rng);
    const key = a.join('');
    orders.set(key, (orders.get(key) ?? 0) + 1);
  }
  assert.equal(orders.size, 6, JSON.stringify([...orders]));
  for (const [order, count] of orders) assert.ok(count >= 930 && count <= 1070, `${order}: ${count}`);
});

test('with seed 3, five of the numbers 1 to 100 are 73, 49, 86, 15, 52', () => {
  const numbers = Array.from({ length: 100 }, (_, i) => i + 1);
  assert.deepEqual(reservoir(numbers, 5, new MinStd(3)), [73, 49, 86, 15, 52]);
});

test('reads a generator once: four squares of 0 to 999 with seed 17', () => {
  function* squares() {
    for (let x = 0; x < 1000; x++) yield x * x;
  }
  assert.deepEqual(reservoir(squares(), 4, new MinStd(17)), [6241, 197136, 962361, 841]);
});

test('a stream shorter than k gives all its items in order; k = 0 gives nothing', () => {
  assert.deepEqual(reservoir('abc', 5, new MinStd(1)), ['a', 'b', 'c']);
  assert.deepEqual(reservoir('abc', 0, new MinStd(1)), []);
  assert.deepEqual(reservoir([], 3, new MinStd(1)), []);
});

test('20,000 draws of 2 of 10 items: every item is drawn between 3,800 and 4,200 times', () => {
  const rng = new MinStd(9);
  const drawn = new Array(10).fill(0);
  for (let n = 0; n < 20000; n++) for (const item of reservoir([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], 2, rng)) drawn[item]++;
  for (let item = 0; item < 10; item++) assert.ok(drawn[item] >= 3800 && drawn[item] <= 4200, `item ${item}: ${drawn[item]}`);
});

Interview questions

What is the difference between a Las Vegas and a Monte Carlo algorithm?

A Las Vegas algorithm always returns a correct answer, and only its running time is random; quicksort with random pivots is one, and its expected running time is what you analyse. A Monte Carlo algorithm has a running time you fix in advance, and its answer may be wrong with a bounded probability; the Miller–Rabin test is one, since it can call a composite number prime. The error of a Monte Carlo algorithm can be driven down by repeating it with fresh random choices: each Miller–Rabin round cuts the bound to a quarter. The two convert into each other. A Las Vegas algorithm stopped at a time limit becomes Monte Carlo, and a Monte Carlo algorithm whose answers can be checked becomes Las Vegas when it repeats until a check passes.

Why is the shuffle that swaps every position with a random position biased?

For n items it makes n random choices of n positions each, so there are nⁿ equally likely sequences of choices, and each sequence produces one order. A fair shuffle needs every one of the n! orders to receive the same number of sequences, which is possible only when nⁿ is a multiple of n!. For three items, 27 sequences cannot split evenly into 6 orders: three orders get 4 sequences and three get 5, so their chances are 4/27 and 5/27 instead of 1/6. Fisher–Yates avoids this by drawing from a shrinking range, 0 to i at step i, which gives exactly n! sequences, one for each order.

Key takeaways

• A Las Vegas algorithm is always right and has a random running time; a Monte Carlo algorithm has a fixed running time and a bounded chance of a wrong answer.
• Expected time averages over the algorithm’s coin flips for every input. Quicksort with random pivots made 9,787 to 13,390 comparisons on 1,000 sorted numbers, with a mean of 10,984 against the predicted 10,986.
• Each Miller–Rabin round lets a composite through with probability at most 1/4; one round called 48 of 182,015 composites prime, two rounds none.
• Fisher–Yates draws j from 0 to i and is exactly uniform. Drawing from the whole list is biased, and drawing from 0 to i − 1 yields only single cycles.
• Reservoir sampling keeps k items of a stream of any length, each with probability k/n, in one pass and O(k) memory.
• Seed a generator you pass in, or bring your own, to repeat randomised tests; never use a seeded generator for tokens, codes or keys.