مختبر الرياضيات

"""Cantor's diagonal trick, run on a finite list.

An infinite binary sequence is a rule giving a digit 0 or 1 at every position
n = 0, 1, 2, ... Here each sequence is a Python function of n. Given a list of
sequences s_0, s_1, s_2, ..., Cantor builds a new sequence d by flipping the
diagonal: d(n) = 1 - s_n(n). Then d differs from s_n at position n, so d is
not anywhere in the list.

We run the mechanism on a finite list of 8 sequences. The finite run only
illustrates the trick; the proof that no list of infinite sequences can be
complete is the general argument, which a computer cannot run to the end.

Runs unchanged in the browser (Pyodide) using only the standard library.
"""

from itertools import product
from math import isqrt
import textwrap

def _is_prime(n):
    """Trial-division primality test (n >= 0)."""
    if n < 2:
        return False
    return all(n % d for d in range(2, isqrt(n) + 1))

def _sqrt2_digit(n):
    """Digit n (counting from 0) after the binary point of sqrt(2) = 1.0110101..."""
    return isqrt(2 << (2 * (n + 1))) & 1

#: A fixed, deterministic list of infinite binary sequences: (name, rule).
SEQUENCES = [
    ("all zeros", lambda n: 0),
    ("all ones", lambda n: 1),
    ("alternating", lambda n: n % 2),
    ("1 at squares", lambda n: int(isqrt(n) ** 2 == n)),
    ("1 at primes", lambda n: int(_is_prime(n))),
    ("Thue-Morse", lambda n: bin(n).count("1") % 2),
    ("1 at multiples of 3", lambda n: int(n % 3 == 0)),
    ("binary digits of sqrt 2", _sqrt2_digit),
]

def prefix(seq, length):
    """The first ``length`` digits of a sequence, as a list."""
    return [seq(n) for n in range(length)]

def diagonal_flip(seqs):
    """Return Cantor's diagonal sequence d with d(n) = 1 - seqs[n](n).

    ``seqs`` is a list of sequences (functions of n). For a finite list, d(n) is
    only defined for n < len(seqs); in the infinite argument the list never
    runs out, so d is defined everywhere.
    """
    seqs = list(seqs)

    def d(n):
        if not 0 <= n < len(seqs):
            raise IndexError("the diagonal is only defined at positions 0..len(list)-1")
        return 1 - seqs[n](n)

    return d

def differs_from_every_row(seqs, d):
    """True when d(i) != seqs[i](i) for every i - so d equals none of the listed sequences."""
    return all(d(i) != s(i) for i, s in enumerate(seqs))

def first_difference(a, b, length):
    """The first position n < length where sequences a and b differ, or None."""
    for n in range(length):
        if a(n) != b(n):
            return n
    return None

def diagonal_flip_rows(rows):
    """Diagonal trick on finite rows of digits: return [1 - rows[i][i] for each i].

    Raises ValueError when some row is shorter than the number of rows - the
    trick needs one position per row.
    """
    if any(len(row) < len(rows) for row in rows):
        raise ValueError("every row needs at least as many digits as there are rows")
    return [1 - row[i] for i, row in enumerate(rows)]

def show_row(digits, mark=None):
    """Digits separated by spaces, with the digit at position ``mark`` in [brackets]."""
    return "".join(f"[{x}]" if i == mark else f" {x} " for i, x in enumerate(digits))

def main():
    seqs = [rule for _, rule in SEQUENCES]
    width = 12
    print(f"A list of {len(seqs)} infinite binary sequences (first {width} digits shown).")
    print("The digit in [brackets] is row n, position n - the diagonal.")
    print()
    for i, (name, rule) in enumerate(SEQUENCES):
        print(f"  s{i} {name:<24}{show_row(prefix(rule, width), mark=i)} ...")
    d = diagonal_flip(seqs)
    diag = [d(n) for n in range(len(seqs))]
    print(f"  d  {'flip each [digit]':<24}{show_row(diag)}   (positions 0..{len(seqs) - 1})")
    print()
    print(f"d differs from s_n at position n, for every n in the list: {differs_from_every_row(seqs, d)}")
    firsts = [first_difference(d, s, len(seqs)) for s in seqs]
    print(f"First position where d and s0, s1, ... differ: {firsts}")
    print("(d may differ earlier than the diagonal; the diagonal guarantees a difference.)")
    reordered = list(reversed(seqs))
    d2 = diagonal_flip(reordered)
    print(f"Reverse the order of the list: new d = "
          f"{''.join(str(d2(n)) for n in range(len(seqs)))}, "
          f"still in no row: {differs_from_every_row(reordered, d2)}")
    print()
    rows = [list(bits) for bits in product((0, 1), repeat=3)]
    flipped = diagonal_flip_rows(rows[:3])
    where = rows.index(flipped)
    print("Finite caution: all 2^3 = 8 binary strings of length 3 CAN be listed.")
    print(f"Flipping the diagonal of the first 3 rows gives {''.join(map(str, flipped))}"
          f" - row {where} of that list:")
    print("the trick only beats one row per position, so it needs infinite sequences.")
    print()
    print(textwrap.fill(
        "Epistemic status: these checks are exact, but they only illustrate the "
        "mechanism on 8 sequences. The proof that no infinite list s0, s1, s2, ... "
        "contains every infinite binary sequence is Cantor's general argument, "
        "which applies to every such list at once; no program can check them all.",
        width=78,
    ))

if __name__ == "__main__":
    main()

"""Clock arithmetic: tables mod k and which residues can be undone.

On a k-hour clock we only keep remainders 0, 1, ..., k-1. We print the
addition and multiplication tables mod 6, find which residues have a
multiplicative inverse (some b with a * b = 1 mod k) by searching, and compare
the answer with the rule "gcd(a, k) = 1" for every modulus k from 2 to 30.

Runs unchanged in the browser (Pyodide) using only the standard library.
"""

from math import gcd
import textwrap

def add_table(k):
    """The addition table mod k as a list of rows: row a, column b holds (a + b) % k."""
    if k < 1:
        raise ValueError("the modulus k must be at least 1")
    return [[(a + b) % k for b in range(k)] for a in range(k)]

def mul_table(k):
    """The multiplication table mod k: row a, column b holds (a * b) % k."""
    if k < 1:
        raise ValueError("the modulus k must be at least 1")
    return [[(a * b) % k for b in range(k)] for a in range(k)]

def inverse_by_search(a, k):
    """Return the b in 0..k-1 with a * b = 1 (mod k), found by trying every b; None if none.

    Note that mod 1 every number is 0 and 0 * 0 = 0 = 1 (mod 1), so 0 counts as invertible.
    """
    if k < 1:
        raise ValueError("the modulus k must be at least 1")
    for b in range(k):
        if (a * b) % k == 1 % k:
            return b
    return None

def invertible_residues(k):
    """The residues a in 0..k-1 that have an inverse mod k, found by search."""
    return [a for a in range(k) if inverse_by_search(a, k) is not None]

def coprime_residues(k):
    """The residues a in 0..k-1 with gcd(a, k) = 1."""
    return [a for a in range(k) if gcd(a, k) == 1]

def inverse_by_pow(a, k):
    """Inverse via Python's pow(a, -1, k) (the extended Euclidean algorithm), or None."""
    try:
        return pow(a, -1, k)
    except ValueError:
        return None

def compare_up_to(k_max, k_min=2):
    """For k = k_min..k_max compare search, gcd rule and pow(a, -1, k).

    Returns a list of moduli where the three methods disagree (empty if they always agree).
    """
    disagreements = []
    for k in range(k_min, k_max + 1):
        searched = invertible_residues(k)
        by_pow = [a for a in range(k) if inverse_by_pow(a, k) is not None]
        same_inverses = all(inverse_by_search(a, k) == inverse_by_pow(a, k) for a in searched)
        if not (searched == coprime_residues(k) == by_pow and same_inverses):
            disagreements.append(k)
    return disagreements

def format_side_by_side(k):
    """Return the lines of the + and x tables mod k printed next to each other."""
    head = " ".join(f"{b:>2}" for b in range(k))
    lines = [f" + |{head}      x |{head}", "---+" + "-" * (3 * k) + "    ---+" + "-" * (3 * k)]
    for a, (add_row, mul_row) in enumerate(zip(add_table(k), mul_table(k))):
        add_s = " ".join(f"{v:>2}" for v in add_row)
        mul_s = " ".join(f"{v:>2}" for v in mul_row)
        lines.append(f"{a:>2} |{add_s}     {a:>2} |{mul_s}")
    return lines

def main():
    print("Addition and multiplication on a 6-hour clock (mod 6):")
    for line in format_side_by_side(6):
        print("  " + line)
    print()
    for k in (6, 7):
        units = invertible_residues(k)
        pairs = " ".join(f"{a}*{inverse_by_search(a, k)}" for a in units if a <= inverse_by_search(a, k))
        print(f"Mod {k}: invertible residues {str(units):<20} inverse pairs: {pairs}")
    print("Mod 6, the residue 2 has no inverse: row 2 of the x table never contains a 1.")
    print()
    k_max = 30
    bad = compare_up_to(k_max)
    print(f"For every k = 2..{k_max}, compare three answers to 'which a have inverses?':")
    print("  search every b,  the rule gcd(a, k) = 1,  and pow(a, -1, k)")
    print(f"  moduli where they disagree: {bad if bad else 'none'}")
    counts = " ".join(str(len(coprime_residues(k))) for k in range(2, 17))
    print(f"  how many residues are invertible, k = 2..16: {counts}")
    print()
    print(textwrap.fill(
        f"Epistemic status: every residue for every k up to {k_max} was checked, so "
        "the gcd rule is verified for those moduli only. That it holds for every "
        "k is a theorem (it follows from Bezout's identity) proved elsewhere, not "
        "by this table.",
        width=78,
    ))

if __name__ == "__main__":
    main()

"""Sorting and counting a market basket.

A *collection* here is a basket of items. Two mangoes are two different items
that happen to share a name, so the basket is a list (order and repeats are
recorded), not a set. We sort the basket by a property, count each group, and
check that the counts do not depend on the order in which the items were
picked.

Runs unchanged in the browser (Pyodide) using only the standard library.
"""

from collections import Counter
from itertools import permutations
from math import factorial
import textwrap

#: The basket, in the order the items were picked. Each item is
#: (name, kind, colour). Repeated names are separate items.
BASKET = (
    ("mango", "fruit", "yellow"),
    ("okra", "vegetable", "green"),
    ("banana", "fruit", "yellow"),
    ("lime", "fruit", "green"),
    ("carrot", "vegetable", "orange"),
    ("mango", "fruit", "yellow"),
    ("aubergine", "vegetable", "purple"),
    ("orange", "fruit", "orange"),
    ("lime", "fruit", "green"),
    ("yam", "vegetable", "brown"),
)

#: Positions of each property inside an item tuple.
PROPERTIES = {"name": 0, "kind": 1, "colour": 2}

def value(item, prop):
    """Return the value of property ``prop`` ("name", "kind" or "colour") of an item."""
    try:
        return item[PROPERTIES[prop]]
    except KeyError:
        raise ValueError(f"unknown property {prop!r}; expected one of {sorted(PROPERTIES)}") from None

def sort_by(items, *props):
    """Return a new list of the items sorted by the given properties, in priority order.

    The original collection is not changed. With no properties, items are sorted
    by (name, kind, colour) as plain tuples.
    """
    if not props:
        return sorted(items)
    return sorted(items, key=lambda item: tuple(value(item, p) for p in props))

def group_by(items, prop):
    """Group items by a property: return {value: [items with that value]} with keys sorted."""
    groups = {}
    for item in items:
        groups.setdefault(value(item, prop), []).append(item)
    return {key: groups[key] for key in sorted(groups)}

def count_by(items, prop):
    """Count how many items have each value of ``prop``: {value: count}, keys sorted."""
    counts = Counter(value(item, prop) for item in items)
    return {key: counts[key] for key in sorted(counts)}

def distinct_values(items, prop):
    """Return the set of different values of ``prop`` that occur in the collection."""
    return {value(item, prop) for item in items}

def counts_same_in_every_order(items, prop):
    """Check every ordering of ``items``: do the group counts by ``prop`` ever change?

    Returns ``(orders_checked, all_same)``. The number of orders is n!, so this
    is only practical for small collections (n <= 8 or so).
    """
    items = list(items)
    reference = count_by(items, prop)
    checked = 0
    for order in permutations(items):
        checked += 1
        if count_by(order, prop) != reference:
            return checked, False
    return checked, True

def format_counts(counts):
    """Format {value: count} as 'a 2, b 3' with the total in brackets."""
    parts = ", ".join(f"{key} {n}" for key, n in counts.items())
    return f"{parts}   (total {sum(counts.values())})"

def main():
    basket = list(BASKET)
    print(f"A market basket of {len(basket)} items, in the order they were picked:")
    print("  " + ", ".join(value(item, "name") for item in basket))
    print()
    print("Sorted by kind, then colour, then name:")
    print(f"  {'kind':<10} {'colour':<8} name")
    for item in sort_by(basket, "kind", "colour", "name"):
        print(f"  {value(item, 'kind'):<10} {value(item, 'colour'):<8} {value(item, 'name')}")
    print()
    by_kind = count_by(basket, "kind")
    by_colour = count_by(basket, "colour")
    print("Counts by kind:   " + format_counts(by_kind))
    print("Counts by colour: " + format_counts(by_colour))
    sums = " + ".join(str(n) for n in by_kind.values())
    print(f"Every item lands in exactly one group, so the groups add up: {sums} = {len(basket)}.")
    names = distinct_values(basket, "name")
    print(f"Different names: {len(names)} (two mangoes are two items with one name).")
    print()
    small = basket[:6]
    n_orders = factorial(len(small))
    print(f"Does order matter? Take the first {len(small)} items and try all {len(small)}! = {n_orders} orders:")
    for prop in ("kind", "colour"):
        checked, same = counts_same_in_every_order(small, prop)
        print(f"  counts by {prop:<6} identical in all {checked} orders: {same}")
    print()
    print(textwrap.fill(
        "Epistemic status: the counts above are exact, and all "
        f"{n_orders} orders of the 6-item basket were checked, which proves "
        "order-independence for that basket only. Why it holds for every finite "
        "collection (counting pairs items with 1, 2, 3, ...) is explained in the "
        "lesson, not by this program.",
        width=78,
    ))

if __name__ == "__main__":
    main()

"""Composing finite functions: associativity, identities and order.

A function on a finite set X is stored as a dict {x: f(x)}. The composite
g o f ("g after f") sends x to g(f(x)). On the 3-element set X = {0, 1, 2} we
check associativity h o (g o f) = (h o g) o f for every triple of functions,
check the identity laws, and search for functions with g o f != f o g.

Runs unchanged in the browser (Pyodide) using only the standard library.
"""

from itertools import product
import textwrap

X = (0, 1, 2)

def all_functions(domain, codomain=None):
    """Return every function domain -> codomain (default: domain -> domain) as a dict."""
    domain = list(domain)
    codomain = domain if codomain is None else list(codomain)
    return [dict(zip(domain, values)) for values in product(codomain, repeat=len(domain))]

def identity(domain):
    """The identity function on ``domain``: every x goes to itself."""
    return {x: x for x in domain}

def compose(g, f):
    """Return g o f, the function x -> g(f(x)) ("first f, then g").

    Raises ValueError if some output of f is not an input of g, because then
    g o f is not defined.
    """
    missing = set(f.values()) - set(g)
    if missing:
        raise ValueError(f"g o f is undefined: g has no value at {sorted(missing)}")
    return {x: g[f[x]] for x in f}

def show(f):
    """Write a function on 0..n-1 as its list of values, e.g. [1, 0, 2] means 0->1, 1->0, 2->2."""
    return "[" + ", ".join(str(f[x]) for x in sorted(f)) + "]"

def check_associativity(functions):
    """Check h o (g o f) == (h o g) o f for all triples; return (triples, failures)."""
    triples = 0
    failures = []
    for f, g, h in product(functions, repeat=3):
        triples += 1
        if compose(h, compose(g, f)) != compose(compose(h, g), f):
            failures.append((f, g, h))
    return triples, failures

def check_identity_laws(functions, domain):
    """Check id o f == f == f o id for every f; return (functions_checked, failures)."""
    ident = identity(domain)
    failures = [f for f in functions if not (compose(ident, f) == f == compose(f, ident))]
    return len(functions), failures

def commuting_pairs(functions):
    """Return the list of ordered pairs (f, g) with g o f == f o g."""
    return [(f, g) for f, g in product(functions, repeat=2) if compose(g, f) == compose(f, g)]

def first_noncommuting_pair(functions):
    """Return the first ordered pair (f, g) with g o f != f o g, or None if all commute."""
    for f, g in product(functions, repeat=2):
        if compose(g, f) != compose(f, g):
            return f, g
    return None

def main():
    functions = all_functions(X)
    n = len(functions)
    print(f"X = {{0, 1, 2}} has 3^3 = {n} functions X -> X.")
    print("A function is written as its list of values: [1, 0, 2] means 0->1, 1->0, 2->2.")
    print("g o f means 'first f, then g': (g o f)(x) = g(f(x)).")
    print()
    triples, failures = check_associativity(functions)
    print(f"Associativity h o (g o f) = (h o g) o f: {triples} triples checked, {len(failures)} failures")
    checked, failures = check_identity_laws(functions, X)
    print(f"Identity laws id o f = f = f o id:     {checked} functions checked, {len(failures)} failures")
    print()
    f, g = first_noncommuting_pair(functions)
    print("But order matters. First pair found with g o f != f o g:")
    print(f"  f = {show(f)}   g = {show(g)}")
    print(f"  g o f = {show(compose(g, f))}   f o g = {show(compose(f, g))}")
    swap, cycle = {0: 1, 1: 0, 2: 2}, {0: 1, 1: 2, 2: 0}
    print("A friendlier pair: s = swap 0 and 1, r = rotate 0->1->2->0:")
    print(f"  s = {show(swap)}   r = {show(cycle)}")
    print(f"  r o s = {show(compose(cycle, swap))}   s o r = {show(compose(swap, cycle))}")
    pairs = commuting_pairs(functions)
    print(f"Only {len(pairs)} of the {n * n} ordered pairs commute.")
    print()
    print("Why associativity always holds: both sides send x to h(g(f(x))).")
    print()
    print(textwrap.fill(
        f"Epistemic status: all {triples} triples and all {n} functions on this "
        "3-element set were checked, a proof by exhaustion for this set only; the "
        "one-line argument above proves associativity for all functions. A single "
        "non-commuting pair is a genuine counterexample: composition is not "
        "commutative in general.",
        width=78,
    ))

if __name__ == "__main__":
    main()

"""Counting functions between small finite sets.

A function f: A -> B picks one output in B for every input in A. If A has m
elements and B has n, we can list every function as a tuple
(f(a_1), ..., f(a_m)) using ``itertools.product``. We count how many are
injective (one-to-one), surjective (onto) and bijective, and compare the
counts with the standard formulas:

    all functions   n^m
    injective       n (n-1) ... (n-m+1)          (the falling factorial)
    surjective      sum_k (-1)^k C(n,k) (n-k)^m  (inclusion-exclusion)
    bijective       m! when m = n, otherwise 0

Runs unchanged in the browser (Pyodide) using only the standard library.
"""

from itertools import product
from math import comb, factorial
import textwrap

def all_functions(domain, codomain):
    """Return every function domain -> codomain as a dict {input: output}.

    There are len(codomain) ** len(domain) of them; for an empty domain there is
    exactly one (the empty function).
    """
    domain, codomain = list(domain), list(codomain)
    return [dict(zip(domain, outputs)) for outputs in product(codomain, repeat=len(domain))]

def is_injective(f):
    """True when different inputs always give different outputs."""
    outputs = list(f.values())
    return len(set(outputs)) == len(outputs)

def is_surjective(f, codomain):
    """True when every element of the codomain is an output of f."""
    return set(f.values()) == set(codomain)

def is_bijective(f, codomain):
    """True when f is both injective and surjective."""
    return is_injective(f) and is_surjective(f, codomain)

def count_by_enumeration(m, n):
    """Enumerate all functions {0..m-1} -> {0..n-1}; return a dict of counts."""
    domain, codomain = range(m), range(n)
    functions = all_functions(domain, codomain)
    return {
        "all": len(functions),
        "injective": sum(is_injective(f) for f in functions),
        "surjective": sum(is_surjective(f, codomain) for f in functions),
        "bijective": sum(is_bijective(f, codomain) for f in functions),
    }

def falling_factorial(n, k):
    """n (n-1) ... (n-k+1), the number of ways to choose k distinct items in order.

    Equals 1 when k = 0 and 0 when k > n.
    """
    if n < 0 or k < 0:
        raise ValueError("n and k must be non-negative")
    result = 1
    for i in range(k):
        result *= n - i
    return result

def surjection_formula(m, n):
    """Number of surjections from an m-set onto an n-set, by inclusion-exclusion.

    sum over k = 0..n of (-1)^k C(n, k) (n - k)^m, using Python's 0 ** 0 == 1.
    """
    if m < 0 or n < 0:
        raise ValueError("m and n must be non-negative")
    return sum((-1) ** k * comb(n, k) * (n - k) ** m for k in range(n + 1))

def formula_counts(m, n):
    """The counts predicted by the formulas, in the same shape as :func:`count_by_enumeration`."""
    return {
        "all": n ** m,
        "injective": falling_factorial(n, m),
        "surjective": surjection_formula(m, n),
        "bijective": factorial(m) if m == n else 0,
    }

def main():
    sizes = [(m, n) for m in range(1, 5) for n in range(1, 5)]
    print("Every function A -> B, |A| = m and |B| = n, listed with itertools.product.")
    print("Each cell shows: counted by enumeration / predicted by the formula.")
    print()
    print(f"{'m':>2} {'n':>2} | {'all  (n^m)':>11} | {'injective':>9} | {'surjective':>10} | {'bijective':>9}")
    print("-" * 57)
    agree = 0
    listed = 0
    for m, n in sizes:
        counted, predicted = count_by_enumeration(m, n), formula_counts(m, n)
        listed += counted["all"]
        agree += counted == predicted
        cells = [f"{counted[k]}/{predicted[k]}" for k in ("all", "injective", "surjective", "bijective")]
        print(f"{m:>2} {n:>2} | {cells[0]:>11} | {cells[1]:>9} | {cells[2]:>10} | {cells[3]:>9}")
    print()
    print(f"Listed {listed} functions; enumeration and formulas agree for {agree} of {len(sizes)} size pairs.")
    empty = [count_by_enumeration(0, n)["all"] for n in range(5)]
    print("Edge cases: from the empty set to a set of size 0, 1, 2, 3, 4 there are")
    print(f"{empty} functions (exactly one: the empty function);")
    print(f"from a 2-element set to the empty set there are {count_by_enumeration(2, 0)['all']}.")
    print()
    print(textwrap.fill(
        f"Epistemic status: exhaustive enumeration proves these counts for the "
        f"{len(sizes)} size pairs shown, and the formulas match them exactly. That the "
        "formulas hold for every m and n is a theorem proved by counting "
        "arguments; this table is consistent with it but does not prove it.",
        width=78,
    ))

if __name__ == "__main__":
    main()

"""Euler's prime-generating polynomial: forty successes, then a failure.

p(n) = n^2 + n + 41 is prime for n = 0, 1, 2, ..., 39 - forty cases in a row.
It is tempting to conclude "always prime". n = 40 breaks the pattern:
p(40) = 1681 = 41 * 41. Evidence, however plentiful, is not proof.

Runs unchanged in the browser (Pyodide) using only the standard library.
"""

import textwrap

def euler(n):
    """Euler's polynomial n^2 + n + 41."""
    return n * n + n + 41

def is_prime(n):
    """Exact primality test by trial division (fine for the small numbers used here)."""
    if n < 2:
        return False
    if n % 2 == 0:
        return n == 2
    d = 3
    while d * d <= n:
        if n % d == 0:
            return False
        d += 2
    return True

def smallest_factor(n):
    """The smallest prime factor of an integer n >= 2."""
    if n < 2:
        raise ValueError("n must be at least 2")
    if n % 2 == 0:
        return 2
    d = 3
    while d * d <= n:
        if n % d == 0:
            return d
        d += 2
    return n

def first_failure(poly, start=0, stop=10_000):
    """Return the first n in start..stop-1 where poly(n) is not prime, or None."""
    for n in range(start, stop):
        if not is_prime(poly(n)):
            return n
    return None

def prime_count(poly, stop):
    """How many of poly(0), ..., poly(stop - 1) are prime."""
    return sum(is_prime(poly(n)) for n in range(stop))

def main():
    print("Euler's polynomial p(n) = n^2 + n + 41. Is p(n) always prime?")
    print()
    values = [euler(n) for n in range(40)]
    for row in range(0, 40, 8):
        cells = " ".join(f"{v:>5}" for v in values[row:row + 8])
        print(f"  n = {row:>2}..{row + 7:<2}: {cells}")
    streak = sum(is_prime(v) for v in values)
    print(f"All {streak} of these 40 values are prime (checked exactly by trial division).")
    print()
    n = first_failure(euler)
    v = euler(n)
    f = smallest_factor(v)
    print(f"First failure: n = {n}, p({n}) = {v} = {f} * {v // f}.")
    print(f"Why n = 40 had to fail: p(40) = 40*40 + 40 + 41 = 40*41 + 41 = 41*41.")
    print(f"And p(41) = 41*41 + 41 + 41 = 41 * 43 = {euler(41)} fails too.")
    print()
    for stop in (100, 1000):
        print(f"Primes among p(0), ..., p({stop - 1}): {prime_count(euler, stop)} of {stop}")
    print("Still remarkably many - but a pattern that is 'usually true' is not a theorem.")
    print()
    print(textwrap.fill(
        "Epistemic status: forty confirming cases were experimental evidence, "
        "and they were not a proof - the single counterexample at n = 40 proves "
        "that n^2+n+41 is not always prime. (No non-constant polynomial with "
        "integer coefficients gives primes for every n; that is a theorem, not "
        "something this program checked.)",
        width=78,
    ))

if __name__ == "__main__":
    main()

"""Exact computation versus numerical approximation.

Floating-point numbers ("floats") store a number in binary with about 16
significant decimal digits, so most decimals - even 0.1 - are rounded. That is
usually fine, but errors can show up and accumulate. Exact rational arithmetic
(``fractions.Fraction``, ``sympy.Rational``) and exact symbolic objects
(``sympy.sqrt(2)``) never round, at the cost of speed and memory.

Runs in the browser (Pyodide) with the ``sympy`` and ``numpy`` packages loaded.
"""

from fractions import Fraction
import math
import textwrap

import numpy as np
import sympy

def float_sum_check():
    """Return (0.1 + 0.2 as a float, whether it equals the float 0.3)."""
    s = 0.1 + 0.2
    return s, s == 0.3

def fraction_sum_check():
    """Return (1/10 + 2/10 as an exact Fraction, whether it equals 3/10)."""
    s = Fraction(1, 10) + Fraction(2, 10)
    return s, s == Fraction(3, 10)

def sympy_sum_check():
    """Return (1/10 + 2/10 as a sympy Rational, whether it equals 3/10)."""
    s = sympy.Rational(1, 10) + sympy.Rational(2, 10)
    return s, s == sympy.Rational(3, 10)

def repeated_sum_loop(x, times):
    """Add the float x to 0.0 one at a time, ``times`` times (rounding at every step)."""
    total = 0.0
    for _ in range(times):
        total += x
    return total

def repeated_sums(times=10_000):
    """Return {method: result} for adding 0.1 ``times`` times in several ways."""
    f64 = np.full(times, 0.1, dtype=np.float64)
    f32 = np.full(times, 0.1, dtype=np.float32)
    return {
        "Python loop, float64": repeated_sum_loop(0.1, times),
        "numpy.sum, float64 (pairwise)": float(np.sum(f64)),
        "numpy running sum, float32": float(np.cumsum(f32)[-1]),
        "math.fsum (correctly rounded)": math.fsum([0.1] * times),
        "Fraction(1, 10) added exactly": sum([Fraction(1, 10)] * times, Fraction(0)),
    }

def newton_sqrt2_float(steps, x0=1.0):
    """Newton's method x -> (x + 2/x)/2 for sqrt(2) in floats; return the list of iterates."""
    xs = [float(x0)]
    for _ in range(steps):
        x = xs[-1]
        xs.append((x + 2 / x) / 2)
    return xs

def newton_sqrt2_exact(steps, x0=Fraction(1)):
    """The same Newton iteration in exact fractions; every iterate is a rational number."""
    xs = [Fraction(x0)]
    for _ in range(steps):
        x = xs[-1]
        xs.append((x + 2 / x) / 2)
    return xs

def main():
    print("1. Is 0.1 + 0.2 equal to 0.3?")
    s, eq = float_sum_check()
    print(f"   float (rounded):   0.1 + 0.2   = {s!r:<21}  equals 0.3?  {eq}")
    s, eq = fraction_sum_check()
    print(f"   Fraction (exact):  1/10 + 2/10 = {str(s):<21}  equals 3/10? {eq}")
    s, eq = sympy_sum_check()
    print(f"   sympy.Rational:    1/10 + 2/10 = {str(s):<21}  equals 3/10? {eq}")
    print(f"   The float written 0.1 is exactly {Fraction(0.1)}.")
    print()
    times = 10_000
    print(f"2. Add 0.1 to itself {times} times (exact answer: 1000)")
    print(f"   {'method':<32}{'result':<22}error")
    for method, result in repeated_sums(times).items():
        err = abs(Fraction(result) - 1000)
        err_s = "0 (exact)" if err == 0 else f"{float(err):.1e}"
        print(f"   {method:<32}{str(result):<22}{err_s}")
    hidden = times * (Fraction(0.1) - Fraction(1, 10))
    print("   fsum adds the stored floats exactly, then rounds once: their hidden excess")
    print(f"   ({float(hidden):.2e}) is just under half the float gap at 1000"
          f" ({math.ulp(1000.0) / 2:.2e}), so 1000.0.")
    print()
    print("3. The square root of 2")
    exact = sympy.sqrt(2)
    print(f"   exact symbol:  sympy.sqrt(2)**2 == 2 is {exact ** 2 == 2}")
    print(f"   float:         math.sqrt(2)**2 = {math.sqrt(2) ** 2!r}")
    print(f"   30 digits from sympy: {sympy.N(exact, 30)}")
    print(f"   {'step':<6}{'Newton in floats':<22}{'error vs sympy':<16}Newton in fractions")
    for i, (x, q) in enumerate(zip(newton_sqrt2_float(5), newton_sqrt2_exact(5))):
        err = abs(sympy.Float(x, 30) - sympy.N(exact, 30))
        err_s = "<1e-16" if err < 1e-16 else f"{float(err):.1e}"
        print(f"   {i:<6}{x!r:<22}{err_s:<16}{q}")
    print("   Floats stall at about 16 digits. Each fraction is exact, but none squares")
    print("   to exactly 2 - no fraction can, because sqrt(2) is irrational.")
    print()
    print(textwrap.fill(
        "Epistemic status: float results are numerical approximations whose errors "
        "we measured against exact values; Fraction and sympy.Rational results are "
        "exact; sympy.sqrt(2) is an exact symbolic object. None of these "
        "computations proves sqrt(2) is irrational - that needs the classical proof.",
        width=78,
    ))

if __name__ == "__main__":
    main()

"""Searching for inverses among the integers, and among the rationals.

For each integer n we search a window of integers m for an *additive inverse*
(n + m = 0) and a *multiplicative inverse* (n * m = 1). Additive inverses are
always found; multiplicative ones only for n = 1 and n = -1. Then we move to
the rational numbers, using ``fractions.Fraction`` for exact arithmetic, where
every non-zero number has a multiplicative inverse.

Runs unchanged in the browser (Pyodide) using only the standard library.
"""

from fractions import Fraction
import textwrap

def window(radius):
    """Return the integers -radius..radius (inclusive) as a range."""
    if radius < 0:
        raise ValueError("radius must be non-negative")
    return range(-radius, radius + 1)

def find_additive_inverse(n, candidates):
    """Return the first m in ``candidates`` with n + m == 0, or None if none is found."""
    for m in candidates:
        if n + m == 0:
            return m
    return None

def find_multiplicative_inverse(n, candidates):
    """Return the first m in ``candidates`` with n * m == 1, or None if none is found.

    ``None`` only means "not found among the candidates": on its own a search
    cannot prove that no inverse exists further away.
    """
    for m in candidates:
        if n * m == 1:
            return m
    return None

def rational_inverse(q):
    """Return the exact multiplicative inverse 1/q of a rational q, or None when q == 0."""
    q = Fraction(q)
    if q == 0:
        return None
    return 1 / q

def search_summary(n_radius, m_radius):
    """Search every n in -n_radius..n_radius over candidates m in -m_radius..m_radius.

    Returns ``(additive_found, multiplicative_found)``: the lists of n for which
    an additive / multiplicative inverse was found.
    """
    candidates = window(m_radius)
    additive = [n for n in window(n_radius) if find_additive_inverse(n, candidates) is not None]
    multiplicative = [n for n in window(n_radius) if find_multiplicative_inverse(n, candidates) is not None]
    return additive, multiplicative

def main():
    candidates = window(50)
    print("For each n, search m in -50..50 for n + m = 0 and for n * m = 1.")
    print("Then, in the rationals Q, compute 1/n exactly with fractions.Fraction.")
    print()
    print(f"{'n':>4} | {'n + m = 0':>9} | {'n * m = 1':>9} | {'1/n in Q':>8} | n * (1/n)")
    print("-" * 52)
    for n in window(6):
        add = find_additive_inverse(n, candidates)
        mul = find_multiplicative_inverse(n, candidates)
        inv = rational_inverse(n)
        mul_s = "none" if mul is None else str(mul)
        inv_s = "none" if inv is None else str(inv)
        check = "-" if inv is None else str(n * inv)
        print(f"{n:>4} | {add:>9} | {mul_s:>9} | {inv_s:>8} | {check}")
    print()
    radius = 30
    additive, multiplicative = search_summary(radius, radius)
    total = 2 * radius + 1
    print(f"Wider search, n and m both in -{radius}..{radius}:")
    print(f"  additive inverse found for {len(additive)} of {total} integers")
    print(f"  multiplicative inverse found for {len(multiplicative)} of {total}: n = {multiplicative}")
    q = Fraction(-7, 12)
    print(f"  in Q, e.g. q = {q}: 1/q = {rational_inverse(q)} and q * (1/q) = {q * rational_inverse(q)} exactly")
    print()
    print("Why the integer search can stop: if |n| >= 2 and m != 0 then |n * m| >= 2,")
    print("and if m = 0 then n * m = 0. So n * m = 1 is impossible however far we look.")
    print()
    print(textwrap.fill(
        "Epistemic status: the searches are experiments over finite windows (a "
        "failed search alone proves nothing about numbers outside the window); "
        "the rational arithmetic is exact. The two-line size argument above, "
        "not the search, proves that only 1 and -1 have integer inverses.",
        width=78,
    ))

if __name__ == "__main__":
    main()

"""A random walk on the integers: what can statistics tell us?

Start at 0. At each step flip a fair coin and move +1 or -1. How likely is
the walker to come back to 0 within h steps? We estimate this by simulating
many walks (with a fixed seed, so the run is reproducible), attach a 95%
confidence interval, and compare with the exact answer, which is known:

    P(first return within 2m steps) = 1 - C(2m, m) / 4^m.

Runs unchanged in the browser (Pyodide) using only the standard library.
"""

from fractions import Fraction
from math import comb, sqrt
import random
import textwrap

SEED = 20260101
HORIZONS = (2, 10, 100, 1000)

def first_return_time(rng, max_steps):
    """Walk from 0 with fair +1/-1 steps; return the first time back at 0, or None if not within max_steps."""
    position = 0
    for t in range(1, max_steps + 1):
        position += 1 if rng.getrandbits(1) else -1
        if position == 0:
            return t
    return None

def simulate(n_walks, max_steps, seed):
    """Run ``n_walks`` independent walks with a seeded generator; return their first-return times."""
    rng = random.Random(seed)
    return [first_return_time(rng, max_steps) for _ in range(n_walks)]

def returned_within(times, horizon):
    """How many walks (from :func:`simulate`) returned to 0 within ``horizon`` steps."""
    return sum(1 for t in times if t is not None and t <= horizon)

def exact_return_probability(horizon):
    """Exact P(first return to 0 within ``horizon`` steps) = 1 - C(2m, m)/4^m with m = horizon // 2.

    Returns a Fraction. Returns are only possible at even times.
    """
    if horizon < 0:
        raise ValueError("horizon must be non-negative")
    m = horizon // 2
    return 1 - Fraction(comb(2 * m, m), 4 ** m)

def wilson_interval(successes, trials, z=1.96):
    """Wilson score confidence interval for a proportion (95% for z = 1.96).

    Returns (low, high). It behaves better than "p +/- 1.96 * standard error"
    when p is close to 0 or 1.
    """
    if trials <= 0:
        raise ValueError("trials must be positive")
    if not 0 <= successes <= trials:
        raise ValueError("successes must be between 0 and trials")
    p = successes / trials
    denom = 1 + z * z / trials
    centre = (p + z * z / (2 * trials)) / denom
    half = z * sqrt(p * (1 - p) / trials + z * z / (4 * trials * trials)) / denom
    return max(0.0, centre - half), min(1.0, centre + half)

def main():
    n_walks = 4000
    times = simulate(n_walks, max(HORIZONS), SEED)
    print(f"{n_walks} random walks from 0 (fair coin, random.Random(seed={SEED})).")
    print("Did the walk come back to 0 within h steps?")
    print()
    print(f"  {'h':>5} {'returned':>9} {'estimate':>9}   {'95% interval':<16} {'exact':>8}  inside?")
    inside_count = 0
    for h in HORIZONS:
        k = returned_within(times, h)
        low, high = wilson_interval(k, n_walks)
        exact = float(exact_return_probability(h))
        inside = low <= exact <= high
        inside_count += inside
        print(f"  {h:>5} {k:>9} {k / n_walks:>9.4f}   [{low:.4f}, {high:.4f}] {exact:>8.4f}  {'yes' if inside else 'no'}")
    print()
    print("Same question (h = 100), five different seeds, 1000 walks each:")
    estimates = [returned_within(simulate(1000, 100, seed), 100) / 1000 for seed in range(1, 6)]
    print("  " + "  ".join(f"{e:.3f}" for e in estimates)
          + f"   (exact {float(exact_return_probability(100)):.3f})")
    print("  Different seeds give different estimates: that spread is sampling error.")
    print()
    never = n_walks - returned_within(times, max(HORIZONS))
    print(f"{never} walks had not returned after {max(HORIZONS)} steps. Will they ever?")
    print("Polya's theorem (1921) says yes, with probability 1 - but no finite simulation")
    print("can show that: every finite run leaves some walks still out, and a true")
    print("probability of 0.9999 would look just like 1 in a sample of this size.")
    print()
    print(textwrap.fill(
        f"Epistemic status: statistical evidence. The exact value fell inside the "
        f"95% interval in {inside_count} of {len(HORIZONS)} rows, which supports the "
        "formula but cannot prove it (about 1 interval in 20 misses by chance). "
        "The formula and Polya's theorem are proved mathematically; this "
        "simulation only illustrates them.",
        width=78,
    ))

if __name__ == "__main__":
    main()

"""Remainder families: splitting the integers into residue classes.

Two integers belong to the same *family* (residue class) mod k when they leave
the same remainder on division by k. We split -12..12 into families, check
that the families never overlap and together contain every number, and check
that the family of a sum (or product) depends only on the families of the
numbers being added - so "family + family" makes sense. Finally we show an
operation that does NOT respect families.

Runs unchanged in the browser (Pyodide) using only the standard library.
"""

from itertools import product
import textwrap

NUMBERS = range(-12, 13)

def residue(n, k):
    """The remainder of n on division by k, always in 0..k-1 (Python's % does this even for n < 0)."""
    if k < 1:
        raise ValueError("the modulus k must be at least 1")
    return n % k

def residue_classes(numbers, k):
    """Split ``numbers`` into families mod k: {remainder: [numbers with that remainder]}.

    Every remainder 0..k-1 appears as a key, even if its family is empty.
    """
    classes = {r: [] for r in range(k)}
    for n in numbers:
        classes[residue(n, k)].append(n)
    return classes

def is_partition(classes, numbers):
    """True when the families are pairwise disjoint and together contain exactly ``numbers``."""
    seen = []
    for members in classes.values():
        seen.extend(members)
    return len(seen) == len(set(seen)) and set(seen) == set(numbers)

def depends_only_on_classes(op, numbers, k, second=None):
    """Check that residue(op(a, b)) depends only on residue(a) and residue(b).

    ``a`` runs over ``numbers`` and ``b`` over ``second`` (default: ``numbers``).
    Returns ``(cases, counterexample)`` where counterexample is None, or a pair
    of input pairs ((a, b), (a2, b2)) in the same classes whose results land in
    different classes.
    """
    second = numbers if second is None else second
    first_seen = {}
    cases = 0
    for a, b in product(numbers, second):
        cases += 1
        key = (residue(a, k), residue(b, k))
        result = residue(op(a, b), k)
        if key not in first_seen:
            first_seen[key] = ((a, b), result)
        elif first_seen[key][1] != result:
            return cases, (first_seen[key][0], (a, b))
    return cases, None

def class_table(op, k):
    """The table of op on classes mod k, computed from the representatives 0..k-1."""
    return [[residue(op(a, b), k) for b in range(k)] for a in range(k)]

def add(a, b):
    """a + b"""
    return a + b

def mul(a, b):
    """a * b"""
    return a * b

def power(a, b):
    """a ** b (used with non-negative exponents b only)"""
    return a ** b

def main():
    k = 3
    print(f"Families mod {k} for the integers -12..12 (Python's % gives remainders 0..{k - 1}):")
    for r, members in residue_classes(NUMBERS, k).items():
        print(f"  [{r}]: " + " ".join(f"{n:>3}" for n in members))
    print(f"  For example -7 = 3 * (-3) + 2, so -7 is in family [2].")
    print()
    moduli = range(1, 13)
    partitions = all(is_partition(residue_classes(NUMBERS, m), NUMBERS) for m in moduli)
    print(f"k = 1..12: families are disjoint and cover -12..12 for every k: {partitions}")
    for name, op in [("a + b", add), ("a * b", mul)]:
        total, broken = 0, []
        for m in moduli:
            cases, bad = depends_only_on_classes(op, NUMBERS, m)
            total += cases
            if bad:
                broken.append(m)
        print(f"  family of {name} depends only on the families of a and b: "
              f"{'yes' if not broken else 'NO for k = ' + str(broken)} ({total} cases)")
    print()
    print(f"So families can be added. The addition table of families mod {k}:")
    print("   + | " + "  ".join(f"[{b}]" for b in range(k)))
    for a, row in enumerate(class_table(add, k)):
        print(f"  [{a}]| " + "  ".join(f"[{v}]" for v in row))
    print()
    _, bad = depends_only_on_classes(power, [2], k, second=range(0, 7))
    (a1, b1), (a2, b2) = bad
    print(f"Not every operation respects families. Try powers of 2 mod {k}, exponents 0..6:")
    print(f"  exponents {b1} and {b2} are both in family [{residue(b1, k)}], but")
    print(f"  {a1}^{b1} = {a1 ** b1} is in [{residue(a1 ** b1, k)}] while "
          f"{a2}^{b2} = {a2 ** b2} is in [{residue(a2 ** b2, k)}].")
    print()
    print(textwrap.fill(
        "Epistemic status: every pair of numbers in -12..12 was checked for each k "
        "from 1 to 12. That is strong evidence, but it is not a proof that family "
        "addition works for all integers; the lesson's argument (if a - a2 and "
        "b - b2 are multiples of k, so is (a + b) - (a2 + b2)) is. The power "
        "example is a genuine counterexample.",
        width=78,
    ))

if __name__ == "__main__":
    main()

"""Checking set identities exhaustively.

Every subset of a small universe U is represented as a Python ``frozenset``
(and, equivalently, as a bitmask: one bit per element). We check De Morgan's
laws and the distributive laws for *every* pair or triple of subsets, and we
show that a plausible-looking "law", A - (B - C) = (A - B) - C, is false by
finding its smallest counterexample.

Runs unchanged in the browser (Pyodide) using only the standard library.
"""

from itertools import combinations, product
import textwrap

UNIVERSE = frozenset({1, 2, 3, 4})

def all_subsets(universe):
    """Return all 2^n subsets of ``universe`` as frozensets, smallest first.

    Within one size, subsets are ordered by their sorted elements, so the order
    is deterministic.
    """
    elements = sorted(universe)
    return [frozenset(c) for r in range(len(elements) + 1) for c in combinations(elements, r)]

def complement(s, universe):
    """Return U - s, the elements of the universe that are not in s."""
    return frozenset(universe) - s

def to_bitmask(s, universe):
    """Encode subset s as an integer: bit i is 1 when the i-th smallest element is in s."""
    elements = sorted(universe)
    if not set(s) <= set(elements):
        raise ValueError("s must be a subset of the universe")
    return sum(1 << i for i, x in enumerate(elements) if x in s)

def from_bitmask(mask, universe):
    """Decode an integer bitmask back into a frozenset (inverse of :func:`to_bitmask`)."""
    elements = sorted(universe)
    if not 0 <= mask < (1 << len(elements)):
        raise ValueError("mask out of range for this universe")
    return frozenset(x for i, x in enumerate(elements) if mask >> i & 1)

def bit_operations_agree(universe):
    """Check union, intersection and complement against |, & and ~ on bitmasks.

    Returns the number of (A, B) pairs checked. Raises AssertionError on a mismatch.
    """
    full = (1 << len(universe)) - 1
    subsets = all_subsets(universe)
    for a, b in product(subsets, repeat=2):
        ma, mb = to_bitmask(a, universe), to_bitmask(b, universe)
        assert to_bitmask(a | b, universe) == ma | mb
        assert to_bitmask(a & b, universe) == ma & mb
        assert to_bitmask(complement(a, universe), universe) == ~ma & full
    return len(subsets) ** 2

#: (description, arity, lhs, rhs). lhs and rhs take (U, A, B[, C]).
IDENTITIES = [
    ("U - (A ∪ B) = (U - A) ∩ (U - B)", 2,
     lambda U, A, B: U - (A | B), lambda U, A, B: (U - A) & (U - B)),
    ("U - (A ∩ B) = (U - A) ∪ (U - B)", 2,
     lambda U, A, B: U - (A & B), lambda U, A, B: (U - A) | (U - B)),
    ("A ∩ (B ∪ C) = (A ∩ B) ∪ (A ∩ C)", 3,
     lambda U, A, B, C: A & (B | C), lambda U, A, B, C: (A & B) | (A & C)),
    ("A ∪ (B ∩ C) = (A ∪ B) ∩ (A ∪ C)", 3,
     lambda U, A, B, C: A | (B & C), lambda U, A, B, C: (A | B) & (A | C)),
]

NON_IDENTITY = ("A - (B - C) = (A - B) - C", 3,
                lambda U, A, B, C: A - (B - C), lambda U, A, B, C: (A - B) - C)

def check_identity(lhs, rhs, universe, arity):
    """Compare lhs and rhs on every tuple of ``arity`` subsets of ``universe``.

    Returns ``(cases_checked, counterexamples)`` where counterexamples is a list
    of the tuples of subsets on which the two sides differ.
    """
    U = frozenset(universe)
    subsets = all_subsets(U)
    failures = []
    cases = 0
    for sets in product(subsets, repeat=arity):
        cases += 1
        if lhs(U, *sets) != rhs(U, *sets):
            failures.append(sets)
    return cases, failures

def smallest_counterexample(lhs, rhs, universe, arity):
    """Return the counterexample with the fewest elements in total, or None.

    Ties are broken by comparing the sorted element lists, so the answer is
    deterministic.
    """
    _, failures = check_identity(lhs, rhs, universe, arity)
    if not failures:
        return None
    return min(failures, key=lambda sets: (sum(len(s) for s in sets), [sorted(s) for s in sets]))

def fmt(s):
    """Format a set as {1, 2} or ∅."""
    return "{" + ", ".join(str(x) for x in sorted(s)) + "}" if s else "∅"

def main():
    U = UNIVERSE
    n = len(U)
    print(f"Universe U = {fmt(U)} has 2^{n} = {2 ** n} subsets.")
    print(f"A subset is also a {n}-bit number: {fmt({1, 3})} -> {to_bitmask({1, 3}, U):0{n}b} (bit i = element i+1).")
    pairs = bit_operations_agree(U)
    print(f"∪, ∩ and complement match bitwise |, &, ~ on all {pairs} pairs: True")
    print()
    print(f"{'statement':<34}{'cases':>7}{'failures':>10}")
    for name, arity, lhs, rhs in IDENTITIES + [NON_IDENTITY]:
        cases, failures = check_identity(lhs, rhs, U, arity)
        print(f"{name:<34}{cases:>7}{len(failures):>10}")
    print()
    name, arity, lhs, rhs = NON_IDENTITY
    A, B, C = smallest_counterexample(lhs, rhs, U, arity)
    print(f"The last line is NOT an identity. Smallest counterexample (fewest elements):")
    print(f"  A = {fmt(A)}, B = {fmt(B)}, C = {fmt(C)}")
    print(f"  A - (B - C) = {fmt(A)} - {fmt(B - C)} = {fmt(A - (B - C))}")
    print(f"  (A - B) - C = {fmt(A - B)} - {fmt(C)} = {fmt((A - B) - C)}")
    cases, fails = check_identity(lambda U, A, B, C: (A - B) - C <= A - (B - C),
                                  lambda U, A, B, C: True, U, 3)
    print(f"What does hold: (A - B) - C ⊆ A - (B - C) in all {cases} cases ({len(fails)} failures).")
    print()
    print(textwrap.fill(
        "Epistemic status: every pair and triple of subsets of this 4-element "
        "universe was checked, so the four laws are proved for this universe "
        "only, not for all sets (the lesson's element-by-element argument does "
        "that). One counterexample, by contrast, is enough to show that "
        "A - (B - C) = (A - B) - C fails for sets in general.",
        width=78,
    ))

if __name__ == "__main__":
    main()

"""Why is a negative times a negative positive? Three kinds of evidence.

1. Experimental: check (-a)(-b) = ab on a grid of integers.
2. Symbolic: let SymPy manipulate the letters a and b themselves.
3. Explanation: the lesson's proof from the distributive law, whose steps we
   display (SymPy cannot tell you *why* - it already has the rules built in).

Runs in the browser (Pyodide) with the ``sympy`` package loaded.
"""

import textwrap

import sympy

def grid(radius):
    """Return all integer pairs (a, b) with -radius <= a, b <= radius."""
    if radius < 0:
        raise ValueError("radius must be non-negative")
    values = range(-radius, radius + 1)
    return [(a, b) for a in values for b in values]

def check_on_grid(rule, pairs):
    """Apply ``rule(a, b) -> bool`` to every pair; return (cases, list of failing pairs)."""
    failures = [(a, b) for a, b in pairs if not rule(a, b)]
    return len(pairs), failures

def neg_times_neg(a, b):
    """The rule (-a)(-b) = ab for one pair of numbers."""
    return (-a) * (-b) == a * b

def neg_times_pos(a, b):
    """The rule (-a)b = -(ab) for one pair of numbers."""
    return (-a) * b == -(a * b)

def symbolic_difference(lhs, rhs):
    """Return ``sympy.expand(lhs - rhs)``; it is 0 exactly when SymPy proves lhs = rhs as polynomials."""
    return sympy.expand(lhs - rhs)

def symbolic_examples():
    """Return a list of (description, SymPy result) rows used in the printout."""
    a, b = sympy.symbols("a b")
    A, B = sympy.symbols("A B", commutative=False)
    return [
        ("(-a)*(-b)", (-a) * (-b)),
        ("(-a)*(-b) - a*b, expanded", symbolic_difference((-a) * (-b), a * b)),
        ("(-a)*b + a*b, expanded", symbolic_difference((-a) * b, -(a * b))),
        ("(a - b)*(a + b), expanded", sympy.expand((a - b) * (a + b))),
        ("how SymPy stores -a", sympy.srepr(-a)),
        ("(-A)*(-B), non-commuting A, B", (-A) * (-B)),
        ("A*B - B*A, non-commuting A, B", sympy.expand(A * B - B * A)),
    ]

#: Steps of the proof in the lesson. Each uses one ring law.
PROOF_STEPS = [
    ("0 = (-a) * 0", "x * 0 = 0, itself proved from distributivity"),
    ("  = (-a) * (b + (-b))", "b + (-b) = 0"),
    ("  = (-a)b + (-a)(-b)", "distributive law"),
    ("  = -(ab) + (-a)(-b)", "(-a)b = -(ab), proved the same way"),
    ("so (-a)(-b) = ab", "add ab to both sides"),
]

def main():
    print(f"1. Experiment: exact integer arithmetic on a grid (SymPy {sympy.__version__} used below)")
    pairs = grid(10)
    for label, rule in [("(-a)(-b) = ab   ", neg_times_neg), ("(-a)b  = -(ab)  ", neg_times_pos)]:
        cases, failures = check_on_grid(rule, pairs)
        print(f"   {label} for all a, b in -10..10: {cases} cases, {len(failures)} failures")
    print("   A grid can never reach every integer: this is evidence, not proof.")
    print()
    print("2. Symbolic algebra: SymPy works with the letters themselves")
    for description, result in symbolic_examples():
        print(f"   {description:<34} -> {result}")
    print(textwrap.fill(
        "SymPy turns (-a)*(-b) into a*b the moment it is typed: it stores -a as "
        "(-1)*a and multiplies the coefficients (-1)*(-1) = 1. A zero difference "
        "means the identity holds for all values in any commutative ring - "
        "provided SymPy's built-in rules are implemented correctly. The last two "
        "rows show the sign rule does not even need commutativity.",
        width=78, initial_indent="   ", subsequent_indent="   ",
    ))
    print()
    print("3. Why it is true (the lesson's proof, using only the ring laws):")
    for step, reason in PROOF_STEPS:
        print(f"   {step:<24} [{reason}]")
    print()
    print(textwrap.fill(
        "Epistemic status: part 1 is experimental evidence for 441 cases; part 2 "
        "is a symbolic verification that is only as trustworthy as SymPy's "
        "implementation of the ring axioms; only the proof in part 3 explains "
        "why (-a)(-b) = ab in every ring.",
        width=78,
    ))

if __name__ == "__main__":
    main()