geoit 0.0.2

use crate::algebra::blade_new::blade_to_mask;
use crate::algebra::blade_new::{grade, BladeMask};
use crate::algebra::mv::Mv;
use crate::algebra::product_new::blade_product;
use crate::algebra::signature::Signature;
use crate::scalar::Scalar;

// === Binary operations ===

/// Geometric product: no grade filter. The full Clifford product.
pub fn geometric(a: &Mv, b: &Mv, sig: &Signature) -> Mv {
    product_filtered(a, b, sig, |_, _, _| true)
}

/// Outer (wedge) product: keep only grade(a) + grade(b).
pub fn outer(a: &Mv, b: &Mv, sig: &Signature) -> Mv {
    product_filtered(a, b, sig, |ga, gb, gr| gr == ga + gb)
}

/// Inner product (Hestenes): |grade(a) - grade(b)|, both nonzero.
pub fn inner(a: &Mv, b: &Mv, sig: &Signature) -> Mv {
    product_filtered(a, b, sig, |ga, gb, gr| {
        ga > 0 && gb > 0 && gr == ga.abs_diff(gb)
    })
}

/// Left contraction: grade(result) = grade(b) - grade(a).
pub fn left_contract(a: &Mv, b: &Mv, sig: &Signature) -> Mv {
    product_filtered(a, b, sig, |ga, gb, gr| gb >= ga && gr == gb - ga)
}

/// Right contraction: grade(result) = grade(a) - grade(b).
pub fn right_contract(a: &Mv, b: &Mv, sig: &Signature) -> Mv {
    product_filtered(a, b, sig, |ga, gb, gr| ga >= gb && gr == ga - gb)
}

/// Scalar product: keep only grade 0.
pub fn scalar_product(a: &Mv, b: &Mv, sig: &Signature) -> Mv {
    product_filtered(a, b, sig, |_, _, gr| gr == 0)
}

/// Extract scalar part as a Scalar value.
pub fn scalar_product_value(a: &Mv, b: &Mv, sig: &Signature) -> Scalar {
    scalar_product(a, b, sig).coefficient(0)
}

/// Commutator product: (ab - ba) / 2.
pub fn commutator(a: &Mv, b: &Mv, sig: &Signature) -> Mv {
    let ab = geometric(a, b, sig);
    let ba = geometric(b, a, sig);
    let diff = ab - ba;
    diff.scale(&Scalar::from(crate::scalar::Rat::new(1, 2)))
}

// === Unary operations (no Signature needed) ===

/// Reverse: negate blades of grade k where k(k-1)/2 is odd.
/// ~(e₁₂₃) = -e₁₂₃ (grade 3: 3*2/2=3, odd)
/// ~(e₁₂) = -e₁₂ (grade 2: 2*1/2=1, odd)
/// ~(e₁) = e₁ (grade 1: 1*0/2=0, even)
/// ~(1) = 1 (grade 0: 0, even)
pub fn reverse(a: &Mv) -> Mv {
    let terms: Vec<_> = a
        .terms
        .iter()
        .map(|(blade, coeff)| {
            let k = blade.grade() as u32;
            let sign = if k < 2 || (k * (k - 1) / 2) & 1 == 0 {
                1i64
            } else {
                -1i64
            };
            let c = coeff.clone() * Scalar::from(sign);
            (blade.clone(), c)
        })
        .filter(|(_, c)| !c.is_zero())
        .collect();
    Mv { terms }
}

/// Grade involution: negate odd-grade blades.
pub fn grade_involution(a: &Mv) -> Mv {
    let terms: Vec<_> = a
        .terms
        .iter()
        .map(|(blade, coeff)| {
            let sign = if blade.grade() & 1 == 0 { 1i64 } else { -1i64 };
            let c = coeff.clone() * Scalar::from(sign);
            (blade.clone(), c)
        })
        .filter(|(_, c)| !c.is_zero())
        .collect();
    Mv { terms }
}

/// Clifford conjugate: reverse composed with grade involution.
pub fn conjugate(a: &Mv) -> Mv {
    grade_involution(&reverse(a))
}

// === Unary with Signature ===

/// Dual: a * pseudoscalar_inverse.
pub fn dual(a: &Mv, sig: &Signature) -> Mv {
    let ps_inv = pseudoscalar_inverse(sig);
    geometric(a, &ps_inv, sig)
}

/// Undual: a * pseudoscalar.
pub fn undual(a: &Mv, sig: &Signature) -> Mv {
    let ps = Mv::term(sig.pseudoscalar_mask(), Scalar::from(1i64));
    geometric(a, &ps, sig)
}

/// Compute the pseudoscalar inverse.
fn pseudoscalar_inverse(sig: &Signature) -> Mv {
    let ps_mask = sig.pseudoscalar_mask();
    let (_, sign) = blade_product(ps_mask, ps_mask, sig);
    assert!(
        sign != 0,
        "Pseudoscalar is not invertible (degenerate algebra)"
    );
    // ps * ps = sign * scalar, so ps⁻¹ = sign * ps (since sign² = 1)
    Mv::term(ps_mask, Scalar::from(sign as i64))
}

// === Derived operations ===

/// Sandwich product: a * b * reverse(a).
pub fn sandwich(a: &Mv, b: &Mv, sig: &Signature) -> Mv {
    let ab = geometric(a, b, sig);
    geometric(&ab, &reverse(a), sig)
}

/// Norm squared: scalar part of a * reverse(a).
pub fn norm_squared(a: &Mv, sig: &Signature) -> Scalar {
    scalar_product_value(a, &reverse(a), sig)
}

/// Grade projection: keep only grade k. (Delegates to Mv method.)
pub fn grade_project(a: &Mv, k: u8) -> Mv {
    a.grade_project(k)
}

/// Regressive product (meet): a ∨ b = undual(dual(a) ∧ dual(b)).
///
/// In a non-degenerate algebra, this computes the intersection
/// of the subspaces represented by a and b.
pub fn regressive(a: &Mv, b: &Mv, sig: &Signature) -> Mv {
    let da = dual(a, sig);
    let db = dual(b, sig);
    let wedge = outer(&da, &db, sig);
    undual(&wedge, sig)
}

// === Internal ===

/// Grade-indexed term groups for sparse product optimization.
struct GradeGroups {
    /// groups[g] = Vec of (mask, index_into_terms) for grade g
    groups: Vec<Vec<(BladeMask, usize)>>,
    /// The actual (mask, coeff) pairs, in order
    terms: Vec<(BladeMask, Scalar)>,
    /// Which grades have any terms
    populated: u64, // bitmask: bit g set iff groups[g] is non-empty
}

impl GradeGroups {
    fn from_mv(mv: &Mv) -> Self {
        let terms: Vec<(BladeMask, Scalar)> = mv.blades().map(|(m, c)| (m, c.clone())).collect();
        let max_grade = terms.iter().map(|(m, _)| grade(*m)).max().unwrap_or(0) as usize;
        let mut groups = vec![Vec::new(); max_grade + 1];
        let mut populated = 0u64;
        for (idx, &(mask, _)) in terms.iter().enumerate() {
            let g = grade(mask) as usize;
            groups[g].push((mask, idx));
            populated |= 1u64 << g;
        }
        GradeGroups {
            groups,
            terms,
            populated,
        }
    }
}

/// Generic filtered product with grade-pair prefiltering.
///
/// Optimization: before iterating term pairs, we check whether the
/// grade combination (ga, gb) can ever pass the filter for any result grade.
/// If not, we skip the entire grade group. For sparse Mvs in high-dimensional
/// algebras, this eliminates most of the work.
fn product_filtered<F>(a: &Mv, b: &Mv, sig: &Signature, filter: F) -> Mv
where
    F: Fn(u8, u8, u8) -> bool,
{
    if a.is_zero() || b.is_zero() {
        return Mv::new();
    }

    let n = sig.n();

    // For small Mvs (few terms), skip the grouping overhead
    if a.len() * b.len() <= 16 {
        return product_filtered_direct(a, b, sig, filter);
    }

    let ga = GradeGroups::from_mv(a);
    let gb = GradeGroups::from_mv(b);

    // Use flat accumulator for small algebras (n ≤ 12 → 4096 blades)
    if n <= 12 {
        return product_filtered_flat(&ga, &gb, sig, n, filter);
    }

    // Standard path with grade-pair prefiltering
    let mut result = Mv::new();
    for grade_a in 0..=n {
        if ga.populated & (1u64 << grade_a) == 0 {
            continue;
        }
        for grade_b in 0..=n {
            if gb.populated & (1u64 << grade_b) == 0 {
                continue;
            }

            // Check if ANY result grade from this pair could pass
            // Result grade is in [|ga-gb|, ga+gb], stepping by 2 (or ga^gb for geometric)
            let min_gr = grade_a.abs_diff(grade_b);
            let max_gr = (grade_a + grade_b).min(n);
            let mut any_pass = false;
            for gr in min_gr..=max_gr {
                if filter(grade_a, grade_b, gr) {
                    any_pass = true;
                    break;
                }
            }
            if !any_pass {
                continue;
            }

            // Process all term pairs in this grade combination
            for &(mask_a, idx_a) in &ga.groups[grade_a as usize] {
                let coeff_a = &ga.terms[idx_a].1;
                for &(mask_b, idx_b) in &gb.groups[grade_b as usize] {
                    let coeff_b = &gb.terms[idx_b].1;
                    let (mask_r, sign) = blade_product(mask_a, mask_b, sig);
                    if sign == 0 {
                        continue;
                    }
                    let gr = grade(mask_r);
                    if !filter(grade_a, grade_b, gr) {
                        continue;
                    }
                    let coeff = coeff_a.clone() * coeff_b.clone() * Scalar::from(sign as i64);
                    result.add_term(mask_r, coeff);
                }
            }
        }
    }
    result
}

/// Direct (no grouping) product for small Mvs.
fn product_filtered_direct<F>(a: &Mv, b: &Mv, sig: &Signature, filter: F) -> Mv
where
    F: Fn(u8, u8, u8) -> bool,
{
    let mut result = Mv::new();
    for (blade_a, coeff_a) in &a.terms {
        let mask_a = blade_to_mask(blade_a);
        for (blade_b, coeff_b) in &b.terms {
            let mask_b = blade_to_mask(blade_b);
            let (mask_r, sign) = blade_product(mask_a, mask_b, sig);
            if sign == 0 {
                continue;
            }
            let gr = grade(mask_r);
            if !filter(grade(mask_a), grade(mask_b), gr) {
                continue;
            }
            let coeff = coeff_a.clone() * coeff_b.clone() * Scalar::from(sign as i64);
            result.add_term(mask_r, coeff);
        }
    }
    result
}

/// Flat-accumulator product for small algebras (n ≤ 12).
/// Uses a pre-allocated array indexed by blade mask instead of BTreeMap.
fn product_filtered_flat<F>(
    ga: &GradeGroups,
    gb: &GradeGroups,
    sig: &Signature,
    n: u8,
    filter: F,
) -> Mv
where
    F: Fn(u8, u8, u8) -> bool,
{
    let dim = 1usize << n;
    let mut accum: Vec<Scalar> = vec![Scalar::from(0i64); dim];
    let mut nonzero_masks: Vec<BladeMask> = Vec::new();

    for grade_a in 0..=n {
        if ga.populated & (1u64 << grade_a) == 0 {
            continue;
        }
        for grade_b in 0..=n {
            if gb.populated & (1u64 << grade_b) == 0 {
                continue;
            }

            let min_gr = grade_a.abs_diff(grade_b);
            let max_gr = (grade_a + grade_b).min(n);
            let mut any_pass = false;
            for gr in min_gr..=max_gr {
                if filter(grade_a, grade_b, gr) {
                    any_pass = true;
                    break;
                }
            }
            if !any_pass {
                continue;
            }

            for &(mask_a, idx_a) in &ga.groups[grade_a as usize] {
                let coeff_a = &ga.terms[idx_a].1;
                for &(mask_b, idx_b) in &gb.groups[grade_b as usize] {
                    let coeff_b = &gb.terms[idx_b].1;
                    let (mask_r, sign) = blade_product(mask_a, mask_b, sig);
                    if sign == 0 {
                        continue;
                    }
                    let gr = grade(mask_r);
                    if !filter(grade_a, grade_b, gr) {
                        continue;
                    }
                    let coeff = coeff_a.clone() * coeff_b.clone() * Scalar::from(sign as i64);
                    let idx = mask_r as usize;
                    let was_zero = accum[idx].is_zero();
                    accum[idx] = accum[idx].clone() + coeff;
                    if was_zero && !accum[idx].is_zero() {
                        nonzero_masks.push(mask_r);
                    }
                }
            }
        }
    }

    // Convert back to sparse Mv
    let mut result = Mv::new();
    for mask in nonzero_masks {
        let idx = mask as usize;
        if !accum[idx].is_zero() {
            result.add_term(mask, std::mem::replace(&mut accum[idx], Scalar::from(0i64)));
        }
    }
    result
}

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

    fn sig_vga3() -> Signature {
        Signature::new(0, 0, 3).unwrap()
    }

    fn e(k: u8) -> Mv {
        Mv::generator(k)
    }

    #[test]
    fn geometric_e0_e1() {
        let sig = sig_vga3();
        let r = geometric(&e(0), &e(1), &sig);
        assert_eq!(r.coefficient(0b11), Scalar::from(1));
        assert_eq!(r.len(), 1);
    }

    #[test]
    fn geometric_e1_e0() {
        let sig = sig_vga3();
        let r = geometric(&e(1), &e(0), &sig);
        assert_eq!(r.coefficient(0b11), Scalar::from(-1));
    }

    #[test]
    fn geometric_associativity() {
        let sig = sig_vga3();
        let a = e(0);
        let b = e(1);
        let c = e(2);
        let ab_c = geometric(&geometric(&a, &b, &sig), &c, &sig);
        let a_bc = geometric(&a, &geometric(&b, &c, &sig), &sig);
        assert_eq!(ab_c, a_bc);
    }

    #[test]
    fn outer_anticommutative() {
        let sig = sig_vga3();
        let ab = outer(&e(0), &e(1), &sig);
        let ba = outer(&e(1), &e(0), &sig);
        assert_eq!(ab, -ba);
    }

    #[test]
    fn outer_self_is_zero() {
        let sig = sig_vga3();
        let r = outer(&e(0), &e(0), &sig);
        assert!(r.is_zero());
    }

    #[test]
    fn outer_produces_bivector() {
        let sig = sig_vga3();
        let r = outer(&e(0), &e(1), &sig);
        assert_eq!(r.coefficient(0b11), Scalar::from(1));
        assert_eq!(r.len(), 1);
        assert_eq!(grade(0b11), 2);
    }

    #[test]
    fn inner_bivector_vector() {
        // e₀₁ · e₁ in VGA(3)
        let sig = sig_vga3();
        let e01 = outer(&e(0), &e(1), &sig);
        let r = inner(&e01, &e(1), &sig);
        // e₀₁ · e₁ = e₀ (grade |2-1|=1)
        assert_eq!(r.len(), 1);
        assert_eq!(r.terms[0].0.grade(), 1);
    }

    #[test]
    fn reverse_grade0() {
        let s = Mv::scalar(Scalar::from(5));
        assert_eq!(reverse(&s), s);
    }

    #[test]
    fn reverse_grade1() {
        assert_eq!(reverse(&e(0)), e(0));
    }

    #[test]
    fn reverse_grade2() {
        let sig = sig_vga3();
        let e01 = outer(&e(0), &e(1), &sig);
        let r = reverse(&e01);
        assert_eq!(r, -e01);
    }

    #[test]
    fn reverse_involution() {
        // reverse(reverse(x)) == x
        let sig = sig_vga3();
        let e01 = outer(&e(0), &e(1), &sig);
        assert_eq!(reverse(&reverse(&e01)), e01);
    }

    #[test]
    fn grade_involution_involution() {
        let x = e(0);
        assert_eq!(grade_involution(&grade_involution(&x)), x);
    }

    #[test]
    fn norm_squared_unit_vector() {
        let sig = sig_vga3();
        let ns = norm_squared(&e(0), &sig);
        assert_eq!(ns, Scalar::from(1));
    }

    #[test]
    fn norm_squared_negative() {
        // In Cl(1,0,0): g₀² = -1, so norm² of g₀ = -1
        let sig = Signature::new(1, 0, 0).unwrap();
        let ns = norm_squared(&e(0), &sig);
        assert_eq!(ns, Scalar::from(-1));
    }

    #[test]
    fn norm_squared_vector() {
        let sig = sig_vga3();
        // v = 3e₀ + 4e₁ + 5e₂, norm² = 9 + 16 + 25 = 50
        let v = Mv::from_rat_terms(&[
            (0b001, Rat::from(3)),
            (0b010, Rat::from(4)),
            (0b100, Rat::from(5)),
        ]);
        let ns = norm_squared(&v, &sig);
        assert_eq!(ns, Scalar::from(50));
    }

    #[test]
    fn dual_undual_roundtrip() {
        let sig = sig_vga3();
        let v = e(0);
        // undual(dual(x)) = x * ps⁻¹ * ps = x
        let roundtrip = undual(&dual(&v, &sig), &sig);
        assert_eq!(roundtrip, v);
    }

    #[test]
    fn sandwich_preserves_grade() {
        let sig = sig_vga3();
        // Rotate e₀ by the rotor e^(π/4 e₀₁) ≈ cos(π/4) + sin(π/4) e₀₁
        // For exact test, use a simple reflection: sandwich(e₀, e₁) = e₀ e₁ e₀ = -e₁
        let r = sandwich(&e(0), &e(1), &sig);
        // e₀ * e₁ * e₀ *= e₁ * e₀
        // = (e₀₁) * e₀ = e₀₁₀ → reorder: e₀ cancels with g₀²=+1, leaving -e₁
        assert_eq!(r.coefficient(0b010), Scalar::from(-1));
    }

    #[test]
    fn scalar_product_orthogonal() {
        let sig = sig_vga3();
        let sp = scalar_product_value(&e(0), &e(1), &sig);
        assert_eq!(sp, Scalar::from(0));
    }

    #[test]
    fn scalar_product_parallel() {
        let sig = sig_vga3();
        let sp = scalar_product_value(&e(0), &e(0), &sig);
        assert_eq!(sp, Scalar::from(1));
    }

    #[test]
    fn commutator_of_same() {
        let sig = sig_vga3();
        let c = commutator(&e(0), &e(0), &sig);
        assert!(c.is_zero()); // [a, a] = 0
    }

    #[test]
    fn grade_project_extracts() {
        let v = Mv::from_rat_terms(&[
            (0b000, Rat::from(1)), // scalar
            (0b001, Rat::from(2)), // vector
            (0b011, Rat::from(3)), // bivector
        ]);
        let g0 = grade_project(&v, 0);
        assert_eq!(g0.coefficient(0), Scalar::from(1));
        assert_eq!(g0.len(), 1);

        let g1 = grade_project(&v, 1);
        assert_eq!(g1.coefficient(0b001), Scalar::from(2));
        assert_eq!(g1.len(), 1);
    }

    // ─── Regressive product ───

    #[test]
    fn regressive_identity_vga3() {
        // a ∨ b should equal undual(dual(a) ∧ dual(b))
        // Test with pseudoscalar ∨ anything = that thing (up to scale)
        let sig = sig_vga3();
        let ps = Mv::term(sig.pseudoscalar_mask(), Scalar::from(1));
        let v = e(0);
        let result = regressive(&ps, &v, &sig);
        // pseudoscalar ∨ vector: dual(ps) = ±1, dual(v) = ±bivector
        // 1 ∧ bivector = bivector, undual(bivector) = ±vector
        assert_eq!(result.len(), 1);
        // Result should be proportional to e(0)
        let g = grade(result.blades().next().unwrap().0);
        assert_eq!(g, 1);
    }

    #[test]
    fn regressive_two_bivectors_vga3() {
        // In VGA3: two bivectors meet in a vector (or zero)
        // e₀₁ ∨ e₀₂: share e₀, so meet should be proportional to e₀
        let sig = sig_vga3();
        let a = Mv::term(0b011, Scalar::from(1)); // e₀₁
        let b = Mv::term(0b101, Scalar::from(1)); // e₀₂
        let result = regressive(&a, &b, &sig);
        // In VGA3 (n=3): grade(a∨b) = grade(a) + grade(b) - n = 2+2-3 = 1
        if !result.is_zero() {
            let g = grade(result.blades().next().unwrap().0);
            assert_eq!(g, 1, "regressive of two bivectors should be grade 1");
        }
    }

    #[test]
    fn regressive_dual_identity() {
        // The regressive product identity: a ∨ b = undual(dual(a) ∧ dual(b))
        // We verify this holds by comparing with the direct computation
        let sig = sig_vga3();
        let a = Mv::from_rat_terms(&[(0b011, Rat::from(2))]); // 2e₀₁
        let b = Mv::from_rat_terms(&[(0b110, Rat::from(3))]); // 3e₁₂
        let direct = regressive(&a, &b, &sig);
        let manual = undual(&outer(&dual(&a, &sig), &dual(&b, &sig), &sig), &sig);
        assert_eq!(direct, manual);
    }

    // ─── Sparse optimization correctness ───

    #[test]
    fn sparse_product_matches_for_many_terms() {
        // Build Mvs with enough terms to trigger grade-grouped path (>16 pairs)
        let sig = Signature::new(0, 0, 5).unwrap(); // 32-blade algebra
        let a = Mv::from_rat_terms(&[
            (0b00001, Rat::from(1)),
            (0b00010, Rat::from(2)),
            (0b00100, Rat::from(3)),
            (0b01000, Rat::from(4)),
            (0b10000, Rat::from(5)),
        ]);
        let b = Mv::from_rat_terms(&[
            (0b00001, Rat::from(1)),
            (0b00010, Rat::from(-1)),
            (0b00100, Rat::from(2)),
            (0b01000, Rat::from(-2)),
        ]);
        // Geometric product should be correct regardless of optimization path
        let geo = geometric(&a, &b, &sig);
        // Verify associativity: (a*b)*c == a*(b*c)
        let c = Mv::from_rat_terms(&[(0b00011, Rat::from(1))]);
        let ab_c = geometric(&geo, &c, &sig);
        let bc = geometric(&b, &c, &sig);
        let a_bc = geometric(&a, &bc, &sig);
        assert_eq!(ab_c, a_bc, "associativity broken in sparse path");
    }

    #[test]
    fn outer_product_sparse_5gen() {
        // Outer product of 5 vectors should produce the pseudoscalar
        let sig = Signature::new(0, 0, 5).unwrap();
        let mut result = e(0);
        for k in 1..5u8 {
            result = outer(&result, &e(k), &sig);
        }
        // Result should be the pseudoscalar (mask = 0b11111 = 31)
        assert_eq!(result.len(), 1);
        assert_eq!(result.blades().next().unwrap().0, 0b11111);
    }
}