geoit 0.0.2

//! Expr compilation: compile governance expressions into flat polynomial
//! evaluation plans for the rendering pipeline.
//!
//! The key insight: if the object Mv is fixed (constructed at governance time)
//! and the probe is the only dynamic input, then field evaluation reduces to
//! a polynomial in probe coordinates. This polynomial can be pre-computed and
//! evaluated as a flat sequence of integer multiply-accumulate operations.
//!
//! Pipeline: Expr → PartialExpr (separate static/dynamic) → IntPoly → EvalPlan

use crate::algebra::blade_new::BladeMask;
use crate::algebra::mv::Mv;
use crate::algebra::signature::Signature;
use crate::governance::field::FieldOp;
use std::collections::BTreeMap;

/// Integer polynomial in probe parameters.
/// Coefficients are i64 (pre-scaled from exact Rat).
#[derive(Clone, Debug, PartialEq)]
pub struct IntPoly {
    /// Monomial → coefficient. Key: exponent vector (one per probe param).
    pub terms: Vec<(Vec<u8>, i64)>,
    /// Number of probe parameters.
    pub arity: usize,
}

impl IntPoly {
    pub fn zero(arity: usize) -> Self {
        IntPoly {
            terms: Vec::new(),
            arity,
        }
    }

    pub fn constant(value: i64, arity: usize) -> Self {
        if value == 0 {
            return Self::zero(arity);
        }
        IntPoly {
            terms: vec![(vec![0; arity], value)],
            arity,
        }
    }

    pub fn variable(index: usize, arity: usize) -> Self {
        let mut exp = vec![0u8; arity];
        exp[index] = 1;
        IntPoly {
            terms: vec![(exp, 1)],
            arity,
        }
    }

    pub fn is_zero(&self) -> bool {
        self.terms.is_empty()
    }

    pub fn add(&self, other: &Self) -> Self {
        let mut result = BTreeMap::new();
        for (exp, c) in &self.terms {
            *result.entry(exp.clone()).or_insert(0i64) += c;
        }
        for (exp, c) in &other.terms {
            *result.entry(exp.clone()).or_insert(0i64) += c;
        }
        let terms: Vec<_> = result.into_iter().filter(|(_, c)| *c != 0).collect();
        IntPoly {
            terms,
            arity: self.arity,
        }
    }

    pub fn scale(&self, s: i64) -> Self {
        if s == 0 {
            return Self::zero(self.arity);
        }
        let terms: Vec<_> = self
            .terms
            .iter()
            .map(|(exp, c)| (exp.clone(), c * s))
            .filter(|(_, c)| *c != 0)
            .collect();
        IntPoly {
            terms,
            arity: self.arity,
        }
    }

    pub fn mul(&self, other: &Self) -> Self {
        let mut result: BTreeMap<Vec<u8>, i64> = BTreeMap::new();
        for (ea, ca) in &self.terms {
            for (eb, cb) in &other.terms {
                let exp: Vec<u8> = ea.iter().zip(eb).map(|(&a, &b)| a + b).collect();
                *result.entry(exp).or_insert(0) += ca * cb;
            }
        }
        let terms: Vec<_> = result.into_iter().filter(|(_, c)| *c != 0).collect();
        IntPoly {
            terms,
            arity: self.arity,
        }
    }

    pub fn neg(&self) -> Self {
        self.scale(-1)
    }

    /// Maximum total degree across all terms.
    pub fn degree(&self) -> u16 {
        self.terms
            .iter()
            .map(|(exp, _)| exp.iter().map(|&e| e as u16).sum::<u16>())
            .max()
            .unwrap_or(0)
    }

    /// Evaluate at integer coordinates.
    pub fn eval(&self, coords: &[i64]) -> i64 {
        let mut sum = 0i64;
        for (exp, coeff) in &self.terms {
            let mut term: i64 = *coeff;
            for (i, &e) in exp.iter().enumerate() {
                for _ in 0..e {
                    term = term.wrapping_mul(coords[i]);
                }
            }
            sum = sum.wrapping_add(term);
        }
        sum
    }
}

/// A compiled field evaluation: for each output blade, a polynomial
/// in probe parameters with i64 coefficients.
#[derive(Clone, Debug)]
pub struct CompiledFieldEval {
    /// For each output component, a polynomial in probe parameters.
    pub polys: Vec<(BladeMask, IntPoly)>,
    /// Common denominator (scale factor from Rat → i64 conversion).
    pub denominator: i64,
    /// Number of probe parameters (2 for 2D, 3 for 3D).
    pub probe_arity: usize,
}

impl CompiledFieldEval {
    /// Evaluate the compiled field at integer probe coordinates.
    /// Returns the scalar result (for IPNS field evaluation).
    pub fn eval_scalar(&self, coords: &[i64]) -> i64 {
        self.polys
            .iter()
            .filter(|(mask, _)| *mask == 0) // scalar blade
            .map(|(_, poly)| poly.eval(coords))
            .sum()
    }

    /// Evaluate all components.
    pub fn eval_all(&self, coords: &[i64]) -> Vec<(BladeMask, i64)> {
        self.polys
            .iter()
            .map(|(mask, poly)| (*mask, poly.eval(coords)))
            .collect()
    }
}

/// Compile a field evaluation for a given object Mv and field operation.
///
/// The object is known at compile time; the probe coordinates are the only dynamic input.
/// The result is a set of polynomials in the probe coordinates.
pub fn compile_field_eval(
    object: &Mv,
    field_op: &FieldOp,
    sig: &Signature,
    probe_construction_arity: usize,
    probe_construction_gen_start: u8,
) -> CompiledFieldEval {
    let arity = probe_construction_arity;

    // Build the probe Mv symbolically: probe = Σ x_k · gen(probe_gen_start + k)
    // Each coefficient of the probe is a polynomial (degree 1 variable).
    let mut probe_polys: Vec<(BladeMask, IntPoly)> = Vec::new();
    for k in 0..arity {
        let gen_idx = probe_construction_gen_start + k as u8;
        let mask = 1u64 << gen_idx;
        probe_polys.push((mask, IntPoly::variable(k, arity)));
    }

    // Apply the field operation symbolically
    match field_op {
        FieldOp::ScalarProduct => {
            // ⟨probe, object⟩ = Σ_{matching blades} probe_coeff * object_coeff
            // Since probe is grade-1 and coefficients are linear in x_k,
            // the scalar product is linear in the probe variables.
            compile_scalar_product(&probe_polys, object, sig, arity)
        }
        FieldOp::InnerProduct => compile_inner_product(&probe_polys, object, sig, arity),
        _ => {
            // Fallback: empty plan (field op not yet compiled)
            CompiledFieldEval {
                polys: vec![],
                denominator: 1,
                probe_arity: arity,
            }
        }
    }
}

/// Compile scalar product ⟨probe, object⟩.
fn compile_scalar_product(
    probe_polys: &[(BladeMask, IntPoly)],
    object: &Mv,
    sig: &Signature,
    arity: usize,
) -> CompiledFieldEval {
    use crate::algebra::product_new::blade_product;

    // Find LCM of object coefficient denominators for integer scaling
    let mut lcm_den = 1i64;
    for (_, coeff) in object.blades() {
        if let Some(r) = coeff.try_as_rat() {
            let d = r.den() as i64;
            lcm_den = lcm(lcm_den, d);
        }
    }

    let mut result_poly = IntPoly::zero(arity);

    for (probe_mask, probe_poly) in probe_polys {
        for (obj_mask, obj_coeff) in object.blades() {
            let (result_mask, sign) = blade_product(*probe_mask, obj_mask, sig);
            if sign == 0 || result_mask != 0 {
                continue;
            } // only keep scalar part

            // Object coeff as scaled integer
            if let Some(r) = obj_coeff.try_as_rat() {
                let scaled = (r.num() * (lcm_den / r.den() as i64) as i128) as i64;
                let term = probe_poly.scale(scaled * sign as i64);
                result_poly = result_poly.add(&term);
            }
        }
    }

    CompiledFieldEval {
        polys: vec![(0, result_poly)],
        denominator: lcm_den,
        probe_arity: arity,
    }
}

/// Compile inner product ⟨probe, object⟩ (Hestenes).
fn compile_inner_product(
    probe_polys: &[(BladeMask, IntPoly)],
    object: &Mv,
    sig: &Signature,
    arity: usize,
) -> CompiledFieldEval {
    use crate::algebra::blade_new::grade;
    use crate::algebra::product_new::blade_product;

    let mut lcm_den = 1i64;
    for (_, coeff) in object.blades() {
        if let Some(r) = coeff.try_as_rat() {
            lcm_den = lcm(lcm_den, r.den() as i64);
        }
    }

    let mut result_polys: BTreeMap<BladeMask, IntPoly> = BTreeMap::new();

    for (probe_mask, probe_poly) in probe_polys {
        let ga = grade(*probe_mask);
        for (obj_mask, obj_coeff) in object.blades() {
            let gb = grade(obj_mask);
            if ga == 0 || gb == 0 {
                continue;
            }
            let (result_mask, sign) = blade_product(*probe_mask, obj_mask, sig);
            if sign == 0 {
                continue;
            }
            let gr = grade(result_mask);
            if gr != ga.abs_diff(gb) {
                continue;
            }

            if let Some(r) = obj_coeff.try_as_rat() {
                let scaled = (r.num() * (lcm_den / r.den() as i64) as i128) as i64;
                let term = probe_poly.scale(scaled * sign as i64);
                let entry = result_polys
                    .entry(result_mask)
                    .or_insert_with(|| IntPoly::zero(arity));
                *entry = entry.add(&term);
            }
        }
    }

    let polys: Vec<_> = result_polys
        .into_iter()
        .filter(|(_, p)| !p.is_zero())
        .collect();

    CompiledFieldEval {
        polys,
        denominator: lcm_den,
        probe_arity: arity,
    }
}

fn gcd_i64(mut a: i64, mut b: i64) -> i64 {
    a = a.abs();
    b = b.abs();
    while b != 0 {
        let t = b;
        b = a % b;
        a = t;
    }
    a
}

fn lcm(a: i64, b: i64) -> i64 {
    if a == 0 || b == 0 {
        return 0;
    }
    (a / gcd_i64(a, b)) * b
}

// ═══════════════════════════════════════════════════════════
// B1: PARTIAL MV — Mv with polynomial coefficients
// ═══════════════════════════════════════════════════════════

/// A multivector whose blade coefficients are integer polynomials in probe parameters.
/// This is the intermediate representation between Expr evaluation and field compilation.
#[derive(Clone, Debug)]
pub struct PartialMv {
    /// (blade_mask, polynomial_coefficient)
    pub terms: Vec<(BladeMask, IntPoly)>,
    pub arity: usize,
}

impl PartialMv {
    pub fn zero(arity: usize) -> Self {
        PartialMv {
            terms: Vec::new(),
            arity,
        }
    }

    /// Constant scalar PartialMv.
    pub fn scalar(value: i64, arity: usize) -> Self {
        if value == 0 {
            return Self::zero(arity);
        }
        PartialMv {
            terms: vec![(0, IntPoly::constant(value, arity))],
            arity,
        }
    }

    /// Single basis blade with constant coefficient.
    pub fn blade(mask: BladeMask, coeff: i64, arity: usize) -> Self {
        if coeff == 0 {
            return Self::zero(arity);
        }
        PartialMv {
            terms: vec![(mask, IntPoly::constant(coeff, arity))],
            arity,
        }
    }

    /// Single basis blade with polynomial coefficient.
    pub fn blade_poly(mask: BladeMask, poly: IntPoly) -> Self {
        let arity = poly.arity;
        if poly.is_zero() {
            return Self::zero(arity);
        }
        PartialMv {
            terms: vec![(mask, poly)],
            arity,
        }
    }

    /// Embed a concrete Mv as a PartialMv with constant polynomial coefficients.
    /// Returns (partial_mv, denominator) where all coefficients are scaled to integers.
    pub fn from_mv(mv: &Mv, arity: usize) -> (Self, i64) {
        let mut lcm_den = 1i64;
        for (_, coeff) in mv.blades() {
            if let Some(r) = coeff.try_as_rat() {
                lcm_den = lcm(lcm_den, r.den() as i64);
            }
        }
        let terms: Vec<(BladeMask, IntPoly)> = mv
            .blades()
            .filter_map(|(mask, coeff)| {
                if coeff.is_zero() {
                    return None;
                }
                if let Some(r) = coeff.try_as_rat() {
                    let scaled = (r.num() * (lcm_den / r.den() as i64) as i128) as i64;
                    Some((mask, IntPoly::constant(scaled, arity)))
                } else {
                    None
                }
            })
            .collect();
        (PartialMv { terms, arity }, lcm_den)
    }

    /// Add two PartialMvs.
    pub fn add(&self, other: &Self) -> Self {
        let mut result: BTreeMap<BladeMask, IntPoly> = BTreeMap::new();
        for (mask, poly) in &self.terms {
            let entry = result
                .entry(*mask)
                .or_insert_with(|| IntPoly::zero(self.arity));
            *entry = entry.add(poly);
        }
        for (mask, poly) in &other.terms {
            let entry = result
                .entry(*mask)
                .or_insert_with(|| IntPoly::zero(self.arity));
            *entry = entry.add(poly);
        }
        let terms: Vec<_> = result.into_iter().filter(|(_, p)| !p.is_zero()).collect();
        PartialMv {
            terms,
            arity: self.arity,
        }
    }

    /// Scale all coefficients by a constant.
    pub fn scale(&self, s: i64) -> Self {
        if s == 0 {
            return Self::zero(self.arity);
        }
        let terms = self
            .terms
            .iter()
            .map(|(mask, poly)| (*mask, poly.scale(s)))
            .collect();
        PartialMv {
            terms,
            arity: self.arity,
        }
    }

    /// Negate.
    pub fn neg(&self) -> Self {
        self.scale(-1)
    }

    /// Geometric product of two PartialMvs.
    pub fn mul(&self, other: &Self, sig: &Signature) -> Self {
        use crate::algebra::product_new::blade_product;
        let mut result: BTreeMap<BladeMask, IntPoly> = BTreeMap::new();
        for (ma, pa) in &self.terms {
            for (mb, pb) in &other.terms {
                let (result_mask, sign) = blade_product(*ma, *mb, sig);
                if sign == 0 {
                    continue;
                }
                let product = pa.mul(pb).scale(sign as i64);
                let entry = result
                    .entry(result_mask)
                    .or_insert_with(|| IntPoly::zero(self.arity));
                *entry = entry.add(&product);
            }
        }
        let terms: Vec<_> = result.into_iter().filter(|(_, p)| !p.is_zero()).collect();
        PartialMv {
            terms,
            arity: self.arity,
        }
    }

    /// Extract the scalar (grade 0) component polynomial.
    pub fn scalar_part(&self) -> IntPoly {
        self.terms
            .iter()
            .find(|(mask, _)| *mask == 0)
            .map(|(_, p)| p.clone())
            .unwrap_or_else(|| IntPoly::zero(self.arity))
    }

    /// Keep only blades at a specific grade.
    pub fn grade_project(&self, g: u8) -> Self {
        use crate::algebra::blade_new::grade;
        let terms: Vec<_> = self
            .terms
            .iter()
            .filter(|(mask, _)| grade(*mask) == g)
            .cloned()
            .collect();
        PartialMv {
            terms,
            arity: self.arity,
        }
    }
}

// ═══════════════════════════════════════════════════════════
// B2: EXPR PARTIAL EVALUATOR
// ═══════════════════════════════════════════════════════════

/// Partially evaluate an Expr, treating specified parameters as polynomial variables
/// and everything else as concrete values.
///
/// This is the core of the compilation pipeline: it separates the dynamic probe
/// coordinates from the static algebra structure, producing a PartialMv whose
/// coefficients are polynomials in the probe parameters.
pub fn partial_eval_expr(
    expr: &crate::governance::expr::Expr,
    sig: &Signature,
    derived_gens: &[Mv],
    arity: usize,
) -> PartialMv {
    use crate::governance::expr::Expr;
    match expr {
        Expr::Param(k) => {
            if *k < arity {
                // Dynamic probe parameter → polynomial variable on the scalar blade
                PartialMv::blade_poly(0, IntPoly::variable(*k, arity))
            } else {
                PartialMv::zero(arity)
            }
        }
        Expr::Generator(g) => PartialMv::blade(1u64 << g, 1, arity),
        Expr::DerivedGen(k) => {
            if *k < derived_gens.len() {
                let (pmv, _den) = PartialMv::from_mv(&derived_gens[*k], arity);
                pmv
            } else {
                PartialMv::zero(arity)
            }
        }
        Expr::Literal(s) => {
            if let Some(r) = s.try_as_rat() {
                // Scale to integer: store numerator, track denominator separately
                // For simplicity, store as-is (numerator only, denominator tracked at top level)
                PartialMv::scalar(r.num() as i64, arity)
            } else {
                PartialMv::zero(arity)
            }
        }
        Expr::Add(a, b) => {
            let pa = partial_eval_expr(a, sig, derived_gens, arity);
            let pb = partial_eval_expr(b, sig, derived_gens, arity);
            pa.add(&pb)
        }
        Expr::Mul(a, b) => {
            let pa = partial_eval_expr(a, sig, derived_gens, arity);
            let pb = partial_eval_expr(b, sig, derived_gens, arity);
            pa.mul(&pb, sig)
        }
        Expr::Neg(a) => partial_eval_expr(a, sig, derived_gens, arity).neg(),
        Expr::Pow(a, n) => {
            let base = partial_eval_expr(a, sig, derived_gens, arity);
            let mut result = PartialMv::scalar(1, arity);
            for _ in 0..*n {
                result = result.mul(&base, sig);
            }
            result
        }
        _ => {
            // FieldEval, ValueRef, ConstructionRef — not handled in partial eval
            PartialMv::zero(arity)
        }
    }
}

// ═══════════════════════════════════════════════════════════
// B3: CONSTRUCTION-BASED FIELD COMPILATION
// ═══════════════════════════════════════════════════════════

/// Compile a field evaluation using the full probe construction expression.
///
/// Unlike `compile_field_eval` which assumes a linear probe (Σ xₖ·gen(k)),
/// this uses the actual Expr body of the probe construction, handling
/// nonlinear terms like CGA's -½Σxᵢ² on e∞.
pub fn compile_field_eval_from_construction(
    probe_construction_body: &crate::governance::expr::Expr,
    object: &Mv,
    field_op: &FieldOp,
    sig: &Signature,
    derived_gens: &[Mv],
    arity: usize,
) -> CompiledFieldEval {
    let probe = partial_eval_expr(probe_construction_body, sig, derived_gens, arity);
    let (obj_partial, obj_den) = PartialMv::from_mv(object, arity);

    match field_op {
        FieldOp::ScalarProduct => {
            compile_partial_scalar_product(&probe, &obj_partial, sig, arity, obj_den)
        }
        FieldOp::InnerProduct => {
            compile_partial_inner_product(&probe, &obj_partial, sig, arity, obj_den)
        }
        FieldOp::OuterProduct => {
            compile_partial_outer_product(&probe, &obj_partial, sig, arity, obj_den)
        }
        FieldOp::LeftContraction => {
            compile_partial_left_contraction(&probe, &obj_partial, sig, arity, obj_den)
        }
        FieldOp::GeometricProduct => {
            compile_partial_geometric(&probe, &obj_partial, sig, arity, obj_den)
        }
        FieldOp::GradeProduct(k) => {
            compile_partial_grade_product(&probe, &obj_partial, sig, arity, obj_den, *k)
        }
    }
}

// ═══════════════════════════════════════════════════════════
// B4: PARTIAL MV FIELD COMPILERS
// ═══════════════════════════════════════════════════════════

/// Compile ⟨probe, object⟩ (scalar product) from PartialMv probe.
fn compile_partial_scalar_product(
    probe: &PartialMv,
    object: &PartialMv,
    sig: &Signature,
    arity: usize,
    denominator: i64,
) -> CompiledFieldEval {
    let product = probe.mul(object, sig);
    let scalar_poly = product.scalar_part();
    CompiledFieldEval {
        polys: vec![(0, scalar_poly)],
        denominator,
        probe_arity: arity,
    }
}

/// Compile probe · object (Hestenes inner product) from PartialMv probe.
fn compile_partial_inner_product(
    probe: &PartialMv,
    object: &PartialMv,
    sig: &Signature,
    arity: usize,
    denominator: i64,
) -> CompiledFieldEval {
    use crate::algebra::blade_new::grade;
    use crate::algebra::product_new::blade_product;
    // Inner product: for each (a, b) blade pair, keep if |grade(a) - grade(b)| == grade(result)
    // and grade(a) != 0 and grade(b) != 0
    let mut result: BTreeMap<BladeMask, IntPoly> = BTreeMap::new();
    for (ma, pa) in &probe.terms {
        let ga = grade(*ma);
        if ga == 0 {
            continue;
        }
        for (mb, pb) in &object.terms {
            let gb = grade(*mb);
            if gb == 0 {
                continue;
            }
            let (rm, sign) = blade_product(*ma, *mb, sig);
            if sign == 0 {
                continue;
            }
            let gr = grade(rm);
            let expected = ga.abs_diff(gb);
            if gr != expected {
                continue;
            }
            let product = pa.mul(pb).scale(sign as i64);
            let entry = result.entry(rm).or_insert_with(|| IntPoly::zero(arity));
            *entry = entry.add(&product);
        }
    }
    let polys: Vec<_> = result.into_iter().filter(|(_, p)| !p.is_zero()).collect();
    CompiledFieldEval {
        polys,
        denominator,
        probe_arity: arity,
    }
}

/// Compile probe ∧ object (outer product) from PartialMv probe.
fn compile_partial_outer_product(
    probe: &PartialMv,
    object: &PartialMv,
    sig: &Signature,
    arity: usize,
    denominator: i64,
) -> CompiledFieldEval {
    use crate::algebra::blade_new::grade;
    use crate::algebra::product_new::blade_product;
    let mut result: BTreeMap<BladeMask, IntPoly> = BTreeMap::new();
    for (ma, pa) in &probe.terms {
        let ga = grade(*ma);
        for (mb, pb) in &object.terms {
            let gb = grade(*mb);
            let (rm, sign) = blade_product(*ma, *mb, sig);
            if sign == 0 {
                continue;
            }
            if grade(rm) != ga + gb {
                continue;
            } // outer product: grades must add
            let product = pa.mul(pb).scale(sign as i64);
            let entry = result.entry(rm).or_insert_with(|| IntPoly::zero(arity));
            *entry = entry.add(&product);
        }
    }
    let polys: Vec<_> = result.into_iter().filter(|(_, p)| !p.is_zero()).collect();
    CompiledFieldEval {
        polys,
        denominator,
        probe_arity: arity,
    }
}

/// Compile probe ⌋ object (left contraction) from PartialMv probe.
fn compile_partial_left_contraction(
    probe: &PartialMv,
    object: &PartialMv,
    sig: &Signature,
    arity: usize,
    denominator: i64,
) -> CompiledFieldEval {
    use crate::algebra::blade_new::grade;
    use crate::algebra::product_new::blade_product;
    let mut result: BTreeMap<BladeMask, IntPoly> = BTreeMap::new();
    for (ma, pa) in &probe.terms {
        let ga = grade(*ma);
        for (mb, pb) in &object.terms {
            let gb = grade(*mb);
            if ga > gb {
                continue;
            } // left contraction: grade(a) ≤ grade(b)
            let (rm, sign) = blade_product(*ma, *mb, sig);
            if sign == 0 {
                continue;
            }
            if grade(rm) != gb - ga {
                continue;
            } // result grade = grade(b) - grade(a)
            let product = pa.mul(pb).scale(sign as i64);
            let entry = result.entry(rm).or_insert_with(|| IntPoly::zero(arity));
            *entry = entry.add(&product);
        }
    }
    let polys: Vec<_> = result.into_iter().filter(|(_, p)| !p.is_zero()).collect();
    CompiledFieldEval {
        polys,
        denominator,
        probe_arity: arity,
    }
}

/// Compile probe * object (full geometric product) from PartialMv probe.
fn compile_partial_geometric(
    probe: &PartialMv,
    object: &PartialMv,
    sig: &Signature,
    arity: usize,
    denominator: i64,
) -> CompiledFieldEval {
    let product = probe.mul(object, sig);
    let polys: Vec<_> = product
        .terms
        .into_iter()
        .filter(|(_, p)| !p.is_zero())
        .collect();
    CompiledFieldEval {
        polys,
        denominator,
        probe_arity: arity,
    }
}

/// Compile ⟨probe * object⟩_k (grade-k component of geometric product).
fn compile_partial_grade_product(
    probe: &PartialMv,
    object: &PartialMv,
    sig: &Signature,
    arity: usize,
    denominator: i64,
    k: u8,
) -> CompiledFieldEval {
    let product = probe.mul(object, sig).grade_project(k);
    let polys: Vec<_> = product
        .terms
        .into_iter()
        .filter(|(_, p)| !p.is_zero())
        .collect();
    CompiledFieldEval {
        polys,
        denominator,
        probe_arity: arity,
    }
}

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

    #[test]
    fn intpoly_constant() {
        let p = IntPoly::constant(5, 2);
        assert_eq!(p.eval(&[0, 0]), 5);
        assert_eq!(p.eval(&[100, 200]), 5);
    }

    #[test]
    fn intpoly_variable() {
        let x = IntPoly::variable(0, 2);
        let y = IntPoly::variable(1, 2);
        assert_eq!(x.eval(&[3, 7]), 3);
        assert_eq!(y.eval(&[3, 7]), 7);
    }

    #[test]
    fn intpoly_add() {
        let x = IntPoly::variable(0, 2);
        let y = IntPoly::variable(1, 2);
        let sum = x.add(&y);
        assert_eq!(sum.eval(&[3, 7]), 10);
    }

    #[test]
    fn intpoly_mul() {
        let x = IntPoly::variable(0, 2);
        let y = IntPoly::variable(1, 2);
        let product = x.mul(&y);
        assert_eq!(product.eval(&[3, 7]), 21);
    }

    #[test]
    fn intpoly_quadratic() {
        // x² + 2xy + y² = (x+y)²
        let x = IntPoly::variable(0, 2);
        let y = IntPoly::variable(1, 2);
        let sum = x.add(&y);
        let sq = sum.mul(&sum);
        assert_eq!(sq.eval(&[3, 4]), 49); // (3+4)²
        assert_eq!(sq.eval(&[0, 5]), 25);
    }

    #[test]
    fn compile_vga_scalar_product() {
        // VGA(3): probe is a vector, object is a vector
        // ⟨probe, object⟩ = Σ probe_i * object_i
        let sig = Signature::new(0, 0, 3).unwrap();
        let object = Mv::from_rat_terms(&[
            (0b001, Rat::from(3)),
            (0b010, Rat::from(4)),
            (0b100, Rat::from(5)),
        ]);

        let compiled = compile_field_eval(
            &object,
            &FieldOp::ScalarProduct,
            &sig,
            3, // 3 probe params
            0, // probe gens start at 0
        );

        // At probe = (1, 0, 0): ⟨probe, object⟩ = 3
        assert_eq!(compiled.eval_scalar(&[1, 0, 0]), 3);
        // At probe = (0, 1, 0): ⟨probe, object⟩ = 4
        assert_eq!(compiled.eval_scalar(&[0, 1, 0]), 4);
        // At probe = (1, 1, 1): ⟨probe, object⟩ = 3+4+5 = 12
        assert_eq!(compiled.eval_scalar(&[1, 1, 1]), 12);
    }

    #[test]
    fn compiled_degree() {
        // Scalar product with linear probe → degree 1 polynomial
        let sig = Signature::new(0, 0, 2).unwrap();
        let object = Mv::from_rat_terms(&[(0b01, Rat::from(1)), (0b10, Rat::from(1))]);
        let compiled = compile_field_eval(&object, &FieldOp::ScalarProduct, &sig, 2, 0);
        assert!(compiled.polys.len() >= 1);
        assert!(compiled.polys[0].1.degree() <= 1);
    }

    #[test]
    fn partial_mv_add() {
        let a = PartialMv::blade(0b01, 3, 2);
        let b = PartialMv::blade(0b10, 4, 2);
        let sum = a.add(&b);
        assert_eq!(sum.terms.len(), 2);
    }

    #[test]
    fn partial_mv_mul_geometric() {
        // e1 * e2 = e12 in VGA(3)
        let sig = Signature::new(0, 0, 3).unwrap();
        let a = PartialMv::blade(0b001, 1, 2);
        let b = PartialMv::blade(0b010, 1, 2);
        let product = a.mul(&b, &sig);
        assert_eq!(product.terms.len(), 1);
        assert_eq!(product.terms[0].0, 0b011); // e12
    }

    #[test]
    fn partial_eval_linear_probe() {
        // VGA(3): probe = p0*e0 + p1*e1 + p2*e2
        use crate::governance::expr::Expr;
        let sig = Signature::new(0, 0, 3).unwrap();
        let body = Expr::Add(
            Expr::add(
                Expr::mul(Expr::param(0), Expr::gen(0)),
                Expr::mul(Expr::param(1), Expr::gen(1)),
            ),
            Expr::mul(Expr::param(2), Expr::gen(2)),
        );
        let probe = partial_eval_expr(&body, &sig, &[], 3);
        // Should have 3 blade terms, each with a degree-1 polynomial
        assert_eq!(probe.terms.len(), 3);
        for (_, poly) in &probe.terms {
            assert_eq!(poly.degree(), 1);
        }
    }

    #[test]
    fn compile_from_construction_vga_scalar() {
        // Same test as compile_vga_scalar_product but through the construction path
        use crate::governance::expr::Expr;
        let sig = Signature::new(0, 0, 3).unwrap();
        let object = Mv::from_rat_terms(&[
            (0b001, Rat::from(3)),
            (0b010, Rat::from(4)),
            (0b100, Rat::from(5)),
        ]);
        let probe_body = Expr::Add(
            Expr::add(
                Expr::mul(Expr::param(0), Expr::gen(0)),
                Expr::mul(Expr::param(1), Expr::gen(1)),
            ),
            Expr::mul(Expr::param(2), Expr::gen(2)),
        );
        let compiled = compile_field_eval_from_construction(
            &probe_body,
            &object,
            &FieldOp::ScalarProduct,
            &sig,
            &[],
            3,
        );
        // ⟨(1,0,0), (3,4,5)⟩ = 3
        assert_eq!(compiled.eval_scalar(&[1, 0, 0]), 3);
        // ⟨(1,1,1), (3,4,5)⟩ = 12
        assert_eq!(compiled.eval_scalar(&[1, 1, 1]), 12);
    }

    #[test]
    fn compile_outer_product_vanishes() {
        // Outer product of parallel vectors = 0
        use crate::governance::expr::Expr;
        let sig = Signature::new(0, 0, 2).unwrap();
        // Object is e1
        let object = Mv::from_rat_terms(&[(0b01, Rat::from(1))]);
        // Probe is x*e1 (parallel to object)
        let probe_body = *Expr::mul(Expr::param(0), Expr::gen(0));
        let compiled = compile_field_eval_from_construction(
            &probe_body,
            &object,
            &FieldOp::OuterProduct,
            &sig,
            &[],
            1,
        );
        // e1 ∧ e1 = 0 for any x
        assert_eq!(compiled.eval_scalar(&[5]), 0);
        assert!(compiled.polys.is_empty() || compiled.polys.iter().all(|(_, p)| p.is_zero()));
    }

    #[test]
    fn compile_geometric_product_all_grades() {
        // Geometric product of two vectors gives scalar + bivector
        use crate::governance::expr::Expr;
        let sig = Signature::new(0, 0, 2).unwrap();
        let object = Mv::from_rat_terms(&[(0b01, Rat::from(1)), (0b10, Rat::from(1))]);
        let probe_body = Expr::Add(
            Expr::mul(Expr::param(0), Expr::gen(0)),
            Expr::mul(Expr::param(1), Expr::gen(1)),
        );
        let compiled = compile_field_eval_from_construction(
            &probe_body,
            &object,
            &FieldOp::GeometricProduct,
            &sig,
            &[],
            2,
        );
        // Should have scalar and bivector components
        let masks: Vec<u64> = compiled.polys.iter().map(|(m, _)| *m).collect();
        assert!(masks.contains(&0b00), "should have scalar component");
        assert!(masks.contains(&0b11), "should have bivector component");
    }
}