symplex 0.22.3

//! Gröbner basis computation via Buchberger's algorithm with Gebauer-Möller criteria.
//!
//! Computes reduced Gröbner bases for polynomial ideals over ℚ,
//! with FGLM order conversion for back-substitution solving.

use crate::poly::multipoly::*;
use num_bigint::BigInt;
use num_rational::Ratio;
use num_traits::{One, Zero};
use std::collections::{BTreeSet, HashMap, HashSet, VecDeque};

// ═══════════════════════════════════════════════════════════════════════════
// Wave G3: Buchberger's Algorithm
// ═══════════════════════════════════════════════════════════════════════════

/// Compute a reduced Gröbner basis for the ideal generated by `polys`.
///
/// Uses Buchberger's algorithm with Gebauer-Möller criteria for pair elimination.
/// The result is a reduced Gröbner basis: monic, inter-reduced, with no leading
/// term dividing another.
///
/// # Performance
/// For systems with > 6 variables or degree > 6, computation may be slow.
/// Consider using grevlex ordering for computation and FGLM for conversion to lex.
pub fn groebner_basis<O: MonomialOrd>(polys: &[MultiPoly<O>]) -> Vec<MultiPoly<O>> {
    // Filter zero polynomials
    let nonzero: Vec<MultiPoly<O>> = polys.iter().filter(|p| !p.is_zero()).cloned().collect();
    if nonzero.is_empty() {
        return vec![];
    }

    // Check if any input is a nonzero constant — ideal is the whole ring
    for p in &nonzero {
        if p.total_degree() == Some(0) {
            let n = p.num_vars();
            return vec![MultiPoly::from_int(n, 1)];
        }
    }

    buchberger_core(&nonzero)
}

/// Core Buchberger algorithm with Gebauer-Möller pair elimination.
fn buchberger_core<O: MonomialOrd>(input: &[MultiPoly<O>]) -> Vec<MultiPoly<O>> {
    // Pre-reduce and make monic
    let mut basis: Vec<MultiPoly<O>> = pre_reduce(input);

    if basis.is_empty() {
        return vec![];
    }

    // Check again for constant after pre-reduction
    for p in &basis {
        if p.total_degree() == Some(0) {
            let n = p.num_vars();
            return vec![MultiPoly::from_int(n, 1)];
        }
    }

    // Critical pairs: set of (i, j) indices into basis, i < j
    let mut pairs: BTreeSet<(usize, usize)> = BTreeSet::new();
    for i in 0..basis.len() {
        for j in (i + 1)..basis.len() {
            pairs.insert((i, j));
        }
    }

    // Track which pairs have been processed (for chain criterion)
    let mut processed: HashSet<(usize, usize)> = HashSet::new();

    while let Some((i, j)) = select_pair(&pairs, &basis) {
        pairs.remove(&(i, j));

        // Safety: indices must be valid
        if i >= basis.len() || j >= basis.len() {
            continue;
        }

        let lm_i = match basis[i].leading_monomial() {
            Some(lm) => lm.to_vec(),
            None => {
                processed.insert((i, j));
                continue;
            }
        };
        let lm_j = match basis[j].leading_monomial() {
            Some(lm) => lm.to_vec(),
            None => {
                processed.insert((i, j));
                continue;
            }
        };

        // Product criterion: if LMs are coprime, S-poly reduces to 0
        if monomial_coprime(&lm_i, &lm_j) {
            processed.insert((i, j));
            continue;
        }

        // Chain criterion: skip if exists k where LM(k) divides LCM(LM(i), LM(j))
        // and both (i,k) and (j,k) have been processed
        let lcm_ij = monomial_lcm(&lm_i, &lm_j);
        let chain_skip = (0..basis.len()).any(|k| {
            if k == i || k == j {
                return false;
            }
            if let Some(lm_k) = basis[k].leading_monomial() {
                if monomial_divides(lm_k, &lcm_ij) {
                    let pair_ik = if i < k { (i, k) } else { (k, i) };
                    let pair_jk = if j < k { (j, k) } else { (k, j) };
                    processed.contains(&pair_ik) && processed.contains(&pair_jk)
                } else {
                    false
                }
            } else {
                false
            }
        });

        if chain_skip {
            processed.insert((i, j));
            continue;
        }

        // Compute S-polynomial and reduce
        let s = s_polynomial(&basis[i], &basis[j]);
        let basis_refs: Vec<&MultiPoly<O>> = basis.iter().collect();
        let remainder = s.reduce(&basis_refs);

        processed.insert((i, j));

        if !remainder.is_zero() {
            let h = remainder.monic().primitive_part_q();

            // Check for constant — if so, ideal is the whole ring
            if h.total_degree() == Some(0) {
                let n = h.num_vars();
                return vec![MultiPoly::from_int(n, 1)];
            }

            let h_idx = basis.len();

            // Remove obsolete pairs where new element makes them redundant
            update_pairs(&mut pairs, &basis, &h);

            basis.push(h);

            // Add new pairs with all existing basis elements
            for k in 0..h_idx {
                pairs.insert((k, h_idx));
            }
        }
    }

    // Inter-reduce to get fully reduced basis
    inter_reduce(&mut basis);

    // Sort for deterministic output (by leading monomial, descending)
    basis.sort_by(|a, b| {
        let lm_a = a.leading_monomial().unwrap_or(&[]);
        let lm_b = b.leading_monomial().unwrap_or(&[]);
        O::cmp_exponents(lm_b, lm_a)
    });

    basis
}

/// Select the critical pair with minimum LCM of leading monomials (normal selection strategy).
fn select_pair<O: MonomialOrd>(
    pairs: &BTreeSet<(usize, usize)>,
    basis: &[MultiPoly<O>],
) -> Option<(usize, usize)> {
    pairs
        .iter()
        .filter_map(|&(i, j)| {
            // Pairs whose members are out of range or zero are skipped.
            let lm_i = basis.get(i)?.leading_monomial()?;
            let lm_j = basis.get(j)?.leading_monomial()?;
            Some(((i, j), monomial_lcm(lm_i, lm_j)))
        })
        .min_by(|(_, lcm1), (_, lcm2)| O::cmp_exponents(lcm1, lcm2))
        .map(|(pair, _)| pair)
}

/// Pre-reduce a set of polynomials: reduce each against predecessors, repeat until stable.
fn pre_reduce<O: MonomialOrd>(input: &[MultiPoly<O>]) -> Vec<MultiPoly<O>> {
    let mut result: Vec<MultiPoly<O>> = input
        .iter()
        .filter(|p| !p.is_zero())
        .map(|p| p.monic())
        .collect();

    loop {
        let mut changed = false;
        let mut new_result: Vec<MultiPoly<O>> = Vec::new();
        for item in &result {
            // Reduce against predecessors already in new_result
            let others: Vec<&MultiPoly<O>> = new_result.iter().collect();
            let reduced = if others.is_empty() {
                item.clone()
            } else {
                item.reduce(&others)
            };
            if reduced.is_zero() {
                changed = true; // element reduced to zero, removed
            } else {
                let r = reduced.monic();
                if r != *item {
                    changed = true;
                }
                new_result.push(r);
            }
        }
        result = new_result;
        if !changed {
            break;
        }
    }
    result
}

/// Inter-reduce a basis: remove elements whose LM is divisible by another's LM,
/// then reduce each element modulo the rest.
fn inter_reduce<O: MonomialOrd>(basis: &mut Vec<MultiPoly<O>>) {
    // Remove zero polynomials
    basis.retain(|p| !p.is_zero());

    // Remove redundant elements: those whose LM is divisible by another element's LM
    let mut i = 0;
    while i < basis.len() {
        let lm_i = match basis[i].leading_monomial() {
            Some(lm) => lm.to_vec(),
            None => {
                basis.remove(i);
                continue;
            }
        };
        let redundant = (0..basis.len()).any(|j| {
            if j == i {
                return false;
            }
            match basis[j].leading_monomial() {
                Some(lm_j) => monomial_divides(lm_j, &lm_i) && lm_j != lm_i,
                None => false,
            }
        });
        if redundant {
            basis.remove(i);
        } else {
            i += 1;
        }
    }

    // If duplicate LMs exist, keep only the first
    let mut seen_lms: HashSet<Vec<u32>> = HashSet::new();
    basis.retain(|p| {
        if let Some(lm) = p.leading_monomial() {
            seen_lms.insert(lm.to_vec())
        } else {
            false
        }
    });

    // Reduce each element modulo the others
    for i in 0..basis.len() {
        let others: Vec<MultiPoly<O>> = basis
            .iter()
            .enumerate()
            .filter(|&(j, _)| j != i)
            .map(|(_, p)| p.clone())
            .collect();
        let refs: Vec<&MultiPoly<O>> = others.iter().collect();
        if !refs.is_empty() {
            let reduced = basis[i].reduce(&refs);
            if !reduced.is_zero() {
                basis[i] = reduced.monic();
            }
        }
    }

    // Remove any polynomials that became zero after reduction
    basis.retain(|p| !p.is_zero());
}

/// Update pairs when a new polynomial h is about to be added (Gebauer-Möller pair update).
///
/// Removes old pairs (g1, g2) where LM(h) divides LCM(LM(g1), LM(g2))
/// and the new pairs involving h would make the old pair redundant.
fn update_pairs<O: MonomialOrd>(
    pairs: &mut BTreeSet<(usize, usize)>,
    basis: &[MultiPoly<O>],
    h: &MultiPoly<O>,
) {
    let lm_h = match h.leading_monomial() {
        Some(lm) => lm,
        None => return,
    };

    let pairs_to_remove: Vec<(usize, usize)> = pairs
        .iter()
        .filter(|&&(i, j)| {
            if i >= basis.len() || j >= basis.len() {
                return false;
            }
            let lm_i = match basis[i].leading_monomial() {
                Some(lm) => lm,
                None => return false,
            };
            let lm_j = match basis[j].leading_monomial() {
                Some(lm) => lm,
                None => return false,
            };
            let lcm_ij = monomial_lcm(lm_i, lm_j);
            if monomial_divides(lm_h, &lcm_ij) {
                // Check that the new pairs (h, g_i) and (h, g_j)
                // have strictly smaller LCMs than the old pair (g_i, g_j)
                let lcm_ih = monomial_lcm(lm_i, lm_h);
                let lcm_jh = monomial_lcm(lm_j, lm_h);
                lcm_ih != lcm_ij && lcm_jh != lcm_ij
            } else {
                false
            }
        })
        .copied()
        .collect();

    for pair in pairs_to_remove {
        pairs.remove(&pair);
    }
}

/// Verify that a set of polynomials forms a Gröbner basis for the ideal they generate.
///
/// Checks that all S-polynomials reduce to zero modulo the basis (Buchberger's criterion).
pub fn is_groebner_basis<O: MonomialOrd>(basis: &[MultiPoly<O>]) -> bool {
    if basis.is_empty() {
        return true;
    }
    let refs: Vec<&MultiPoly<O>> = basis.iter().collect();
    for i in 0..basis.len() {
        for j in (i + 1)..basis.len() {
            if basis[i].is_zero() || basis[j].is_zero() {
                continue;
            }
            let s = s_polynomial(&basis[i], &basis[j]);
            let r = s.reduce(&refs);
            if !r.is_zero() {
                return false;
            }
        }
    }
    true
}

// ═══════════════════════════════════════════════════════════════════════════
// Wave G4: FGLM Order Conversion
// ═══════════════════════════════════════════════════════════════════════════

/// Check if an ideal (represented by its Gröbner basis) is zero-dimensional.
///
/// An ideal is zero-dimensional iff for each variable, there exists a basis element
/// whose leading monomial is a pure power of that variable.
pub fn is_zero_dimensional<O: MonomialOrd>(basis: &[MultiPoly<O>]) -> bool {
    if basis.is_empty() {
        return false;
    }
    let n = basis[0].num_vars();

    // For each variable, check if some basis element has LM = x_i^k
    for var in 0..n {
        let has_pure_power = basis.iter().any(|p| {
            if let Some(lm) = p.leading_monomial() {
                // Check: only variable `var` has nonzero exponent, and it's positive
                lm.iter().enumerate().all(|(i, &e)| i == var || e == 0) && lm[var] > 0
            } else {
                false
            }
        });
        if !has_pure_power {
            return false;
        }
    }
    true
}

/// Enumerate the standard monomials (staircase) of a zero-dimensional ideal.
///
/// These are monomials NOT divisible by any leading monomial of the basis.
/// Returns them in an arbitrary order.
fn standard_monomials<O: MonomialOrd>(basis: &[MultiPoly<O>], num_vars: usize) -> Vec<Vec<u32>> {
    let leading_monomials: Vec<Vec<u32>> = basis
        .iter()
        .filter_map(|p| p.leading_monomial().map(|m| m.to_vec()))
        .collect();

    let mut staircase = Vec::new();
    let mut queue: VecDeque<Vec<u32>> = VecDeque::new();
    queue.push_back(vec![0u32; num_vars]); // Start with the constant monomial
    let mut visited: HashSet<Vec<u32>> = HashSet::new();

    while let Some(mono) = queue.pop_front() {
        if visited.contains(&mono) {
            continue;
        }
        visited.insert(mono.clone());

        // Check if this monomial is divisible by any leading monomial
        let is_divisible = leading_monomials
            .iter()
            .any(|lm| monomial_divides(lm, &mono));

        if !is_divisible {
            staircase.push(mono.clone());

            // Generate successors: multiply by each variable
            for var in 0..num_vars {
                let mut next = mono.clone();
                next[var] += 1;
                if !visited.contains(&next) {
                    queue.push_back(next);
                }
            }
        }
        // If divisible, don't add successors (all multiples are also divisible)
    }

    staircase
}

/// Convert a reduced Gröbner basis from ordering `From` to ordering `To` using FGLM.
///
/// Only works for zero-dimensional ideals (finitely many solutions).
/// Returns `None` if the ideal is not zero-dimensional.
///
/// # Algorithm
/// 1. Enumerate standard monomials in `From` ordering
/// 2. Build multiplication matrices for the quotient algebra
/// 3. Iterate monomials in `To` order, checking for linear dependence
/// 4. Dependence relations give the `To`-ordering basis elements
pub fn fglm<From: MonomialOrd, To: MonomialOrd>(
    basis: &[MultiPoly<From>],
) -> Option<Vec<MultiPoly<To>>> {
    if basis.is_empty() {
        return Some(vec![]);
    }
    let n = basis[0].num_vars();

    if !is_zero_dimensional(basis) {
        return None;
    }

    // Step 1: Enumerate standard monomials in the `From` ordering
    let staircase = standard_monomials(basis, n);
    let d = staircase.len(); // Quotient space dimension

    if d == 0 {
        // Ideal is the whole ring
        return Some(vec![MultiPoly::from_int(n, 1)]);
    }

    // Create index map: monomial → position in staircase
    let staircase_idx: HashMap<Vec<u32>, usize> = staircase
        .iter()
        .enumerate()
        .map(|(i, m)| (m.clone(), i))
        .collect();

    // Step 2: Build multiplication matrices
    // For each variable x_i, compute the D×D matrix M_i where
    // column j = normal_form(x_i * staircase[j]) expressed in the staircase basis
    let basis_refs: Vec<&MultiPoly<From>> = basis.iter().collect();

    let mut mult_matrices: Vec<Vec<Vec<Ratio<BigInt>>>> = Vec::with_capacity(n);
    for var in 0..n {
        let mut matrix = vec![vec![Ratio::<BigInt>::zero(); d]; d]; // D×D
        for j in 0..d {
            // Compute x_var * staircase[j]
            let mut product_exp = staircase[j].clone();
            product_exp[var] += 1;
            let product = MultiPoly::<From>::monomial(Ratio::one(), product_exp);

            // Reduce modulo basis to get normal form
            let nf = product.reduce(&basis_refs);

            // Express normal form in staircase basis
            for (key, coeff) in nf.terms() {
                if let Some(&idx) = staircase_idx.get(key) {
                    matrix[idx][j] = coeff.clone();
                }
                // If key is not in staircase, the reduction should have eliminated it.
                // If not, it means the basis isn't a proper GB, but we trust the input.
            }
        }
        mult_matrices.push(matrix);
    }

    // Step 3: FGLM iteration
    // Process monomials in ascending `To` order, compute their normal forms
    // via multiplication matrices, and check for linear dependence.
    let mut new_basis: Vec<MultiPoly<To>> = Vec::new();
    let mut to_staircase: Vec<Vec<u32>> = Vec::new();
    let mut echelon = IncrementalEchelon::new(d);

    // Cache: monomial exponent → normal form vector (coordinates in `From` staircase)
    let mut nf_cache: HashMap<Vec<u32>, Vec<Ratio<BigInt>>> = HashMap::new();

    // Normal form of the constant monomial 1
    let zero_mono = vec![0u32; n];
    let mut nf_one = vec![Ratio::<BigInt>::zero(); d];
    if let Some(&idx) = staircase_idx.get(&zero_mono) {
        nf_one[idx] = Ratio::one();
    }
    nf_cache.insert(zero_mono.clone(), nf_one);

    // Priority queue of candidates in `To` order (ascending: smallest first via BTreeSet)
    let mut candidates: BTreeSet<MonoKey<To>> = BTreeSet::new();
    candidates.insert(MonoKey::new(zero_mono));

    let mut visited_to: HashSet<Vec<u32>> = HashSet::new();

    while let Some(mono_key) = candidates.iter().next().cloned() {
        candidates.remove(&mono_key);
        let mono = mono_key.exponents.clone();

        if visited_to.contains(&mono) {
            continue;
        }
        visited_to.insert(mono.clone());

        // Check if mono is divisible by any LT of new_basis — if so, skip
        let is_in_ideal = new_basis.iter().any(|p| {
            if let Some(lm) = p.leading_monomial() {
                monomial_divides(lm, &mono)
            } else {
                false
            }
        });
        if is_in_ideal {
            continue; // All multiples are also in the ideal
        }

        // Compute normal form of mono via multiplication matrices
        let nf = compute_nf_via_matrices(&mono, &mult_matrices, &nf_cache, n, d);
        nf_cache.insert(mono.clone(), nf.clone());

        // Check linear independence with existing To-staircase NFs
        match echelon.add(&nf) {
            EchelonResult::Independent => {
                // Linearly independent: add to To-staircase
                to_staircase.push(mono.clone());

                // Add successors (multiply by each variable)
                for var in 0..n {
                    let mut next = mono.clone();
                    next[var] += 1;
                    if !visited_to.contains(&next) {
                        candidates.insert(MonoKey::<To>::new(next));
                    }
                }
            }
            EchelonResult::Dependent(coeffs) => {
                // Linear dependency found: construct new basis element
                // mono = Σ coeffs[i] * to_staircase[i] in the quotient
                // So: mono - Σ coeffs[i] * to_staircase[i] is in the ideal
                let mut new_poly = MultiPoly::<To>::monomial(Ratio::one(), mono.clone());
                for (i, c) in coeffs.iter().enumerate() {
                    if !c.is_zero() {
                        let term = MultiPoly::<To>::monomial(-c.clone(), to_staircase[i].clone());
                        new_poly = new_poly.add(&term);
                    }
                }
                new_basis.push(new_poly.monic());
                // Don't add successors — they're all in the ideal
            }
        }

        // Early termination: if we've found enough basis elements and staircase
        // has reached full dimension D
        if to_staircase.len() >= d {
            // We have D independent monomials; all further monomials will be dependent.
            // But we still need to find the remaining basis elements.
            // Process remaining candidates to collect all basis elements,
            // but only those not already divisible by known LTs.
            // Continue the loop naturally.
        }
    }

    // Final inter-reduction of the new basis
    inter_reduce(&mut new_basis);

    // Sort for deterministic output
    new_basis.sort_by(|a, b| {
        let lm_a = a.leading_monomial().unwrap_or(&[]);
        let lm_b = b.leading_monomial().unwrap_or(&[]);
        To::cmp_exponents(lm_b, lm_a)
    });

    Some(new_basis)
}

/// Compute the normal form of a monomial using multiplication matrices and cached NFs.
fn compute_nf_via_matrices(
    mono: &[u32],
    mult_matrices: &[Vec<Vec<Ratio<BigInt>>>],
    nf_cache: &HashMap<Vec<u32>, Vec<Ratio<BigInt>>>,
    n: usize,
    d: usize,
) -> Vec<Ratio<BigInt>> {
    // If already cached, return it
    if let Some(nf) = nf_cache.get(mono) {
        return nf.clone();
    }

    // Find a parent: some variable var such that mono / x_var is in the cache
    for var in 0..n {
        if mono[var] > 0 {
            let mut parent = mono.to_vec();
            parent[var] -= 1;
            if let Some(parent_nf) = nf_cache.get(&parent) {
                // NF(mono) = M_var * NF(parent)
                return matrix_vector_mul(&mult_matrices[var], parent_nf, d);
            }
        }
    }

    // Fallback: compute by applying multiplication matrices one variable at a time
    // Start from the constant monomial NF
    let zero_mono = vec![0u32; n];
    let base_nf = match nf_cache.get(&zero_mono) {
        Some(nf) => nf.clone(),
        None => {
            // This shouldn't happen — the constant monomial should always be cached
            let mut nf = vec![Ratio::<BigInt>::zero(); d];
            if d > 0 {
                nf[0] = Ratio::one();
            }
            nf
        }
    };

    let mut current_nf = base_nf;
    for var in 0..n {
        for _ in 0..mono[var] {
            current_nf = matrix_vector_mul(&mult_matrices[var], &current_nf, d);
        }
    }

    current_nf
}

/// Multiply a D×D matrix by a D-vector.
fn matrix_vector_mul(
    matrix: &[Vec<Ratio<BigInt>>],
    vec: &[Ratio<BigInt>],
    d: usize,
) -> Vec<Ratio<BigInt>> {
    let mut result = vec![Ratio::<BigInt>::zero(); d];
    for i in 0..d {
        for j in 0..d {
            if !vec[j].is_zero() && !matrix[i][j].is_zero() {
                result[i] = &result[i] + &(&matrix[i][j] * &vec[j]);
            }
        }
    }
    result
}

// ═══════════════════════════════════════════════════════════════════════════
// Incremental Gaussian Elimination for FGLM
// ═══════════════════════════════════════════════════════════════════════════

/// Result of adding a vector to the incremental echelon form.
enum EchelonResult {
    /// The vector is linearly independent of existing vectors.
    Independent,
    /// The vector is linearly dependent; coefficients express it as
    /// a linear combination of the previously added independent vectors.
    Dependent(Vec<Ratio<BigInt>>),
}

/// Maintains a set of vectors in row echelon form and can incrementally
/// check whether a new vector is in their span.
///
/// Each echelon row stores:
///   - `pivot`: the pivot column in the vector part
///   - `vec_part`: the reduced vector (length d)
///   - `coeff_part`: expresses the echelon row as Σ `coeff_part[i] * v_i`
///     where v_i are the original independent input vectors
///
/// When a new vector is checked:
///   - If dependent: returns coefficients c such that nf = Σ `c[i] * v_i`
///   - If independent: adds it as a new echelon row
type EchelonRow = (usize, Vec<Ratio<BigInt>>, Vec<Ratio<BigInt>>);

struct IncrementalEchelon {
    d: usize,
    /// (pivot_column, vec_part, coeff_part)
    rows: Vec<EchelonRow>,
    /// Number of independent vectors added so far.
    num_independent: usize,
}

impl IncrementalEchelon {
    fn new(d: usize) -> Self {
        Self {
            d,
            rows: Vec::new(),
            num_independent: 0,
        }
    }

    /// Try to add a new vector. Returns whether it's independent or dependent.
    ///
    /// For a dependent vector, returns coefficients c_0, ..., c_{m-1} such that
    /// nf = c_0 * v_0 + c_1 * v_1 + ... + c_{m-1} * v_{m-1}
    /// where v_i are the previously added independent vectors.
    fn add(&mut self, nf: &[Ratio<BigInt>]) -> EchelonResult {
        let m = self.num_independent;
        let d = self.d;

        let mut v = nf.to_vec();
        // trail tracks: nf = Σ trail[i] * v_i + v_remaining
        // When we reduce v by subtracting factor * echelon_row, we accumulate
        // trail += factor * row.coeff_part (since echelon_row = Σ coeff_part[i] * v_i)
        let mut trail = vec![Ratio::<BigInt>::zero(); m];

        // Reduce v using existing echelon rows
        for (pivot_col, row_vec, row_coeff) in &self.rows {
            if !v[*pivot_col].is_zero() {
                let factor = v[*pivot_col].clone() / &row_vec[*pivot_col];
                // v -= factor * row_vec
                for k in 0..d {
                    if !row_vec[k].is_zero() {
                        let sub = &factor * &row_vec[k];
                        v[k] = &v[k] - &sub;
                    }
                }
                // trail += factor * row_coeff
                for k in 0..row_coeff.len().min(trail.len()) {
                    if !row_coeff[k].is_zero() {
                        let add = &factor * &row_coeff[k];
                        trail[k] = &trail[k] + &add;
                    }
                }
            }
        }

        // v reduced to zero: dependent, nf = Σ trail[i] * v_i.
        // Otherwise its first non-zero entry is the pivot of a new echelon row.
        match v.iter().position(|c| !c.is_zero()) {
            None => EchelonResult::Dependent(trail),
            Some(pivot_col) => {
                // Independent: this is the m-th independent vector

                // The echelon row = nf - Σ trail[i] * v_i = v_remaining
                // Express it as: echelon_row = (-trail[0])*v_0 + ... + (-trail[m-1])*v_{m-1} + 1*v_m
                // where v_m is this new independent vector (nf itself)
                let mut coeff: Vec<Ratio<BigInt>> = trail.iter().map(|t| -t.clone()).collect();
                coeff.push(Ratio::one()); // coefficient for v_m = this vector

                self.rows.push((pivot_col, v, coeff));
                self.num_independent += 1;
                EchelonResult::Independent
            }
        }
    }
}

// ═══════════════════════════════════════════════════════════════════════════
// Convenience: compute in grevlex, convert to lex
// ═══════════════════════════════════════════════════════════════════════════

/// Compute a lex-order Gröbner basis by first computing in grevlex and
/// then converting via FGLM (if zero-dimensional) or direct Buchberger (otherwise).
///
/// This is the recommended approach for solving polynomial systems:
/// grevlex computation is fast, and FGLM conversion to lex yields
/// a triangular system suitable for back-substitution.
pub fn groebner_basis_lex(polys: &[MultiPoly<GrevLex>]) -> Vec<MultiPoly<Lex>> {
    let grevlex_gb = groebner_basis(polys);
    if grevlex_gb.is_empty() {
        return vec![];
    }

    // Try FGLM first
    if let Some(lex_gb) = fglm::<GrevLex, Lex>(&grevlex_gb) {
        return lex_gb;
    }

    // Fallback: convert to lex and compute directly
    let lex_polys: Vec<MultiPoly<Lex>> = polys.iter().map(|p| p.convert_order()).collect();
    groebner_basis(&lex_polys)
}

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

    fn r(n: i64) -> Ratio<BigInt> {
        Ratio::from_integer(BigInt::from(n))
    }

    #[test]
    fn test_pre_reduce_simple() {
        // {x + y, x} should pre-reduce to {y, x} or {x, y}
        let x = MultiPoly::<GrevLex>::var(2, 0);
        let y = MultiPoly::<GrevLex>::var(2, 1);
        let _one = MultiPoly::<GrevLex>::from_int(2, 1);
        let p1 = &x + &y;
        let p2 = x.clone();
        let result = pre_reduce(&[p1, p2]);
        assert!(!result.is_empty());
        // After pre-reduction, we should have basis for same ideal
    }

    #[test]
    fn test_inter_reduce_removes_redundant() {
        // {x^2, x} should reduce to {x}
        let x = MultiPoly::<GrevLex>::var(2, 0);
        let x2 = &x * &x;
        let mut basis = vec![x2, x];
        inter_reduce(&mut basis);
        assert_eq!(basis.len(), 1);
    }

    #[test]
    fn test_standard_monomials_simple() {
        // Ideal (x^2, y^2) has standard monomials {1, x, y, xy}
        let x2 = MultiPoly::<GrevLex>::monomial(Ratio::one(), vec![2, 0]);
        let y2 = MultiPoly::<GrevLex>::monomial(Ratio::one(), vec![0, 2]);
        let basis = vec![x2, y2];
        let sm = standard_monomials(&basis, 2);
        assert_eq!(sm.len(), 4);
    }

    #[test]
    fn test_echelon_independent() {
        let mut ech = IncrementalEchelon::new(3);
        let v1 = vec![r(1), r(0), r(0)];
        match ech.add(&v1) {
            EchelonResult::Independent => {}
            _ => panic!("should be independent"),
        }
        let v2 = vec![r(0), r(1), r(0)];
        match ech.add(&v2) {
            EchelonResult::Independent => {}
            _ => panic!("should be independent"),
        }
    }

    #[test]
    fn test_echelon_dependent() {
        let mut ech = IncrementalEchelon::new(2);
        let v1 = vec![r(1), r(0)];
        match ech.add(&v1) {
            EchelonResult::Independent => {}
            _ => panic!("should be independent"),
        }
        let v2 = vec![r(0), r(1)];
        match ech.add(&v2) {
            EchelonResult::Independent => {}
            _ => panic!("should be independent"),
        }
        // v3 = 3*v1 + 2*v2 = [3, 2]
        let v3 = vec![r(3), r(2)];
        match ech.add(&v3) {
            EchelonResult::Dependent(coeffs) => {
                assert_eq!(coeffs.len(), 2);
                assert_eq!(coeffs[0], r(3));
                assert_eq!(coeffs[1], r(2));
            }
            _ => panic!("should be dependent"),
        }
    }
}