geoit 0.0.2

#![allow(clippy::needless_range_loop)]
use crate::governance::poly::Poly;

/// Error from Gröbner basis computation.
#[derive(Clone, Debug)]
pub enum GroebnerError {
    /// Buchberger's algorithm did not terminate within the pair limit.
    IterationLimit {
        pairs_processed: usize,
        basis_size: usize,
    },
}

impl std::fmt::Display for GroebnerError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            GroebnerError::IterationLimit {
                pairs_processed,
                basis_size,
            } => write!(
                f,
                "Gröbner basis did not converge after {} pairs (basis size {})",
                pairs_processed, basis_size
            ),
        }
    }
}

/// Compute a Gröbner basis for the ideal generated by `polys`
/// using Buchberger's algorithm with grlex monomial ordering.
///
/// Returns `Err(GroebnerError::IterationLimit)` if the algorithm
/// does not converge within a generous bound.
pub fn groebner_basis(polys: Vec<Poly>) -> Result<Vec<Poly>, GroebnerError> {
    if polys.is_empty() {
        return Ok(vec![]);
    }

    let mut basis: Vec<Poly> = polys.into_iter().filter(|p| !p.is_zero()).collect();
    if basis.is_empty() {
        return Ok(vec![]);
    }

    let max_pairs = basis.len().max(1).pow(2) * 64;
    let mut pair_count: usize = 0;

    // Buchberger's algorithm: process all pairs
    let mut i = 0;
    while i < basis.len() {
        let mut j = i + 1;
        while j < basis.len() {
            pair_count += 1;
            if pair_count > max_pairs {
                return Err(GroebnerError::IterationLimit {
                    pairs_processed: pair_count,
                    basis_size: basis.len(),
                });
            }
            let s = basis[i].s_polynomial(&basis[j]);
            let remainder = reduce_by_set(&s, &basis);
            if !remainder.is_zero() {
                basis.push(remainder);
            }
            j += 1;
        }
        i += 1;
    }

    minimize_basis(&mut basis);
    interreduce_basis(&mut basis);
    Ok(basis)
}

/// Reduce a polynomial by a set of divisors.
/// Repeatedly tries to reduce by each divisor until no reduction is possible.
pub fn reduce_by_set(poly: &Poly, divisors: &[Poly]) -> Poly {
    let mut current = poly.clone();
    let mut changed = true;
    while changed {
        changed = false;
        for d in divisors {
            if d.is_zero() {
                continue;
            }
            let reduced = current.reduce_by(d);
            if reduced != current {
                current = reduced;
                changed = true;
                break; // restart
            }
        }
    }
    current
}

/// Minimize: remove any basis element whose leading monomial is
/// divisible by another basis element's leading monomial.
fn minimize_basis(basis: &mut Vec<Poly>) {
    let mut i = 0;
    while i < basis.len() {
        let lm_i = match basis[i].leading_monomial() {
            Some(m) => m.clone(),
            None => {
                let _ = basis.remove(i);
                continue;
            }
        };
        let mut redundant = false;
        for j in 0..basis.len() {
            if i == j {
                continue;
            }
            if let Some(lm_j) = basis[j].leading_monomial() {
                if Poly::monomial_divides(lm_j, &lm_i) && lm_j != &lm_i {
                    redundant = true;
                    break;
                }
            }
        }
        if redundant {
            let _ = basis.remove(i);
        } else {
            i += 1;
        }
    }

    // Make each element monic
    for p in basis.iter_mut() {
        p.make_monic();
    }
}

/// Inter-reduce: reduce each basis element by all others.
fn interreduce_basis(basis: &mut Vec<Poly>) {
    let n = basis.len();
    for i in 0..n {
        let mut others: Vec<Poly> = Vec::with_capacity(n - 1);
        for j in 0..n {
            if i != j {
                others.push(basis[j].clone());
            }
        }
        basis[i] = reduce_by_set(&basis[i], &others);
        if basis[i].is_zero() {
            continue;
        }
        basis[i].make_monic();
    }
    basis.retain(|p| !p.is_zero());
}

/// Identify free variables from a Gröbner basis.
///
/// A variable is "determined" if there exists a basis polynomial that
/// can be solved for that variable — meaning the variable appears at
/// degree 1 in some term, and the polynomial expresses it as a function
/// of other variables.
///
/// For linear polynomials: the leading variable is determined.
/// For mixed polynomials (quadratic + linear terms): the variable
/// appearing at degree 1 (not in the leading monomial) is determined,
/// because we can solve for it.
///
/// Variables not determined by any basis polynomial are free.
pub fn free_variables(basis: &[Poly], num_vars: usize) -> Vec<usize> {
    let mut determined = vec![false; num_vars];

    for p in basis {
        if p.is_zero() {
            continue;
        }

        // Strategy 1: if the polynomial is univariate (only one variable
        // appears across all terms), that variable is determined.
        let mut all_vars = vec![false; num_vars];
        for (exp, coeff) in &p.terms {
            if coeff.is_zero() {
                continue;
            }
            for (i, &e) in exp.iter().enumerate() {
                if e > 0 {
                    all_vars[i] = true;
                }
            }
        }
        let vars_present: Vec<usize> = all_vars
            .iter()
            .enumerate()
            .filter(|(_, &v)| v)
            .map(|(i, _)| i)
            .collect();
        if vars_present.len() == 1 {
            determined[vars_present[0]] = true;
            continue;
        }

        // Strategy 2: find a variable that appears as a pure linear term
        // (degree 1, only variable in that monomial). That variable can
        // be solved for in terms of the others.
        let mut solvable_var: Option<usize> = None;
        for (exp, coeff) in &p.terms {
            if coeff.is_zero() {
                continue;
            }
            let vars_in_mono: Vec<(usize, u8)> = exp
                .iter()
                .enumerate()
                .filter(|(_, &e)| e > 0)
                .map(|(i, &e)| (i, e))
                .collect();
            if vars_in_mono.len() == 1 && vars_in_mono[0].1 == 1 {
                let var = vars_in_mono[0].0;
                if solvable_var.map_or(true, |sv| var < sv) {
                    solvable_var = Some(var);
                }
            }
        }
        if let Some(var) = solvable_var {
            determined[var] = true;
        }
    }

    (0..num_vars).filter(|&i| !determined[i]).collect()
}

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

    #[test]
    fn linear_system() {
        // {x - 1, y - 2}: already a Gröbner basis
        let n = 2;
        let p1 = {
            let mut p = Poly::zero(n);
            p.terms.insert(vec![1, 0], Rat::ONE);
            p.terms.insert(vec![0, 0], -Rat::ONE);
            p
        };
        let p2 = {
            let mut p = Poly::zero(n);
            p.terms.insert(vec![0, 1], Rat::ONE);
            p.terms.insert(vec![0, 0], Rat::from(-2));
            p
        };
        let basis = groebner_basis(vec![p1, p2]).unwrap();
        assert_eq!(basis.len(), 2);
        // Free variables: none (both determined)
        assert!(free_variables(&basis, n).is_empty());
    }

    #[test]
    fn linear_system_2() {
        // {x + y - 1, x - y - 1} → x = 1, y = 0
        let n = 2;
        let p1 = {
            let mut p = Poly::zero(n);
            p.terms.insert(vec![1, 0], Rat::ONE);
            p.terms.insert(vec![0, 1], Rat::ONE);
            p.terms.insert(vec![0, 0], -Rat::ONE);
            p
        };
        let p2 = {
            let mut p = Poly::zero(n);
            p.terms.insert(vec![1, 0], Rat::ONE);
            p.terms.insert(vec![0, 1], -Rat::ONE);
            p.terms.insert(vec![0, 0], -Rat::ONE);
            p
        };
        let basis = groebner_basis(vec![p1, p2]).unwrap();
        assert!(free_variables(&basis, n).is_empty());
        // Verify: substituting x=1, y=0 gives zero for all basis elements
        let vals = vec![Rat::ONE, Rat::ZERO];
        for p in &basis {
            assert!(p.eval(&vals).is_zero(), "basis poly not zero at solution");
        }
    }

    #[test]
    fn quadratic_system() {
        // {x² + y² - 1, x - y} → should produce univariate in y
        let n = 2;
        let p1 = {
            let mut p = Poly::zero(n);
            p.terms.insert(vec![2, 0], Rat::ONE);
            p.terms.insert(vec![0, 2], Rat::ONE);
            p.terms.insert(vec![0, 0], -Rat::ONE);
            p
        };
        let p2 = {
            let mut p = Poly::zero(n);
            p.terms.insert(vec![1, 0], Rat::ONE);
            p.terms.insert(vec![0, 1], -Rat::ONE);
            p
        };
        let basis = groebner_basis(vec![p1, p2]).unwrap();
        // Should have a univariate polynomial in y (2y² - 1)
        // and x - y
        assert!(basis.len() >= 2);
        // Both variables should be determined
        assert!(free_variables(&basis, n).is_empty());
    }

    #[test]
    fn vga_constraints() {
        // VGA Vector: grade = 1 in Cl(0,0,3), 8 blades total.
        // Variables: c0 (scalar), c1 (g0), c2 (g1), c4 (g2),
        //           c3 (g01), c5 (g02), c6 (g12), c7 (g012)
        // GradeEquals(1) means: c0 = 0, c3 = 0, c5 = 0, c6 = 0, c7 = 0
        // Free: c1, c2, c4 (the grade-1 coefficients)
        let n = 8;
        let mut polys = Vec::new();
        // c0 = 0
        polys.push(Poly::variable(0, n));
        // c3 = 0
        polys.push(Poly::variable(3, n));
        // c5 = 0
        polys.push(Poly::variable(5, n));
        // c6 = 0
        polys.push(Poly::variable(6, n));
        // c7 = 0
        polys.push(Poly::variable(7, n));

        let basis = groebner_basis(polys).unwrap();
        let free = free_variables(&basis, n);
        // Free: indices 1, 2, 4 (the grade-1 blade coefficients)
        assert_eq!(free, vec![1, 2, 4]);
    }

    #[test]
    fn cga_like_constraints() {
        // Simplified CGA: 5 variables (c0..c4), two constraints:
        //   c0 + c1 = 1  (linear, from normalization)
        //   -c0² + c1² + c2² + c3² + c4² = 0  (quadratic, from null condition)
        let n = 5;
        let p_linear = {
            let mut p = Poly::zero(n);
            p.terms.insert(vec![1, 0, 0, 0, 0], Rat::ONE); // c0
            p.terms.insert(vec![0, 1, 0, 0, 0], Rat::ONE); // c1
            p.terms.insert(vec![0, 0, 0, 0, 0], -Rat::ONE); // -1
            p
        };
        let p_quadratic = {
            let mut p = Poly::zero(n);
            p.terms.insert(vec![2, 0, 0, 0, 0], -Rat::ONE); // -c0²
            p.terms.insert(vec![0, 2, 0, 0, 0], Rat::ONE); // c1²
            p.terms.insert(vec![0, 0, 2, 0, 0], Rat::ONE); // c2²
            p.terms.insert(vec![0, 0, 0, 2, 0], Rat::ONE); // c3²
            p.terms.insert(vec![0, 0, 0, 0, 2], Rat::ONE); // c4²
            p
        };
        let basis = groebner_basis(vec![p_linear, p_quadratic]).unwrap();

        // Free variables should be c2, c3, c4 (the Euclidean coordinates)
        let free = free_variables(&basis, n);
        assert!(free.contains(&2), "c2 should be free");
        assert!(free.contains(&3), "c3 should be free");
        assert!(free.contains(&4), "c4 should be free");
        assert!(!free.contains(&0), "c0 should be determined");
        assert!(!free.contains(&1), "c1 should be determined");

        // Verify: at c0=(1+50)/2=51/2, c1=(1-50)/2=-49/2, c2=3, c3=4, c4=5
        // both constraints should be zero
        let vals = vec![
            Rat::new(51, 2),
            Rat::new(-49, 2),
            Rat::from(3),
            Rat::from(4),
            Rat::from(5),
        ];
        for p in &basis {
            assert!(
                p.eval(&vals).is_zero(),
                "basis poly not zero at CGA point (3,4,5)"
            );
        }
    }

    #[test]
    fn empty_input() {
        let basis = groebner_basis(vec![]).unwrap();
        assert!(basis.is_empty());
    }

    #[test]
    fn all_zero_input() {
        let basis = groebner_basis(vec![Poly::zero(3), Poly::zero(3)]).unwrap();
        assert!(basis.is_empty());
    }

    #[test]
    fn single_polynomial() {
        let n = 2;
        let p = {
            let mut p = Poly::zero(n);
            p.terms.insert(vec![1, 0], Rat::from(2)); // 2x
            p.terms.insert(vec![0, 0], Rat::from(4)); // + 4
            p
        };
        let basis = groebner_basis(vec![p]).unwrap();
        assert_eq!(basis.len(), 1);
        // Should be monic: x + 2
        assert_eq!(basis[0].leading_coefficient(), Rat::ONE);
    }

    #[test]
    fn free_vars_underdetermined() {
        // One equation in 3 unknowns: x + y + z = 0
        let n = 3;
        let p = {
            let mut p = Poly::zero(n);
            p.terms.insert(vec![1, 0, 0], Rat::ONE);
            p.terms.insert(vec![0, 1, 0], Rat::ONE);
            p.terms.insert(vec![0, 0, 1], Rat::ONE);
            p
        };
        let basis = groebner_basis(vec![p]).unwrap();
        let free = free_variables(&basis, n);
        // x is determined (leading variable), y and z are free
        assert_eq!(free.len(), 2);
    }
}