Compute Goodstein's sequences — https://lipsum.dev/2023-03-1-logique-arithmetiques-completude/

goodstein.py

from typing import Iterator


HereditaryRepresentation = list[tuple[int, "HereditaryRepresentation"]]

def to_base(n: int, base: int, prefix: bool = False) -> list[int]:
    """
    Computes the representation of n in the specified base.
    """
    if n == 0:
        return [] if prefix else [0]
    q, r = divmod(n, base)
    return to_base(q, base, True) + [r]

def hereditary(n: int, base: int) -> HereditaryRepresentation:
    """
    Returns the hereditary representation of n in the specified base.

    Note that HereditaryRepresentation doesn't make a reference to the
    specific base in this implementation.
    """
    base_repr = to_base(n, base)
    if len(base_repr) == 1:
        return [(base_repr[0], [])]
    her_repr = []
    for pos, i in enumerate(base_repr):
        her_repr.append((i, hereditary(len(base_repr) - pos - 1, base)))
    return her_repr

def evaluate(her_repr: HereditaryRepresentation, base: int) -> int:
    """
    Evaluates an hereditary representation in the specified base.
    """
    return sum(a * base**evaluate(b, base) for a, b in her_repr)

def to_string(hereditary_repr: HereditaryRepresentation, base: int | str) -> str:
    """
    Returns a latex-friendly representation of the hereditary representation,
    with the specified base displayed.
    Pass base="ω" to obtain a representation of the corresponding ordinal.
    """
    parts = []
    for a, b in hereditary_repr:
        if a == 0:
            continue
        exp = to_string(b, base)
        if exp == "0":
            parts.append(str(a))
        elif a == 1:
            parts.append(f"{base}^{{{exp}}}")
        elif a > 1:
            parts.append(f"{base}^{{{exp}}} \\times {a}")
    if not parts:
        return "0"
    return " + ".join(parts)

def g_seq(n: int) -> Iterator[tuple[int, int, HereditaryRepresentation]]:
    """
    Iterates over the Goodstein sequence with seed n, yielding tuples
    corresponding to the base, the value and the hereditary representation.
    """
    base = 2
    while True:
        her_repr = hereditary(n, base)
        yield base, n, her_repr
        if n == 0:
            continue
        base += 1
        n = evaluate(her_repr, base) - 1

def explain_g_seq(m: int, k: int, base_symbol="b") -> None:
    """
    Displays the detail of the k first values of the Goodstein sequence of seed m,
    in a Markdown table format (without headers).
    """
    seq = g_seq(m)
    for n in range(k):
        base, value, her_repr = next(seq)
        print(f"|$G({m}, {n})$|b={base}|${to_string(her_repr, base_symbol)}$|{value}|")

# See https://oeis.org/A056193
seq4 = g_seq(4)
expected = [4, 26, 41, 60, 83, 109, 139, 173, 211, 253, 299, 348]
assert [next(seq4)[1] for _ in range(len(expected))] == expected

seq19 = g_seq(19)
assert next(seq19)[1] == 19
assert next(seq19)[1] == 7625597484990

seq3 = g_seq(3)
expected = [3, 3, 3, 2, 1, 0, 0]
assert [next(seq3)[1] for _ in range(len(expected))] == expected

# explain_g_seq(3, 6)
# explain_g_seq(7, 7)
# explain_g_seq(7, 7, "ω")