ogdoad 1.0.0

//! The **Brown invariant** `β ∈ ℤ/8` of a `ℤ/4`-valued quadratic refinement — the
//! char-2 cell of the mod-8 spine (Bridge M).
//!
//! The mod-8 spine otherwise lives entirely on the char-0 side: the exact rational
//! signature, the genus oddity ([`genus_signature_mod8`](crate::forms::genus_signature_mod8)),
//! the Milgram Gauss-sum phase
//! ([`milgram_signature_mod8`](crate::forms::DiscriminantForm::milgram_signature_mod8)),
//! and the Weil `S` prefactor are four routes to `σ mod 8`. The char-2 side carries
//! only the `ℤ/2` Arf bit. The classical object filling the char-2 mod-8 cell is the
//! Brown invariant.
//!
//! ## The `ℤ/4`-quadratic category *(standard math)*
//!
//! A **`ℤ/4`-quadratic form** on an `F₂`-space `V` is a map `q : V → ℤ/4` with
//!
//! ```text
//! q(x + y) = q(x) + q(y) + 2·b(x,y),
//! ```
//!
//! where `b : V×V → F₂` is symmetric bilinear and `2· : F₂ ↪ ℤ/4`. Setting `y = x`
//! forces `b(x,x) = q(x) mod 2`, so `b` is symmetric **but not alternating**. For
//! `b` nondegenerate the Gauss sum is a `ℤ[i]`-integer of absolute value `2^{n/2}`,
//!
//! ```text
//! Σ_{x ∈ V} i^{q(x)} = 2^{n/2} · ζ₈^β,    ζ₈ = e^{2πi/8},
//! ```
//!
//! and `β ∈ ℤ/8` is the **Brown invariant** (E. H. Brown, *Generalizations of the
//! Kervaire invariant*, Ann. of Math. **95** (1972); C. T. C. Wall, *Quadratic
//! forms on finite groups*, Topology **2** (1963)): additive under `⊥`, and a
//! complete invariant up to split planes, making the Witt group of the category
//! cyclic of order 8 generated by `⟨1⟩` (`q(x) = 1`).
//!
//! ## Category trap (load-bearing)
//!
//! This `b` is **not** the engine's polar form. The crate's char-2
//! [`Metric`](crate::clifford::Metric) carries an **alternating** `b` (`b_ii = 0`)
//! with `q` valued in the field; Brown's category has `ℤ/4`-valued `q` with
//! `b_ii = q_i mod 2`. Hard rule 2 ("don't collapse `q` and `b`") has a cousin
//! here: don't identify the two categories. The doubling map [`double_f2`] (a
//! classical `F₂` form `q'` ↦ `2q' : V → ℤ/4`) is the only bridge between them, and
//! it lands the shipped Arf bit as the 2-torsion of `β`:
//!
//! ```text
//! β(2q') = 4 · Arf(q')   ∈ {0, 4} ⊂ ℤ/8.
//! ```
//!
//! ## Implementation: the reduction route
//!
//! Primary route is **orthogonal reduction**, with no floating point and no
//! `2^rank` public budget. The input mirrors [`arf_f2`](crate::forms::arf_f2)
//! field-for-field: `q4` (the values mod 4 on the basis) replaces the `F₂`
//! diagonal, and `bmat` carries the **off-diagonal** symmetric polar `b` (the
//! diagonal `b_ii = q4[i] mod 2` is forced, not input).
//!
//! Splitting off the radical of `b`, `q` restricted to `rad(b)` is **linear** into
//! `{0, 2}` (since `b = 0` there forces `2·q(x) ≡ 0 mod 4`), so `V = core ⊥ rad`
//! is an orthogonal sum on which `q` is additive and the Gauss sum **factors**:
//!
//! ```text
//! Σ_V i^{q(x)} = (Σ_core i^{q}) · (Σ_rad i^{q}),
//!   Σ_rad = 2^{radical_dim}  if q|rad ≡ 0,  else 0  (radical_anisotropic).
//! ```
//!
//! On the nonsingular core, the algorithm repeatedly splits off either an odd
//! line (`q(v) = 1` or `3`) or an even symplectic plane; the known block phases
//! (`1`, `7`, `0`, or `4`) add in `ℤ/8`. So `β` is reduced block-by-block and
//! reported alongside the radical data exactly as [`ArfInvariants`](crate::forms::ArfInvariants)
//! reports its radical. The old direct enumeration route remains only as a test
//! oracle. The lattice tie ([`DiscriminantForm::brown_invariant`]) and the third
//! identification (`β ≡ sign(L) mod 8` on 2-elementary discriminant forms) live in
//! `integral/discriminant.rs`; this module is the pure-`F₂` core.
//!
//! **Claim level:** standard math (Brown 1972; Wall 1963) made computational; the
//! bridge is the wiring to Arf (shipped) and Milgram (Bridge A), no new theorem.

/// The Brown invariant of a `ℤ/4`-quadratic form, mirroring
/// [`ArfInvariants`](crate::forms::ArfInvariants) field-for-field.
#[derive(Debug, Clone, PartialEq, Eq)]
pub struct BrownInvariants {
    /// Brown invariant `β ∈ ℤ/8` of the **nonsingular core** (the part where the
    /// symmetric polar `b` is nondegenerate).
    pub beta: u128,
    /// Rank of the nonsingular part of `b` (`= n − radical_dim`).
    pub rank: usize,
    /// Dimension of the polar-form radical (where `b` vanishes).
    pub radical_dim: usize,
    /// Whether `q` takes the value `2` somewhere on the radical (a "defective"
    /// direction). When true the *full* Gauss sum vanishes; `beta` still reports
    /// the core. Data, not a panic — exactly as `ArfInvariants::radical_anisotropic`.
    pub radical_anisotropic: bool,
}

impl BrownInvariants {
    /// `display()` alias kept for Python callers.
    pub fn display(&self) -> String {
        self.to_string()
    }
}

impl std::fmt::Display for BrownInvariants {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        let rad = if self.radical_dim > 0 {
            let aniso = if self.radical_anisotropic {
                " aniso"
            } else {
                ""
            };
            format!(" ⊗̂ Λ({}^{}{})", "F", self.radical_dim, aniso)
        } else {
            String::new()
        };
        write!(f, "β={} rank={}{}", self.beta, self.rank, rad)
    }
}

/// Bits of a mask strictly above position `i`.
fn above(i: usize) -> u128 {
    if i >= 127 {
        0
    } else {
        (!0u128) << (i + 1)
    }
}

/// `q(x) ∈ ℤ/4` for a bitmask vector `x`, from the basis values `q4` and the
/// **off-diagonal** symmetric polar `bmat`:
/// `q(x) = Σ_{i∈x} q4[i] + 2·#{ i<j∈x : b_{ij}=1 }  (mod 4)`.
fn q_of4(x: u128, q4: &[u128], bmat: &[u128]) -> u128 {
    let mut lin = 0u128;
    let mut cross = 0u32; // parity of crossing pairs
    let mut vv = x;
    while vv != 0 {
        let i = vv.trailing_zeros() as usize;
        vv &= vv - 1;
        lin += q4[i] % 4;
        let inter = bmat[i] & x & above(i);
        cross ^= inter.count_ones() & 1;
    }
    (lin + 2 * cross as u128) % 4
}

/// `b(x,y)` for the full symmetric polar matrix, diagonal included.
fn b_pair(x: u128, y: u128, full_b: &[u128]) -> bool {
    let mut parity = 0u32;
    let mut xx = x;
    while xx != 0 {
        let i = xx.trailing_zeros() as usize;
        xx &= xx - 1;
        parity ^= (full_b[i] & y).count_ones() & 1;
    }
    parity == 1
}

/// Insert `v` into an `F₂` XOR-basis keyed by lowest set bit. Returns `true` iff
/// `v` was independent of the basis so far (and was added).
fn xor_insert(basis: &mut [Option<u128>; 128], mut v: u128) -> bool {
    while v != 0 {
        let c = v.trailing_zeros() as usize;
        match basis[c] {
            None => {
                basis[c] = Some(v);
                return true;
            }
            Some(b) => v ^= b,
        }
    }
    false
}

/// A basis of `{ x : Σ_{i∈x} full_b[i] = 0 }` — the linear dependencies among the
/// rows of the **full** symmetric polar matrix (diagonal included), i.e. the
/// radical `{ x : b(x, ·) = 0 }`. Returned as coefficient bitmasks over the basis.
fn radical_basis(full_b: &[u128], n: usize) -> Vec<u128> {
    let mut piv: [Option<(u128, u128)>; 128] = [None; 128]; // col -> (reduced row, provenance)
    let mut null = Vec::new();
    for i in 0..n {
        let mut row = full_b[i];
        let mut prov = 1u128 << i;
        loop {
            if row == 0 {
                null.push(prov);
                break;
            }
            let c = row.trailing_zeros() as usize;
            match piv[c] {
                None => {
                    piv[c] = Some((row, prov));
                    break;
                }
                Some((r, p)) => {
                    row ^= r;
                    prov ^= p;
                }
            }
        }
    }
    null
}

/// A coordinate-axis complement to the radical. Since the radical is orthogonal to
/// all of `V`, any linear complement carries the nonsingular core.
fn core_complement_basis(radical: &[u128], n: usize) -> Vec<u128> {
    let mut lin: [Option<u128>; 128] = [None; 128];
    for &x in radical {
        xor_insert(&mut lin, x);
    }
    let mut basis = Vec::new();
    for i in 0..n {
        if xor_insert(&mut lin, 1u128 << i) {
            basis.push(1u128 << i);
        }
    }
    basis
}

/// Reduce a nonsingular Brown core into odd lines and even planes, adding the
/// block Brown phases in `ℤ/8`.
fn reduce_brown_core(mut basis: Vec<u128>, q4: &[u128], bmat: &[u128], full_b: &[u128]) -> u128 {
    let mut beta = 0u128;
    while !basis.is_empty() {
        if let Some(p) = basis.iter().position(|&v| b_pair(v, v, full_b)) {
            let v = basis.swap_remove(p);
            let qv = q_of4(v, q4, bmat);
            debug_assert!(qv == 1 || qv == 3);
            beta = (beta + if qv == 1 { 1 } else { 7 }) % 8;
            for w in &mut basis {
                if b_pair(*w, v, full_b) {
                    *w ^= v;
                }
            }
            continue;
        }

        let v = basis
            .pop()
            .expect("nonempty basis already checked before even-plane reduction");
        let p = basis
            .iter()
            .position(|&w| b_pair(v, w, full_b))
            .expect("a nonsingular alternating core has a symplectic partner");
        let w = basis.swap_remove(p);
        let qv = q_of4(v, q4, bmat);
        let qw = q_of4(w, q4, bmat);
        debug_assert!(qv == 0 || qv == 2);
        debug_assert!(qw == 0 || qw == 2);
        if qv == 2 && qw == 2 {
            beta = (beta + 4) % 8;
        }

        for x in &mut basis {
            let xv = b_pair(*x, v, full_b);
            let xw = b_pair(*x, w, full_b);
            if xw {
                *x ^= v;
            }
            if xv {
                *x ^= w;
            }
        }
    }
    beta
}

/// Recover `β ∈ ℤ/8` from the **exact** core Gauss sum `G = re + im·i`, where
/// `G = 2^{rank/2}·ζ₈^β`. The eight cases split into axis values (`β` even, one of
/// `re`/`im` zero) and diagonal values (`β` odd, `|re| = |im|`), so `β` is read off
/// by integer sign and magnitude comparison alone — no floating point. `None` only
/// if `G` is not of eighth-root shape (a singular core whose sum vanished, or
/// malformed input).
pub(crate) fn beta_from_gauss(re: i128, im: i128) -> Option<u128> {
    match (re.signum(), im.signum()) {
        (0, 0) => None, // vanished — caller must enumerate the core, not the whole space
        (s, 0) => Some(if s > 0 { 0 } else { 4 }),
        (0, s) => Some(if s > 0 { 2 } else { 6 }),
        (sr, si) => {
            if re.unsigned_abs() != im.unsigned_abs() {
                return None;
            }
            Some(match (sr, si) {
                (1, 1) => 1,
                (-1, 1) => 3,
                (-1, -1) => 5,
                (1, -1) => 7,
                _ => unreachable!("nonzero signums already matched"),
            })
        }
    }
}

/// The Brown invariant `β ∈ ℤ/8` of the `ℤ/4`-quadratic form given by basis values
/// `q4` (mod 4) and the **off-diagonal** symmetric polar `bmat` (`bmat[i]` bit `j`
/// ⇔ `b_{ij}=1`, `i ≠ j`; the diagonal `b_ii = q4[i] mod 2` is derived, so any
/// diagonal bits in `bmat` are ignored). Mirrors [`arf_f2`](crate::forms::arf_f2).
///
/// Reduction route: split off `rad(b)`, then reduce the nonsingular core into odd
/// lines and even planes, adding their known Brown phases in `ℤ/8`.
pub fn brown_f2(n: usize, q4: &[u128], bmat: &[u128]) -> BrownInvariants {
    assert!(
        n <= 128,
        "brown_f2 uses u128 bitmasks, so n must be at most 128"
    );
    assert!(
        q4.len() >= n && bmat.len() >= n,
        "brown_f2 needs q4 and bmat entries for every basis vector"
    );

    // Full symmetric polar including the forced diagonal b_ii = q4[i] mod 2.
    let full_b: Vec<u128> = (0..n)
        .map(|i| {
            let off = bmat[i] & !(1u128 << i);
            off | (if q4[i] % 2 == 1 { 1u128 << i } else { 0 })
        })
        .collect();

    let radical = radical_basis(&full_b, n);
    let radical_dim = radical.len();
    // q|rad is linear into {0,2}; anisotropic iff some radical basis vector has q=2.
    let radical_anisotropic = radical.iter().any(|&x| q_of4(x, q4, bmat) == 2);

    let core = core_complement_basis(&radical, n);
    let rank = core.len(); // = n − radical_dim
    let beta = reduce_brown_core(core, q4, bmat, &full_b);
    BrownInvariants {
        beta,
        rank,
        radical_dim,
        radical_anisotropic,
    }
}

/// The **doubling bridge**: a classical nonsingular `F₂` quadratic form `q'`
/// (alternating polar, given as `arf_f2` data) maps to `2q' : V → ℤ/4`, whose Brown
/// invariant is `β(2q') = 4·Arf(q')`. Builds `q4 = 2·q'` (values in `{0,2}`, so the
/// derived `ℤ/4` polar is the original alternating `bmat`) and runs [`brown_f2`].
pub fn double_f2(qd: &[bool], bmat: &[u128]) -> BrownInvariants {
    let q4: Vec<u128> = qd.iter().map(|&b| if b { 2 } else { 0 }).collect();
    brown_f2(qd.len(), &q4, bmat)
}

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

    /// The old exact Gauss-sum route, retained as an oracle for the reduction path.
    fn brown_f2_by_enumeration(n: usize, q4: &[u128], bmat: &[u128]) -> BrownInvariants {
        let full_b: Vec<u128> = (0..n)
            .map(|i| {
                let off = bmat[i] & !(1u128 << i);
                off | (if q4[i] % 2 == 1 { 1u128 << i } else { 0 })
            })
            .collect();

        let radical = radical_basis(&full_b, n);
        let radical_dim = radical.len();
        let radical_anisotropic = radical.iter().any(|&x| q_of4(x, q4, bmat) == 2);
        let core = core_complement_basis(&radical, n);
        let rank = core.len();
        assert!(
            rank < 128,
            "test-only enumeration needs a representable 2^rank loop bound"
        );

        let mut counts = [0i128; 4];
        for mask in 0u128..(1u128 << rank) {
            let mut x = 0u128;
            for (b, &v) in core.iter().enumerate() {
                if (mask >> b) & 1 == 1 {
                    x ^= v;
                }
            }
            counts[q_of4(x, q4, bmat) as usize] += 1;
        }

        let re = counts[0] - counts[2];
        let im = counts[1] - counts[3];
        BrownInvariants {
            beta: beta_from_gauss(re, im).expect("a nonsingular core has an eighth-root Gauss sum"),
            rank,
            radical_dim,
            radical_anisotropic,
        }
    }

    /// Block-diagonal `⊥` of two Brown forms (disjoint generators, no cross polar).
    fn ortho_sum(
        (n1, q1, b1): (usize, &[u128], &[u128]),
        (n2, q2, b2): (usize, &[u128], &[u128]),
    ) -> (usize, Vec<u128>, Vec<u128>) {
        let n = n1 + n2;
        let mut q4 = q1.to_vec();
        q4.extend_from_slice(q2);
        let mut bmat = vec![0u128; n];
        bmat[..n1].copy_from_slice(b1);
        for i in 0..n2 {
            bmat[n1 + i] = b2[i] << n1; // shift block-2 polar into the high coordinates
        }
        (n, q4, bmat)
    }

    // --- the ℤ/8 generator and its inverse ---

    #[test]
    fn one_dimensional_generators() {
        // ⟨1⟩: q(e)=1 ⇒ β=1 (the ℤ/8 generator).  ⟨−1⟩: q(e)=3 ⇒ β=7.
        assert_eq!(brown_f2(1, &[1], &[0]).beta, 1);
        assert_eq!(brown_f2(1, &[3], &[0]).beta, 7);
    }

    #[test]
    fn order_eight_relation() {
        // ⟨1⟩^{⊥8}: β = 8 ≡ 0 (the Witt group is cyclic of order 8).
        let r = brown_f2(8, &[1; 8], &[0; 8]);
        assert_eq!(r.beta, 0);
        assert_eq!((r.rank, r.radical_dim), (8, 0));
    }

    // --- the split objects (β = 0) ---

    #[test]
    fn split_objects_vanish() {
        // ⟨1⟩ ⊥ ⟨−1⟩: β = 1 + 7 = 0.
        assert_eq!(brown_f2(2, &[1, 3], &[0, 0]).beta, 0);
        // the even hyperbolic plane [q(e)=0, q(f)=0, b(e,f)=1]: β = 0.
        assert_eq!(brown_f2(2, &[0, 0], &[0b10, 0b01]).beta, 0);
    }

    // --- the doubling bridge: β(2q′) = 4·Arf(q′) ---

    #[test]
    #[allow(clippy::type_complexity)] // compact table of (q-diagonal, polar pairs)
    fn double_is_four_times_arf() {
        let cases: &[(&[bool], &[(usize, usize)])] = &[
            (&[false, false], &[(0, 1)]),                       // hyperbolic, Arf 0
            (&[true, true], &[(0, 1)]),                         // anisotropic, Arf 1
            (&[false, false, false, false], &[(0, 1), (2, 3)]), // H⊕H, Arf 0
            (&[false, false, true, true], &[(0, 1), (2, 3)]),   // H⊕A, Arf 1
            (&[true, true, true, true], &[(0, 1), (2, 3)]),     // A⊕A, Arf 0
        ];
        for (qd, pairs) in cases {
            let n = qd.len();
            let mut bmat = vec![0u128; n];
            for &(i, j) in *pairs {
                bmat[i] |= 1 << j;
                bmat[j] |= 1 << i;
            }
            let arf = arf_f2(n, qd, &bmat).arf;
            let beta = double_f2(qd, &bmat).beta;
            assert_eq!(beta, 4 * arf, "β(2q′) ≠ 4·Arf for q={qd:?}");
            // doubled forms land in the 2-torsion {0,4} ⊂ ℤ/8.
            assert!(beta == 0 || beta == 4);
        }
    }

    // --- additivity under ⊥ (the defining group law) ---

    #[test]
    fn beta_is_additive_under_orthogonal_sum() {
        // β(A ⊥ B) = β(A) + β(B) mod 8, across a spread of components.
        let comps: &[(usize, &[u128], &[u128])] = &[
            (1, &[1], &[0]),             // β 1
            (1, &[3], &[0]),             // β 7
            (1, &[2], &[0]),             // ⟨2⟩: q(e)=2
            (2, &[0, 0], &[0b10, 0b01]), // even plane, β 0
            (2, &[1, 1], &[0b10, 0b01]), // a genuine ℤ/4 plane
        ];
        for a in comps {
            for b in comps {
                let ba = brown_f2(a.0, a.1, a.2).beta;
                let bb = brown_f2(b.0, b.1, b.2).beta;
                let (n, q4, bmat) = ortho_sum(*a, *b);
                let bab = brown_f2(n, &q4, &bmat).beta;
                assert_eq!(bab, (ba + bb) % 8, "additivity failed for {a:?} ⊥ {b:?}");
            }
        }
    }

    // --- the radical: q=2 direction kills the full sum, β reports the core ---

    #[test]
    fn anisotropic_radical_is_detected() {
        // [q(e)=0, q(f)=0, b(e,f)=1] ⊥ ⟨q(r)=2⟩: core is the even plane (β=0), the
        // radical direction r has q=2 (anisotropic). The full Gauss sum vanishes;
        // brown_f2 reports the core β and flags the radical.
        let r = brown_f2(3, &[0, 0, 2], &[0b10, 0b01, 0]);
        assert_eq!(
            (r.beta, r.rank, r.radical_dim, r.radical_anisotropic),
            (0, 2, 1, true)
        );
        // a benign (q=0) radical direction: same core, not anisotropic.
        let r0 = brown_f2(3, &[0, 0, 0], &[0b10, 0b01, 0]);
        assert_eq!(
            (r0.beta, r0.rank, r0.radical_dim, r0.radical_anisotropic),
            (0, 2, 1, false)
        );
    }

    #[test]
    fn reduction_matches_enumeration_on_all_four_dimensional_inputs() {
        for qmask in 0u128..(1u128 << 8) {
            let q4: Vec<u128> = (0..4).map(|i| (qmask >> (2 * i)) & 0b11).collect();
            for edges in 0u128..(1u128 << 6) {
                let mut bmat = vec![0u128; 4];
                let mut bit = 0;
                for i in 0..4 {
                    for j in (i + 1)..4 {
                        if (edges >> bit) & 1 == 1 {
                            bmat[i] |= 1u128 << j;
                            bmat[j] |= 1u128 << i;
                        }
                        bit += 1;
                    }
                }
                assert_eq!(
                    brown_f2(4, &q4, &bmat),
                    brown_f2_by_enumeration(4, &q4, &bmat),
                    "reduction/enumeration mismatch for q={q4:?}, b={bmat:?}"
                );
            }
        }
    }

    #[test]
    fn reduction_matches_old_budget_edge() {
        // Twenty-six odd lines hit the old public enumeration ceiling exactly.
        let q4 = vec![1u128; 26];
        let bmat = vec![0u128; 26];
        assert_eq!(
            brown_f2(26, &q4, &bmat),
            brown_f2_by_enumeration(26, &q4, &bmat)
        );
    }

    #[test]
    fn brown_f2_reduces_past_the_old_enumeration_budget() {
        // Forty odd lines would have panicked under the old rank <= 26 budget.
        let q4 = vec![1u128; 40];
        let bmat = vec![0u128; 40];
        let r = brown_f2(40, &q4, &bmat);
        assert_eq!(
            (r.beta, r.rank, r.radical_dim, r.radical_anisotropic),
            (0, 40, 0, false)
        );
    }
}