Algorithm design and analysis

"""How many comparisons does sorting need? Mergesort against the lower bound.

We count the comparisons made by mergesort and by insertion sort. For small n
we try every one of the n! orderings, so the worst case is computed exactly
and can be compared with the decision-tree lower bound ceil(log2 n!) and with
the recurrence C(n) = C(floor(n/2)) + C(ceil(n/2)) + n - 1 for mergesort. For
larger n we sort seeded random permutations and compare the average counts
with n log2 n and with n^2/4, the average for insertion sort.
"""

import itertools
import math
import random

def merge_sort(items):
    """Return (sorted list, number of comparisons) for mergesort."""
    if len(items) <= 1:
        return list(items), 0
    mid = len(items) // 2
    left, c1 = merge_sort(items[:mid])
    right, c2 = merge_sort(items[mid:])
    out, i, j, count = [], 0, 0, c1 + c2
    while i < len(left) and j < len(right):
        count += 1
        if left[i] <= right[j]:
            out.append(left[i])
            i += 1
        else:
            out.append(right[j])
            j += 1
    out.extend(left[i:])
    out.extend(right[j:])
    return out, count

def insertion_sort(items):
    """Return (sorted list, number of comparisons) for insertion sort."""
    out, count = list(items), 0
    for k in range(1, len(out)):
        x, j = out[k], k - 1
        while j >= 0:
            count += 1
            if out[j] <= x:
                break
            out[j + 1] = out[j]
            j -= 1
        out[j + 1] = x
    return out, count

def merge_recurrence(n):
    """Worst-case mergesort comparisons from C(n) = C(n//2) + C(n - n//2) + n - 1."""
    if n <= 1:
        return 0
    return merge_recurrence(n // 2) + merge_recurrence(n - n // 2) + n - 1

def lower_bound(n):
    """The decision-tree bound: any comparison sort needs ceil(log2 n!) in the worst case."""
    return math.ceil(math.log2(math.factorial(n))) if n > 1 else 0

def worst_case(sort, n):
    """Largest comparison count of `sort` over all n! orderings of range(n)."""
    return max(sort(list(p))[1] for p in itertools.permutations(range(n)))

def average_case(sort, n, trials, rng):
    """Mean comparison count over `trials` random permutations."""
    total = 0
    for _ in range(trials):
        items = list(range(n))
        rng.shuffle(items)
        result, count = sort(items)
        if result != sorted(items):
            raise AssertionError("sort failed")
        total += count
    return total / trials

def main():
    print("Worst cases over all n! inputs (exact):")
    print(" n   n!   lower bound  mergesort  recurrence  insertion  n(n-1)/2")
    for n in range(1, 8):
        print(f"{n:2d} {math.factorial(n):5d} {lower_bound(n):12d}"
              f" {worst_case(merge_sort, n):10d} {merge_recurrence(n):11d}"
              f" {worst_case(insertion_sort, n):10d} {n * (n - 1) // 2:9d}")
    print()
    rng = random.Random(2024)
    print("Averages over seeded random inputs (experimental):")
    print("    n   mergesort  n log2 n   insertion    n^2/4")
    for n, trials in [(16, 200), (128, 20), (1024, 2)]:
        m = average_case(merge_sort, n, trials, rng)
        s = average_case(insertion_sort, n, trials, rng)
        print(f"{n:5d} {m:11.1f} {n * math.log2(n):9.0f} {s:11.1f}"
              f" {n * n / 4:8.0f}")
    print()
    print("Epistemic status: the first table is exact for n <= 7 (every input")
    print("tried); mergesort's worst case matches its recurrence and stays")
    print("within one or two comparisons of the lower bound, which is a")
    print("theorem for every n. The averages are experiments on a few seeded")
    print("random inputs: evidence about typical behaviour, not a proof.")

if __name__ == "__main__":
    main()

"""Greedy algorithms against brute force: when they work and when they fail.

1. Kruskal's algorithm is compared with an exhaustive search over all
   spanning trees of seeded random weighted graphs on 5 vertices.
2. Dijkstra's algorithm is compared with Bellman-Ford on random directed
   acyclic graphs, first with non-negative weights and then allowing negative
   ones, where Dijkstra's greedy choice can be wrong.
3. Greedy change-making is compared with dynamic programming for every coin
   system {1, a, b} with b <= 12. By a theorem of Kozen and Zaks (1994), if
   greedy is ever suboptimal, it already fails for some amount below a + b,
   so checking those amounts decides each system completely.
"""

import heapq
import itertools
import random

def find(parent, x):
    """Union-find root of x, with path halving."""
    while parent[x] != x:
        parent[x] = parent[parent[x]]
        x = parent[x]
    return x

def kruskal(n, edges):
    """Total weight of a minimum spanning tree; edges are (w, u, v)."""
    parent, total, used = list(range(n)), 0, 0
    for w, u, v in sorted(edges):
        ru, rv = find(parent, u), find(parent, v)
        if ru != rv:
            parent[ru] = rv
            total, used = total + w, used + 1
    if used != n - 1:
        raise ValueError("graph is not connected")
    return total

def brute_force_mst(n, edges):
    """Minimum weight over all (n-1)-edge subsets that form a spanning tree."""
    best = None
    for subset in itertools.combinations(edges, n - 1):
        parent, ok = list(range(n)), True
        for _, u, v in subset:
            ru, rv = find(parent, u), find(parent, v)
            if ru == rv:
                ok = False
                break
            parent[ru] = rv
        if ok:
            weight = sum(w for w, _, _ in subset)
            best = weight if best is None else min(best, weight)
    return best

def dijkstra(n, adj, source=0):
    """Distances by Dijkstra's algorithm; a vertex is final once popped."""
    dist, done = [None] * n, [False] * n
    dist[source], heap = 0, [(0, source)]
    while heap:
        d, u = heapq.heappop(heap)
        if done[u]:
            continue
        done[u] = True
        for v, w in adj[u]:
            if not done[v] and (dist[v] is None or d + w < dist[v]):
                dist[v] = d + w
                heapq.heappush(heap, (dist[v], v))
    return dist

def bellman_ford(n, adj, source=0):
    """Distances by Bellman-Ford (correct with negative weights, no negative cycles)."""
    dist = [None] * n
    dist[source] = 0
    for _ in range(n - 1):
        for u in range(n):
            if dist[u] is None:
                continue
            for v, w in adj[u]:
                if dist[v] is None or dist[u] + w < dist[v]:
                    dist[v] = dist[u] + w
    return dist

def random_dag(n, rng, low, high):
    """Random DAG on 0..n-1 (edges i -> j for i < j), weights in [low, high]."""
    adj = [[] for _ in range(n)]
    for i, j in itertools.combinations(range(n), 2):
        if rng.random() < 0.5:
            adj[i].append((j, rng.randint(low, high)))
    return adj

def greedy_coins(coins, amount):
    """Number of coins used by the greedy method (coins include 1)."""
    count = 0
    for c in sorted(coins, reverse=True):
        count += amount // c
        amount %= c
    return count

def optimal_coins(coins, limit):
    """Fewest coins for every amount 0..limit, by dynamic programming."""
    best = [0] + [None] * limit
    for x in range(1, limit + 1):
        best[x] = 1 + min(best[x - c] for c in coins if c <= x)
    return best

def first_greedy_failure(coins):
    """Smallest amount below a + b (the two largest coins) where greedy is worse."""
    a, b = sorted(coins)[-2:]
    best = optimal_coins(coins, a + b)
    for x in range(a + b):
        if greedy_coins(coins, x) > best[x]:
            return x
    return None

def main():
    rng = random.Random(7)
    agree = 0
    for _ in range(300):
        edges = [(rng.randint(1, 9), u, v)
                 for u, v in itertools.combinations(range(5), 2)]
        agree += kruskal(5, edges) == brute_force_mst(5, edges)
    print(f"1. Kruskal = exhaustive minimum over spanning trees: {agree}/300")
    for low, label in [(0, "weights 0..9"), (-4, "weights -4..9")]:
        same = sum(dijkstra(6, g) == bellman_ford(6, g)
                   for g in (random_dag(6, rng, low, 9) for _ in range(300)))
        print(f"2. Dijkstra = Bellman-Ford, {label}: {same}/300 graphs")
    adj = [[(1, 1), (2, 2)], [], [(1, -2)]]
    print(f"   s->a (1), s->b (2), b->a (-2): Dijkstra gives d(a) ="
          f" {dijkstra(3, adj)[1]}, true d(a) = {bellman_ford(3, adj)[1]}")
    systems = [(1, a, b) for a in range(2, 12) for b in range(a + 1, 13)]
    fails = {s: first_greedy_failure(s) for s in systems}
    good = sum(1 for x in fails.values() if x is None)
    print(f"3. Coin systems {{1, a, b}}, b <= 12: greedy optimal for {good}"
          f" of {len(systems)}")
    print(f"   {{1, 3, 4}}: first failure at {fails[(1, 3, 4)]}"
          f" (greedy {greedy_coins((1, 3, 4), 6)} coins, best"
          f" {optimal_coins((1, 3, 4), 6)[6]})")
    print(f"   {{1, 5, 10}}: failures below 15: {fails[(1, 5, 10)]}")
    print()
    print("Epistemic status: parts 1 and 2 test random instances - evidence")
    print("that agrees with the correctness proofs of Kruskal (cut property)")
    print("and Dijkstra (non-negative weights), and counterexamples showing")
    print("the weight hypothesis is needed. Part 3 is an exhaustive check of")
    print("finitely many coin systems, conclusive only because of the")
    print("Kozen-Zaks bound.")

if __name__ == "__main__":
    main()

"""Longest common subsequence: the dynamic-programming table, and its savings.

The length L(i, j) of a longest common subsequence of the first i letters of
x and the first j letters of y satisfies L(i, j) = L(i-1, j-1) + 1 when the
i-th and j-th letters agree, and max(L(i-1, j), L(i, j-1)) otherwise. We fill
the (m+1) x (n+1) table for the classic example ABCBDAB / BDCABA, recover a
longest common subsequence by walking back through the table, and compare
the number of table cells with the number of calls made by the same
recursion without memory, on seeded random DNA strings.
"""

import random

def lcs_table(x, y):
    """The (len(x)+1) x (len(y)+1) table of LCS lengths of prefixes."""
    m, n = len(x), len(y)
    table = [[0] * (n + 1) for _ in range(m + 1)]
    for i in range(1, m + 1):
        for j in range(1, n + 1):
            if x[i - 1] == y[j - 1]:
                table[i][j] = table[i - 1][j - 1] + 1
            else:
                table[i][j] = max(table[i - 1][j], table[i][j - 1])
    return table

def traceback(x, y, table):
    """One longest common subsequence, read off the table from the end."""
    i, j, out = len(x), len(y), []
    while i > 0 and j > 0:
        if x[i - 1] == y[j - 1]:
            out.append(x[i - 1])
            i, j = i - 1, j - 1
        elif table[i - 1][j] >= table[i][j - 1]:
            i -= 1
        else:
            j -= 1
    return "".join(reversed(out))

def is_subsequence(s, t):
    """True if s can be obtained from t by deleting letters."""
    it = iter(t)
    return all(ch in it for ch in s)

def naive_calls(x, y):
    """(LCS length, number of calls) for the recursion without memory."""
    calls = 0

    def rec(i, j):
        nonlocal calls
        calls += 1
        if i == 0 or j == 0:
            return 0
        if x[i - 1] == y[j - 1]:
            return rec(i - 1, j - 1) + 1
        return max(rec(i - 1, j), rec(i, j - 1))

    return rec(len(x), len(y)), calls

def main():
    x, y = "ABCBDAB", "BDCABA"
    table = lcs_table(x, y)
    print(" " * 6 + "".join(f"{c:>3}" for c in y))
    for i, row in enumerate(table):
        label = " " if i == 0 else x[i - 1]
        print(f"  {label} " + " ".join(str(v).rjust(2) for v in row))
    z = traceback(x, y, table)
    print(f"LCS length {table[-1][-1]}; one LCS: {z}; subsequence of both:"
          f" {is_subsequence(z, x) and is_subsequence(z, y)}")
    print()
    rng = random.Random(42)
    print(" length  LCS  table cells  calls without memory")
    for n in (4, 6, 8, 10, 12):
        a = "".join(rng.choice("ACGT") for _ in range(n))
        b = "".join(rng.choice("ACGT") for _ in range(n))
        length, calls = naive_calls(a, b)
        assert length == lcs_table(a, b)[-1][-1]
        print(f"{n:7d} {length:4d} {(n + 1) ** 2:12d} {calls:21d}")
    print()
    print("Epistemic status: the table and traceback are exact computations")
    print("for one example; the lesson proves the recurrence correct for all")
    print("inputs (and Lean checks it). The call counts are measurements on a")
    print("few seeded random strings: they illustrate exponential growth of")
    print("the plain recursion but are not a worst-case analysis.")

if __name__ == "__main__":
    main()