ogdoad 1.0.0

//! The transfinite nim-multiplication tower: ordinals `< ω^(ω^ω)` (Conway's
//! algebraically closed field `\bar{F_2}`) as monomials in the prime-power
//! generators `χ_{u^n}`, multiplied by generator-power-vector addition with the
//! per-prime carry relations.
//!
//! This generalizes the degree-3 cube-root tower — the prime-3, exponent-place-`ω^0`
//! special case — to all primes and places (Conway *ONAG* ch. 6; Lenstra *Nim
//! multiplication* 1978; DiMuro *On Onp* arXiv:1108.0962, Thm 3.1 + Table 1).
//!
//! ## The uniform representation
//!
//! Every ordinal `< ω^(ω^ω)` is `Σ ω^E · c` in CNF, with each exponent `E < ω^ω`
//! and `c` a finite nimber. An `E < ω^ω` in turn has CNF `Σ ω^m · V_m` with **finite**
//! exponents `m` and finite coefficients `V_m`. The generator picture:
//!
//! - the exponent-place `ω^m` of `E` is governed by the prime `p(m)` = the
//!   `(m+2)`-th prime (`p(0)=3`, `p(1)=5`, `p(2)=7`, …);
//! - reading `V_m` in base `p(m)` gives the digit `d_{m,k}` = the power of the
//!   generator `χ_{p(m)^{k+1}}` (each `d_{m,k} < p(m)`);
//! - so a pure monomial `ω^E` is `⊗_{m,k} χ_{p(m)^{k+1}}^{⊗ d_{m,k}}`, and a general
//!   element is a finite nim-sum of such monomials with finite-nimber coefficients.
//!
//! We therefore key a monomial by its exponent's decomposition,
//! [`GenKey`] = `place m ↦ base-p(m) digit vector`.
//!
//! ## Multiplication
//!
//! `ω^{E1} ⊗ ω^{E2}`: the generator powers **add** (`χ^i ⊗ χ^j = χ^{i+j}`), i.e. the
//! two keys' digit vectors add componentwise per `(m,k)`; then each place reduces by
//! the tower carries (high→low):
//!
//! - `χ_{u^{k+1}}^u = χ_{u^k}` (`k ≥ 1`): a digit `≥ u` at level `k` carries one down
//!   to level `k-1` — exact, keeps a single monomial;
//! - `χ_u^u = α_u` (`k = 0`, the **Kummer** relation): a digit `≥ u` at level 0
//!   consumes a factor of the *excess* `α_u`.
//!
//! ## The branching reduction (non-scalar `α_u`)
//!
//! The excess `α_u` is a finite nimber for some primes (`α_3=2`, `α_5=4`, `α_17=16`)
//! and a genuine transfinite ordinal for others (`α_7=ω+1`, `α_11=ω^ω+1`, `α_13=ω+4`,
//! …). A scalar `α_u` keeps a level-0 carry inside the coefficient — the product stays
//! one monomial. A non-scalar `α_u` is a *sum*, so the carry **branches** the monomial:
//! `χ_u^u = α_u` replaces a generator power by `α_u`, and the (reduced) monomial must be
//! nim-multiplied by that sum, mixing across exponent places.
//!
//! This recursion **descends by place**: every `α_{p(m)}` is built from generators at
//! places strictly `< m` (`α_7 = ω+1` uses `ω = χ_3`, place 0 < 2; `α_11 = ω^ω+1` uses
//! `χ_5`, place 1 < 3; verified from DiMuro Table 1). So `base ⊗ excess` can never
//! re-trigger a carry at the place that produced it, and the recursion bottoms out at
//! `α_3 = 2` in the finite field — the crate's "recurse only on strictly-simpler
//! exponents" discipline. Termination depth is bounded by the largest place index.
//!
//! ## Staging (honest scope)
//!
//! We carry the finite Lenstra excess integers `m_u` from **OEIS A380496** ("Lenstra
//! excess of the n-th odd prime"): the b-file's 126 known rows, the odd primes
//! `3..=709` (see [`finite_excess`] and `docs/OPEN.md`). The ordinal expression is
//! assembled in code: compute `f(u)=ord_u(2)`, compute DiMuro's `Q(f(u))`, form the
//! `χ`-sum, then nim-add `m_u`. The product of any two ordinals `< ω^(ω^ω)` is
//! therefore exact whenever every Kummer carry it triggers is at a prime `≤ 709`; a
//! carry needing `m_719` (the first OEIS-unknown row) or beyond returns `None` — the
//! honest operational boundary. (Anything `≥ ω^(ω^ω)`, an infinite exponent place, is
//! out of range outright.)

use super::Ordinal;
use crate::scalar::{is_prime_u128, nim_mul};
use std::collections::{BTreeMap, BTreeSet};

/// A monomial's exponent `E < ω^ω`, decomposed per `ω`-place and per prime:
/// `place m ↦ base-p(m) digit vector`, where digit `k` is the power of the generator
/// `χ_{p(m)^{k+1}}` (`0 ≤ digit < p(m)`). An absent place / empty vector is all-zero;
/// the empty map is the exponent `0` (the monomial `1`).
type GenKey = BTreeMap<u128, Vec<u128>>;

/// The prime governing exponent-place `ω^m`: `p(m)` = the `(m+2)`-th prime
/// (`p(0)=3`, `p(1)=5`, `p(2)=7`, …). Prime 2 is excluded — the prime-2 (Fermat)
/// tower is the finite nimber field, handled by [`crate::scalar::nim_mul`].
pub(super) fn place_prime(m: u128) -> u128 {
    let mut count = 0u128;
    let mut n = 2u128; // skip the prime 2
    loop {
        n += 1;
        if is_prime_u128(n) {
            count += 1;
            if count == m + 1 {
                return n;
            }
        }
    }
}

/// The excess `α_u` (`χ_u^u = α_u`, the Kummer relation) as an ordinal, or `None` for
/// primes beyond the verified table (`u > 709` — the staged boundary). The only
/// hardcoded row datum is the finite excess integer `m_u` (see [`finite_excess`]); the
/// ordinal expression is reconstructed from `f(u)=ord_u(2)` and DiMuro's recursive
/// `Q(f(u))`. Every `α_u` is built from generators at strictly-lower places than
/// `χ_u`'s own (the prime factors of `f(u) | u-1` are all `< u`), which is what makes
/// the branching reduction descend and terminate.
pub(super) fn alpha_ordinal(u: u128) -> Option<Ordinal> {
    let f = multiplicative_order_two_mod_prime(u)?;
    let mut val = chi_sum(&q_set(f)?)?;
    val = val.nim_add(&Ordinal::from_u128(finite_excess(u)?));
    Some(val)
}

/// The finite Lenstra excess `m_p` for an odd prime `p`, from **OEIS A380496**
/// ("Lenstra excess of the n-th odd prime"). The b-file's 126 known rows cover the odd
/// primes `3..=709`, indexed here by the 0-based odd-prime place (place 0 = prime 3),
/// so `EXCESS[odd_prime_place(p)] = m_p`. `None` past prime 709 — the staged boundary,
/// where the first OEIS-unknown row is `p = 719`. The place ↦ prime map is the tower's
/// own coordinate (`place_prime`, pinned by `place_primes_are_the_odd_primes`); the
/// values are pinned against the b-file in `excess_table_matches_oeis_a380496`.
/// (Provenance: `a(1–3)` Conway, `a(4–13)` Lenstra, `a(14–18)` Le Bruyn, `a(19–59)`
/// Siegel, `a(60–126)` Peeters, via CGSuite's Lenstra-excess calculator. The excess is
/// `0` or `1` except `m_19 = m_163 = 4`. See `docs/OPEN.md` for the open boundary.)
fn finite_excess(u: u128) -> Option<u128> {
    // OEIS A380496 b-file column `a(n)`, `n = 1..126`, i.e. the odd primes `3..=709` in
    // order (`a(n)` at index `n-1` = the (n-1)-th place). Diff against `b380496.txt`.
    const EXCESS: [u128; 126] = [
        0, 0, 1, 1, 0, 0, 4, 1, 0, 1, 0, 1, 1, 1, // 3 .. 47   (DiMuro Table 1 + m_47)
        1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, // 53 .. 109
        0, 1, 0, 1, 0, 1, 0, 1, 4, 1, 0, 1, 0, 0, // 113 .. 191  (m_163 = 4)
        0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 1, // 193 .. 269
        0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 1, 1, 1, 0, // 271 .. 353
        1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, // 359 .. 439
        0, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, // 443 .. 523
        0, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, // 541 .. 617
        1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, // 619 .. 709
    ];
    EXCESS.get(odd_prime_place(u)? as usize).copied()
}

fn multiplicative_order_two_mod_prime(u: u128) -> Option<u128> {
    if u <= 2 || !is_prime_u128(u) {
        return None;
    }
    let mut x = 2 % u;
    let mut k = 1u128;
    while x != 1 {
        x = x.checked_mul(2)? % u;
        k = k.checked_add(1)?;
        if k > u - 1 {
            return None;
        }
    }
    Some(k)
}

fn q_set(h: u128) -> Option<Vec<u128>> {
    if h == 1 {
        return Some(Vec::new());
    }
    let (r, g) = smallest_prime_power_factor(h)?;
    let mut qs = q_set(g)?;
    let chi_g = chi_sum(&qs)?;
    let degree = chi_g.finite_subfield_degree()?;
    if !degree.is_multiple_of(r) {
        qs.push(r);
        qs.sort_unstable();
    }
    Some(qs)
}

fn smallest_prime_power_factor(h: u128) -> Option<(u128, u128)> {
    if h <= 1 {
        return None;
    }
    let u = smallest_prime_factor(h)?;
    let mut r = 1u128;
    let mut g = h;
    while g.is_multiple_of(u) {
        r = r.checked_mul(u)?;
        g /= u;
    }
    Some((r, g))
}

fn smallest_prime_factor(n: u128) -> Option<u128> {
    if n < 2 {
        return None;
    }
    if n.is_multiple_of(2) {
        return Some(2);
    }
    let mut d = 3u128;
    while d <= n / d {
        if n.is_multiple_of(d) {
            return Some(d);
        }
        d += 2;
    }
    Some(n)
}

fn chi_sum(qs: &[u128]) -> Option<Ordinal> {
    qs.iter()
        .try_fold(Ordinal::zero(), |acc, &q| Some(acc.nim_add(&chi(q)?)))
}

fn chi(q: u128) -> Option<Ordinal> {
    let (p, n) = prime_power(q)?;
    if p == 2 {
        let shift = q / 2;
        if shift >= u128::BITS as u128 {
            return None;
        }
        return Some(Ordinal::from_u128(1u128 << shift));
    }
    let place = odd_prime_place(p)?;
    let coeff = checked_pow(p, n - 1)?;
    let exp = Ordinal::monomial(Ordinal::from_u128(place), coeff);
    Some(Ordinal::omega_pow(exp))
}

fn prime_power(q: u128) -> Option<(u128, u128)> {
    if q < 2 {
        return None;
    }
    let p = smallest_prime_factor(q)?;
    let mut n = 0u128;
    let mut rest = q;
    while rest.is_multiple_of(p) {
        n += 1;
        rest /= p;
    }
    (rest == 1).then_some((p, n))
}

fn odd_prime_place(p: u128) -> Option<u128> {
    if p == 2 || !is_prime_u128(p) {
        return None;
    }
    let mut place = 0u128;
    loop {
        let q = place_prime(place);
        if q == p {
            return Some(place);
        }
        if q > p {
            return None;
        }
        place += 1;
    }
}

/// Shared with [`subfield`](super::subfield) (`pub(super)`, not duplicated there).
pub(super) fn checked_pow(base: u128, exp: u128) -> Option<u128> {
    let mut acc = 1u128;
    for _ in 0..exp {
        acc = acc.checked_mul(base)?;
    }
    Some(acc)
}

/// Base-`base` digit vector of `v` (least-significant first, no trailing zeros).
/// Shared with [`subfield`](super::subfield) (`pub(super)`, not duplicated there).
pub(super) fn base_digits(mut v: u128, base: u128) -> Vec<u128> {
    let mut d = Vec::new();
    while v > 0 {
        d.push(v % base);
        v /= base;
    }
    d
}

/// Decompose an exponent `E` into its [`GenKey`], or `None` if `E ≥ ω^ω` (some CNF
/// exponent of `E` is itself infinite — the whole ordinal is then `≥ ω^(ω^ω)`, beyond
/// the algebraically-closed segment this tower represents).
fn decompose_exp(e: &Ordinal) -> Option<GenKey> {
    let mut key = GenKey::new();
    for (exp, c) in e.terms() {
        let m = exp.as_finite()?; // infinite place ⇒ E ≥ ω^ω ⇒ out of range
        key.insert(m, base_digits(*c, place_prime(m)));
    }
    Some(key)
}

/// Rebuild the exponent ordinal `E = Σ ω^m · V_m` from a [`GenKey`] (`V_m` read back
/// from its base-`p(m)` digits).
fn recompose_exp(key: &GenKey) -> Option<Ordinal> {
    key.iter().try_fold(Ordinal::zero(), |acc, (&m, digits)| {
        let u = place_prime(m);
        let mut v: u128 = 0;
        let mut pw: u128 = 1;
        for &d in digits {
            let term = d.checked_mul(pw)?;
            v = v.checked_add(term)?;
            pw = pw.checked_mul(u)?;
        }
        Some(if v == 0 {
            acc
        } else {
            acc.nim_add(&Ordinal::monomial(Ordinal::from_u128(m), v))
        })
    })
}

/// Reduce one place's raw (post-addition) generator-power digits to canonical digits
/// `< u`, returning the canonical digits and the number of **level-0 (Kummer) carries**
/// `q` (each owes one factor of the excess `α_u`, resolved by the caller). Processes
/// high→low: a digit `≥ u` at level `k ≥ 1` carries one to level `k-1`
/// (`χ_{u^{k+1}}^u = χ_{u^k}`); at level 0 it is removed and counted (`χ_u^u = α_u`).
fn reduce_place(raw: &[u128], u: u128) -> Option<(Vec<u128>, u128)> {
    let mut d = raw.to_vec();
    for k in (0..d.len()).rev() {
        let carry = d[k] / u;
        d[k] %= u;
        if carry == 0 {
            continue;
        }
        if k == 0 {
            while d.last() == Some(&0) {
                d.pop();
            }
            return Some((d, carry));
        }
        d[k - 1] = d[k - 1].checked_add(carry)?;
    }
    let mut digits: Vec<u128> = d;
    while digits.last() == Some(&0) {
        digits.pop();
    }
    Some((digits, 0))
}

/// Add two generator monomials' digit vectors per `(m,k)` and reduce each place,
/// returning the canonical base [`GenKey`] and the per-place count of level-0 Kummer
/// carries (the excess `α_{p(m)}` owed). Pure digit bookkeeping — no `α` resolution.
fn reduce_keys(a: &GenKey, b: &GenKey) -> Option<(GenKey, BTreeMap<u128, u128>)> {
    let mut base = GenKey::new();
    let mut overflow: BTreeMap<u128, u128> = BTreeMap::new();
    let places: BTreeSet<u128> = a.keys().chain(b.keys()).copied().collect();
    for m in places {
        let da = a.get(&m).map(Vec::as_slice).unwrap_or(&[]);
        let db = b.get(&m).map(Vec::as_slice).unwrap_or(&[]);
        let len = da.len().max(db.len());
        let raw: Vec<u128> = (0..len)
            .map(|i| {
                da.get(i)
                    .copied()
                    .unwrap_or(0)
                    .checked_add(db.get(i).copied().unwrap_or(0))
            })
            .collect::<Option<_>>()?;
        let (red, q) = reduce_place(&raw, place_prime(m))?;
        if q > 0 {
            overflow.insert(m, q);
        }
        if !red.is_empty() {
            base.insert(m, red);
        }
    }
    Some((base, overflow))
}

/// The product of two generator monomials `ω^{E_a}·c_a` and `ω^{E_b}·c_b`, as a full
/// ordinal (a *sum*, once a non-scalar Kummer carry branches it). Adds the generator
/// powers, reduces, then nim-multiplies in the excess `α` factors the level-0 carries
/// owe — recursively, since `α_u` is itself a (strictly-lower-place) ordinal. `None` if
/// some owed `α_u` is past the verified table (`u > 709`).
fn mul_mono(ka: &GenKey, ca: u128, kb: &GenKey, cb: u128) -> Option<Ordinal> {
    let (base_key, overflow) = reduce_keys(ka, kb)?;
    let coeff = nim_mul(ca, cb);
    let base = if base_key.is_empty() {
        Ordinal::from_u128(coeff)
    } else {
        Ordinal::monomial(recompose_exp(&base_key)?, coeff)
    };
    if overflow.is_empty() {
        return Some(base);
    }
    // Excess factor `∏_m α_{p(m)}^{⊗ q_m}`. Each `α` lives at places `< m` (DiMuro), so
    // both this fold and `base ⊗ excess` descend in place and terminate.
    let mut excess = Ordinal::from_u128(1);
    for (&m, &q) in &overflow {
        let alpha = alpha_ordinal(place_prime(m))?;
        for _ in 0..q {
            excess = mul(&excess, &alpha)?;
        }
    }
    mul(&base, &excess)
}

/// Nim-multiply two ordinals `< ω^(ω^ω)`, or `None` outside that range / when a Kummer
/// carry needs an excess `α_u` past the verified table (`u > 709`). Distributes over CNF
/// (char-2 field addition = nim-add); each monomial pair is handled by [`mul_mono`].
pub(super) fn mul(a: &Ordinal, b: &Ordinal) -> Option<Ordinal> {
    let mut acc = Ordinal::zero();
    for (ea, ca) in a.terms() {
        let ka = decompose_exp(ea)?;
        for (eb, cb) in b.terms() {
            let kb = decompose_exp(eb)?;
            let term = mul_mono(&ka, *ca, &kb, *cb)?;
            acc = acc.nim_add(&term);
        }
    }
    Some(acc)
}

#[cfg(test)]
mod tests {
    use super::*;

    fn fin(n: u128) -> Ordinal {
        Ordinal::from_u128(n)
    }
    fn w() -> Ordinal {
        Ordinal::omega()
    }
    fn ww() -> Ordinal {
        Ordinal::omega_pow(Ordinal::omega()) // ω^ω
    }
    fn chi7() -> Ordinal {
        Ordinal::omega_pow(Ordinal::omega_pow(fin(2))) // ω^(ω²) = χ_7
    }
    /// `χ_7^⊗n` by repeated nim-multiplication.
    fn chi7_pow(n: u128) -> Ordinal {
        let mut p = fin(1);
        for _ in 0..n {
            p = mul(&p, &chi7()).unwrap();
        }
        p
    }

    fn expected_alpha(u: u128) -> Ordinal {
        match u {
            3 => fin(2),
            5 => fin(4),
            7 => w().nim_add(&fin(1)),
            11 => Ordinal::omega_pow(w()).nim_add(&fin(1)),
            13 => w().nim_add(&fin(4)),
            17 => fin(16),
            19 => Ordinal::omega_pow(fin(3)).nim_add(&fin(4)),
            23 => Ordinal::omega_pow(Ordinal::omega_pow(fin(3))).nim_add(&fin(1)),
            29 => Ordinal::omega_pow(Ordinal::omega_pow(fin(2))).nim_add(&fin(4)),
            31 => Ordinal::omega_pow(w()).nim_add(&fin(1)),
            37 => Ordinal::omega_pow(fin(3)).nim_add(&fin(4)),
            41 => Ordinal::omega_pow(w()).nim_add(&fin(1)),
            43 => Ordinal::omega_pow(Ordinal::omega_pow(fin(2))).nim_add(&fin(1)),
            47 => Ordinal::omega_pow(Ordinal::omega_pow(fin(7))).nim_add(&fin(1)),
            _ => panic!("unexpected test row"),
        }
    }

    #[test]
    fn place_primes_are_the_odd_primes() {
        for (m, p) in [
            (0, 3),
            (1, 5),
            (2, 7),
            (3, 11),
            (4, 13),
            (5, 17),
            (6, 19),
            (7, 23),
            (8, 29),
            (9, 31),
            (10, 37),
            (11, 41),
            (12, 43),
            (13, 47),
            (14, 53),
        ] {
            assert_eq!(place_prime(m), p);
        }
    }

    #[test]
    fn dimuro_rows_are_assembled_from_order_qset_and_finite_excess() {
        // Test-only transcription of DiMuro Table 1 plus the 47 row (= OEIS A380496
        // a(1)..a(14)). The production code above hardcodes only `m_u`; these checks pin the
        // computed `f(u)`, recursive Q-set, and final χ-sum separately.
        for (u, f, qs, m) in [
            (3, 2, &[2][..], 0),
            (5, 4, &[4][..], 0),
            (7, 3, &[3][..], 1),
            (11, 10, &[5][..], 1),
            (13, 12, &[3, 4][..], 0),
            (17, 8, &[8][..], 0),
            (19, 18, &[9][..], 4),
            (23, 11, &[11][..], 1),
            (29, 28, &[4, 7][..], 0),
            (31, 5, &[5][..], 1),
            (37, 36, &[4, 9][..], 0),
            (41, 20, &[5][..], 1),
            (43, 14, &[7][..], 1),
            (47, 23, &[23][..], 1),
        ] {
            assert_eq!(multiplicative_order_two_mod_prime(u), Some(f));
            assert_eq!(q_set(f).as_deref(), Some(qs));
            assert_eq!(finite_excess(u), Some(m));
            assert_eq!(alpha_ordinal(u), Some(expected_alpha(u)));
        }
    }

    #[test]
    fn beyond_dimuro_table_rows_are_assembled_from_order_qset_and_finite_excess() {
        // Rows past DiMuro Table 1 (u=47 already lives above): `f(u)` here is computed
        // independently (plain modular-order arithmetic, not run through this module),
        // `qs` is cross-checked against `experiments/ordinal_excess_probe.py`'s
        // hand-curated `Q_SET` (itself sourced from DiMuro/CGSuite, not from this file),
        // and `m` against the vendored OEIS A380496 b-file. `expected_alpha` is then a
        // literal CNF built from `chi_sum(qs) + m` by hand, cross-checked against
        // structurally related rows already pinned above:
        //   - u=73: same q_set {9} as the already-pinned u=19, but excess 1 not 4 — the
        //     `Q={9}` exception-column contrast OPEN.md documents (`m_19=4`, `m_73=1`).
        //   - u=89: f(89)=11 and q_set(11)=[11] happen to coincide exactly with u=23's
        //     row, so alpha_89 == alpha_23 by construction (same chi_sum, same excess) —
        //     an internal-consistency cross-check, not a coincidence bug.
        // u=179 (f=178=2·89, q_set={89}, m=1, probe-certified) is deliberately NOT pinned
        // here: `alpha_ordinal(179)` recurses into `finite_subfield_degree` on
        // `chi(89)=ω^(ω^22)`, whose Frobenius minimization is exactly the "huge component
        // field" cost this module's docs already warn about (confirmed by hand — it does
        // not finish in a unit-test budget). `finite_excess(179)` alone (the OEIS integer,
        // not the ordinal reconstruction) is still exercised by the full b-file diff above.
        for (u, f, qs, m, expected) in [
            (
                73,
                9,
                &[9][..],
                1,
                Ordinal::omega_pow(fin(3)).nim_add(&fin(1)),
            ),
            (
                89,
                11,
                &[11][..],
                1,
                Ordinal::omega_pow(Ordinal::omega_pow(fin(3))).nim_add(&fin(1)),
            ),
        ] {
            assert_eq!(multiplicative_order_two_mod_prime(u), Some(f), "f({u})");
            assert_eq!(q_set(f).as_deref(), Some(qs), "q_set(f({u}))");
            assert_eq!(finite_excess(u), Some(m), "m_{u}");
            assert_eq!(alpha_ordinal(u), Some(expected), "alpha_{u}");
        }
    }

    #[test]
    fn alpha_excesses_descend_in_place() {
        // The termination invariant: every verified α_u is built from generators at
        // places strictly below χ_u's own place (the (u)-index in the odd primes). This
        // is what makes `base ⊗ excess` descend; a typo that violated it would loop.
        for m in 0..=13u128 {
            let u = place_prime(m);
            let alpha = alpha_ordinal(u).unwrap();
            // the highest place appearing anywhere in α_u must be < m.
            let mut hi: Option<u128> = None;
            for (exp, _) in alpha.terms() {
                if let Some(sub) = decompose_exp(exp) {
                    if let Some(&mx) = sub.keys().last() {
                        hi = Some(hi.map_or(mx, |h| h.max(mx)));
                    }
                }
            }
            if let Some(h) = hi {
                assert!(h < m, "α_{u} reaches place {h} ≥ its own place {m}");
            }
        }
    }

    #[test]
    fn reproduces_cube_tower_below_omega_omega() {
        // ω ⊗ ω = ω², ω⊗³ = 2, ω⊗⁴ = ω·2 — the prime-3, place-0 behavior.
        let wsq = mul(&w(), &w()).unwrap();
        assert_eq!(wsq, Ordinal::omega_pow(fin(2)));
        assert_eq!(mul(&wsq, &w()).unwrap(), fin(2)); // ω³ = 2
        assert_eq!(mul(&wsq, &wsq).unwrap(), Ordinal::monomial(fin(1), 2)); // ω⁴ = ω·2
    }

    #[test]
    fn quintic_landmarks_from_dimuro() {
        // ω^ω = χ_5, the degree-5 generator: free powers are ordinary ordinal powers,
        // then the Kummer reduction (ω^ω)⊗⁵ = α_5 = 4.
        let w2 = mul(&ww(), &ww()).unwrap();
        assert_eq!(w2, Ordinal::omega_pow(Ordinal::monomial(fin(1), 2))); // ω^{ω·2}
        let w3 = mul(&w2, &ww()).unwrap();
        assert_eq!(w3, Ordinal::omega_pow(Ordinal::monomial(fin(1), 3))); // ω^{ω·3}
        let w4 = mul(&w3, &ww()).unwrap();
        assert_eq!(w4, Ordinal::omega_pow(Ordinal::monomial(fin(1), 4))); // ω^{ω·4}
        let w5 = mul(&w4, &ww()).unwrap();
        assert_eq!(w5, fin(4)); // (ω^ω)⊗⁵ = α_5 = 4
    }

    #[test]
    fn cross_place_products() {
        // ω^ω ⊗ ω = ω^{ω+1} (exponents ω and 1 add, no carry).
        assert_eq!(
            mul(&ww(), &w()).unwrap(),
            Ordinal::omega_pow(Ordinal::omega().nim_add(&fin(1)))
        );
        // ω^ω ⊗ 2 = ω^ω·2 (a finite-nimber coefficient).
        assert_eq!(
            mul(&ww(), &fin(2)).unwrap(),
            Ordinal::monomial(Ordinal::omega(), 2)
        );
    }

    #[test]
    fn septic_kummer_landmark() {
        // THE Stage-2 headline, from DiMuro Table 1 (NOT from the engine — non-circular):
        // χ_7 = ω^(ω²), and χ_7^⊗7 = α_7 = ω + 1. The 7th power is the first non-scalar
        // Kummer carry; it branches the monomial into the sum ω + 1.
        assert_eq!(chi7_pow(7), w().nim_add(&fin(1))); // ω + 1

        // χ_7^⊗8 = α_7 ⊗ χ_7 = (ω+1)⊗ω^(ω²) = ω^(ω²+1) + ω^(ω²).
        let e_w2_1 = Ordinal::omega_pow(fin(2)).nim_add(&fin(1)); // ω² + 1
        let w2 = Ordinal::omega_pow(fin(2)); // ω²
        assert_eq!(
            chi7_pow(8),
            Ordinal::omega_pow(e_w2_1).nim_add(&Ordinal::omega_pow(w2))
        );

        // χ_7^⊗9 = α_7 ⊗ χ_7^⊗2 = (ω+1)⊗ω^(ω²·2) = ω^(ω²·2+1) + ω^(ω²·2). Hand-verified
        // both ways (= χ_7^⊗8 ⊗ χ_7), so it also pins associativity through the carry.
        let w2_2 = Ordinal::monomial(fin(2), 2); // ω²·2  (exponent)
        let w2_2_1 = w2_2.nim_add(&fin(1)); // ω²·2 + 1
        assert_eq!(
            chi7_pow(9),
            Ordinal::omega_pow(w2_2_1).nim_add(&Ordinal::omega_pow(w2_2))
        );
        assert_eq!(chi7_pow(9), mul(&chi7_pow(8), &chi7()).unwrap());
    }

    #[test]
    fn locally_verified_alpha_47_landmark() {
        // `experiments/ordinal_excess_probe.py` independently verifies Lenstra excess
        // m_47 = 1 by a fixed-base finite-field power test. Since f(47)=23, this gives
        // α_47 = κ_23 + 1 = ω^(ω^7) + 1.
        let chi47 = Ordinal::omega_pow(Ordinal::omega_pow(fin(13)));
        let mut pow = fin(1);
        for _ in 0..47 {
            pow = mul(&pow, &chi47).unwrap();
        }
        assert_eq!(
            pow,
            Ordinal::omega_pow(Ordinal::omega_pow(fin(7))).nim_add(&fin(1))
        );
    }

    #[test]
    fn quintic_stage_field_axioms() {
        // The prime-3/prime-5 (scalar-α) commutative-ring sweep on a sample of ordinals
        // < ω^(ω²) spanning both the place-ω⁰ and place-ω¹ towers, coeffs in F_4. Every
        // product is defined here; associativity is what a digit-carry bug would break.
        let w = Ordinal::omega();
        let wn = |n| Ordinal::monomial(fin(1), n); // ω·n = ω^1·n (finite exponent 1)
        let elems: Vec<Ordinal> = vec![
            fin(1),
            fin(2),
            fin(3),                                                         // F_4 scalars
            w.clone(),                                                      // ω = χ_3
            Ordinal::omega_pow(fin(2)),                                     // ω²
            Ordinal::omega_pow(fin(3)),                                     // ω³ = χ_9 (g_1)
            ww(),                                                           // ω^ω = χ_5
            Ordinal::omega_pow(wn(2)),                                      // ω^(ω·2)
            Ordinal::omega_pow(wn(5)),                                      // ω^(ω·5) = χ_25
            Ordinal::omega_pow(w.nim_add(&fin(1))),                         // ω^(ω+1)
            ww().nim_add(&w).nim_add(&fin(1)),                              // ω^ω + ω + 1
            wn(3).nim_add(&fin(2)),                                         // ω·3 + 2
            Ordinal::omega_pow(wn(2)).nim_add(&Ordinal::omega_pow(fin(3))), // ω^(ω·2)+ω³
        ];
        check_field_axioms(&elems);
    }

    #[test]
    fn septic_stage_field_axioms() {
        // The decisive Stage-2 check: the commutative-ring axioms on a sample built from
        // χ_7 = ω^(ω²) (the first non-scalar-α generator), its powers 1..6, ω (= χ_3,
        // which α_7 = ω+1 drags in), F_4 scalars, and mixed sums. Every pairwise product
        // stays within primes {3, 7} (well inside the verified range ⇒ all `Some`), and associativity /
        // distributivity through the α_7 branching is exactly what a mis-mixed carry
        // would break. The α_7 = ω+1 *value* is source-pinned in `septic_kummer_landmark`.
        let mut elems: Vec<Ordinal> = vec![fin(1), fin(2), fin(3), w(), w().nim_add(&fin(1))];
        for n in 1..=6u128 {
            elems.push(chi7_pow(n));
        }
        elems.push(chi7().nim_add(&w())); // χ_7 + ω
        elems.push(Ordinal::monomial(Ordinal::omega_pow(fin(2)), 2).nim_add(&fin(1))); // χ_7·2 + 1
        elems.push(chi7_pow(3).nim_add(&w()).nim_add(&fin(1))); // χ_7³ + ω + 1
        check_field_axioms(&elems);
    }

    /// Commutativity, identity, associativity, and distributivity over `⊕`, on every
    /// triple of a sample whose pairwise products are all defined.
    fn check_field_axioms(elems: &[Ordinal]) {
        let one = fin(1);
        for a in elems {
            for b in elems {
                let ab = a.nim_mul(b).expect("sample is closed under ⊗");
                assert_eq!(ab, b.nim_mul(a).unwrap(), "non-commutative");
                assert_eq!(a.nim_mul(&one).unwrap(), *a, "identity");
                for c in elems {
                    let l = ab.nim_mul(c).unwrap();
                    let r = a.nim_mul(&b.nim_mul(c).unwrap()).unwrap();
                    assert_eq!(l, r, "× not associative");
                    let l = a.nim_mul(&b.nim_add(c)).unwrap();
                    let r = ab.nim_add(&a.nim_mul(c).unwrap());
                    assert_eq!(l, r, "× not distributive over ⊕");
                }
            }
        }
    }

    #[test]
    fn excess_table_matches_oeis_a380496() {
        // Spot-check named landmark rows against OEIS A380496. The first 14 rows are
        // DiMuro Table 1 + m_47 — a cross-check that the OEIS import agrees with the old
        // hardcoded values — then higher landmarks. The full 126-row diff against the
        // vendored b-file lives in `excess_table_matches_vendored_b380496_in_full` below.
        for (u, m) in [
            (3, 0),
            (5, 0),
            (7, 1),
            (11, 1),
            (13, 0),
            (17, 0),
            (19, 4),
            (23, 1),
            (29, 0),
            (31, 1),
            (37, 0),
            (41, 1),
            (43, 1),
            (47, 1),
            (53, 1),
            (61, 0),
            (73, 1),
            (109, 0),
            (163, 4),
            (263, 1),
            (521, 0),
            (659, 0),
            (701, 0),
            (709, 1),
        ] {
            assert_eq!(finite_excess(u), Some(m), "m_{u}");
        }
        // Every in-range row is 0 or 1 except m_19 = m_163 = 4 (the f(p)=2·3^k column).
        for place in 0..126u128 {
            let u = place_prime(place);
            let m = finite_excess(u).expect("all 126 b-file rows are defined");
            assert!(
                m == 0 || m == 1 || ((u == 19 || u == 163) && m == 4),
                "m_{u} = {m} is outside the A380496 value set"
            );
        }
    }

    #[test]
    fn excess_table_matches_vendored_b380496_in_full() {
        // OEIS A380496 b-file (a380496.txt / "Lenstra excess of the n-th odd prime"),
        // fetched 2026-07-02 and vendored at `b380496.txt` beside this module. Each line
        // is `n a(n)`, 1-indexed (`a(1)` = the excess at the odd-prime place 0, i.e.
        // prime 3). Diffs all 126 rows against `EXCESS`, closing the gap the old
        // landmark-only test left (it never actually read the file).
        let b_file = include_str!("b380496.txt");
        let mut rows = 0u128;
        for line in b_file.lines() {
            let line = line.trim();
            if line.is_empty() {
                continue;
            }
            let mut parts = line.split_whitespace();
            let n: u128 = parts
                .next()
                .expect("row has an index")
                .parse()
                .expect("n is u128");
            let a_n: u128 = parts
                .next()
                .expect("row has a value")
                .parse()
                .expect("a(n) is u128");
            assert!(parts.next().is_none(), "unexpected extra column in row {n}");
            let place = n - 1;
            let u = place_prime(place);
            assert_eq!(finite_excess(u), Some(a_n), "OEIS A380496 a({n}), u={u}");
            rows += 1;
        }
        assert_eq!(rows, 126, "expected all 126 known A380496 b-file rows");
    }

    #[test]
    fn boundary_returns_none_past_prime_709() {
        // The verified table now reaches prime 709 (OEIS A380496, 126 b-file rows). The
        // old boundary case — a Kummer carry at place 14 (prime 53) — is now defined:
        // ω^(ω^14·50) = χ_53^⊗50; squaring drives the place-14 digit to 100 ≥ 53, owing
        // α_53 = κ_{f(53)}+1, which is now in range. End-to-end carry resolution, so this
        // pins the table extension and not just the lookup. (f(53)=52=4·13 keeps the
        // reconstruction cheap; large primes like 709 are *in* the table but their
        // q_set/finite-subfield reconstruction over the huge component field is costly —
        // we don't pay for that in a unit test, so the boundary side uses table lookups.)
        let at53 = Ordinal::omega_pow(Ordinal::monomial(fin(14), 50)); // ω^(ω^14·50)
        assert!(mul(&at53, &at53).is_some());

        // The boundary now sits at place 126 = prime 719, the first OEIS-unknown row.
        assert_eq!(place_prime(125), 709);
        assert_eq!(place_prime(126), 719);
        assert_eq!(finite_excess(709), Some(1)); // last known row (table lookup)
        assert_eq!(finite_excess(719), None);
        // A Kummer carry at place 126 (prime 719) owes the unverified α_719 ⇒ None.
        // ω^(ω^126·360) squared drives the place-126 digit to 720 ≥ 719.
        let at719 = Ordinal::omega_pow(Ordinal::monomial(fin(126), 360));
        assert_eq!(mul(&at719, &at719), None);

        // And anything ≥ ω^(ω^ω) (an infinite exponent place) is out of range outright.
        let w_ww = Ordinal::omega_pow(ww()); // ω^(ω^ω)
        assert_eq!(mul(&w_ww, &w()), None);
    }

    #[test]
    fn identity_preserves_large_valid_prime_digits() {
        // p(53)=257, so digit 256 is legal. This used to truncate through `u128`
        // storage and collapse the monomial to 1 even when multiplying by 1.
        assert_eq!(place_prime(53), 257);
        let exp = Ordinal::monomial(fin(53), 256);
        let x = Ordinal::omega_pow(exp);
        assert_eq!(mul(&x, &fin(1)), Some(x.clone()));
        assert_eq!(mul(&fin(1), &x), Some(x));
    }
}