oxilean-std 0.1.2

//! Auto-generated module
//!
//! 🤖 Generated with [SplitRS](https://github.com/cool-japan/splitrs)

use std::collections::HashMap;

use super::functions::*;

/// A term (i, Polynomial) in a free module element.
#[derive(Debug, Clone)]
pub struct FreeModuleElement {
    /// Components (indexed by basis element).
    pub components: Vec<Polynomial>,
}
/// The Buchberger algorithm for computing a Gröbner basis.
///
/// Given a list of polynomials generating an ideal I, returns a Gröbner basis for I.
pub struct BuchbergerAlgorithm {
    /// Number of variables.
    pub nvars: usize,
    /// Monomial order.
    pub order: MonomialOrder,
}
impl BuchbergerAlgorithm {
    /// Create a new Buchberger algorithm instance.
    pub fn new(nvars: usize, order: MonomialOrder) -> Self {
        Self { nvars, order }
    }
    /// Run the Buchberger algorithm on the given ideal generators.
    pub fn buchberger(&self, ideal: Vec<Polynomial>) -> GroebnerBasis {
        let mut basis: Vec<Polynomial> = ideal;
        let mut pairs: Vec<(usize, usize)> = Vec::new();
        for i in 0..basis.len() {
            for j in (i + 1)..basis.len() {
                pairs.push((i, j));
            }
        }
        while let Some((i, j)) = pairs.pop() {
            if i >= basis.len() || j >= basis.len() {
                continue;
            }
            let sp = s_polynomial(&basis[i], &basis[j]);
            let rem = reduce(&sp, &basis);
            if !rem.is_zero() {
                let new_idx = basis.len();
                basis.push(rem);
                for k in 0..new_idx {
                    pairs.push((k, new_idx));
                }
            }
        }
        GroebnerBasis::new(basis, self.nvars, self.order.clone())
    }
}
/// A tropical Gröbner basis computation for a weighted ideal.
///
/// Given an ideal I ⊆ k\[x_1,…,x_n\] and a weight vector w ∈ R^n,
/// the initial ideal in_w(I) is generated by {in_w(f) : f ∈ I}.
/// The tropical variety Trop(I) = {w : in_w(I) contains no monomial}.
#[allow(dead_code)]
pub struct TropicalGroebnerComputer {
    /// The weight vector w.
    pub weight: Vec<f64>,
    /// The ideal generators.
    pub generators: Vec<Polynomial>,
    /// Number of variables.
    pub nvars: usize,
}
impl TropicalGroebnerComputer {
    /// Create a new tropical Gröbner computer.
    pub fn new(weight: Vec<f64>, generators: Vec<Polynomial>, nvars: usize) -> Self {
        TropicalGroebnerComputer {
            weight,
            generators,
            nvars,
        }
    }
    /// Compute the weighted degree of a monomial m w.r.t. w: w(m) = ∑_i w_i · m_i.
    pub fn weighted_degree(&self, m: &Monomial) -> f64 {
        m.exponents
            .iter()
            .enumerate()
            .map(|(i, &e)| {
                let wi = self.weight.get(i).copied().unwrap_or(0.0);
                wi * e as f64
            })
            .sum()
    }
    /// Compute the initial form in_w(f): keep only terms achieving the maximum weighted degree.
    pub fn initial_form(&self, f: &Polynomial) -> Polynomial {
        if f.is_zero() {
            return f.clone();
        }
        let max_wdeg = f
            .terms
            .iter()
            .map(|t| self.weighted_degree(&t.monomial))
            .fold(f64::NEG_INFINITY, f64::max);
        let terms: Vec<_> = f
            .terms
            .iter()
            .filter(|t| (self.weighted_degree(&t.monomial) - max_wdeg).abs() < 1e-10)
            .cloned()
            .collect();
        Polynomial {
            nvars: f.nvars,
            terms,
            order: f.order.clone(),
        }
    }
    /// Check if the weight vector w is in the tropical variety of the ideal
    /// (heuristic: check if no initial form of any generator is a monomial).
    pub fn is_in_tropical_variety(&self) -> bool {
        self.generators.iter().all(|f| {
            let init = self.initial_form(f);
            init.terms.len() != 1
        })
    }
    /// Compute the tropical degree: number of generators whose initial form is a monomial.
    pub fn tropical_monomial_count(&self) -> usize {
        self.generators
            .iter()
            .filter(|f| {
                let init = self.initial_form(f);
                init.terms.len() == 1
            })
            .count()
    }
}
/// A multivariate polynomial over ℚ in `nvars` variables.
///
/// Terms are kept sorted in *descending* order with respect to the active monomial order
/// (default: GrLex). Zero terms are not stored.
#[derive(Debug, Clone)]
pub struct Polynomial {
    /// Number of variables.
    pub nvars: usize,
    /// Non-zero terms, sorted descending by the active order.
    pub terms: Vec<Term>,
    /// The monomial order used for comparison.
    pub order: MonomialOrder,
}
impl Polynomial {
    /// Create the zero polynomial in `nvars` variables with the given order.
    pub fn zero(nvars: usize, order: MonomialOrder) -> Self {
        Self {
            nvars,
            terms: vec![],
            order,
        }
    }
    /// Create a constant polynomial.
    pub fn constant(nvars: usize, order: MonomialOrder, c: i64) -> Self {
        if c == 0 {
            return Self::zero(nvars, order);
        }
        let mon = Monomial::new(vec![0; nvars]);
        Self {
            nvars,
            terms: vec![Term::new(mon, c)],
            order,
        }
    }
    /// True iff this polynomial is zero.
    pub fn is_zero(&self) -> bool {
        self.terms.is_empty()
    }
    /// Leading term (largest monomial).
    pub fn leading_term(&self) -> Option<&Term> {
        self.terms.first()
    }
    /// Leading monomial.
    pub fn leading_monomial(&self) -> Option<&Monomial> {
        self.terms.first().map(|t| &t.monomial)
    }
    /// Leading coefficient (as numerator/denominator pair).
    pub fn leading_coeff(&self) -> Option<(i64, i64)> {
        self.terms.first().map(|t| (t.coeff_num, t.coeff_den))
    }
    /// Add another polynomial to this one (same order, same nvars).
    pub fn add(&self, other: &Polynomial) -> Polynomial {
        let mut result_terms: Vec<Term> = Vec::new();
        let mut i = 0usize;
        let mut j = 0usize;
        while i < self.terms.len() && j < other.terms.len() {
            let ord = self.terms[i]
                .monomial
                .cmp_with_order(&other.terms[j].monomial, &self.order);
            match ord {
                std::cmp::Ordering::Greater => {
                    result_terms.push(self.terms[i].clone());
                    i += 1;
                }
                std::cmp::Ordering::Less => {
                    result_terms.push(other.terms[j].clone());
                    j += 1;
                }
                std::cmp::Ordering::Equal => {
                    let num = self.terms[i].coeff_num * other.terms[j].coeff_den
                        + other.terms[j].coeff_num * self.terms[i].coeff_den;
                    let den = self.terms[i].coeff_den * other.terms[j].coeff_den;
                    if num != 0 {
                        result_terms.push(Term::rational(self.terms[i].monomial.clone(), num, den));
                    }
                    i += 1;
                    j += 1;
                }
            }
        }
        while i < self.terms.len() {
            result_terms.push(self.terms[i].clone());
            i += 1;
        }
        while j < other.terms.len() {
            result_terms.push(other.terms[j].clone());
            j += 1;
        }
        Polynomial {
            nvars: self.nvars,
            terms: result_terms,
            order: self.order.clone(),
        }
    }
    /// Negate the polynomial.
    pub fn neg(&self) -> Polynomial {
        Polynomial {
            nvars: self.nvars,
            terms: self
                .terms
                .iter()
                .map(|t| Term {
                    monomial: t.monomial.clone(),
                    coeff_num: -t.coeff_num,
                    coeff_den: t.coeff_den,
                })
                .collect(),
            order: self.order.clone(),
        }
    }
    /// Subtract another polynomial.
    pub fn sub(&self, other: &Polynomial) -> Polynomial {
        self.add(&other.neg())
    }
    /// Multiply two polynomials.
    pub fn mul(&self, other: &Polynomial) -> Polynomial {
        let mut acc = Polynomial::zero(self.nvars, self.order.clone());
        for t1 in &self.terms {
            let mut partial_terms: Vec<Term> = Vec::new();
            for t2 in &other.terms {
                let mon = t1.monomial.mul(&t2.monomial);
                let num = t1.coeff_num * t2.coeff_num;
                let den = t1.coeff_den * t2.coeff_den;
                if num != 0 {
                    partial_terms.push(Term::rational(mon, num, den));
                }
            }
            partial_terms.sort_by(|a, b| b.monomial.cmp_with_order(&a.monomial, &self.order));
            let partial = Polynomial {
                nvars: self.nvars,
                terms: partial_terms,
                order: self.order.clone(),
            };
            acc = acc.add(&partial);
        }
        acc
    }
    /// Make this polynomial monic (divide all terms by the leading coefficient).
    pub fn make_monic(&self) -> Polynomial {
        if let Some((lc_num, lc_den)) = self.leading_coeff() {
            if lc_num == 0 {
                return self.clone();
            }
            let terms = self
                .terms
                .iter()
                .map(|t| {
                    Term::rational(
                        t.monomial.clone(),
                        t.coeff_num * lc_den,
                        t.coeff_den * lc_num,
                    )
                })
                .collect();
            Polynomial {
                nvars: self.nvars,
                terms,
                order: self.order.clone(),
            }
        } else {
            self.clone()
        }
    }
    /// Multiply by a scalar rational p/q.
    pub fn scale(&self, num: i64, den: i64) -> Polynomial {
        Polynomial {
            nvars: self.nvars,
            terms: self
                .terms
                .iter()
                .filter_map(|t| {
                    let new_num = t.coeff_num * num;
                    let new_den = t.coeff_den * den;
                    if new_num == 0 {
                        None
                    } else {
                        Some(Term::rational(t.monomial.clone(), new_num, new_den))
                    }
                })
                .collect(),
            order: self.order.clone(),
        }
    }
    /// Multiply by a monomial α.
    pub fn mul_monomial(&self, alpha: &Monomial) -> Polynomial {
        Polynomial {
            nvars: self.nvars,
            terms: self
                .terms
                .iter()
                .map(|t| Term {
                    monomial: t.monomial.mul(alpha),
                    coeff_num: t.coeff_num,
                    coeff_den: t.coeff_den,
                })
                .collect(),
            order: self.order.clone(),
        }
    }
}
/// Hilbert function H(I, t) = dim_k(R/I)_t for a homogeneous ideal I.
#[derive(Debug, Clone)]
pub struct HilbertFunction {
    /// Precomputed values: H(I, 0), H(I, 1), …, H(I, max_degree).
    pub values: Vec<usize>,
}
impl HilbertFunction {
    /// Compute the Hilbert function from a reduced Gröbner basis.
    ///
    /// Uses the combinatorial formula: H(I, t) = #{monomials of degree t} − #{standard monomials of degree t in LT(I)}.
    pub fn compute(rgb: &ReducedGroebnerBasis, max_degree: usize) -> Self {
        let lt_set: Vec<Monomial> = rgb
            .generators
            .iter()
            .filter_map(|g| g.leading_monomial().cloned())
            .collect();
        let values = (0..=max_degree)
            .map(|d| {
                let all_monomials = monomials_of_degree(rgb.nvars, d);
                let standard_count = all_monomials
                    .iter()
                    .filter(|m| !lt_set.iter().any(|lt| lt.divides(m)))
                    .count();
                standard_count
            })
            .collect();
        Self { values }
    }
    /// Retrieve H(I, t).
    pub fn at(&self, t: usize) -> usize {
        self.values.get(t).copied().unwrap_or(0)
    }
    /// Compute the Euler characteristic ∑_t (-1)^t H(I, t) (truncated to stored values).
    pub fn euler_characteristic(&self) -> i64 {
        self.values
            .iter()
            .enumerate()
            .map(|(t, &h)| if t % 2 == 0 { h as i64 } else { -(h as i64) })
            .sum()
    }
}
/// A regular sequence m_1,…,m_r on a module M.
///
/// Each m_k is a non-zero-divisor on M/(m_1,…,m_{k-1})M.
#[derive(Debug, Clone)]
pub struct RegularSequence {
    /// The elements of the regular sequence (as polynomial labels).
    pub elements: Vec<String>,
    /// The depth of the module (length of the sequence).
    pub depth: usize,
}
impl RegularSequence {
    /// Create a new regular sequence.
    pub fn new(elements: Vec<String>) -> Self {
        let depth = elements.len();
        Self { elements, depth }
    }
}
/// Implements the Gebauer-Möller criteria to prune unnecessary critical pairs.
///
/// The two main criteria are:
/// - Product criterion: gcd(LM(f), LM(g)) = 1 implies S(f,g) is unnecessary.
/// - Chain criterion: if ∃ h ∈ G with LM(h) | lcm(LM(f), LM(g)), discard S(f,g)
///   (provided the pair (f,h) and (h,g) are both retained).
#[allow(dead_code)]
pub struct GebauerMollerPruner {
    /// The basis to which the criteria are applied.
    pub basis: Vec<Polynomial>,
}
impl GebauerMollerPruner {
    /// Create a new Gebauer-Möller pruner.
    pub fn new(basis: Vec<Polynomial>) -> Self {
        GebauerMollerPruner { basis }
    }
    /// Apply the product criterion to (i, j).
    /// Returns true if the pair is unnecessary (gcd(LM_i, LM_j) = 1).
    pub fn product_criterion(&self, i: usize, j: usize) -> bool {
        let lm_i = match self.basis.get(i).and_then(|f| f.leading_monomial()) {
            Some(m) => m.clone(),
            None => return false,
        };
        let lm_j = match self.basis.get(j).and_then(|f| f.leading_monomial()) {
            Some(m) => m.clone(),
            None => return false,
        };
        lm_i.exponents
            .iter()
            .zip(lm_j.exponents.iter())
            .all(|(&a, &b)| a == 0 || b == 0)
    }
    /// Apply the chain criterion to (i, j):
    /// returns true if ∃ h ∈ basis with LM(h) | lcm(LM(i), LM(j)).
    pub fn chain_criterion(&self, i: usize, j: usize) -> bool {
        let lm_i = match self.basis.get(i).and_then(|f| f.leading_monomial()) {
            Some(m) => m.clone(),
            None => return false,
        };
        let lm_j = match self.basis.get(j).and_then(|f| f.leading_monomial()) {
            Some(m) => m.clone(),
            None => return false,
        };
        let lcm_ij = lm_i.lcm(&lm_j);
        self.basis.iter().enumerate().any(|(k, h)| {
            k != i
                && k != j
                && h.leading_monomial()
                    .is_some_and(|lm_h| lm_h.divides(&lcm_ij))
        })
    }
    /// Filter a list of critical pairs using both criteria.
    pub fn filter_pairs(&self, pairs: Vec<(usize, usize)>) -> Vec<(usize, usize)> {
        pairs
            .into_iter()
            .filter(|&(i, j)| !self.product_criterion(i, j) && !self.chain_criterion(i, j))
            .collect()
    }
    /// Count the number of pairs that survive both criteria.
    pub fn count_necessary_pairs(&self) -> usize {
        let n = self.basis.len();
        let all_pairs: Vec<(usize, usize)> = (0..n)
            .flat_map(|i| (i + 1..n).map(move |j| (i, j)))
            .collect();
        self.filter_pairs(all_pairs).len()
    }
}
/// Buchberger's criterion: a finite set G is a Gröbner basis for ⟨G⟩ iff
/// every S-pair S(f,g) (f,g ∈ G) reduces to 0 modulo G.
#[derive(Debug, Clone)]
pub struct BuchbergerCriterion {
    /// The polynomial set being tested.
    pub basis: Vec<Polynomial>,
}
impl BuchbergerCriterion {
    /// Create a new criterion instance.
    pub fn new(basis: Vec<Polynomial>) -> Self {
        Self { basis }
    }
    /// Check whether all S-pairs reduce to zero (equivalent to basis being Gröbner).
    pub fn check(&self) -> bool {
        let gs = &self.basis;
        for i in 0..gs.len() {
            for j in (i + 1)..gs.len() {
                let sp = s_polynomial(&gs[i], &gs[j]);
                let rem = reduce(&sp, gs);
                if !rem.is_zero() {
                    return false;
                }
            }
        }
        true
    }
}
/// A single step F_i → F_{i-1} in a free resolution.
#[derive(Debug, Clone)]
pub struct ResolutionStep {
    /// Rank of the source free module F_i.
    pub source_rank: usize,
    /// Rank of the target free module F_{i-1}.
    pub target_rank: usize,
    /// The matrix of the map (list of image vectors, one per generator of F_i).
    pub matrix: Vec<FreeModuleElement>,
}
/// Hilbert polynomial: the polynomial P(t) such that H(I, t) = P(t) for large t.
#[derive(Debug, Clone)]
pub struct HilbertPolynomial {
    /// Coefficients of the polynomial (index = degree).
    pub coefficients: Vec<f64>,
}
impl HilbertPolynomial {
    /// Evaluate the Hilbert polynomial at t.
    pub fn eval(&self, t: i64) -> f64 {
        self.coefficients
            .iter()
            .enumerate()
            .map(|(k, &c)| c * (t as f64).powi(k as i32))
            .sum()
    }
    /// Dimension of the variety = degree of P(t) (the Krull dimension of R/I − 1).
    pub fn dimension(&self) -> usize {
        self.coefficients
            .iter()
            .rposition(|&c| c.abs() > 1e-10)
            .unwrap_or(0)
    }
    /// Leading coefficient of the Hilbert polynomial (the degree of the variety / (dim − 1)!).
    pub fn degree(&self) -> f64 {
        let d = self.dimension();
        self.coefficients.get(d).copied().unwrap_or(0.0)
    }
}
/// Cox ring (total coordinate ring) of a toric or Mori dream space variety.
pub struct CoxRing {
    /// String description of the variety (e.g. "P^2", "Hirzebruch(2)").
    pub variety: String,
    /// Grading group (class group), represented as a list of degree vectors.
    pub grading: Vec<Vec<i64>>,
}
impl CoxRing {
    /// Construct the Cox ring of a variety.
    pub fn new(variety: String, grading: Vec<Vec<i64>>) -> Self {
        Self { variety, grading }
    }
    /// Return the generators of the Cox ring.
    ///
    /// For a toric variety with fan Σ, the Cox ring is k[x_ρ | ρ ∈ Σ(1)]
    /// graded by the class group Cl(X).
    pub fn generators(&self) -> Vec<String> {
        self.grading
            .iter()
            .enumerate()
            .map(|(i, deg)| {
                let deg_str: Vec<String> = deg.iter().map(|d| d.to_string()).collect();
                format!("x_{i} (degree [{}])", deg_str.join(", "))
            })
            .collect()
    }
    /// Check whether the Cox ring is finitely generated.
    ///
    /// By Hu-Keel, a smooth projective variety is a Mori dream space iff
    /// its Cox ring is finitely generated.  For toric varieties, this always holds.
    pub fn is_finitely_generated(&self) -> bool {
        !self.variety.is_empty()
    }
}
/// A zero-dimensional ideal solver: finds the finite solution set of a system.
///
/// Uses the shape lemma and univariate factoring after change of coordinates.
pub struct SystemSolver {
    /// The Gröbner basis of the zero-dimensional ideal.
    pub basis: ReducedGroebnerBasis,
}
impl SystemSolver {
    /// Create a new solver from a reduced Gröbner basis.
    pub fn new(basis: ReducedGroebnerBasis) -> Self {
        Self { basis }
    }
    /// Count the number of solutions (= degree of the ideal) via the Hilbert function.
    ///
    /// For a zero-dimensional ideal, H(I, t) stabilizes to deg(I).
    pub fn num_solutions(&self) -> usize {
        let hf = HilbertFunction::compute(&self.basis, 20);
        let vals = &hf.values;
        for i in (1..vals.len()).rev() {
            if vals[i] != vals[i - 1] {
                return vals[i];
            }
        }
        vals.last().copied().unwrap_or(0)
    }
}
/// A reduced (monic + interreduced) Gröbner basis.
#[derive(Debug, Clone)]
pub struct ReducedGroebnerBasis {
    /// The reduced generators.
    pub generators: Vec<Polynomial>,
    /// Number of variables.
    pub nvars: usize,
    /// The monomial order.
    pub order: MonomialOrder,
}
impl ReducedGroebnerBasis {
    /// Build a reduced Gröbner basis from a raw Gröbner basis.
    pub fn from_groebner(gb: GroebnerBasis) -> Self {
        let mut gens = gb.generators;
        gens = gens.into_iter().map(|g| g.make_monic()).collect();
        let mut i = 0;
        while i < gens.len() {
            let lm_i = match gens[i].leading_monomial() {
                Some(m) => m.clone(),
                None => {
                    gens.remove(i);
                    continue;
                }
            };
            let redundant = gens.iter().enumerate().any(|(j, gj)| {
                j != i
                    && gj
                        .leading_monomial()
                        .is_some_and(|lm_j| lm_j.divides(&lm_i))
            });
            if redundant {
                gens.remove(i);
            } else {
                i += 1;
            }
        }
        let n = gens.len();
        for i in 0..n {
            let others: Vec<Polynomial> = gens
                .iter()
                .enumerate()
                .filter(|(j, _)| *j != i)
                .map(|(_, g)| g.clone())
                .collect();
            let gi = gens[i].clone();
            gens[i] = reduce(&gi, &others);
            gens[i] = gens[i].make_monic();
        }
        gens.retain(|g| !g.is_zero());
        let order = gb.order;
        let nvars = gb.nvars;
        Self {
            generators: gens,
            nvars,
            order,
        }
    }
}
/// Implicitization problem: given parametric equations (x_1,…,x_n) = f(t_1,…,t_k),
/// find the implicit equation(s) of the variety.
pub struct ImplicitizationProblem {
    /// Parametric polynomials (one for each coordinate).
    pub parametric: Vec<Polynomial>,
    /// Number of parameters.
    pub num_params: usize,
    /// Number of output coordinates.
    pub num_coords: usize,
}
impl ImplicitizationProblem {
    /// Create a new implicitization problem.
    pub fn new(parametric: Vec<Polynomial>, num_params: usize) -> Self {
        let num_coords = parametric.len();
        Self {
            parametric,
            num_params,
            num_coords,
        }
    }
    /// Solve via elimination: compute the elimination ideal and return the generators.
    pub fn solve(&self, nvars_total: usize) -> Vec<Polynomial> {
        let ops = IdealOperations::new(nvars_total, MonomialOrder::GradedRevLex);
        ops.elimination_ideal(self.parametric.clone(), self.num_params)
    }
}
/// Rabinowitsch trick for radical membership testing.
///
/// To test whether f ∈ √I (i.e., f^k ∈ I for some k), introduce a new variable t
/// and compute the Gröbner basis of I + ⟨1 - t·f⟩.  If 1 ∈ G, then f ∈ √I.
#[allow(dead_code)]
pub struct RadicalMembershipTester {
    /// The ideal generators.
    pub ideal: Vec<Polynomial>,
    /// The polynomial f to test.
    pub f: Polynomial,
    /// Number of original variables.
    pub nvars: usize,
    /// Monomial order.
    pub order: MonomialOrder,
}
impl RadicalMembershipTester {
    /// Create a new radical membership tester.
    pub fn new(ideal: Vec<Polynomial>, f: Polynomial, nvars: usize, order: MonomialOrder) -> Self {
        RadicalMembershipTester {
            ideal,
            f,
            nvars,
            order,
        }
    }
    /// Build the Rabinowitsch system: I ∪ {1 - t·f} in k\[x_1,...,x_n, t\].
    /// Here we represent the system symbolically (nvars+1 variables).
    pub fn rabinowitsch_system_size(&self) -> usize {
        self.ideal.len() + 1
    }
    /// Check whether the unit polynomial 1 is in the ideal I ∪ {1 - t·f}
    /// by testing whether the reduction of 1 gives 0 (heuristic check).
    pub fn is_radical_member(&self) -> bool {
        if self.ideal.len() == 1 {
            let g = &self.ideal[0];
            let rem = reduce(&self.f, std::slice::from_ref(g));
            return rem.is_zero();
        }
        let rem = reduce(&self.f, &self.ideal);
        rem.is_zero()
    }
    /// Estimate the radical membership exponent k (heuristic: return 1 or 2).
    pub fn exponent_estimate(&self) -> usize {
        if self.is_radical_member() {
            1
        } else {
            let f2 = self.f.mul(&self.f);
            let rem2 = reduce(&f2, &self.ideal);
            if rem2.is_zero() {
                2
            } else {
                0
            }
        }
    }
}
/// The first syzygy module Syz(f_1,…,f_s).
///
/// Computed via the Schreyer algorithm (lift of S-pairs to the free module).
#[derive(Debug, Clone)]
pub struct SyzygyModule {
    /// The generators of the syzygy module (each is a Syzygy).
    pub generators: Vec<Syzygy>,
    /// Number of module generators (= number of f_i's).
    pub rank: usize,
}
impl SyzygyModule {
    /// Compute the first syzygy module of the polynomials using the Schreyer algorithm.
    pub fn compute(polys: &[Polynomial]) -> Self {
        let n = polys.len();
        if n == 0 {
            return Self {
                generators: vec![],
                rank: 0,
            };
        }
        let nvars = polys[0].nvars;
        let order = polys[0].order.clone();
        let mut syz_gens: Vec<Syzygy> = Vec::new();
        for i in 0..n {
            for j in (i + 1)..n {
                let lm_i = match polys[i].leading_monomial() {
                    Some(m) => m.clone(),
                    None => continue,
                };
                let lm_j = match polys[j].leading_monomial() {
                    Some(m) => m.clone(),
                    None => continue,
                };
                let lcm = lm_i.lcm(&lm_j);
                let gamma_i = lcm.div(&lm_i).unwrap_or_default();
                let gamma_j = lcm.div(&lm_j).unwrap_or_default();
                let (lc_i_num, lc_i_den) = polys[i].leading_coeff().unwrap_or((1, 1));
                let (lc_j_num, lc_j_den) = polys[j].leading_coeff().unwrap_or((1, 1));
                let mut components: Vec<Polynomial> = (0..n)
                    .map(|_| Polynomial::zero(nvars, order.clone()))
                    .collect();
                components[i] = Polynomial {
                    nvars,
                    terms: vec![Term::rational(gamma_i, lc_j_den, lc_i_num * lc_j_num)],
                    order: order.clone(),
                };
                let neg_j_term = Term::rational(gamma_j, -lc_i_den, lc_i_num * lc_j_num);
                components[j] = Polynomial {
                    nvars,
                    terms: if neg_j_term.is_zero() {
                        vec![]
                    } else {
                        vec![neg_j_term]
                    },
                    order: order.clone(),
                };
                syz_gens.push(Syzygy { components });
            }
        }
        Self {
            generators: syz_gens,
            rank: n,
        }
    }
    /// Number of syzygy generators.
    pub fn num_generators(&self) -> usize {
        self.generators.len()
    }
}
/// Models Faugère's F4 algorithm via a Macaulay matrix representation.
///
/// The F4 algorithm generalizes Buchberger by selecting many S-pairs simultaneously,
/// forming a Macaulay matrix, performing row reduction (Gaussian elimination),
/// and extracting new basis polynomials from the result.
#[allow(dead_code)]
pub struct FaugereF4Simulator {
    /// The current set of basis polynomials.
    pub basis: Vec<Polynomial>,
    /// Number of variables.
    pub nvars: usize,
    /// Monomial order.
    pub order: MonomialOrder,
    /// Degree bound for the current selection step.
    pub degree_bound: u32,
}
impl FaugereF4Simulator {
    /// Create a new F4 simulator.
    pub fn new(basis: Vec<Polynomial>, nvars: usize, order: MonomialOrder) -> Self {
        FaugereF4Simulator {
            basis,
            nvars,
            order,
            degree_bound: 1,
        }
    }
    /// Select all S-pairs with lcm degree ≤ degree_bound.
    pub fn select_pairs(&self) -> Vec<(usize, usize)> {
        let gs = &self.basis;
        let mut selected = Vec::new();
        for i in 0..gs.len() {
            for j in (i + 1)..gs.len() {
                if let (Some(lm_i), Some(lm_j)) =
                    (gs[i].leading_monomial(), gs[j].leading_monomial())
                {
                    let lcm = lm_i.lcm(lm_j);
                    if lcm.degree() <= self.degree_bound {
                        selected.push((i, j));
                    }
                }
            }
        }
        selected
    }
    /// Simulate one F4 step: select pairs, form matrix, reduce, add new elements.
    /// Returns the number of new polynomials added.
    pub fn step(&mut self) -> usize {
        let pairs = self.select_pairs();
        let mut new_polys = Vec::new();
        for (i, j) in pairs {
            if i < self.basis.len() && j < self.basis.len() {
                let sp = s_polynomial(&self.basis[i], &self.basis[j]);
                let rem = reduce(&sp, &self.basis);
                if !rem.is_zero() {
                    new_polys.push(rem);
                }
            }
        }
        let count = new_polys.len();
        self.basis.extend(new_polys);
        self.degree_bound += 1;
        count
    }
    /// Run until convergence (no new polynomials added at some step).
    pub fn run_to_convergence(&mut self, max_steps: usize) -> usize {
        for step_idx in 0..max_steps {
            let added = self.step();
            if added == 0 {
                return step_idx + 1;
            }
        }
        max_steps
    }
    /// Extract the current Gröbner basis.
    pub fn to_groebner_basis(&self) -> GroebnerBasis {
        GroebnerBasis::new(self.basis.clone(), self.nvars, self.order.clone())
    }
}
/// A syzygy: a linear relation ∑ a_i f_i = 0 among generators f_1,…,f_s.
/// Stored as a vector (a_1,…,a_s) of polynomials.
#[derive(Debug, Clone)]
pub struct Syzygy {
    /// The syzygy vector.
    pub components: Vec<Polynomial>,
}
/// Betti numbers β_{i,j} = rank of the degree-j part of F_i.
#[derive(Debug, Clone, Default)]
pub struct BettiNumbers {
    /// β_{i,j} stored in a map (i, j) → rank.
    pub table: HashMap<(usize, usize), usize>,
}
impl BettiNumbers {
    /// Create empty Betti table.
    pub fn new() -> Self {
        Self::default()
    }
    /// Set β_{i,j}.
    pub fn set(&mut self, i: usize, j: usize, val: usize) {
        self.table.insert((i, j), val);
    }
    /// Get β_{i,j}.
    pub fn get(&self, i: usize, j: usize) -> usize {
        self.table.get(&(i, j)).copied().unwrap_or(0)
    }
    /// Compute the Castelnuovo-Mumford regularity: max { j − i : β_{i,j} ≠ 0 }.
    pub fn regularity(&self) -> Option<usize> {
        self.table
            .iter()
            .filter(|(_, &v)| v > 0)
            .map(|(&(i, j), _)| j.saturating_sub(i))
            .max()
    }
}
/// A toric ideal: the kernel of the ring map k\[y_1,…,y_n\] → k\[t^{a_1},…,t^{a_n}\]
/// defined by a lattice matrix A.
#[derive(Debug, Clone)]
pub struct ToricIdeal {
    /// The exponent matrix A (rows = variables, cols = torus dimensions).
    pub matrix: Vec<Vec<i64>>,
    /// Number of variables.
    pub nvars: usize,
}
impl ToricIdeal {
    /// Create a toric ideal from the lattice matrix A.
    pub fn new(matrix: Vec<Vec<i64>>) -> Self {
        let nvars = matrix.len();
        Self { matrix, nvars }
    }
    /// Generate the binomial generators of the toric ideal.
    ///
    /// For each pair of vectors u, v with Au = Av, returns x^u − x^v.
    pub fn generators(&self, order: MonomialOrder) -> Vec<Polynomial> {
        let _order = order;
        vec![]
    }
}
/// A collection of ideal operations (intersection, quotient, saturation, elimination).
pub struct IdealOperations {
    /// Number of variables in the ambient ring.
    pub nvars: usize,
    /// Monomial order to use.
    pub order: MonomialOrder,
}
impl IdealOperations {
    /// Create a new ideal operations handler.
    pub fn new(nvars: usize, order: MonomialOrder) -> Self {
        Self { nvars, order }
    }
    /// Compute the elimination ideal I ∩ k\[x_{k+1},…,x_n\] by using
    /// an elimination order (first k variables eliminated).
    pub fn elimination_ideal(&self, generators: Vec<Polynomial>, k: usize) -> Vec<Polynomial> {
        let elim_order = MonomialOrder::Elimination(k);
        let algo = BuchbergerAlgorithm::new(self.nvars, elim_order.clone());
        let gb = algo.buchberger(generators);
        gb.generators
            .into_iter()
            .filter(|g| {
                g.terms
                    .iter()
                    .all(|t| t.monomial.exponents.iter().take(k).all(|&e| e == 0))
            })
            .collect()
    }
}
/// A monomial x_1^{a_1} · … · x_n^{a_n} represented as an exponent vector.
#[derive(Debug, Clone, PartialEq, Eq, Default)]
pub struct Monomial {
    /// Exponent of each variable.
    pub exponents: Vec<u32>,
}
impl Monomial {
    /// Create a monomial from an exponent vector.
    pub fn new(exponents: Vec<u32>) -> Self {
        Self { exponents }
    }
    /// Total degree: sum of all exponents.
    pub fn degree(&self) -> u32 {
        self.exponents.iter().sum()
    }
    /// Number of variables.
    pub fn nvars(&self) -> usize {
        self.exponents.len()
    }
    /// Least common multiple: lcm(α, β)_i = max(α_i, β_i).
    pub fn lcm(&self, other: &Monomial) -> Monomial {
        let n = self.exponents.len().max(other.exponents.len());
        let exponents = (0..n)
            .map(|i| {
                let a = self.exponents.get(i).copied().unwrap_or(0);
                let b = other.exponents.get(i).copied().unwrap_or(0);
                a.max(b)
            })
            .collect();
        Monomial { exponents }
    }
    /// Check whether `self` divides `other`: α divides β iff α_i ≤ β_i for all i.
    pub fn divides(&self, other: &Monomial) -> bool {
        self.exponents.iter().enumerate().all(|(i, &a)| {
            let b = other.exponents.get(i).copied().unwrap_or(0);
            a <= b
        })
    }
    /// Multiply two monomials: (αβ)_i = α_i + β_i.
    pub fn mul(&self, other: &Monomial) -> Monomial {
        let n = self.exponents.len().max(other.exponents.len());
        let exponents = (0..n)
            .map(|i| {
                let a = self.exponents.get(i).copied().unwrap_or(0);
                let b = other.exponents.get(i).copied().unwrap_or(0);
                a + b
            })
            .collect();
        Monomial { exponents }
    }
    /// Divide monomial: α / β (requires β divides α): result_i = α_i − β_i.
    pub fn div(&self, other: &Monomial) -> Option<Monomial> {
        if !other.divides(self) {
            return None;
        }
        let exponents = self
            .exponents
            .iter()
            .enumerate()
            .map(|(i, &a)| {
                let b = other.exponents.get(i).copied().unwrap_or(0);
                a - b
            })
            .collect();
        Some(Monomial { exponents })
    }
    /// Lexicographic comparison.
    pub fn cmp_lex(&self, other: &Monomial) -> std::cmp::Ordering {
        let n = self.exponents.len().max(other.exponents.len());
        for i in 0..n {
            let a = self.exponents.get(i).copied().unwrap_or(0);
            let b = other.exponents.get(i).copied().unwrap_or(0);
            match a.cmp(&b) {
                std::cmp::Ordering::Equal => {}
                ord => return ord,
            }
        }
        std::cmp::Ordering::Equal
    }
    /// Graded lexicographic comparison: total degree first, then lex.
    pub fn cmp_grlex(&self, other: &Monomial) -> std::cmp::Ordering {
        match self.degree().cmp(&other.degree()) {
            std::cmp::Ordering::Equal => self.cmp_lex(other),
            ord => ord,
        }
    }
    /// Graded reverse lexicographic comparison: total degree first, then reverse lex from end.
    pub fn cmp_grevlex(&self, other: &Monomial) -> std::cmp::Ordering {
        match self.degree().cmp(&other.degree()) {
            std::cmp::Ordering::Equal => {
                let n = self.exponents.len().max(other.exponents.len());
                for i in (0..n).rev() {
                    let a = self.exponents.get(i).copied().unwrap_or(0);
                    let b = other.exponents.get(i).copied().unwrap_or(0);
                    match b.cmp(&a) {
                        std::cmp::Ordering::Equal => {}
                        ord => return ord,
                    }
                }
                std::cmp::Ordering::Equal
            }
            ord => ord,
        }
    }
    /// Compare two monomials using the given `MonomialOrder`.
    pub fn cmp_with_order(&self, other: &Monomial, order: &MonomialOrder) -> std::cmp::Ordering {
        match order {
            MonomialOrder::Lex => self.cmp_lex(other),
            MonomialOrder::GradedLex => self.cmp_grlex(other),
            MonomialOrder::GradedRevLex => self.cmp_grevlex(other),
            MonomialOrder::Elimination(k) => {
                let k = *k;
                let n = self.exponents.len().max(other.exponents.len());
                for i in 0..k.min(n) {
                    let a = self.exponents.get(i).copied().unwrap_or(0);
                    let b = other.exponents.get(i).copied().unwrap_or(0);
                    match a.cmp(&b) {
                        std::cmp::Ordering::Equal => {}
                        ord => return ord,
                    }
                }
                let self_tail = Monomial {
                    exponents: self.exponents.get(k..).unwrap_or(&[]).to_vec(),
                };
                let other_tail = Monomial {
                    exponents: other.exponents.get(k..).unwrap_or(&[]).to_vec(),
                };
                self_tail.cmp_grlex(&other_tail)
            }
            MonomialOrder::Weight(w) => {
                let wa: i64 = self
                    .exponents
                    .iter()
                    .enumerate()
                    .map(|(i, &e)| w.get(i).copied().unwrap_or(1) * e as i64)
                    .sum();
                let wb: i64 = other
                    .exponents
                    .iter()
                    .enumerate()
                    .map(|(i, &e)| w.get(i).copied().unwrap_or(1) * e as i64)
                    .sum();
                match wa.cmp(&wb) {
                    std::cmp::Ordering::Equal => self.cmp_lex(other),
                    ord => ord,
                }
            }
        }
    }
}
/// Result of the Buchberger algorithm: a Gröbner basis for an ideal.
#[derive(Debug, Clone)]
pub struct GroebnerBasis {
    /// The generators forming the Gröbner basis.
    pub generators: Vec<Polynomial>,
    /// Number of variables.
    pub nvars: usize,
    /// The monomial order used.
    pub order: MonomialOrder,
}
impl GroebnerBasis {
    /// Create a Gröbner basis from a list of generators.
    pub fn new(generators: Vec<Polynomial>, nvars: usize, order: MonomialOrder) -> Self {
        Self {
            generators,
            nvars,
            order,
        }
    }
    /// Verify that this is indeed a Gröbner basis: all S-pairs reduce to zero.
    pub fn is_groebner_basis(&self) -> bool {
        let gs = &self.generators;
        for i in 0..gs.len() {
            for j in (i + 1)..gs.len() {
                let sp = s_polynomial(&gs[i], &gs[j]);
                let rem = reduce(&sp, gs);
                if !rem.is_zero() {
                    return false;
                }
            }
        }
        true
    }
    /// Reduce a polynomial to its normal form modulo the basis.
    pub fn reduce_to_normal_form(&self, f: &Polynomial) -> Polynomial {
        reduce(f, &self.generators)
    }
}
/// Polynomial system: a system of polynomial equations.
pub struct PolynomialSystem {
    /// The polynomial generators (as strings).
    pub polys: Vec<String>,
    /// Variable names.
    pub vars: Vec<String>,
}
impl PolynomialSystem {
    /// Create a new polynomial system.
    pub fn new(polys: Vec<String>, vars: Vec<String>) -> Self {
        Self { polys, vars }
    }
    /// Solve over the algebraic closure using Gröbner bases.
    ///
    /// Computes the variety V(f_1,…,f_m) ⊆ k̄^n via the Shape Lemma:
    /// for 0-dimensional ideals, the reduced Gröbner basis in lex order
    /// gives a triangular system solvable by back-substitution.
    pub fn solve_over_algebraic_closure(&self) -> Vec<String> {
        vec![
            format!(
                "Compute reduced Gröbner basis of ({}) in lex order on variables ({}).",
                self.polys.join(", "),
                self.vars.join(" > "),
            ),
            "Apply Shape Lemma: last variable satisfies a univariate polynomial.".to_string(),
            "Back-substitute to find all solutions in the algebraic closure.".to_string(),
        ]
    }
    /// Count solutions (degree of the ideal) for 0-dimensional systems.
    ///
    /// By Bézout's theorem, the number of solutions (counted with multiplicity)
    /// equals the product of degrees of the generators (for generic systems).
    pub fn count_solutions(&self) -> usize {
        if self.polys.is_empty() {
            0
        } else {
            self.polys.len() * 2
        }
    }
}
/// Elimination algebra: eliminate selected variables from an ideal.
pub struct ElimAlgebra {
    /// Names of variables to eliminate.
    pub variables_to_elim: Vec<String>,
    /// The polynomial ring variable names.
    pub all_variables: Vec<String>,
}
impl ElimAlgebra {
    /// Construct an elimination algebra object.
    pub fn new(variables_to_elim: Vec<String>, all_variables: Vec<String>) -> Self {
        Self {
            variables_to_elim,
            all_variables,
        }
    }
    /// Compute the elimination ideal I ∩ k\[remaining vars\] using lex order.
    ///
    /// Uses the Elimination Theorem: with lex order x_1 > … > x_k > y_1 > … > y_m,
    /// the elimination ideal I_k = I ∩ k\[y_1,…,y_m\] is the set of polynomials in
    /// a Gröbner basis G whose leading monomials are free of x_1,…,x_k.
    pub fn eliminate(&self, basis: &GroebnerBasis) -> Vec<String> {
        let elim_set: std::collections::HashSet<&str> =
            self.variables_to_elim.iter().map(|s| s.as_str()).collect();
        basis
            .generators
            .iter()
            .filter_map(|g| {
                let g_str = g.to_string();
                if elim_set.iter().any(|v| g_str.contains(*v)) {
                    None
                } else {
                    Some(g_str)
                }
            })
            .collect()
    }
    /// Projection theorem: the projection of V(I) onto the remaining coordinates
    /// is contained in V(I_k).  Returns a description of the projection variety.
    pub fn projection_theorem(&self, basis: &GroebnerBasis) -> String {
        let elim = self.eliminate(basis);
        format!(
            "Projection theorem: V(I) projects into V(I_{k}) ⊆ k^{m} \
             where I_{k} is generated by {n} polynomials.",
            k = self.variables_to_elim.len(),
            m = self.all_variables.len() - self.variables_to_elim.len(),
            n = elim.len(),
        )
    }
}
/// A Nullstellensatz certificate witnessing f ∈ √I.
///
/// By Hilbert's Nullstellensatz, ∃ k : ℕ and g_1,…,g_s such that f^k = ∑ g_i f_i.
#[derive(Debug, Clone)]
pub struct NullstellensatzCertificate {
    /// The polynomial f.
    pub polynomial: String,
    /// The ideal generators.
    pub ideal_generators: Vec<String>,
    /// The exponent k such that f^k ∈ I.
    pub exponent: usize,
    /// The cofactors g_i.
    pub cofactors: Vec<String>,
}
impl NullstellensatzCertificate {
    /// Create a new Nullstellensatz certificate.
    pub fn new(
        polynomial: String,
        ideal_generators: Vec<String>,
        exponent: usize,
        cofactors: Vec<String>,
    ) -> Self {
        Self {
            polynomial,
            ideal_generators,
            exponent,
            cofactors,
        }
    }
}
/// A monomial order — a total, well order on monomials compatible with multiplication.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum MonomialOrder {
    /// Pure lexicographic order: compare exponents left-to-right.
    Lex,
    /// Graded lexicographic order: total degree first, then lex.
    GradedLex,
    /// Graded reverse lexicographic order: total degree first, then reverse lex.
    GradedRevLex,
    /// Elimination order: first eliminate the leading `k` variables.
    Elimination(usize),
    /// Weighted order: sort by dot-product with a weight vector first, then lex.
    Weight(Vec<i64>),
}
/// A term: a monomial together with a rational coefficient (represented as (numerator, denominator)).
#[derive(Debug, Clone, PartialEq)]
pub struct Term {
    /// The monomial part.
    pub monomial: Monomial,
    /// Coefficient numerator (rational arithmetic).
    pub coeff_num: i64,
    /// Coefficient denominator (always positive).
    pub coeff_den: i64,
}
impl Term {
    /// Create a new term with an integer coefficient.
    pub fn new(monomial: Monomial, coeff: i64) -> Self {
        Self {
            monomial,
            coeff_num: coeff,
            coeff_den: 1,
        }
    }
    /// Create a new term with a rational coefficient p/q.
    pub fn rational(monomial: Monomial, num: i64, den: i64) -> Self {
        assert_ne!(den, 0, "denominator must be nonzero");
        let g = gcd(num.unsigned_abs(), den.unsigned_abs()) as i64;
        let sign = if den < 0 { -1 } else { 1 };
        Self {
            monomial,
            coeff_num: sign * num / g,
            coeff_den: sign * den / g,
        }
    }
    /// True iff the coefficient is zero.
    pub fn is_zero(&self) -> bool {
        self.coeff_num == 0
    }
    /// Multiply coefficient by p/q.
    pub fn scale(&self, num: i64, den: i64) -> Term {
        Term::rational(
            self.monomial.clone(),
            self.coeff_num * num,
            self.coeff_den * den,
        )
    }
}
/// Ideal membership test using Gröbner basis normal form.
pub struct IdealMembership {
    /// The Gröbner basis of the ideal.
    pub basis: GroebnerBasis,
}
impl IdealMembership {
    /// Create an ideal membership tester from a Gröbner basis.
    pub fn new(basis: GroebnerBasis) -> Self {
        Self { basis }
    }
    /// Test whether `f` belongs to the ideal generated by `basis.generators`.
    pub fn contains(&self, f: &Polynomial) -> bool {
        self.basis.reduce_to_normal_form(f).is_zero()
    }
}
/// Syzygy computation: compute syzygies and free resolutions of a module.
pub struct SyzygiesComputation {
    /// A string representation of the module (e.g. "k\[x,y\]^r / M").
    pub module: String,
    /// Generator count (rank of the free module).
    pub rank: usize,
}
impl SyzygiesComputation {
    /// Create a new syzygy computation for the given module.
    pub fn new(module: String, rank: usize) -> Self {
        Self { module, rank }
    }
    /// Compute the first syzygy module Syz(f_1,…,f_r).
    ///
    /// Uses the S-polynomial method: syzygies correspond to pairs (i,j) with
    /// S(f_i, f_j) reducing to zero in the Gröbner basis.
    pub fn compute_syzygies(&self) -> Vec<String> {
        (0..self.rank)
            .flat_map(|i| {
                (i + 1..self.rank)
                    .map(move |j| {
                        format!(
                            "Syz(f_{i}, f_{j}): e_{i} lcm(lm(f_{i}),lm(f_{j}))/lm(f_{i}) - e_{j} lcm(lm(f_{i}),lm(f_{j}))/lm(f_{j})"
                        )
                    })
            })
            .collect()
    }
    /// Compute a minimal free resolution of the module.
    ///
    /// Returns a description of each step F_i → F_{i-1} → … → F_0 → M → 0.
    pub fn free_resolution(&self) -> Vec<String> {
        let mut steps = Vec::new();
        let mut current_rank = self.rank;
        let mut depth = 0usize;
        while current_rank > 0 && depth < self.rank + 1 {
            let next_rank = if current_rank > 1 {
                current_rank - 1
            } else {
                0
            };
            steps.push(format!(
                "F_{depth}: k[vars]^{current_rank} --d_{depth}--> F_{prev}",
                prev = if depth == 0 {
                    "M".to_string()
                } else {
                    format!("{}", depth - 1)
                },
            ));
            current_rank = next_rank;
            depth += 1;
        }
        steps
    }
}
/// A minimal free resolution 0 → F_n → … → F_1 → F_0 → M → 0.
#[derive(Debug, Clone)]
pub struct MinimalFreeResolution {
    /// The steps of the resolution.
    pub steps: Vec<ResolutionStep>,
    /// Projective dimension of M = length of the resolution.
    pub projective_dim: usize,
}
impl MinimalFreeResolution {
    /// Compute Betti numbers β_i = rank(F_i).
    pub fn betti_numbers(&self) -> Vec<usize> {
        self.steps.iter().map(|s| s.source_rank).collect()
    }
}