Ordinals and cardinals

"""Ordinal arithmetic below epsilon_0, and Goodstein sequences.

Every ordinal below epsilon_0 has a Cantor normal form
w^e1*c1 + ... + w^ek*ck with e1 > ... > ek (themselves such ordinals) and
positive integer coefficients. We represent these forms as nested tuples and
implement comparison, addition and multiplication exactly, following the
recursive definitions in the lesson. The program shows that addition and
multiplication are not commutative, tests the associative and distributive
laws on random triples (left distributivity holds, right distributivity
fails), and computes Goodstein sequences, whose terms are matched with a
strictly decreasing sequence of ordinals. Evidence: exact symbolic
computation; the random tests are experiments, not proofs.
"""

import random

ZERO = ()

def nat(n):
    """The finite ordinal n."""
    if n < 0:
        raise ValueError("ordinals are non-negative")
    return () if n == 0 else ((ZERO, n),)

ONE = nat(1)
OMEGA = ((ONE, 1),)

def cmp(a, b):
    """Return -1, 0 or 1 as a < b, a = b or a > b (lexicographic on terms)."""
    for (ea, ca), (eb, cb) in zip(a, b):
        c = cmp(ea, eb)
        if c:
            return c
        if ca != cb:
            return -1 if ca < cb else 1
    return (len(a) > len(b)) - (len(a) < len(b))

def add(a, b):
    """Ordinal sum a + b: terms of a below the leading exponent of b vanish."""
    if not b:
        return a
    e, c = b[0]
    higher = tuple(t for t in a if cmp(t[0], e) > 0)
    same = [t for t in a if cmp(t[0], e) == 0]
    if same:
        return higher + ((e, same[0][1] + c),) + b[1:]
    return higher + b

def mul(a, b):
    """Ordinal product a * b, term by term in b (left distributivity)."""
    if not a or not b:
        return ZERO
    lead_e, lead_c = a[0]
    result = ZERO
    for e, c in b:
        if e == ZERO:
            term = ((lead_e, lead_c * c),) + a[1:]
        else:
            term = ((add(lead_e, e), c),)
        result = add(result, term)
    return result

def fmt(a):
    """Write a Cantor normal form, e.g. w^2*3 + w + 5."""
    if not a:
        return "0"
    parts = []
    for e, c in a:
        if e == ZERO:
            parts.append(str(c))
            continue
        if e == ONE:
            base = "w"
        elif (len(e) == 1 and e[0][0] == ZERO) or e == OMEGA:
            base = f"w^{fmt(e)}"
        else:
            base = f"w^({fmt(e)})"
        parts.append(base if c == 1 else f"{base}*{c}")
    return " + ".join(parts)

def random_ordinal(rng, max_exp=2, max_coeff=3):
    """A random ordinal below w^(max_exp + 1) with small coefficients."""
    terms = []
    for e in range(max_exp, -1, -1):
        c = rng.randint(0, max_coeff)
        if c:
            terms.append((nat(e), c))
    return tuple(terms)

def law_failures(rng, trials, law):
    """Count random triples (a, b, c) for which ``law(a, b, c)`` is False."""
    bad = []
    for _ in range(trials):
        a, b, c = (random_ordinal(rng) for _ in range(3))
        if not law(a, b, c):
            bad.append((a, b, c))
    return bad

def hereditary(n, base):
    """The ordinal obtained from n in hereditary base `base` by base -> w."""
    terms, exp = [], 0
    while n:
        n, digit = divmod(n, base)
        if digit:
            terms.append((hereditary(exp, base), digit))
        exp += 1
    return tuple(reversed(terms))

def bump(n, base):
    """Replace the base by base + 1 everywhere in hereditary notation."""
    total, exp = 0, 0
    while n:
        n, digit = divmod(n, base)
        if digit:
            total += digit * (base + 1) ** bump(exp, base)
        exp += 1
    return total

def goodstein(n, steps):
    """First terms (base, value) of the Goodstein sequence starting at n."""
    out, base = [], 2
    for _ in range(steps):
        out.append((base, n))
        if n == 0:
            break
        n = bump(n, base) - 1
        base += 1
    return out

def main():
    w = OMEGA
    two = nat(2)
    print("1. Arithmetic in Cantor normal form (w = omega):")
    rows = [
        ("1 + w", add(ONE, w), "w + 1", add(w, ONE)),
        ("2 * w", mul(two, w), "w * 2", mul(w, two)),
        ("(w + 1) * 2", mul(add(w, ONE), two), "2 * (w + 1)", mul(two, add(w, ONE))),
        ("(1 + 1) * w", mul(add(ONE, ONE), w), "1*w + 1*w", add(w, w)),
    ]
    for left, x, right, y in rows:
        print(f"   {left} = {fmt(x):<11} {right} = {fmt(y):<11} equal: {x == y}")
    rng = random.Random(2024)
    laws = [
        ("(a + b) + c = a + (b + c)",
         lambda a, b, c: add(add(a, b), c) == add(a, add(b, c))),
        ("(ab)c = a(bc)", lambda a, b, c: mul(mul(a, b), c) == mul(a, mul(b, c))),
        ("a(b + c) = ab + ac",
         lambda a, b, c: mul(a, add(b, c)) == add(mul(a, b), mul(a, c))),
        ("(a + b)c = ac + bc",
         lambda a, b, c: mul(add(a, b), c) == add(mul(a, c), mul(b, c))),
        ("a + b = b + a", lambda a, b, c: add(a, b) == add(b, a)),
    ]
    print("2. 3000 random triples below w^3 for each law:")
    examples = []
    for name, law in laws:
        bad = law_failures(rng, 3000, law)
        print(f"   {name:<27}{len(bad):>5} failures")
        if bad:
            examples.append((name, bad[0]))
    name, (a, b, c) = examples[0]
    print(f"   e.g. {name} fails for a = {fmt(a)}, b = {fmt(b)}, c = {fmt(c)}:")
    print(f"   {fmt(mul(add(a, b), c))}  versus  {fmt(add(mul(a, c), mul(b, c)))}")
    print("3. Goodstein: write n in hereditary base b, change b to b + 1, subtract 1.")
    print("   n = 3:", [v for _, v in goodstein(3, 10)])
    print("   n = 4:  base  value  ordinal (base -> w)")
    seq = goodstein(4, 7)
    for base, value in seq:
        print(f"           {base:>3} {value:>6}  {fmt(hereditary(value, base))}")
    ords = [hereditary(v, b) for b, v in seq]
    falling = all(cmp(x, y) > 0 for x, y in zip(ords, ords[1:]))
    print(f"   ordinals strictly decreasing: {falling}")
    print()
    print("Epistemic status: the normal-form arithmetic is exact, and the")
    print("non-commutativity examples are genuine counterexamples. The 3000-triple")
    print("tests are only experimental evidence for associativity and left")
    print("distributivity, which are proved by transfinite induction. Goodstein's")
    print("theorem (every sequence reaches 0) follows because ordinals are")
    print("well-ordered; n = 4 takes more steps than any computer can follow.")

if __name__ == "__main__":
    main()

"""The Schroeder-Bernstein bijection, computed by following chains.

Given injections f : A -> B and g : B -> A, follow any a in A backwards:
a = g(b)?  b = f(a')?  a' = g(b')? ... The backward chain either stops in A
(at a point not of the form g(b)), stops in B (at a point not of the form
f(a)), or never stops. The proof in the lesson sets h(a) = f(a) when the chain
stops in A and h(a) = g^(-1)(a) otherwise; these are exactly the points of the
least fixed point S of the map S -> A minus g[B minus f[S]]. Here we compute h
for two examples: the closed and half-open unit intervals (on rational points,
with exact fractions), where h turns out to be the famous "halve the powers of
1/2" map, and the natural numbers and the integers. Evidence: exact
computation on finitely many points; bijectivity on the whole set is the
theorem.
"""

from fractions import Fraction

def chain_type(a, f_inv, g_inv, limit=200):
    """Follow a in A backwards; return 'A', 'B' (where it stops) or 'endless'."""
    seen = set()
    x, side = a, "A"
    for _ in range(limit):
        if (x, side) in seen:
            return "endless"
        seen.add((x, side))
        prev = g_inv(x) if side == "A" else f_inv(x)
        if prev is None:
            return side
        x, side = prev, ("B" if side == "A" else "A")
    return "endless"

def sb_map(a, f, f_inv, g_inv):
    """The Schroeder-Bernstein bijection: f on A-stopping chains, else g^-1."""
    if chain_type(a, f_inv, g_inv) == "A":
        return f(a)
    return g_inv(a)

# Example 1: A = [0, 1], B = [0, 1), f(x) = x/2 (into B), g(y) = y (into A).
def f1(x):
    return x / 2

def f1_inv(y):
    return 2 * y if y <= Fraction(1, 2) else None

def g1_inv(x):
    return x if x < 1 else None

# Example 2: A = N, B = Z, f(n) = n, g(z) = 3z (z >= 0), -3z - 1 (z < 0).
def f2(n):
    return n

def f2_inv(z):
    return z if z >= 0 else None

def g2_inv(n):
    if n % 3 == 0:
        return n // 3
    if n % 3 == 2:
        return -(n + 1) // 3
    return None

def rational_points(max_den):
    """All fractions p/q in [0, 1] with q <= max_den, sorted."""
    pts = {Fraction(p, q) for q in range(1, max_den + 1) for p in range(q + 1)}
    return sorted(pts)

def main():
    print("Example 1: A = [0, 1], B = [0, 1), f(x) = x/2, g(y) = y.")
    for x in [Fraction(1), Fraction(1, 2), Fraction(1, 8), Fraction(3, 4),
              Fraction(1, 3), Fraction(0)]:
        kind = chain_type(x, f1_inv, g1_inv)
        label = "never stops" if kind == "endless" else f"stops in {kind}"
        print(f"   x = {str(x):>4}: backward chain {label:<11}  h(x) = "
              f"{sb_map(x, f1, f1_inv, g1_inv)}")
    pts = rational_points(30)
    images = [sb_map(x, f1, f1_inv, g1_inv) for x in pts]
    in_b = all(0 <= y < 1 for y in images)
    moved = sorted(x for x, y in zip(pts, images) if y != x)
    print(f"   {len(pts)} rational points with denominator <= 30:")
    print(f"   images in [0, 1): {in_b}; injective: "
          f"{len(set(images)) == len(images)}")
    print(f"   points moved: {[str(x) for x in reversed(moved)]}")
    print("   h halves 1, 1/2, 1/4, ... and fixes everything else.")
    print()
    print("Example 2: A = N, B = Z, f(n) = n, g(z) = 3z or -3z - 1.")
    h = [sb_map(n, f2, f2_inv, g2_inv) for n in range(3000)]
    print("   h(0..14) =", h[:15])
    window = set(range(-200, 201))
    hit = window <= set(h)
    print(f"   h injective on 0..2999: {len(set(h)) == len(h)}; "
          f"hits every z in -200..200: {hit}")
    kinds = [chain_type(n, f2_inv, g2_inv) for n in range(30)]
    print("   chains of 0..29 stop in A / B / never:",
          [kinds.count(k) for k in ("A", "B", "endless")])
    print()
    print("Epistemic status: exact arithmetic, checked on finitely many points")
    print(f"({len(pts)} rationals; 3000 natural numbers). That h is a bijection of the")
    print("whole sets is the Schroeder-Bernstein theorem, proved in the lesson;")
    print("the chain classification uses no choice, so the theorem holds in ZF.")

if __name__ == "__main__":
    main()