Monte Carlo / Las Vegas Algorithms

Why Monte Carlo / Las Vegas Algorithms Matters

Imagine checking a large matrix product produced by another program. Recomputing every entry gives a direct answer, but it may cost almost as much as producing the product. Randomly testing a few carefully chosen projections can detect an incorrect answer much faster, at the price of occasionally missing an error. Now imagine finding a median: random pivots can make the search fast while preserving the exact answer for every pivot choice.

These are two different uses of randomness. Monte Carlo algorithms allow a controlled probability of an incorrect answer, usually within a specified work budget. Las Vegas algorithms return only correct answers, while their running time depends on the random choices and is analyzed in expectation. These names classify guarantees; they do not describe one common function or one particular implementation.

This lesson implements Freivalds' matrix-product test as the Monte Carlo example and randomized three-way Quickselect as the Las Vegas example. You will trace an actual false acceptance, follow an exact selection through several pivots, prove the relevant guarantees, and learn why repeating a seeded run is different from amplifying an error bound.

Core Concepts and Terminology

For an input $x$ and random choices $R$, write the output as $A(x;R)$ and the work as $T(x;R)$. A Monte Carlo guarantee might say $\Pr_R[A(x;R)\text{ is wrong}]\le\delta$ for every fixed valid $x$, together with a bounded amount of work. A Las Vegas guarantee says the returned answer is correct for every allowed choice sequence, with a finite expected running time over $R$ and termination with probability $1$ in the stated model.

The Monte Carlo task here is to decide whether $AB=C$ for square matrices over the field $\mathbb F_q$, with $q=1000003$. Every entry is represented by a residue in $[0,q)$. A round draws a vector $r\in\{0,1\}^n$ and compares $A(Br)$ with $Cr$, with all arithmetic modulo $q$. A mismatch proves inequality. A match can happen even when the product is wrong.

The Las Vegas task is to return the element at zero-based rank $k$ in the sorted order of an integer array, counting duplicates. Quickselect draws a pivot, partitions into smaller, equal, and larger values, and retains only the region containing rank $k$. It does not need to guess whether the selected answer is correct: the partition invariant establishes that fact.

“Las Vegas” does not mean that every execution must have an unbounded worst case. Quickselect's partition work has a quadratic worst-case bound even though its expected work is linear. Conversely, a Monte Carlo algorithm can stop early on a decisive witness; bounded work need not mean identical work on every run.

How Monte Carlo / Las Vegas Algorithms Works, Step by Step

For Freivalds' test:

1. Validate that all three matrices are equally sized, square, and contain valid field residues. Choose a positive round budget $t$.
2. Draw $n$ bits for the current vector $r$. The probability proof assumes these bits are independent and fair.
3. Compute $u=Br$, then $v=Au$, and separately $w=Cr$. Parenthesizing as matrix–vector products avoids forming the full matrix product.
4. If $v\ne w$, return REJECT and retain $r$ as a witness. No later round can make an unequal product equal.
5. If all $t$ rounds match, return ACCEPT. The algorithm has accepted the equality claim, but a false acceptance remains possible.

For Quickselect, work on a copy so the caller's input is unchanged. Maintain inclusive active bounds lo,hi and the absolute target index $k$. Choose a pivot from that interval and partition with indices lt,i,gt. During the scan, positions before lt are smaller than the pivot, positions from lt to i-1 are equal, positions after gt are larger, and positions from i to gt are still unknown.

If the current value is smaller, swap it into the left region and advance both lt and i. If it is larger, swap it with position gt, decrease gt, and examine the replacement at i next. If equal, advance i. Once the unknown region is empty, keep the left region when $k\lt\mathrm{lt}$, keep the right region when $k\gt\mathrm{gt}$, or return when $\mathrm{lt}\le k\le\mathrm{gt}$. The equality block is nonempty because the pivot was selected from the active interval.

flowchart TD
    A["Choose the problem and guarantee"] --> B{"Matrix verification or selection?"}
    B -->|Verification| C["Validate matrices; set round budget"]
    C --> D["Draw a vector; compute A(Br) and Cr"]
    D --> E{"Vectors differ?"}
    E -->|Yes| F["REJECT with an exact witness"]
    E -->|No| G{"Round budget exhausted?"}
    G -->|No| D
    G -->|Yes| H["ACCEPT; one-sided error remains"]
    B -->|Selection| I["Validate rank; copy input"]
    I --> J["Draw a pivot; partition three ways"]
    J --> K{"Target in equality block?"}
    K -->|Yes| L["Return exact selected value"]
    K -->|No| M["Keep the smaller active interval containing k"]
    M --> J

The programs deliberately keep these stopping rules separate. A matching projection is sufficient for one Monte Carlo round to pass. Quickselect stops only when an ordering invariant identifies the answer exactly. Random choices alone are not a certificate of either equality or balance.

Worked Example with a Full Trace

Use the following matrices over $\mathbb F_{1000003}$:

$$A=\begin{pmatrix}1&2\\3&4\end{pmatrix},\qquad B=\begin{pmatrix}2&0\\1&2\end{pmatrix},\qquad AB=\begin{pmatrix}4&4\\10&8\end{pmatrix}.$$

Give the verifier the incorrect claimed product

$$C=\begin{pmatrix}3&5\\10&8\end{pmatrix},\qquad D=AB-C=\begin{pmatrix}1&-1\\0&0\end{pmatrix}.$$

The displayed $-1$ denotes its field residue $q-1$. Here $Dr=0$ precisely when the two coordinates of the binary vector are equal. Enumerating all four possibilities gives:

Under uniform independent bits, the one-round false-acceptance probability is exactly $1/2$ for this input. The runnable programs use the recurrence $s_{j+1}=16807s_j\bmod2147483647$ and take each output's parity as a bit. Starting at seed $7$, the first outputs are $117649,1977326743,621132276,452154665$, yielding first $r=(1,1)^{\mathsf T}$ and then $r=(0,1)^{\mathsf T}$. A one-round call falsely accepts; a four-round-budget call rejects at round $2$, using four raw outputs. Each call starts a new generator with its specified seed.

For the correct product, a separate four-round call with seed $7$ tests vectors $(1,1)$, $(0,1)$, $(1,1)$, and $(0,0)$ and accepts after eight raw draws. These seeded observations demonstrate the mechanics; the ideal probability bound is not a statistical guarantee about this particular deterministic stream.

Now select zero-based rank $k=4$ from $[9,1,8,2,7,3,6,4,5]$, using seed $11$. Bounded sampling uses the raw value minus one and rejects an uneven suffix of the generator's output domain before reducing modulo the active length. The first three raw outputs are $184877,959744092,669281627$, with no rejected draws here.

Every tested less-than comparison counts once. When it fails, the subsequent greater-than comparison counts once more. Thus the final singleton partition contributes two comparisons. The copied work array is not fully sorted; finding one order statistic does not require sorting either discarded region.

Implementation in Python, C++, and Java

Each panel contains both complete examples, independent inputs, seeded generators, input validation, diagnostics, and the same sample output. Matrix dimensions are limited to $1\le n\le128$, the round budget to $1\le t\le1000$, and entries to residues in $[0,1000003)$. Quickselect accepts $1$ through $10000$ signed $32$-bit integer values and a valid zero-based rank. Both routines leave caller-owned arrays unchanged.

The matrix task is equality over the stated field. The code does not certify arbitrary integer-matrix equality or approximate floating-point equality. Reducing arbitrary integer inputs modulo one prime without a separate argument could hide a nonzero integer error that is a multiple of that prime.

For safe native arithmetic, the generator uses signed $64$-bit state updates. Each matrix–vector product also uses a $64$-bit multiplication and reduces the accumulator after every term. Its intermediate sum is at most $(q-1)^2+(q-1)\lt2^{63}$; an overflowing $32$-bit product would invalidate the computation before the modular reduction could help.

For a bound $m$, below sets $L=(P-1)-((P-1)\bmod m)$, draws $z=s-1$, rejects $z\ge L$, and returns $z\bmod m$. Uniform independent raw values would make this uniform on $[0,m)$. The recurrence is a reproducibility tool with correlated outputs, so neither independence nor the ideal error or expected-time bound follows merely from using a seed.

P = 2147483647
MOD = 1000003


class RNG:
    def __init__(self, seed):
        if not 1 <= seed < P:
            raise ValueError("invalid seed")
        self.state, self.draws = seed, 0

    def next(self):
        self.state = 16807 * self.state % P
        self.draws += 1
        return self.state

    def below(self, bound):
        if not 1 <= bound <= P - 1:
            raise ValueError("invalid random bound")
        limit = (P - 1) - (P - 1) % bound
        while True:
            value = self.next() - 1
            if value < limit:
                return value % bound


def matvec(matrix, vector):
    result = []
    for row in matrix:
        total = 0
        for a, b in zip(row, vector):
            total = (total + a * b) % MOD
        result.append(total)
    return result


def freivalds(a, b, c, rounds, seed):
    rng = RNG(seed)
    n = len(a)
    if not 1 <= n <= 128 or not 1 <= rounds <= 1000:
        raise ValueError("invalid size or round count")
    for matrix in (a, b, c):
        if len(matrix) != n or any(len(row) != n for row in matrix):
            raise ValueError("matrices must be square and equally sized")
        if any(not 0 <= x < MOD for row in matrix for x in row):
            raise ValueError("entries must be field residues")
    for trial in range(1, rounds + 1):
        r = [rng.next() % 2 for _ in range(n)]
        if matvec(a, matvec(b, r)) != matvec(c, r):
            return False, trial, rng.draws, r
    return True, rounds, rng.draws, []


def quickselect(values, k, seed):
    rng = RNG(seed)
    if not 1 <= len(values) <= 10000 or not 0 <= k < len(values):
        raise ValueError("invalid input length or rank")
    work = list(values)
    lo, hi, comparisons, trace = 0, len(work) - 1, 0, []
    while lo <= hi:
        pivot_index = lo + rng.below(hi - lo + 1)
        pivot = work[pivot_index]
        lt, i, gt = lo, lo, hi
        while i <= gt:
            comparisons += 1
            if work[i] < pivot:
                work[lt], work[i] = work[i], work[lt]
                lt += 1
                i += 1
            else:
                comparisons += 1
                if work[i] > pivot:
                    work[i], work[gt] = work[gt], work[i]
                    gt -= 1
                else:
                    i += 1
        trace.append((lo, hi, pivot_index, pivot, lt, gt))
        if k < lt:
            hi = lt - 1
        elif k > gt:
            lo = gt + 1
        else:
            return work[k], comparisons, rng.draws, trace, work
    raise AssertionError("unreachable active interval")


def fmt(values):
    return "[" + ",".join(map(str, values)) + "]"


def main():
    a, b = [[1, 2], [3, 4]], [[2, 0], [1, 2]]
    good, bad = [[4, 4], [10, 8]], [[3, 5], [10, 8]]
    results = []
    for name, c, rounds in [("good / 4 rounds", good, 4), ("bad / 1 round", bad, 1), ("bad / 4 rounds", bad, 4)]:
        result = freivalds(a, b, c, rounds, 7)
        results.append(result)
        passes, used, draws, witness = result
        print(f"{name} = {'ACCEPT' if passes else 'REJECT'}, rounds={used}, draws={draws}")
    print("rejecting witness = " + fmt(results[2][3]))
    values, k = [9, 1, 8, 2, 7, 3, 6, 4, 5], 4
    answer, comparisons, draws, trace, work = quickselect(values, k, 11)
    print(f"quickselect k={k}: value={answer}, partitions={len(trace)}, comparisons={comparisons}, draws={draws}")
    for index, (lo, hi, pivot_index, pivot, lt, gt) in enumerate(trace, 1):
        print(f"partition {index}: range=[{lo},{hi}], pivot index={pivot_index}, pivot={pivot}, equal=[{lt},{gt}]")
    print("work = " + fmt(work))
    unchanged = values == [9, 1, 8, 2, 7, 3, 6, 4, 5]
    print("input unchanged = " + str(int(unchanged)))
    ok = results[0][0] and results[1][0] and not results[2][0] and results[2][3] == [0, 1]
    ok = ok and answer == sorted(values)[k] and unchanged
    print("checks = " + ("PASS" if ok else "FAIL"))


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cstdint>
#include <iostream>
#include <stdexcept>
#include <string>
#include <vector>
const std::int64_t P = 2147483647;
const int MOD = 1000003;
using Matrix = std::vector<std::vector<int>>;
struct RNG {
    std::int64_t state, draws = 0;
    explicit RNG(std::int64_t seed) : state(seed) {
        if (seed < 1 || seed >= P) throw std::invalid_argument("invalid seed");
    }
    std::int64_t next() { state = 16807 * state % P; ++draws; return state; }
    int below(int bound) {
        if (bound < 1 || bound > P - 1) throw std::invalid_argument("invalid random bound");
        std::int64_t limit = (P - 1) - (P - 1) % bound;
        for (;;) {
            std::int64_t value = next() - 1;
            if (value < limit) return static_cast<int>(value % bound);
        }
    }
};
std::vector<int> matvec(const Matrix& matrix, const std::vector<int>& vector) {
    std::vector<int> result;
    for (const auto& row : matrix) {
        std::int64_t total = 0;
        for (std::size_t j = 0; j < row.size(); ++j) total = (total + static_cast<std::int64_t>(row[j]) * vector[j]) % MOD;
        result.push_back(static_cast<int>(total));
    }
    return result;
}
struct MCResult { bool passes; int rounds; std::int64_t draws; std::vector<int> witness; };
MCResult freivalds(const Matrix& a, const Matrix& b, const Matrix& c, int rounds, std::int64_t seed) {
    RNG rng(seed); std::size_t n = a.size();
    if (n < 1 || n > 128 || rounds < 1 || rounds > 1000) throw std::invalid_argument("invalid size or round count");
    for (const Matrix* matrix : {&a, &b, &c}) {
        if (matrix->size() != n) throw std::invalid_argument("matrix size mismatch");
        for (const auto& row : *matrix) {
            if (row.size() != n) throw std::invalid_argument("matrix must be square");
            for (int x : row) if (x < 0 || x >= MOD) throw std::invalid_argument("invalid field residue");
        }
    }
    for (int trial = 1; trial <= rounds; ++trial) {
        std::vector<int> r(n);
        for (int& x : r) x = static_cast<int>(rng.next() % 2);
        if (matvec(a, matvec(b, r)) != matvec(c, r)) return {false, trial, rng.draws, r};
    }
    return {true, rounds, rng.draws, {} };
}
struct Part { int lo, hi, pivot_index, pivot, lt, gt; };
struct QResult { int value; std::int64_t comparisons, draws; std::vector<Part> trace; std::vector<int> work; };
QResult quickselect(const std::vector<int>& values, int k, std::int64_t seed) {
    RNG rng(seed);
    if (values.empty() || values.size() > 10000 || k < 0 || static_cast<std::size_t>(k) >= values.size())
        throw std::invalid_argument("invalid input length or rank");
    std::vector<int> work = values; std::vector<Part> trace;
    int lo = 0, hi = static_cast<int>(work.size()) - 1; std::int64_t comparisons = 0;
    while (lo <= hi) {
        int pivot_index = lo + rng.below(hi - lo + 1), pivot = work[pivot_index];
        int lt = lo, i = lo, gt = hi;
        while (i <= gt) {
            ++comparisons;
            if (work[i] < pivot) { std::swap(work[lt], work[i]); ++lt; ++i; }
            else {
                ++comparisons;
                if (work[i] > pivot) { std::swap(work[i], work[gt]); --gt; }
                else ++i;
            }
        }
        trace.push_back({lo, hi, pivot_index, pivot, lt, gt});
        if (k < lt) hi = lt - 1;
        else if (k > gt) lo = gt + 1;
        else return {work[k], comparisons, rng.draws, trace, work};
    }
    throw std::logic_error("unreachable active interval");
}
std::string fmt(const std::vector<int>& values) {
    std::string result = "[";
    for (std::size_t i = 0; i < values.size(); ++i) { if (i) result += ","; result += std::to_string(values[i]); }
    return result + "]";
}
int main() {
    Matrix a = { {1, 2}, {3, 4} }, b = { {2, 0}, {1, 2} };
    Matrix good = { {4, 4}, {10, 8} }, bad = { {3, 5}, {10, 8} };
    std::vector<MCResult> results;
    std::vector<std::string> names = {"good / 4 rounds", "bad / 1 round", "bad / 4 rounds"};
    for (int index = 0; index < 3; ++index) {
        auto result = freivalds(a, b, index == 0 ? good : bad, index == 1 ? 1 : 4, 7);
        results.push_back(result);
        std::cout << names[index] << " = " << (result.passes ? "ACCEPT" : "REJECT") << ", rounds=" << result.rounds << ", draws=" << result.draws << '\n';
    }
    std::cout << "rejecting witness = " << fmt(results[2].witness) << '\n';
    std::vector<int> values = {9, 1, 8, 2, 7, 3, 6, 4, 5}; int k = 4;
    QResult result = quickselect(values, k, 11);
    std::cout << "quickselect k=" << k << ": value=" << result.value << ", partitions=" << result.trace.size() << ", comparisons=" << result.comparisons << ", draws=" << result.draws << '\n';
    for (std::size_t i = 0; i < result.trace.size(); ++i) {
        const auto& p = result.trace[i];
        std::cout << "partition " << i + 1 << ": range=[" << p.lo << ',' << p.hi << "], pivot index=" << p.pivot_index << ", pivot=" << p.pivot << ", equal=[" << p.lt << ',' << p.gt << "]\n";
    }
    std::cout << "work = " << fmt(result.work) << '\n';
    bool unchanged = values == std::vector<int>({9, 1, 8, 2, 7, 3, 6, 4, 5});
    std::cout << "input unchanged = " << unchanged << '\n';
    auto sorted = values; std::sort(sorted.begin(), sorted.end());
    bool ok = results[0].passes && results[1].passes && !results[2].passes && results[2].witness == std::vector<int>({0, 1});
    ok = ok && result.value == sorted[k] && unchanged;
    std::cout << "checks = " << (ok ? "PASS" : "FAIL") << '\n';
    return 0;
}

import java.util.ArrayList;
import java.util.Arrays;

public class Main {
    static final long P = 2147483647L;
    static final int MOD = 1000003;
    static class RNG {
        long state, draws;
        RNG(long seed) {
            if (seed < 1 || seed >= P) throw new IllegalArgumentException("invalid seed");
            state = seed;
        }
        long next() { state = 16807L * state % P; ++draws; return state; }
        int below(int bound) {
            if (bound < 1 || bound > P - 1) throw new IllegalArgumentException("invalid random bound");
            long limit = (P - 1) - (P - 1) % bound;
            while (true) {
                long value = next() - 1;
                if (value < limit) return (int)(value % bound);
            }
        }
    }
    static int[] matvec(int[][] matrix, int[] vector) {
        int[] result = new int[matrix.length];
        for (int i = 0; i < matrix.length; ++i) {
            long total = 0;
            for (int j = 0; j < matrix[i].length; ++j) total = (total + (long)matrix[i][j] * vector[j]) % MOD;
            result[i] = (int)total;
        }
        return result;
    }
    static class MCResult {
        final boolean passes; final int rounds; final long draws; final int[] witness;
        MCResult(boolean passes, int rounds, long draws, int[] witness) {
            this.passes = passes; this.rounds = rounds; this.draws = draws; this.witness = witness;
        }
    }
    static MCResult freivalds(int[][] a, int[][] b, int[][] c, int rounds, long seed) {
        RNG rng = new RNG(seed); int n = a.length;
        if (n < 1 || n > 128 || rounds < 1 || rounds > 1000) throw new IllegalArgumentException("invalid size or round count");
        for (int[][] matrix : new int[][][] {a, b, c}) {
            if (matrix.length != n) throw new IllegalArgumentException("matrix size mismatch");
            for (int[] row : matrix) {
                if (row.length != n) throw new IllegalArgumentException("matrix must be square");
                for (int x : row) if (x < 0 || x >= MOD) throw new IllegalArgumentException("invalid field residue");
            }
        }
        for (int trial = 1; trial <= rounds; ++trial) {
            int[] r = new int[n]; for (int i = 0; i < n; ++i) r[i] = (int)(rng.next() % 2);
            if (!Arrays.equals(matvec(a, matvec(b, r)), matvec(c, r))) return new MCResult(false, trial, rng.draws, r);
        }
        return new MCResult(true, rounds, rng.draws, new int[0]);
    }
    static class Part {
        final int lo, hi, pivotIndex, pivot, lt, gt;
        Part(int lo, int hi, int pivotIndex, int pivot, int lt, int gt) {
            this.lo = lo; this.hi = hi; this.pivotIndex = pivotIndex; this.pivot = pivot; this.lt = lt; this.gt = gt;
        }
    }
    static class QResult {
        final int value; final long comparisons, draws; final ArrayList<Part> trace; final int[] work;
        QResult(int value, long comparisons, long draws, ArrayList<Part> trace, int[] work) {
            this.value = value; this.comparisons = comparisons; this.draws = draws; this.trace = trace; this.work = work;
        }
    }
    static void swap(int[] a, int i, int j) { int tmp = a[i]; a[i] = a[j]; a[j] = tmp; }
    static QResult quickselect(int[] values, int k, long seed) {
        RNG rng = new RNG(seed);
        if (values.length < 1 || values.length > 10000 || k < 0 || k >= values.length)
            throw new IllegalArgumentException("invalid input length or rank");
        int[] work = values.clone(); int lo = 0, hi = work.length - 1; long comparisons = 0;
        ArrayList<Part> trace = new ArrayList<>();
        while (lo <= hi) {
            int pivotIndex = lo + rng.below(hi - lo + 1), pivot = work[pivotIndex];
            int lt = lo, i = lo, gt = hi;
            while (i <= gt) {
                ++comparisons;
                if (work[i] < pivot) { swap(work, lt, i); ++lt; ++i; }
                else {
                    ++comparisons;
                    if (work[i] > pivot) { swap(work, i, gt); --gt; }
                    else ++i;
                }
            }
            trace.add(new Part(lo, hi, pivotIndex, pivot, lt, gt));
            if (k < lt) hi = lt - 1;
            else if (k > gt) lo = gt + 1;
            else return new QResult(work[k], comparisons, rng.draws, trace, work);
        }
        throw new IllegalStateException("unreachable active interval");
    }
    static String fmt(int[] values) {
        StringBuilder b = new StringBuilder("[");
        for (int i = 0; i < values.length; ++i) { if (i > 0) b.append(','); b.append(values[i]); }
        return b.append(']').toString();
    }
    public static void main(String[] args) {
        int[][] a = { {1, 2}, {3, 4} }, b = { {2, 0}, {1, 2} };
        int[][] good = { {4, 4}, {10, 8} }, bad = { {3, 5}, {10, 8} };
        ArrayList<MCResult> results = new ArrayList<>();
        String[] names = {"good / 4 rounds", "bad / 1 round", "bad / 4 rounds"};
        for (int index = 0; index < 3; ++index) {
            MCResult result = freivalds(a, b, index == 0 ? good : bad, index == 1 ? 1 : 4, 7);
            results.add(result);
            System.out.println(names[index] + " = " + (result.passes ? "ACCEPT" : "REJECT") + ", rounds=" + result.rounds + ", draws=" + result.draws);
        }
        System.out.println("rejecting witness = " + fmt(results.get(2).witness));
        int[] values = {9, 1, 8, 2, 7, 3, 6, 4, 5}; int k = 4;
        QResult result = quickselect(values, k, 11);
        System.out.println("quickselect k=" + k + ": value=" + result.value + ", partitions=" + result.trace.size() + ", comparisons=" + result.comparisons + ", draws=" + result.draws);
        for (int i = 0; i < result.trace.size(); ++i) {
            Part p = result.trace.get(i);
            System.out.println("partition " + (i + 1) + ": range=[" + p.lo + "," + p.hi + "], pivot index=" + p.pivotIndex + ", pivot=" + p.pivot + ", equal=[" + p.lt + "," + p.gt + "]");
        }
        System.out.println("work = " + fmt(result.work));
        boolean unchanged = Arrays.equals(values, new int[] {9, 1, 8, 2, 7, 3, 6, 4, 5});
        System.out.println("input unchanged = " + (unchanged ? 1 : 0));
        int[] sorted = values.clone(); Arrays.sort(sorted);
        boolean ok = results.get(0).passes && results.get(1).passes && !results.get(2).passes && Arrays.equals(results.get(2).witness, new int[] {0, 1});
        ok = ok && result.value == sorted[k] && unchanged;
        System.out.println("checks = " + (ok ? "PASS" : "FAIL"));
    }
}

Run Python 3 directly, compile C++ with g++ -std=c++14 -O2 main.cpp -o randomized, or compile Java 17 or later with javac Main.java and run java Main. All three print these eleven lines:

• good / 4 rounds = ACCEPT, rounds=4, draws=8
• bad / 1 round = ACCEPT, rounds=1, draws=2
• bad / 4 rounds = REJECT, rounds=2, draws=4
• rejecting witness = [0,1]
• quickselect k=4: value=5, partitions=3, comparisons=26, draws=3
• partition 1: range=[0,8], pivot index=7, pivot=4, equal=[3,3]
• partition 2: range=[4,8], pivot index=5, pivot=6, equal=[5,5]
• partition 3: range=[4,4], pivot index=4, pivot=5, equal=[4,4]
• work = [1,3,2,4,5,6,8,9,7]
• input unchanged = 1
• checks = PASS

The final check intentionally expects the bad product's one-round acceptance: that fixture demonstrates an allowed Monte Carlo error. General testing must establish the one-sided property and its mathematical bound, rather than asserting that every bad product is rejected by every seed.

Complexity, Correctness, and Invariants

For Freivalds, let $D=AB-C$. If $D=0$, then $Dr=0$ for every vector, so a correct product is never rejected. If $D\ne0$, choose a nonzero entry $D_{ij}$ in one nonzero row. Condition on every coordinate of $r$ except $r_j$. The row equation has the form

$$D_{ij}r_j+c=0.$$

The two choices $r_j=0$ and $r_j=1$ cannot both satisfy this equation because their results differ by the nonzero field element $D_{ij}$. A fresh independent fair bit therefore makes the row vanish with probability at most $1/2$, and the whole vector can vanish only if that row does. Thus $\Pr[Dr=0]\le1/2$ for a fixed incorrect product.

With $t$ independent rounds, acceptance requires every projection to miss the error:

$$\Pr[\text{false acceptance}]\le2^{-t},\qquad t\ge\left\lceil\log_2\frac1\delta\right\rceil\ \Longrightarrow\ \Pr[\text{error}]\le\delta.$$

For example, $t=20$ gives an ideal upper bound $2^{-20}=1/1048576$. This bounds the chance that a fixed incorrect input passes the randomized experiment. It is not the posterior probability that an accepted input is wrong; that would also require a model for how inputs are chosen.

One round performs three dense matrix–vector products, each costing $\Theta(n^2)$ field operations. Validation also costs $\Theta(n^2)$. Correct products execute all rounds, taking $\Theta(tn^2)$ time. Incorrect products can reject in one round; under independent ideal rounds, the expected number executed is at most $\sum_{j=0}^{t-1}2^{-j}\lt2$, so their expected total work is $O(n^2)$ with early stopping. The work bound for any input is $O(tn^2)$, and extra vector storage is $O(n)$ beyond the three input matrices.

Quickselect preserves the input multiset and keeps the target's absolute sorted position inside the active interval. Every retained interval is strictly shorter because the nonempty pivot-equality block is either discarded with the other side or contains the returned answer. Its partition phase terminates for every sequence of valid pivot choices and returns an exact selected value.

With uniform independent pivots, a pivot in the middle half of the active sorted positions occurs with constant probability and shrinks any retained side to at most about three quarters of the previous size. The expected number of partitions before such a reduction is constant. Summing the geometrically decreasing scan lengths gives $O(N)$ expected selection work for input length $N$, including duplicates; the first scan already costs $\Omega(N)$. Repeated extreme pivots can instead produce quadratic work.

The field-operation analysis treats fixed-size modular arithmetic as constant cost. Arbitrarily large integer entries would require bit-complexity analysis instead. Quickselect is iterative and needs no recursive call stack; its copy and six-integer trace records use linear extra space, even on a quadratic-time pivot sequence. Under ideal uniform raw draws, the bounded sampler needs constant expected draws per pivot; the seeded implementation is analyzed separately through its actual draw count.

A useful conversion principle has limits. If a randomized candidate generator succeeds with probability at least $p\gt0$ per fresh trial and a complete exact verifier recognizes correct outputs, independently generate and verify candidates until one passes. The expected number of trials is at most $1/p$, and expected cost is $O((G+V)/p)$ when generation and verification each have per-trial cost bounded by $G$ and $V$. The verifier must recognize the required outputs; Freivalds' inequality witness does not certify its acceptance verdict.

In the other direction, stopping a Las Vegas computation after a fixed budget and returning an arbitrary answer introduces error on timeout. If $\mathbb E[T]\le\mu$, Markov's inequality gives $\Pr[T\gt c\mu]\le1/c$. Returning UNKNOWN instead preserves zero wrong answers but creates a partial-result interface, which is different from always producing an exact answer within that budget.

Variants, Pitfalls, and Optimizations

The most common amplification mistake is replaying the same random experiment. In the sample, restarting seed $7$ before every one-round call tests $(1,1)^{\mathsf T}$ every time, so every call misses the same error. Continuing to new choices exposes $(0,1)^{\mathsf T}$ next and detects it. Fresh choices are necessary for the independent-round proof; a deterministic pseudorandom stream still needs a separate randomness model before assigning that proof's numerical bound.

Other choices change the guarantee or its cost:

• One-sided versus two-sided repetition: Freivalds accepts only when every round passes because correct products always pass. For a two-sided test with success probability uniformly above $1/2$, independent majority voting can amplify success; applying the all-pass rule to such a test can increase the wrong kind of error.
• Many verification jobs: if every job has false-acceptance probability at most $2^{-t}$, a union bound gives at most $M2^{-t}$ for any false acceptance among $M$ jobs. Independence between jobs is not needed for that union bound, but each job must satisfy its own assumptions.
• Adversarial inputs: the bound is for a product fixed independently of the test vectors. Choosing a nonzero difference after observing the vectors can arrange for those projections to miss it. A replay seed is useful for debugging, while unpredictable production sampling is a separate requirement.
• Arithmetic domain: exact field equalities, integer equalities, and floating-point tolerances are different specifications. A tolerance-based comparison does not inherit the finite-field proof unchanged.
• Rejection sampling: direct remainder reduction of an uneven random range can bias pivots. Reject the unused suffix when uniform bounded samples are required by the model, and account for its work.
• Selection duplicates: retain an entire equality block rather than repeatedly peeling off one pivot-equal item. When every input value is equal, one three-way partition already identifies every valid rank.
• Budgeted selection and fallbacks: stopping early must produce UNKNOWN or use a separately correct fallback. A claimed worst-case bound for a hybrid requires bounding the work already spent as well as the fallback's work.
• Empirical tests: seeing no false acceptances in a small experiment does not prove a tiny error probability. Use exact small-domain enumeration and the algebraic argument to establish the bound, then test implementation behavior independently.

Practice Problems and Real-World Use

Use an exact product oracle on small matrices to distinguish equal and unequal instances. Equal products must always pass, while any rejection must carry a vector whose two computed projections differ. Enumerate all binary vectors for small nonzero difference matrices to verify that at most half conceal an error; enumerate short independent vector schedules to check amplification exactly.

For Quickselect, compare the returned value with an independently sorted copy for every rank. Check that the copied work array preserves all multiplicities, the caller's array is unchanged, the equality block contains its stated pivot, and every nonterminal step strictly shrinks the target's interval. Include empty or invalid inputs, singletons, all-equal arrays, negative values, integer boundaries, sorted and reversed arrays, and deliberately extreme pivot sequences.

For further study, see the Boston University notes on fingerprinting and Freivalds' technique and the Princeton selection API. For a deeper analysis that distinguishes key comparisons from bit comparisons, see Fill and Nakama's Quickselect paper. The examples here keep mathematical error, exact output correctness, and random work as separately testable claims.