geoit 0.0.2

use crate::algebra::blade_new::{grade, BladeMask};
use crate::algebra::mv::Mv;
use crate::algebra::signature::Signature;
use crate::governance::field::FieldOp;
use crate::governance::poly::Poly;
use crate::governance::profile::GeneratorProfile;
use crate::scalar::Rat;

/// A geometric class: defines a set of multivectors by grade constraint
/// and polynomial equations/inequalities on blade coefficients.
///
/// The grade_mask determines which grades are permitted (linear subspace).
/// The equations are polynomials that must equal zero (algebraic variety).
/// The inequalities are polynomials that must be nonzero (open conditions).
///
/// Polynomial variables use compact indexing: variable i = the i-th blade
/// mask permitted by the grade_mask, in sorted order. The VariableMap
/// (computed from signature + grade_mask) provides the translation.
#[derive(Clone, Debug, Default)]
pub struct GeomClass {
    /// Bitmask: bit k set means grade k is permitted.
    pub grade_mask: u64,
    /// Polynomial equations on blade coefficients (each must = 0).
    pub equations: Vec<Poly>,
    /// Polynomial inequalities on blade coefficients (each must ≠ 0).
    pub inequalities: Vec<Poly>,
    /// How this class's Mv is evaluated against a probe point for rendering.
    pub field_op: FieldOp,
    /// v0.0.3: Expected generator participation profile.
    /// If set, `govern()` rejects Mvs whose profile uses generators outside this mask.
    pub expected_profile: Option<GeneratorProfile>,
}

impl GeomClass {
    /// Create a GeomClass with only a grade constraint (default field_op = ScalarProduct).
    pub fn grades_only(grades: &[u8]) -> Self {
        let mut mask = 0u64;
        for &g in grades {
            mask |= 1u64 << g;
        }
        GeomClass {
            grade_mask: mask,
            equations: vec![],
            inequalities: vec![],
            field_op: FieldOp::default(),
            expected_profile: None,
        }
    }

    /// Check if a grade is permitted.
    pub fn grade_permitted(&self, g: u8) -> bool {
        self.grade_mask & (1u64 << g) != 0
    }

    /// List of permitted grades.
    pub fn permitted_grades(&self) -> Vec<u8> {
        (0..64u8).filter(|&g| self.grade_permitted(g)).collect()
    }

    /// Number of free parameters (dimension of the solution variety).
    /// Computed via Gröbner basis: total blade variables minus determined variables.
    pub fn dimension(
        &self,
        sig: &Signature,
    ) -> Result<usize, crate::governance::groebner::GroebnerError> {
        let vm = crate::governance::reading::VariableMap::for_grade_mask(sig, self.grade_mask);
        if self.equations.is_empty() {
            return Ok(vm.num_vars);
        }
        let basis = crate::governance::groebner::groebner_basis(self.equations.clone())?;
        Ok(crate::governance::groebner::free_variables(&basis, vm.num_vars).len())
    }

    /// Number of independent equations constraining the class.
    pub fn codimension(
        &self,
        sig: &Signature,
    ) -> Result<usize, crate::governance::groebner::GroebnerError> {
        let vm = crate::governance::reading::VariableMap::for_grade_mask(sig, self.grade_mask);
        Ok(vm.num_vars - self.dimension(sig)?)
    }

    /// Predict the phase of all instances, if determinable from class structure.
    /// Returns Some(phase) if all instances must have that phase, None if instance-dependent.
    pub fn expected_phase(&self, sig: &Signature) -> Option<crate::governance::phase::Phase> {
        use crate::governance::phase::Phase;
        // If any permitted grade involves a degenerate generator, all instances are Evaluation
        let n = sig.n();
        for g in self.permitted_grades() {
            for mask in 0..(1u64 << n) {
                if grade(mask) != g {
                    continue;
                }
                for k in 0..n {
                    if mask & (1u64 << k) != 0 && sig.is_degenerate(k) {
                        return Some(Phase::Evaluation);
                    }
                }
            }
        }
        // If the norm polynomial is in the equation ideal, norm is forced to zero → Evaluation
        let vm = crate::governance::reading::VariableMap::for_grade_mask(sig, self.grade_mask);
        let np = norm_poly(sig, self.grade_mask, vm.num_vars, &vm.mask_to_var);
        if !np.is_zero() && !self.equations.is_empty() {
            let remainder = crate::governance::groebner::reduce_by_set(&np, &self.equations);
            if remainder.is_zero() {
                return Some(Phase::Evaluation);
            }
            // Also check against the Gröbner basis for stronger reduction
            if let Ok(basis) = crate::governance::groebner::groebner_basis(self.equations.clone()) {
                let remainder2 = crate::governance::groebner::reduce_by_set(&np, &basis);
                if remainder2.is_zero() {
                    return Some(Phase::Evaluation);
                }
            }
        }
        // Check if norm is explicitly in the equations list
        for eq in &self.equations {
            if *eq == np {
                return Some(Phase::Evaluation);
            }
        }
        None // instance-dependent
    }

    /// Does every Mv satisfying `other` also satisfy `self`?
    /// other ⊂ self iff: other's grade_mask is a subset of self's grade_mask AND
    /// every equation of self reduces to zero modulo other's equation ideal.
    pub fn contains(
        &self,
        other: &GeomClass,
        _sig: &Signature,
    ) -> Result<bool, crate::governance::groebner::GroebnerError> {
        // Grade check: other's grades must be a subset of self's grades
        if other.grade_mask & !self.grade_mask != 0 {
            return Ok(false);
        }
        // If self has no equations, it contains everything with compatible grades
        if self.equations.is_empty() {
            return Ok(true);
        }
        // Every equation of self must be implied by other's equations
        if other.equations.is_empty() {
            return Ok(false);
        }
        let basis = crate::governance::groebner::groebner_basis(other.equations.clone())?;
        for eq in &self.equations {
            let remainder = crate::governance::groebner::reduce_by_set(eq, &basis);
            if !remainder.is_zero() {
                return Ok(false);
            }
        }
        Ok(true)
    }
}

/// Build the norm-squared polynomial for blades at permitted grades in a given signature.
///
/// norm²(P) = P · rev(P) = Σ_i sign_i * c_i²
/// where sign_i = rev_sign(grade(blade_i)) * blade_i²
pub fn norm_poly(
    sig: &Signature,
    grade_mask: u64,
    num_vars: usize,
    mask_to_var: &std::collections::BTreeMap<BladeMask, usize>,
) -> Poly {
    let n = sig.n();
    let mut p = Poly::zero(num_vars);
    for mask in 0..(1u64 << n) {
        let g = grade(mask);
        if grade_mask & (1u64 << g) == 0 {
            continue;
        }
        if let Some(&vi) = mask_to_var.get(&mask) {
            // Reverse sign for this grade
            let k = g as u32;
            let rev_sign: i64 = if k < 2 || (k * (k - 1) / 2) & 1 == 0 {
                1
            } else {
                -1
            };
            // Blade self-product: product of squares of generators in the blade
            let mut blade_sq: i64 = rev_sign;
            let mut m = mask;
            while m != 0 {
                let gen = m.trailing_zeros() as u8;
                blade_sq *= sig.generator_square(gen) as i64;
                m &= m - 1;
            }
            if blade_sq != 0 {
                let mut exp = vec![0u8; num_vars];
                exp[vi] = 2;
                let entry = p.terms.entry(exp).or_insert(Rat::ZERO);
                *entry += Rat::from(blade_sq);
            }
        }
    }
    p
}

/// Build inner-product polynomial: P · dg = value, as polynomial in P's blade coefficients.
/// For grade-1 blades: inner(P, dg) = Σ_i sig.gen_sq(i) * c_i * dg_i
pub fn inner_product_poly(
    dg: &Mv,
    sig: &Signature,
    grade_mask: u64,
    value: Rat,
    num_vars: usize,
    mask_to_var: &std::collections::BTreeMap<BladeMask, usize>,
) -> Poly {
    let n = sig.n();
    let mut p = Poly::zero(num_vars);
    for mask in 0..(1u64 << n) {
        let g = grade(mask);
        if grade_mask & (1u64 << g) == 0 {
            continue;
        }
        let dg_coeff = dg.coefficient(mask);
        if dg_coeff.is_zero() {
            continue;
        }
        if let Some(&vi) = mask_to_var.get(&mask) {
            if g == 1 {
                let gen = mask.trailing_zeros() as u8;
                let sq = sig.generator_square(gen);
                let dg_rat = dg_coeff.try_as_rat().unwrap_or(Rat::ZERO);
                let coeff = Rat::from(sq as i64) * dg_rat;
                if !coeff.is_zero() {
                    let mut exp = vec![0u8; num_vars];
                    exp[vi] = 1;
                    let entry = p.terms.entry(exp).or_insert(Rat::ZERO);
                    *entry += coeff;
                }
            }
        }
    }
    // Subtract expected value
    if !value.is_zero() {
        let entry = p.terms.entry(vec![0u8; num_vars]).or_insert(Rat::ZERO);
        *entry -= value;
    }
    p
}

/// Build a single-variable polynomial for a blade coefficient: c_mask
pub fn blade_var_poly(
    mask: BladeMask,
    num_vars: usize,
    mask_to_var: &std::collections::BTreeMap<BladeMask, usize>,
) -> Option<Poly> {
    mask_to_var
        .get(&mask)
        .map(|&vi| Poly::variable(vi, num_vars))
}

// ═══════════════════════════════════════════════════════════
// DUALITY
// ═══════════════════════════════════════════════════════════

/// Compute the dual of a GeomClass under the algebra's Hodge dual.
/// Grade k maps to grade n-k. Polynomial variables are remapped.
pub fn dual_class(class: &GeomClass, sig: &Signature) -> GeomClass {
    let n = sig.n();
    let ps_mask = sig.pseudoscalar_mask();

    // Dual grade mask: bit k → bit (n-k)
    let mut dual_mask = 0u64;
    for g in 0..=n {
        if class.grade_mask & (1u64 << g) != 0 {
            dual_mask |= 1u64 << (n - g);
        }
    }

    // Build variable maps for source and target
    let src_vm = crate::governance::reading::VariableMap::for_grade_mask(sig, class.grade_mask);
    let dst_vm = crate::governance::reading::VariableMap::for_grade_mask(sig, dual_mask);

    // Build remap: for each source variable (blade mask m),
    // find the dual blade (ps_mask ^ m) and its sign, and its target variable index.
    let nv_dst = dst_vm.num_vars;

    let remap_equations = |polys: &[Poly]| -> Vec<Poly> {
        polys
            .iter()
            .map(|poly| {
                let mut result = Poly::zero(nv_dst);
                for (exp, &coeff) in &poly.terms {
                    if coeff.is_zero() {
                        continue;
                    }
                    // Map each monomial: variable i with exponent e[i]
                    // becomes variable remap[i] with same exponent
                    let mut new_exp = vec![0u8; nv_dst];
                    let mut valid = true;
                    for (src_vi, &e) in exp.iter().enumerate() {
                        if e == 0 {
                            continue;
                        }
                        let src_mask = src_vm.var_to_mask[src_vi];
                        let dual_blade = ps_mask ^ src_mask;
                        if let Some(&dst_vi) = dst_vm.mask_to_var.get(&dual_blade) {
                            new_exp[dst_vi] = e;
                        } else {
                            valid = false;
                            break;
                        }
                    }
                    if valid {
                        // Sign: determined by blade * complement = ±pseudoscalar
                        // For simplicity, use unsigned remap (signs are algebra-dependent)
                        let entry = result.terms.entry(new_exp).or_insert(Rat::ZERO);
                        *entry += coeff;
                    }
                }
                result
            })
            .collect()
    };

    GeomClass {
        grade_mask: dual_mask,
        equations: remap_equations(&class.equations),
        inequalities: remap_equations(&class.inequalities),
        field_op: FieldOp::default(),
        expected_profile: None,
    }
}

// ═══════════════════════════════════════════════════════════
// INCIDENCE
// ═══════════════════════════════════════════════════════════

/// Build the incidence polynomial: the scalar product P · S as a polynomial
/// in the combined variables of two classes.
///
/// Variables 0..n_a are class_a blade coefficients, n_a..n_a+n_b are class_b's.
/// The result polynomial equals zero iff the two objects are incident.
pub fn incidence_poly(class_a: &GeomClass, class_b: &GeomClass, sig: &Signature) -> Poly {
    let vm_a = crate::governance::reading::VariableMap::for_grade_mask(sig, class_a.grade_mask);
    let vm_b = crate::governance::reading::VariableMap::for_grade_mask(sig, class_b.grade_mask);
    let n_a = vm_a.num_vars;
    let n_total = n_a + vm_b.num_vars;
    let mut poly = Poly::zero(n_total);

    // Scalar product: sum over matching blades of metric_sign * a_coeff * b_coeff
    for (&mask_a, &vi_a) in &vm_a.mask_to_var {
        for (&mask_b, &vi_b) in &vm_b.mask_to_var {
            if mask_a != mask_b {
                continue;
            }
            // scalar_product(e_mask, e_mask) = reverse_sign * blade_square
            let g = grade(mask_a);
            let k = g as u32;
            let rev_sign: i64 = if k < 2 || (k * (k - 1) / 2) & 1 == 0 {
                1
            } else {
                -1
            };
            let mut blade_sq: i64 = rev_sign;
            let mut m = mask_a;
            while m != 0 {
                let gen = m.trailing_zeros() as u8;
                blade_sq *= sig.generator_square(gen) as i64;
                m &= m - 1;
            }
            if blade_sq != 0 {
                let mut exp = vec![0u8; n_total];
                exp[vi_a] = 1;
                exp[n_a + vi_b] = 1;
                let entry = poly.terms.entry(exp).or_insert(Rat::ZERO);
                *entry += Rat::from(blade_sq);
            }
        }
    }
    poly
}

// ═══════════════════════════════════════════════════════════
// OPERATION GRADE PREDICTION
// ═══════════════════════════════════════════════════════════

/// Predict possible result grades for the outer product of two grade masks.
/// Grade j ∧ Grade k → Grade j+k (if j+k ≤ n).
pub fn outer_result_grades(a_mask: u64, b_mask: u64, n: u8) -> u64 {
    let mut result = 0u64;
    for j in 0..=n {
        if a_mask & (1u64 << j) == 0 {
            continue;
        }
        for k in 0..=n {
            if b_mask & (1u64 << k) == 0 {
                continue;
            }
            let s = j + k;
            if s <= n {
                result |= 1u64 << s;
            }
        }
    }
    result
}

/// Predict possible result grades for the inner product.
/// Grade j · Grade k → Grade |j-k|.
pub fn inner_result_grades(a_mask: u64, b_mask: u64, n: u8) -> u64 {
    let mut result = 0u64;
    for j in 0..=n {
        if a_mask & (1u64 << j) == 0 {
            continue;
        }
        for k in 0..=n {
            if b_mask & (1u64 << k) == 0 {
                continue;
            }
            let d = j.abs_diff(k);
            result |= 1u64 << d;
        }
    }
    result
}

/// Predict possible result grades for the geometric product.
/// Grade j * Grade k → Grades |j-k|, |j-k|+2, ..., j+k.
pub fn geometric_result_grades(a_mask: u64, b_mask: u64, n: u8) -> u64 {
    let mut result = 0u64;
    for j in 0..=n {
        if a_mask & (1u64 << j) == 0 {
            continue;
        }
        for k in 0..=n {
            if b_mask & (1u64 << k) == 0 {
                continue;
            }
            let lo = j.abs_diff(k);
            let hi = (j + k).min(n);
            let mut g = lo;
            while g <= hi {
                result |= 1u64 << g;
                g += 2;
            }
        }
    }
    result
}

// ═══════════════════════════════════════════════════════════
// CLASS RELATIONS
// ═══════════════════════════════════════════════════════════

/// A declared relationship between two geometric classes.
#[derive(Clone, Debug)]
pub enum ClassRelation {
    /// Classes are related by Hodge duality.
    Dual(usize, usize),
    /// First class is a subclass of the second.
    Subclass(usize, usize),
    /// Two classes participate in an incidence relation.
    Incidence(usize, usize),
}

impl std::fmt::Display for ClassRelation {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            ClassRelation::Dual(a, b) => write!(f, "Dual({}, {})", a, b),
            ClassRelation::Subclass(sub, sup) => write!(f, "Subclass({} ⊂ {})", sub, sup),
            ClassRelation::Incidence(a, b) => write!(f, "Incidence({}, {})", a, b),
        }
    }
}

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

    #[test]
    fn grades_only_class() {
        let c = GeomClass::grades_only(&[1]);
        assert!(c.grade_permitted(1));
        assert!(!c.grade_permitted(0));
        assert!(!c.grade_permitted(2));
        assert_eq!(c.permitted_grades(), vec![1]);
    }

    #[test]
    fn mixed_grade_class() {
        let c = GeomClass::grades_only(&[0, 2]);
        assert!(c.grade_permitted(0));
        assert!(!c.grade_permitted(1));
        assert!(c.grade_permitted(2));
        assert_eq!(c.permitted_grades(), vec![0, 2]);
    }

    #[test]
    fn empty_equations() {
        let c = GeomClass::grades_only(&[1]);
        assert!(c.equations.is_empty());
        assert!(c.inequalities.is_empty());
    }

    #[test]
    fn vga_vector_dimension() {
        let sig = Signature::new(0, 0, 3).unwrap();
        let c = GeomClass::grades_only(&[1]);
        assert_eq!(c.dimension(&sig).unwrap(), 3); // 3 free blade coefficients
        assert_eq!(c.codimension(&sig).unwrap(), 0);
    }

    #[test]
    fn vga_bivector_dimension() {
        let sig = Signature::new(0, 0, 3).unwrap();
        let c = GeomClass::grades_only(&[2]);
        assert_eq!(c.dimension(&sig).unwrap(), 3); // 3 bivector blades, no equations
    }

    #[test]
    fn cga_point_dimension() {
        let sig = Signature::new(1, 0, 4).unwrap();
        let gm = 0b10u64;
        let vm = crate::governance::reading::VariableMap::for_grade_mask(&sig, gm);
        let einf = Mv::from_rat_terms(&[(0b00001, Rat::from(-1)), (0b00010, Rat::from(1))]);
        let null_eq = norm_poly(&sig, gm, vm.num_vars, &vm.mask_to_var);
        let ip_eq = inner_product_poly(&einf, &sig, gm, Rat::ONE, vm.num_vars, &vm.mask_to_var);
        let c = GeomClass {
            grade_mask: gm,
            equations: vec![null_eq, ip_eq],
            inequalities: vec![],
            field_op: FieldOp::default(),
            expected_profile: None,
        };
        assert_eq!(c.dimension(&sig).unwrap(), 3); // x, y, z
        assert_eq!(c.codimension(&sig).unwrap(), 2); // null + normalization
    }

    #[test]
    fn cga_sphere_dimension() {
        let sig = Signature::new(1, 0, 4).unwrap();
        let gm = 0b10u64;
        let vm = crate::governance::reading::VariableMap::for_grade_mask(&sig, gm);
        let einf = Mv::from_rat_terms(&[(0b00001, Rat::from(-1)), (0b00010, Rat::from(1))]);
        let ip_eq = inner_product_poly(&einf, &sig, gm, Rat::ONE, vm.num_vars, &vm.mask_to_var);
        let c = GeomClass {
            grade_mask: gm,
            equations: vec![ip_eq],
            inequalities: vec![],
            field_op: FieldOp::default(),
            expected_profile: None,
        };
        assert_eq!(c.dimension(&sig).unwrap(), 4); // center + radius
        assert_eq!(c.codimension(&sig).unwrap(), 1); // normalization only
    }

    #[test]
    fn pga_point_dimension() {
        let sig = Signature::new(0, 1, 3).unwrap();
        let c = GeomClass {
            grade_mask: 0b1000,
            equations: vec![],
            inequalities: vec![],
            field_op: FieldOp::default(),
            expected_profile: None,
        };
        assert_eq!(c.dimension(&sig).unwrap(), 4); // 4 grade-3 blades, no equations
    }

    #[test]
    fn vga_vector_phase_unknown() {
        let sig = Signature::new(0, 0, 3).unwrap();
        let c = GeomClass::grades_only(&[1]);
        assert_eq!(c.expected_phase(&sig), None); // could be any norm sign
    }

    #[test]
    fn pga_point_phase_evaluation() {
        let sig = Signature::new(0, 1, 3).unwrap();
        let c = GeomClass {
            grade_mask: 0b1000,
            equations: vec![],
            inequalities: vec![],
            field_op: FieldOp::default(),
            expected_profile: None,
        };
        // PGA has degenerate generator g0, grade-3 blades involve g0
        assert_eq!(
            c.expected_phase(&sig),
            Some(crate::governance::phase::Phase::Evaluation)
        );
    }

    #[test]
    fn cga_point_phase_evaluation() {
        let sig = Signature::new(1, 0, 4).unwrap();
        let gm = 0b10u64;
        let vm = crate::governance::reading::VariableMap::for_grade_mask(&sig, gm);
        let einf = Mv::from_rat_terms(&[(0b00001, Rat::from(-1)), (0b00010, Rat::from(1))]);
        let null_eq = norm_poly(&sig, gm, vm.num_vars, &vm.mask_to_var);
        let ip_eq = inner_product_poly(&einf, &sig, gm, Rat::ONE, vm.num_vars, &vm.mask_to_var);
        let c = GeomClass {
            grade_mask: gm,
            equations: vec![null_eq, ip_eq],
            inequalities: vec![],
            field_op: FieldOp::default(),
            expected_profile: None,
        };
        // CGA Point: norm is forced to zero by equations → Evaluation
        assert_eq!(
            c.expected_phase(&sig),
            Some(crate::governance::phase::Phase::Evaluation)
        );
    }

    #[test]
    fn sphere_contains_point() {
        let sig = Signature::new(1, 0, 4).unwrap();
        let gm = 0b10u64;
        let vm = crate::governance::reading::VariableMap::for_grade_mask(&sig, gm);
        let einf = Mv::from_rat_terms(&[(0b00001, Rat::from(-1)), (0b00010, Rat::from(1))]);
        let null_eq = norm_poly(&sig, gm, vm.num_vars, &vm.mask_to_var);
        let ip_eq = inner_product_poly(&einf, &sig, gm, Rat::ONE, vm.num_vars, &vm.mask_to_var);

        let point_class = GeomClass {
            grade_mask: gm,
            equations: vec![null_eq.clone(), ip_eq.clone()],
            inequalities: vec![],
            field_op: FieldOp::default(),
            expected_profile: None,
        };
        let sphere_class = GeomClass {
            grade_mask: gm,
            equations: vec![ip_eq.clone()],
            inequalities: vec![],
            field_op: FieldOp::default(),
            expected_profile: None,
        };

        // Sphere contains Point: sphere's normalization equation is implied by point's equations
        assert!(sphere_class.contains(&point_class, &sig).unwrap());
        // Point does NOT contain Sphere: sphere doesn't satisfy null equation
        assert!(!point_class.contains(&sphere_class, &sig).unwrap());
    }

    #[test]
    fn different_grades_not_contained() {
        let sig = Signature::new(0, 0, 3).unwrap();
        let vectors = GeomClass::grades_only(&[1]);
        let bivectors = GeomClass::grades_only(&[2]);
        assert!(!vectors.contains(&bivectors, &sig).unwrap());
        assert!(!bivectors.contains(&vectors, &sig).unwrap());
    }

    // ─── Stage 3 tests ───

    #[test]
    fn dual_vga_vector_is_bivector() {
        let sig = Signature::new(0, 0, 3).unwrap();
        let vectors = GeomClass::grades_only(&[1]);
        let dual = dual_class(&vectors, &sig);
        // Grade 1 dual in Cl(0,0,3) with n=3 → grade 2
        assert!(dual.grade_permitted(2));
        assert!(!dual.grade_permitted(1));
        assert_eq!(dual.permitted_grades(), vec![2]);
    }

    #[test]
    fn dual_grade_mask_roundtrip() {
        let sig = Signature::new(0, 0, 3).unwrap();
        let vectors = GeomClass::grades_only(&[1]);
        let bivectors = dual_class(&vectors, &sig);
        let back = dual_class(&bivectors, &sig);
        assert_eq!(back.grade_mask, vectors.grade_mask);
    }

    #[test]
    fn dual_cga_point() {
        let sig = Signature::new(1, 0, 4).unwrap(); // n=5
        let point = GeomClass::grades_only(&[1]);
        let dual = dual_class(&point, &sig);
        // Grade 1 in 5-gen algebra → grade 4
        assert_eq!(dual.permitted_grades(), vec![4]);
    }

    #[test]
    fn outer_vector_vector_is_bivector() {
        // VGA3: Vector ∧ Vector → Bivector
        let result = outer_result_grades(0b10, 0b10, 3); // grade 1 ∧ grade 1
        assert_eq!(result, 0b100); // grade 2
    }

    #[test]
    fn outer_point_point_in_cga() {
        // CGA3: Point ∧ Point → grade 2
        let result = outer_result_grades(0b10, 0b10, 5); // grade 1 ∧ grade 1
        assert_eq!(result, 0b100); // grade 2 (PointPair in CGA)
    }

    #[test]
    fn inner_bivector_vector() {
        // Grade 2 · Grade 1 → Grade 1
        let result = inner_result_grades(0b100, 0b10, 3);
        assert_eq!(result, 0b10); // grade 1
    }

    #[test]
    fn geometric_vector_vector() {
        // Grade 1 * Grade 1 → Grades 0 and 2
        let result = geometric_result_grades(0b10, 0b10, 3);
        assert_eq!(result, 0b101); // grades 0 and 2
    }

    #[test]
    fn incidence_poly_structure() {
        let sig = Signature::new(0, 0, 3).unwrap();
        let vectors = GeomClass::grades_only(&[1]);
        let poly = incidence_poly(&vectors, &vectors, &sig);
        // 3 variables from each class = 6 total
        assert_eq!(poly.num_vars, 6);
        // Polynomial should be bilinear: terms like v_a0 * v_b0 + v_a1 * v_b1 + v_a2 * v_b2
        assert!(!poly.is_zero());
    }

    #[test]
    fn incidence_poly_evaluates() {
        let sig = Signature::new(0, 0, 3).unwrap();
        let vectors = GeomClass::grades_only(&[1]);
        let poly = incidence_poly(&vectors, &vectors, &sig);
        // Two orthogonal vectors: (1,0,0) and (0,1,0) → inner product = 0
        let values = vec![
            Rat::ONE,
            Rat::ZERO,
            Rat::ZERO, // vector A = e1
            Rat::ZERO,
            Rat::ONE,
            Rat::ZERO, // vector B = e2
        ];
        assert!(
            poly.eval(&values).is_zero(),
            "orthogonal vectors should be incident"
        );
        // Two parallel vectors: (1,0,0) and (1,0,0) → inner product = 1 ≠ 0
        let values2 = vec![
            Rat::ONE,
            Rat::ZERO,
            Rat::ZERO,
            Rat::ONE,
            Rat::ZERO,
            Rat::ZERO,
        ];
        assert!(
            !poly.eval(&values2).is_zero(),
            "parallel vectors should not be incident"
        );
    }

    #[test]
    fn class_relation_display() {
        let r = ClassRelation::Dual(0, 1);
        assert_eq!(format!("{}", r), "Dual(0, 1)");
        let r2 = ClassRelation::Subclass(0, 1);
        assert_eq!(format!("{}", r2), "Subclass(0 ⊂ 1)");
    }
}