Randomized Algorithms

import numpy as np

def freivalds_check(A, B, C, k=20):
    """
    Monte Carlo check: does A @ B == C?
    O(k * n^2) time. Error probability <= 2^(-k).

    If AB == C: always returns True.
    If AB != C: returns False with probability >= 1 - 2^(-k).

    Works over any field; here we use integers (exact arithmetic).
    """
    n = len(A)
    for _ in range(k):
        # Random binary vector r in {0,1}^n
        r = np.random.randint(0, 2, size=(n, 1))

        # Compute A(Br) and Cr — two matrix-vector products: O(n^2)
        Br  = B @ r         # n x 1
        ABr = A @ Br        # n x 1
        Cr  = C @ r         # n x 1

        if not np.array_equal(ABr, Cr):
            return False    # definite NO: AB != C
    return True             # probably YES (error <= 2^(-k))


def freivalds_mod2(A, B, C, k=30):
    """
    Freivalds over GF(2) (mod 2 arithmetic) — the original formulation.
    Even faster: bitwise operations.
    """
    n = len(A)
    A, B, C = (np.array(M) % 2 for M in (A, B, C))
    for _ in range(k):
        r   = np.random.randint(0, 2, size=(n, 1))
        Br  = B @ r % 2
        ABr = A @ Br % 2
        Cr  = C @ r % 2
        if not np.array_equal(ABr, Cr):
            return False
    return True


# Example: verify Strassen result
import numpy as np
n = 64
A = np.random.randint(0, 10, (n, n))
B = np.random.randint(0, 10, (n, n))
C_correct = A @ B
C_wrong   = C_correct.copy(); C_wrong[0, 0] += 1

print(freivalds_check(A, B, C_correct))  # True  (always)
print(freivalds_check(A, B, C_wrong))    # False (with prob >= 1 - 2^-20)

import random

# ── RANDOMIZED QUICKSORT ── O(n log n) expected, O(n^2) worst-case ──
def randomized_quicksort(A, lo=0, hi=None):
    """
    In-place randomized quicksort. Las Vegas: always correct.
    Expected O(n log n) for any input (random pivot defeats adversary).
    Worst case O(n^2) with probability 0 for random pivot.
    """
    if hi is None: hi = len(A) - 1
    if lo >= hi: return

    # Random pivot: defeats adversarial sorted/reverse-sorted input
    pivot_idx = random.randint(lo, hi)
    A[pivot_idx], A[hi] = A[hi], A[pivot_idx]   # move pivot to end

    # Partition (Lomuto scheme): O(n)
    pivot = A[hi]
    i = lo - 1
    for j in range(lo, hi):
        if A[j] <= pivot:
            i += 1
            A[i], A[j] = A[j], A[i]
    A[i+1], A[hi] = A[hi], A[i+1]
    q = i + 1

    randomized_quicksort(A, lo, q - 1)
    randomized_quicksort(A, q + 1, hi)


# ── PARANOID QUICKSORT ── guaranteed balanced partition ──────────────
def paranoid_quicksort(A, lo=0, hi=None):
    """
    Retry pivot until both sides <= 3/4 n.
    E[retries per level] <= 2. Total: O(n log n) expected.
    Recurrence: T(n) <= T(n/4) + T(3n/4) + 2cn.
    """
    if hi is None: hi = len(A) - 1
    if lo >= hi: return

    n = hi - lo + 1
    while True:
        pivot_idx = random.randint(lo, hi)
        A[pivot_idx], A[hi] = A[hi], A[pivot_idx]
        pivot = A[hi]
        i = lo - 1
        for j in range(lo, hi):
            if A[j] <= pivot:
                i += 1
                A[i], A[j] = A[j], A[i]
        A[i+1], A[hi] = A[hi], A[i+1]
        q = i + 1
        left_size  = q - lo
        right_size = hi - q
        if left_size <= 3*n//4 and right_size <= 3*n//4:
            break   # good pivot found

    paranoid_quicksort(A, lo, q - 1)
    paranoid_quicksort(A, q + 1, hi)


# ── RANDOMIZED SELECT ── O(n) expected ───────────────────────────────
def randomized_select(A, i, lo=0, hi=None):
    """
    Find i-th smallest (1-indexed) in A[lo..hi].
    Las Vegas: always correct. Expected O(n) time.
    """
    if hi is None: hi = len(A) - 1
    if lo == hi: return A[lo]

    # Random pivot, Lomuto partition
    pivot_idx = random.randint(lo, hi)
    A[pivot_idx], A[hi] = A[hi], A[pivot_idx]
    pivot = A[hi]
    j = lo - 1
    for k in range(lo, hi):
        if A[k] <= pivot:
            j += 1
            A[j], A[k] = A[k], A[j]
    A[j+1], A[hi] = A[hi], A[j+1]
    q = j + 1

    rank = q - lo + 1        # rank of pivot within A[lo..hi]
    if i == rank:
        return A[q]
    elif i < rank:
        return randomized_select(A, i, lo, q - 1)
    else:
        return randomized_select(A, i - rank, q + 1, hi)

import time, random, sys
sys.setrecursionlimit(100000)

def det_quicksort(A, lo=0, hi=None):
    """Deterministic quicksort (last element pivot) — O(n^2) on sorted input."""
    if hi is None: hi = len(A) - 1
    if lo >= hi: return
    pivot = A[hi]                        # FIXED pivot — adversary's target
    i = lo - 1
    for j in range(lo, hi):
        if A[j] <= pivot:
            i += 1
            A[i], A[j] = A[j], A[i]
    A[i+1], A[hi] = A[hi], A[i+1]
    q = i + 1
    det_quicksort(A, lo, q - 1)
    det_quicksort(A, q + 1, hi)

def benchmark(n=5000):
    sorted_input = list(range(n))        # adversarial for det_quicksort

    # Deterministic: O(n^2) on sorted input
    A = sorted_input.copy()
    t0 = time.perf_counter()
    try:
        det_quicksort(A)
        t_det = time.perf_counter() - t0
    except RecursionError:
        t_det = float('inf')

    # Randomized: O(n log n) expected on any input
    A = sorted_input.copy()
    t0 = time.perf_counter()
    randomized_quicksort(A)
    t_rand = time.perf_counter() - t0

    print(f"n={n} sorted input:")
    print(f"  Deterministic: {t_det*1000:.1f}ms  (O(n^2))")
    print(f"  Randomized:    {t_rand*1000:.1f}ms  (O(n log n) expected)")
    if t_det != float('inf'):
        print(f"  Speedup: {t_det/t_rand:.1f}x")

benchmark(3000)
# Typical output: Deterministic: ~80ms, Randomized: ~2ms, Speedup: ~40x