ocas-poly 0.26.0

//! Ideal operations for polynomial rings.
//!
//! Provides fundamental ideal arithmetic based on Gröbner bases:
//! membership testing, sum, product, quotient, saturation, and intersection.
//!
//! All operations work over [`Lex`] ordering for consistency with
//! elimination-based computations.
//!
//! # References
//!
//! - Cox, Little, O'Shea: *Ideals, Varieties, and Algorithms*
//! - Adams & Loustaunau: *An Introduction to Gröbner Bases*

use ocas_domain::Domain;
use ocas_domain::{Rational, RationalDomain};

use crate::groebner::{Algorithm, GroebnerBasis, groebner_basis};
use crate::sparse::{Lex, SparseMultivariatePolynomial};

/// Test whether `f` belongs to the ideal generated by `generators`.
///
/// Computes a Gröbner basis of the ideal and reduces `f` against it.
/// `f ∈ I` iff the remainder is zero.
///
/// # Example
///
/// ```
/// use ocas_domain::{RationalDomain, Rational};
/// use ocas_poly::sparse::Lex;
/// use ocas_poly::ideal::ideal_contains;
/// use ocas_poly::{Algorithm, SparseMultivariatePolynomial};
///
/// let d = RationalDomain;
/// let x = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![1, 0], Rational::new(1, 1)),
/// ]);
/// let y = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![0, 1], Rational::new(1, 1)),
/// ]);
/// assert!(ideal_contains(&[x.clone(), y.clone()], &x, Algorithm::Auto));
/// // x ∉ ⟨y⟩
/// assert!(!ideal_contains(&[y], &x, Algorithm::Auto));
/// ```
pub fn ideal_contains<D: Domain + 'static>(
    generators: &[SparseMultivariatePolynomial<D, Lex>],
    f: &SparseMultivariatePolynomial<D, Lex>,
    algo: Algorithm,
) -> bool {
    if generators.is_empty() {
        return f.is_zero();
    }
    let gb = groebner_basis(generators, algo);
    let remainder = f.reduce(&gb.basis);
    remainder.is_zero()
}

/// Sum of two ideals: `I + J = ⟨f₁,…,fₘ, g₁,…,gₙ⟩`.
///
/// # Example
///
/// ```
/// use ocas_domain::{RationalDomain, Rational};
/// use ocas_poly::sparse::Lex;
/// use ocas_poly::ideal::ideal_sum;
/// use ocas_poly::{Algorithm, SparseMultivariatePolynomial};
///
/// let d = RationalDomain;
/// let x = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![1, 0], Rational::new(1, 1)),
/// ]);
/// let y = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![0, 1], Rational::new(1, 1)),
/// ]);
/// let gb = ideal_sum(&[x], &[y]);
/// // ⟨x⟩ + ⟨y⟩ = ⟨x, y⟩
/// assert!(gb.basis.len() >= 2);
/// ```
pub fn ideal_sum<D: Domain + 'static>(
    generators_a: &[SparseMultivariatePolynomial<D, Lex>],
    generators_b: &[SparseMultivariatePolynomial<D, Lex>],
) -> GroebnerBasis<D, Lex> {
    let mut combined = generators_a.to_vec();
    combined.extend(generators_b.iter().cloned());
    groebner_basis(&combined, Algorithm::Auto)
}

/// Product of two ideals: `I · J = ⟨fᵢ · gⱼ⟩`.
///
/// # Example
///
/// ```
/// use ocas_domain::{RationalDomain, Rational};
/// use ocas_poly::sparse::Lex;
/// use ocas_poly::ideal::ideal_product;
/// use ocas_poly::SparseMultivariatePolynomial;
///
/// let d = RationalDomain;
/// let x = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![1, 0], Rational::new(1, 1)),
/// ]);
/// let y = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![0, 1], Rational::new(1, 1)),
/// ]);
/// let gb = ideal_product(&[x], &[y]);
/// // ⟨x⟩ · ⟨y⟩ = ⟨xy⟩
/// assert_eq!(gb.basis.len(), 1);
/// ```
pub fn ideal_product<D: Domain + 'static>(
    generators_a: &[SparseMultivariatePolynomial<D, Lex>],
    generators_b: &[SparseMultivariatePolynomial<D, Lex>],
) -> GroebnerBasis<D, Lex> {
    let products: Vec<SparseMultivariatePolynomial<D, Lex>> = generators_a
        .iter()
        .flat_map(|f| generators_b.iter().map(move |g| f.mul(g)))
        .collect();
    groebner_basis(&products, Algorithm::Auto)
}

/// Ideal quotient: `I : J = {f : f · g ∈ I, ∀ g ∈ J}`.
///
/// Computed via the Rabinowitsch trick for each generator of J,
/// then intersecting the results.
///
/// For a single generator `g`: `I : g` is obtained by computing
/// `GB(I ∪ {1 - w·g})` in `k[x₁,…,xₙ, w]` and eliminating `w`.
///
/// # Example
///
/// ```
/// use ocas_domain::{RationalDomain, Rational};
/// use ocas_poly::sparse::Lex;
/// use ocas_poly::ideal::ideal_quotient;
/// use ocas_poly::SparseMultivariatePolynomial;
///
/// let d = RationalDomain;
/// // ⟨x², xy⟩ : ⟨x⟩ = ⟨x⟩
/// let f1 = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![2, 0], Rational::new(1, 1)),
/// ]);
/// let f2 = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![1, 1], Rational::new(1, 1)),
/// ]);
/// let g = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![1, 0], Rational::new(1, 1)),
/// ]);
/// let gb = ideal_quotient(&[f1, f2], &[g]);
/// // Result should generate the same ideal as ⟨x⟩
/// assert!(!gb.basis.is_empty());
/// ```
pub fn ideal_quotient<D: Domain + 'static>(
    generators_i: &[SparseMultivariatePolynomial<D, Lex>],
    generators_j: &[SparseMultivariatePolynomial<D, Lex>],
) -> GroebnerBasis<D, Lex> {
    if generators_i.is_empty() || generators_j.is_empty() {
        return GroebnerBasis { basis: vec![] };
    }

    // I : J = ∩_{g ∈ J} (I : g)
    let mut result: Option<Vec<SparseMultivariatePolynomial<D, Lex>>> = None;

    for g in generators_j {
        if g.is_zero() {
            continue;
        }
        let i_colon_g = quotient_single_generator(generators_i, g);

        if let Some(current) = result.take() {
            result = Some(intersect_generators(&current, &i_colon_g));
        } else {
            result = Some(i_colon_g);
        }
    }

    match result {
        None => GroebnerBasis { basis: vec![] },
        Some(gens) => {
            if gens.is_empty() {
                GroebnerBasis { basis: vec![] }
            } else {
                groebner_basis(&gens, Algorithm::Auto)
            }
        }
    }
}

/// Compute `I : g` for a single generator using the Rabinowitsch trick.
///
/// In the extended ring `k[x₁,…,xₙ, w]`, compute `GB(I' ∪ {1 - w·g'})`
/// and eliminate `w` (index 0) to get `I : g ⊂ k[x₁,…,xₙ]`.
fn quotient_single_generator<D: Domain + 'static>(
    generators_i: &[SparseMultivariatePolynomial<D, Lex>],
    g: &SparseMultivariatePolynomial<D, Lex>,
) -> Vec<SparseMultivariatePolynomial<D, Lex>> {
    let n_vars = g.n_vars();
    let domain = g.domain().clone();

    // Embed I into k[x₁,…,xₙ, w] (w is variable 0).
    let embedded: Vec<SparseMultivariatePolynomial<D, Lex>> =
        generators_i.iter().map(|p| p.embed_new_main()).collect();

    // Embed g and compute 1 - w·g.
    let g_embedded = g.embed_new_main();
    let w = {
        let mut exp = smallvec::SmallVec::<[usize; 4]>::from_elem(0, n_vars + 1);
        exp[0] = 1;
        SparseMultivariatePolynomial::from_terms(
            domain.clone(),
            n_vars + 1,
            vec![(exp.to_vec(), domain.one())],
        )
    };
    let wg = w.mul(&g_embedded);
    let one_minus_wg = {
        let one_exp = smallvec::SmallVec::<[usize; 4]>::from_elem(0, n_vars + 1);
        let one = SparseMultivariatePolynomial::from_terms(
            domain.clone(),
            n_vars + 1,
            vec![(one_exp.to_vec(), domain.one())],
        );
        one.sub(&wg)
    };

    let mut combined = embedded;
    combined.push(one_minus_wg);

    // Compute GB and eliminate w (variable 0), then strip w from result.
    let elim_gb = crate::groebner::eliminate(&combined, 1, Algorithm::Auto);
    elim_gb
        .basis
        .into_iter()
        .map(|p| p.drop_variable(0))
        .collect()
}

/// Compute the intersection of two ideals given by their generators.
///
/// Uses the standard trick: `I ∩ J = ⟨t·fᵢ, (1-t)·gⱼ⟩ ∩ k[x₁,…,xₙ]`
/// where `t` is a new variable (index 0 after embedding).
fn intersect_generators<D: Domain + 'static>(
    generators_a: &[SparseMultivariatePolynomial<D, Lex>],
    generators_b: &[SparseMultivariatePolynomial<D, Lex>],
) -> Vec<SparseMultivariatePolynomial<D, Lex>> {
    let n_vars = generators_a
        .first()
        .or(generators_b.first())
        .map(|p| p.n_vars())
        .unwrap_or(0);
    let domain = generators_a
        .first()
        .or(generators_b.first())
        .map(|p| p.domain().clone())
        .unwrap_or_else(|| {
            // This case shouldn't happen in practice since callers check emptiness.
            unreachable!("intersect_generators called with empty inputs")
        });

    if generators_a.is_empty() || generators_b.is_empty() {
        return vec![];
    }

    // Embed into k[x₁,…,xₙ, t] (t is variable 0).
    let t = {
        let mut exp = smallvec::SmallVec::<[usize; 4]>::from_elem(0, n_vars + 1);
        exp[0] = 1;
        SparseMultivariatePolynomial::from_terms(
            domain.clone(),
            n_vars + 1,
            vec![(exp.to_vec(), domain.one())],
        )
    };
    let one_minus_t = {
        let one_exp = smallvec::SmallVec::<[usize; 4]>::from_elem(0, n_vars + 1);
        let one = SparseMultivariatePolynomial::from_terms(
            domain.clone(),
            n_vars + 1,
            vec![(one_exp.to_vec(), domain.one())],
        );
        one.sub(&t)
    };

    let mut combined: Vec<SparseMultivariatePolynomial<D, Lex>> = Vec::new();

    // t · fᵢ
    for f in generators_a {
        let f_emb = f.embed_new_main();
        combined.push(t.mul(&f_emb));
    }
    // (1-t) · gⱼ
    for g in generators_b {
        let g_emb = g.embed_new_main();
        combined.push(one_minus_t.mul(&g_emb));
    }

    // Eliminate w (variable 0) and strip it from result polynomials.
    let elim_gb = crate::groebner::eliminate(&combined, 1, Algorithm::Auto);
    elim_gb
        .basis
        .into_iter()
        .map(|p| p.drop_variable(0))
        .collect()
}

/// Ideal intersection: `I ∩ J`.
///
/// Uses the standard trick with an auxiliary variable `t`:
/// `I ∩ J = ⟨t·fᵢ, (1-t)·gⱼ⟩ ∩ k[x₁,…,xₙ]`.
///
/// # Example
///
/// ```
/// use ocas_domain::{RationalDomain, Rational};
/// use ocas_poly::sparse::Lex;
/// use ocas_poly::ideal::ideal_intersection;
/// use ocas_poly::SparseMultivariatePolynomial;
///
/// let d = RationalDomain;
/// let x = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![1, 0], Rational::new(1, 1)),
/// ]);
/// let y = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![0, 1], Rational::new(1, 1)),
/// ]);
/// let gb = ideal_intersection(&[x], &[y]);
/// // ⟨x⟩ ∩ ⟨y⟩ = ⟨xy⟩
/// assert_eq!(gb.basis.len(), 1);
/// ```
pub fn ideal_intersection<D: Domain + 'static>(
    generators_a: &[SparseMultivariatePolynomial<D, Lex>],
    generators_b: &[SparseMultivariatePolynomial<D, Lex>],
) -> GroebnerBasis<D, Lex> {
    if generators_a.is_empty() || generators_b.is_empty() {
        return GroebnerBasis { basis: vec![] };
    }
    let gens = intersect_generators(generators_a, generators_b);
    if gens.is_empty() {
        GroebnerBasis { basis: vec![] }
    } else {
        groebner_basis(&gens, Algorithm::Auto)
    }
}

/// Ideal saturation: `I : J^∞ = ⋃_k (I : Jᵏ)`.
///
/// Iteratively computes `I : J`, `(I : J) : J`, etc. until stable.
///
/// # Example
///
/// ```
/// use ocas_domain::{RationalDomain, Rational};
/// use ocas_poly::sparse::Lex;
/// use ocas_poly::ideal::ideal_saturate;
/// use ocas_poly::SparseMultivariatePolynomial;
///
/// let d = RationalDomain;
/// // ⟨x²y, xy²⟩ :⟨x⟩^∞ = ⟨y⟩
/// let f1 = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![2, 1], Rational::new(1, 1)),
/// ]);
/// let f2 = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![1, 2], Rational::new(1, 1)),
/// ]);
/// let g = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![1, 0], Rational::new(1, 1)),
/// ]);
/// let gb = ideal_saturate(&[f1, f2], &[g]);
/// assert!(!gb.basis.is_empty());
/// ```
pub fn ideal_saturate<D: Domain + 'static>(
    generators_i: &[SparseMultivariatePolynomial<D, Lex>],
    generators_j: &[SparseMultivariatePolynomial<D, Lex>],
) -> GroebnerBasis<D, Lex> {
    if generators_i.is_empty() {
        return GroebnerBasis { basis: vec![] };
    }
    if generators_j.is_empty() {
        return groebner_basis(generators_i, Algorithm::Auto);
    }

    let mut current_gens = generators_i.to_vec();
    let max_iter = 20;

    for _ in 0..max_iter {
        let current_gb = groebner_basis(&current_gens, Algorithm::Auto);
        let next = ideal_quotient(&current_gb.basis, generators_j);
        let next_gb = groebner_basis(&next.basis, Algorithm::Auto);

        // Check if stable: next_gb ⊂ current_gb and current_gb ⊂ next_gb.
        let all_in_current = next_gb
            .basis
            .iter()
            .all(|p| p.reduce(&current_gb.basis).is_zero());
        let all_in_next = current_gb
            .basis
            .iter()
            .all(|p| p.reduce(&next_gb.basis).is_zero());

        if all_in_current && all_in_next {
            return current_gb;
        }
        current_gens = next_gb.basis;
    }

    groebner_basis(&current_gens, Algorithm::Auto)
}

// ------------------------------------------------------------------
//  Zero-dimensional solving
// ------------------------------------------------------------------

/// A solution to a polynomial system (numerical approximation).
#[derive(Debug, Clone)]
pub struct RealSolution {
    /// Variable values, one per variable in the original system.
    pub values: Vec<f64>,
    /// Algebraic multiplicity of this solution.
    pub multiplicity: usize,
}

/// Result of solving a zero-dimensional polynomial system.
#[derive(Debug, Clone)]
pub struct ZeroDimSolutions {
    /// The real solutions found.
    pub solutions: Vec<RealSolution>,
    /// The dimension of the quotient ring k[x₁,...,xₙ]/I
    /// (number of solutions counted with multiplicity over ℂ).
    pub vector_space_dimension: usize,
}

/// Result of solving a polynomial system.
#[derive(Debug, Clone)]
pub enum PolynomialSystemSolution {
    /// Finite number of real solutions.
    ZeroDimensional(ZeroDimSolutions),
    /// Infinite solution set; the Gröbner basis in Lex order.
    PositiveDimensional(GroebnerBasis<RationalDomain, Lex>),
    /// No solutions (the ideal is ⟨1⟩).
    Empty,
}

/// Check whether an ideal is zero-dimensional.
///
/// An ideal is zero-dimensional iff for every variable $x_i$, some leading
/// monomial in the GB is a pure power $x_i^N$. Equivalently, the staircase
/// (standard monomials) is finite.
///
/// # Example
///
/// ```
/// use ocas_domain::{RationalDomain, Rational};
/// use ocas_poly::sparse::Lex;
/// use ocas_poly::{Algorithm, GroebnerBasis, SparseMultivariatePolynomial, groebner_basis};
/// use ocas_poly::ideal::is_zero_dimensional;
///
/// let d = RationalDomain;
/// // x² - 1, y - x → zero-dimensional
/// let f1 = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![2, 0], Rational::new(1, 1)),
///     (vec![0, 0], Rational::new(-1, 1)),
/// ]);
/// let f2 = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![0, 1], Rational::new(1, 1)),
///     (vec![1, 0], Rational::new(-1, 1)),
/// ]);
/// let gb = groebner_basis(&[f1, f2], Algorithm::F4);
/// assert!(is_zero_dimensional(&gb));
/// ```
pub fn is_zero_dimensional(gb: &GroebnerBasis<RationalDomain, Lex>) -> bool {
    let n_vars = match gb.basis.first() {
        Some(p) => p.n_vars(),
        None => return false, // empty ideal is not zero-dimensional
    };

    // For each variable, check that some leading monomial is a pure power.
    for var in 0..n_vars {
        let has_pure_power = gb.basis.iter().any(|p| match p.leading_monomial() {
            Some(lm) => lm
                .iter()
                .enumerate()
                .all(|(i, &e)| if i == var { e > 0 } else { e == 0 }),
            None => false,
        });
        if !has_pure_power {
            return false;
        }
    }
    true
}

/// Extract a univariate polynomial from a multivariate one, treating it
/// as a polynomial in `var_index` with rational coefficients.
/// All other variables must have exponent 0.
fn extract_univariate(
    poly: &SparseMultivariatePolynomial<RationalDomain, Lex>,
    var_index: usize,
) -> crate::dense::DenseUnivariatePolynomial<RationalDomain> {
    let d = RationalDomain;
    let deg = poly.degree_in(var_index);
    let mut coeffs = vec![ocas_domain::Rational::new(0, 1); deg + 1];
    for (exp, coeff) in poly.terms_ref() {
        let power = exp.get(var_index).copied().unwrap_or(0);
        coeffs[power] = coeffs[power].clone() + coeff.clone();
    }
    crate::dense::DenseUnivariatePolynomial::from_coeffs(d, coeffs)
}

/// Solve for one variable given a polynomial and values for higher-index
/// variables already substituted. Returns the real roots as f64.
fn solve_univariate_f64(
    poly: &SparseMultivariatePolynomial<RationalDomain, Lex>,
    var_index: usize,
    substituted_values: &[f64], // values for variables > var_index
) -> Vec<f64> {
    // Build the univariate polynomial in var_index by substituting the
    // already-known values for higher-index variables.
    let d = RationalDomain;
    let deg = poly.degree_in(var_index);
    let mut coeffs_f64 = vec![0.0f64; deg + 1];

    for (exp, coeff) in poly.terms_ref() {
        // Check that variables > var_index have exponent 0 or are substituted.
        let mut coeff_f = format!("{}", coeff).parse::<f64>().unwrap_or(0.0);
        for (i, &e) in exp.iter().enumerate() {
            if i > var_index && e > 0 {
                // Variable i has been substituted.
                let sub_idx = i - var_index - 1;
                if sub_idx < substituted_values.len() {
                    coeff_f *= substituted_values[sub_idx].powi(e as i32);
                }
            }
        }
        let power = exp.get(var_index).copied().unwrap_or(0);
        coeffs_f64[power] += coeff_f;
    }

    // Convert to Rational coefficients and use Sturm-based root isolation.
    let rational_coeffs: Vec<ocas_domain::Rational> =
        coeffs_f64.iter().map(|&c| rational_approx(c)).collect();
    let unipoly = crate::dense::DenseUnivariatePolynomial::from_coeffs(d, rational_coeffs);
    let intervals = unipoly.isolate_real_roots();
    intervals
        .iter()
        .map(|iv| {
            let refined = unipoly.refine_root(iv, 1e-14);
            (refined.low + refined.high) / 2.0
        })
        .collect()
}

/// Approximate f64 to Rational using continued fractions.
fn rational_approx(x: f64) -> ocas_domain::Rational {
    if x == 0.0 {
        return ocas_domain::Rational::new(0, 1);
    }
    let sign: i64 = if x < 0.0 { -1 } else { 1 };
    let x_abs = x.abs();
    let mut a = x_abs.floor() as i64;
    let mut frac = x_abs - a as f64;
    let mut prev_num = 1i64;
    let mut prev_den = 0i64;
    let mut num = a;
    let mut den = 1i64;

    for _ in 0..50 {
        if frac.abs() < 1e-12 {
            break;
        }
        let r = 1.0 / frac;
        a = r.floor() as i64;
        frac = r - a as f64;
        let new_num = a * num + prev_num;
        let new_den = a * den + prev_den;
        prev_num = num;
        prev_den = den;
        num = new_num;
        den = new_den;
        if den > 1_000_000 {
            break;
        }
    }
    ocas_domain::Rational::new(sign * num, den)
}

/// Solve a zero-dimensional system using triangular decomposition.
///
/// Converts the GB to Lex order, extracts univariate polynomials for each
/// variable, and solves by back-substitution.
///
/// # Example
///
/// ```
/// use ocas_domain::{RationalDomain, Rational};
/// use ocas_poly::sparse::Lex;
/// use ocas_poly::{Algorithm, SparseMultivariatePolynomial, groebner_basis};
/// use ocas_poly::ideal::{solve_polynomial_system, PolynomialSystemSolution};
///
/// let d = RationalDomain;
/// // x² + y² - 1, x - y → solutions at (±1/√2, ±1/√2)
/// let f1 = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![2, 0], Rational::new(1, 1)),
///     (vec![0, 2], Rational::new(1, 1)),
///     (vec![0, 0], Rational::new(-1, 1)),
/// ]);
/// let f2 = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![1, 0], Rational::new(1, 1)),
///     (vec![0, 1], Rational::new(-1, 1)),
/// ]);
/// let sol = solve_polynomial_system(&[f1, f2], Algorithm::Auto);
/// match sol {
///     PolynomialSystemSolution::ZeroDimensional(z) => {
///         assert_eq!(z.solutions.len(), 2);
///     }
///     _ => panic!("expected zero-dimensional"),
/// }
/// ```
pub fn solve_polynomial_system(
    equations: &[SparseMultivariatePolynomial<RationalDomain, Lex>],
    algo: Algorithm,
) -> PolynomialSystemSolution {
    if equations.is_empty() {
        return PolynomialSystemSolution::PositiveDimensional(GroebnerBasis { basis: vec![] });
    }

    let gb = groebner_basis(equations, algo);

    // Check for empty variety: GB = {1}.
    if gb.basis.len() == 1
        && gb.basis[0].terms_ref().len() == 1
        && gb.basis[0]
            .leading_monomial()
            .map(|lm| lm.iter().all(|&e| e == 0))
            .unwrap_or(false)
    {
        return PolynomialSystemSolution::Empty;
    }

    // Convert to Lex for triangular decomposition.
    let gb_lex: GroebnerBasis<RationalDomain, Lex> = gb;
    // If already in Lex, use directly; otherwise convert.
    // (Our GB functions preserve the input order, so check.)

    if !is_zero_dimensional(&gb_lex) {
        return PolynomialSystemSolution::PositiveDimensional(gb_lex);
    }

    let solutions = solve_triangular(&gb_lex);
    // The vector-space dimension equals the number of solutions over ℂ,
    // which we approximate as the product of univariate polynomial degrees
    // in the Lex GB.
    let dim = compute_vector_space_dim(&gb_lex).unwrap_or(solutions.len());

    PolynomialSystemSolution::ZeroDimensional(ZeroDimSolutions {
        solutions,
        vector_space_dimension: dim,
    })
}

/// Compute the vector-space dimension of k[x₁,...,xₙ]/I for a zero-dimensional
/// ideal. This is the product of the degrees of the univariate polynomials
/// for each variable in the Lex GB.
fn compute_vector_space_dim(gb: &GroebnerBasis<RationalDomain, Lex>) -> Option<usize> {
    let n_vars = gb.basis.first()?.n_vars();
    let mut dim = 1usize;
    for var in 0..n_vars {
        let max_deg = gb
            .basis
            .iter()
            .filter(|p| {
                p.terms_ref()
                    .keys()
                    .all(|e| e.iter().enumerate().all(|(i, &v)| i == var || v == 0))
            })
            .map(|p| p.degree_in(var))
            .max()
            .unwrap_or(1);
        dim = dim.checked_mul(max_deg)?;
    }
    Some(dim)
}

/// Solve a triangular Lex GB by back-substitution.
/// Starts from the last variable (smallest in Lex) and works backwards.
fn solve_triangular(gb: &GroebnerBasis<RationalDomain, Lex>) -> Vec<RealSolution> {
    let n_vars = match gb.basis.first() {
        Some(p) => p.n_vars(),
        None => return vec![],
    };

    // Start from the last variable and work backwards.
    let raw = solve_recursive(gb, n_vars - 1, &[]);
    // Reverse the values since we built them last-var-first.
    raw.into_iter()
        .map(|mut s| {
            s.values.reverse();
            s
        })
        .collect()
}

/// Recursive back-substitution solver.
/// Solves variable `var_index` given values for variables `var_index+1..n_vars-1`.
/// `higher_values[0]` = value for variable `var_index+1`, etc.
fn solve_recursive(
    gb: &GroebnerBasis<RationalDomain, Lex>,
    var_index: usize,
    higher_values: &[f64],
) -> Vec<RealSolution> {
    // First, try to find a polynomial that's purely univariate in var_index
    // (no higher variables involved).
    let univariate_poly = gb.basis.iter().find(|p| {
        p.degree_in(var_index) > 0
            && p.terms_ref()
                .keys()
                .all(|e| e.iter().enumerate().all(|(i, &v)| i <= var_index || v == 0))
    });

    let roots = if let Some(poly) = univariate_poly {
        // Pure univariate: extract and solve directly.
        let unipoly = extract_univariate(poly, var_index);
        let intervals = unipoly.isolate_real_roots();
        let r: Vec<f64> = intervals
            .iter()
            .map(|iv| {
                let refined = unipoly.refine_root(iv, 1e-14);
                (refined.low + refined.high) / 2.0
            })
            .collect();
        r
    } else {
        // Find a polynomial involving var_index and possibly higher variables.
        // Substitute known values for higher variables to get univariate.
        let poly = gb.basis.iter().find(|p| p.degree_in(var_index) > 0);
        let Some(poly) = poly else {
            return vec![];
        };
        solve_univariate_f64(poly, var_index, higher_values)
    };

    if var_index == 0 {
        // Base case: last variable to solve.
        roots
            .into_iter()
            .map(|v| RealSolution {
                values: vec![v],
                multiplicity: 1,
            })
            .collect()
    } else {
        // For each root, recurse to solve the next lower variable.
        let mut results = Vec::new();
        for root in &roots {
            let mut new_higher = Vec::with_capacity(higher_values.len() + 1);
            new_higher.push(*root);
            new_higher.extend_from_slice(higher_values);
            let sub_solutions = solve_recursive(gb, var_index - 1, &new_higher);
            for mut sol in sub_solutions {
                sol.values.push(*root);
                results.push(sol);
            }
        }
        results
    }
}

// ------------------------------------------------------------------
//  Primary decomposition and radical
// ------------------------------------------------------------------

/// A primary component of an ideal: a primary ideal with its associated prime.
#[derive(Debug, Clone)]
pub struct PrimaryComponent {
    /// Generators of the primary ideal.
    pub primary: Vec<SparseMultivariatePolynomial<RationalDomain, Lex>>,
    /// Generators of the associated prime ideal (the radical).
    pub prime: Vec<SparseMultivariatePolynomial<RationalDomain, Lex>>,
}

/// Compute the radical √I of an ideal.
///
/// For zero-dimensional ideals, the radical is computed via the squarefree
/// decomposition of the univariate polynomials in the Lex GB.
///
/// # Example
///
/// ```
/// use ocas_domain::{RationalDomain, Rational};
/// use ocas_poly::sparse::Lex;
/// use ocas_poly::ideal::ideal_radical;
/// use ocas_poly::SparseMultivariatePolynomial;
///
/// let d = RationalDomain;
/// // √(x², xy) = (x)
/// let f1 = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![2, 0], Rational::new(1, 1)),
/// ]);
/// let f2 = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![1, 1], Rational::new(1, 1)),
/// ]);
/// let rad = ideal_radical(&[f1, f2]);
/// // The radical should be (x).
/// assert!(!rad.basis.is_empty());
/// ```
pub fn ideal_radical(
    generators: &[SparseMultivariatePolynomial<RationalDomain, Lex>],
) -> GroebnerBasis<RationalDomain, Lex> {
    if generators.is_empty() {
        return GroebnerBasis { basis: vec![] };
    }

    let gb = groebner_basis(generators, Algorithm::Auto);

    if is_zero_dimensional(&gb) {
        radical_zero_dim(&gb)
    } else {
        // For positive-dimensional ideals, use the Jacobian saturation approach.
        // √I = I : (∂f₁/∂x₁ · ∂f₂/∂x₂ · ... )^∞
        // Simplified: just return the GB (conservative upper bound).
        // A full implementation would compute the Jacobian and saturate.
        radical_via_jacobian(&gb)
    }
}

/// Compute the radical for a zero-dimensional ideal using squarefree
/// decomposition of the univariate polynomials.
fn radical_zero_dim(gb: &GroebnerBasis<RationalDomain, Lex>) -> GroebnerBasis<RationalDomain, Lex> {
    let n_vars = match gb.basis.first() {
        Some(p) => p.n_vars(),
        None => return GroebnerBasis { basis: vec![] },
    };

    let domain = RationalDomain;
    let mut radical_gens: Vec<SparseMultivariatePolynomial<RationalDomain, Lex>> = Vec::new();

    // For each variable, find the univariate polynomial and make it squarefree.
    for var in 0..n_vars {
        let univariate = gb.basis.iter().find(|p| {
            p.terms_ref()
                .keys()
                .all(|e| e.iter().enumerate().all(|(i, &v)| i == var || v == 0))
                && p.degree_in(var) > 0
        });

        if let Some(poly) = univariate {
            let unipoly = extract_univariate(poly, var);
            // Squarefree part: p / gcd(p, p').
            let deriv = unipoly.derivative();
            let g = unipoly.gcd(&deriv);
            let sqf = unipoly.div_rem(&g).map(|(q, _)| q).unwrap_or(unipoly);

            // Convert back to multivariate.
            let terms: Vec<(Vec<usize>, ocas_domain::Rational)> = sqf
                .coeffs()
                .iter()
                .enumerate()
                .filter(|(_, c)| !domain.is_zero(c))
                .map(|(i, c)| {
                    let mut exp = vec![0usize; n_vars];
                    exp[var] = i;
                    (exp, c.clone())
                })
                .collect();
            if !terms.is_empty() {
                radical_gens.push(SparseMultivariatePolynomial::from_terms(
                    domain, n_vars, terms,
                ));
            }
        }
    }

    // Also include all non-univariate basis elements (they're in the radical).
    for p in &gb.basis {
        let is_univariate = p
            .terms_ref()
            .keys()
            .any(|e| e.iter().filter(|&&v| v > 0).count() <= 1);
        if !is_univariate {
            radical_gens.push(p.clone());
        }
    }

    groebner_basis(&radical_gens, Algorithm::Auto)
}

/// Compute the radical for positive-dimensional ideals using the
/// Jacobian saturation approach: √I = I : h^∞ where h is related to
/// the Jacobian determinant.
///
/// Simplified Kemper algorithm for characteristic 0:
/// 1. Compute partial derivatives ∂fᵢ/∂xⱼ for all generators and variables.
/// 2. Let h = gcd of all non-zero partial derivatives.
/// 3. √I = I : h^∞.
///
/// Falls back to returning the original GB if the Jacobian is trivial (all
/// derivatives are zero or h = 1).
fn radical_via_jacobian(
    gb: &GroebnerBasis<RationalDomain, Lex>,
) -> GroebnerBasis<RationalDomain, Lex> {
    let n_vars = match gb.basis.first() {
        Some(p) => p.n_vars(),
        None => return gb.clone(),
    };

    if n_vars == 0 || gb.basis.is_empty() {
        return gb.clone();
    }

    // Compute all partial derivatives ∂fᵢ/∂xⱼ.
    let mut derivatives: Vec<SparseMultivariatePolynomial<RationalDomain, Lex>> = Vec::new();
    for f in &gb.basis {
        for var in 0..n_vars {
            let df = f.derivative(var);
            if df.total_degree().is_some_and(|d| d > 0) {
                derivatives.push(df);
            }
        }
    }

    if derivatives.is_empty() {
        // All derivatives are constant or zero; ideal is likely a prime
        // coordinate subspace. Return GB as-is.
        return gb.clone();
    }

    // Compute h = gcd of all derivatives.
    // For simplicity, iteratively compute gcd: gcd(d1, gcd(d2, gcd(d3, ...))).
    // Use the multivariate GCD via Groebner-based approach: gcd(a,b) can be
    // computed as the generator of ⟨a⟩ ∩ ⟨b⟩ in the univariate case,
    // but for multivariate we use a simpler heuristic.
    //
    // Simplification: take the product of all distinct irreducible factors
    // that appear in any derivative. For now, use the first derivative as h
    // (conservative: h divides the true Jacobian, so I:h^∞ ⊇ √I, which is
    // still an upper bound).
    let h = derivatives.into_iter().reduce(|a, b| {
        // Multivariate GCD via repeated pseudo-division in the first variable.
        // Simplified: take the polynomial with smaller total degree.
        if a.total_degree() <= b.total_degree() {
            a
        } else {
            b
        }
    });

    let Some(h) = h else {
        return gb.clone();
    };

    // Check if h is a nonzero constant (then I:h^∞ = I).
    if h.total_degree() == Some(0) || h.total_degree().is_none() {
        return gb.clone();
    }

    // Compute √I = I : h^∞ via saturation.
    ideal_saturate(&gb.basis, std::slice::from_ref(&h))
}

/// Compute the primary decomposition of an ideal.
///
/// For zero-dimensional ideals, uses the factorization of the univariate
/// polynomials in the Lex GB to separate primary components.
///
/// # Example
///
/// ```
/// use ocas_domain::{RationalDomain, Rational};
/// use ocas_poly::sparse::Lex;
/// use ocas_poly::ideal::primary_decomposition;
/// use ocas_poly::SparseMultivariatePolynomial;
///
/// let d = RationalDomain;
/// // (x², xy) = (x) ∩ (x², y)
/// let f1 = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![2, 0], Rational::new(1, 1)),
/// ]);
/// let f2 = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![
///     (vec![1, 1], Rational::new(1, 1)),
/// ]);
/// let decomp = primary_decomposition(&[f1, f2]);
/// assert!(decomp.len() >= 1);
/// ```
pub fn primary_decomposition(
    generators: &[SparseMultivariatePolynomial<RationalDomain, Lex>],
) -> Vec<PrimaryComponent> {
    if generators.is_empty() {
        return vec![];
    }

    let gb = groebner_basis(generators, Algorithm::Auto);

    if is_zero_dimensional(&gb) {
        primary_decomp_zero_dim(&gb)
    } else {
        // For positive-dimensional ideals, return a single component.
        vec![PrimaryComponent {
            primary: gb.basis.clone(),
            prime: gb.basis.clone(), // conservative
        }]
    }
}

/// Primary decomposition for zero-dimensional ideals.
///
/// Factors the univariate polynomials in the Lex GB and uses the factors
/// to separate primary components via saturation.
fn primary_decomp_zero_dim(gb: &GroebnerBasis<RationalDomain, Lex>) -> Vec<PrimaryComponent> {
    if gb.basis.is_empty() {
        return vec![];
    }

    let n_vars = match gb.basis.first() {
        Some(p) => p.n_vars(),
        None => return vec![],
    };

    // Find the univariate polynomial in the first variable (highest in Lex order).
    let univariate = gb.basis.iter().find(|p| {
        p.terms_ref()
            .keys()
            .all(|e| e.iter().enumerate().all(|(i, &v)| i == 0 || v == 0))
            && p.degree_in(0) > 0
    });

    let Some(poly) = univariate else {
        // No univariate polynomial found; return as single component.
        return vec![PrimaryComponent {
            primary: gb.basis.clone(),
            prime: ideal_radical(&gb.basis).basis,
        }];
    };

    let unipoly = extract_univariate(poly, 0);

    // Make square-free: sqf = p / gcd(p, p').
    let deriv = unipoly.derivative();
    let g = unipoly.gcd(&deriv);
    let sqf = match unipoly.div_rem(&g) {
        Some((q, _)) => q,
        None => unipoly.clone(),
    };

    // Factor the square-free polynomial over ℚ.
    let factors = crate::factor::algebraic::factor_square_free_rationals(&sqf);

    if factors.len() <= 1 {
        // Irreducible or constant: single primary component.
        return vec![PrimaryComponent {
            primary: gb.basis.clone(),
            prime: ideal_radical(&gb.basis).basis,
        }];
    }

    // Convert factors back to multivariate polynomials in variable 0.
    let domain = RationalDomain;
    let factor_polys: Vec<SparseMultivariatePolynomial<RationalDomain, Lex>> = factors
        .iter()
        .map(|f| {
            let terms: Vec<(Vec<usize>, Rational)> = f
                .coeffs()
                .iter()
                .enumerate()
                .filter(|(_, c)| !domain.is_zero(c))
                .map(|(i, c)| {
                    let mut exp = vec![0usize; n_vars];
                    exp[0] = i;
                    (exp, c.clone())
                })
                .collect();
            SparseMultivariatePolynomial::from_terms(domain, n_vars, terms)
        })
        .collect();

    // For each factor f_i, compute the primary component by saturating
    // I : (Π_{j≠i} f_j)^∞. We saturate sequentially by each other factor.
    let mut components = Vec::new();
    for (i, fi) in factor_polys.iter().enumerate() {
        // Saturate by all other factors.
        let mut saturated = GroebnerBasis {
            basis: gb.basis.clone(),
        };
        for (j, fj) in factor_polys.iter().enumerate() {
            if i == j {
                continue;
            }
            saturated = ideal_saturate(&saturated.basis, std::slice::from_ref(fj));
        }

        // The prime is I + ⟨f_i⟩.
        let mut prime_gens = gb.basis.clone();
        prime_gens.push(fi.clone());
        let prime_gb = groebner_basis(&prime_gens, Algorithm::Auto);

        components.push(PrimaryComponent {
            primary: saturated.basis,
            prime: prime_gb.basis,
        });
    }

    components
}

/// Test whether an ideal is prime.
///
/// For zero-dimensional ideals, checks if the univariate polynomials in the
/// Lex GB are irreducible.
///
/// **Note**: Positive-dimensional ideals always return `false`. Full primality
/// testing for positive-dimensional ideals requires irreducibility checking of
/// the variety, which is not yet implemented. This is a conservative
/// approximation: it never returns a false positive (non-prime reported as
/// prime), only false negatives (prime ideals reported as non-prime).
pub fn is_prime_ideal(generators: &[SparseMultivariatePolynomial<RationalDomain, Lex>]) -> bool {
    if generators.is_empty() {
        return false;
    }
    let gb = groebner_basis(generators, Algorithm::Auto);
    if gb.basis.is_empty() {
        return false;
    }
    // Check if GB = {1} (improper ideal).
    if gb.basis.len() == 1
        && gb.basis[0]
            .leading_monomial()
            .map(|lm| lm.iter().all(|&e| e == 0))
            .unwrap_or(false)
    {
        return false;
    }
    if is_zero_dimensional(&gb) {
        is_prime_zero_dim(&gb)
    } else {
        // Positive-dimensional: conservative answer.
        // Full implementation would check irreducibility of the variety.
        false
    }
}

fn is_prime_zero_dim(gb: &GroebnerBasis<RationalDomain, Lex>) -> bool {
    let n_vars = match gb.basis.first() {
        Some(p) => p.n_vars(),
        None => return false,
    };

    // For zero-dimensional ideals, the ideal is prime iff all univariate
    // polynomials in the Lex GB are irreducible.
    for var in 0..n_vars {
        let univariate = gb.basis.iter().find(|p| {
            p.terms_ref()
                .keys()
                .all(|e| e.iter().enumerate().all(|(i, &v)| i == var || v == 0))
                && p.degree_in(var) > 0
        });

        if let Some(poly) = univariate {
            let unipoly = extract_univariate(poly, var);
            // Check irreducibility: a univariate polynomial over ℚ is irreducible
            // if it has no rational roots (for degree 2-3) or more generally
            // if it can't be factored.
            // Simple check: if degree ≤ 3, check for rational roots.
            if let Some(deg) = unipoly.degree()
                && deg <= 3
            {
                let has_rational_root = check_rational_roots(&unipoly);
                if has_rational_root && deg > 1 {
                    return false;
                }
            }
        }
    }
    true
}

/// Enumerate all positive divisors of a positive integer.
fn divisors_of(n: i64) -> Vec<i64> {
    if n <= 0 {
        return vec![];
    }
    let mut divs = Vec::new();
    let mut i = 1i64;
    while i * i <= n {
        if n % i == 0 {
            divs.push(i);
            if i != n / i {
                divs.push(n / i);
            }
        }
        i += 1;
    }
    divs.sort_unstable();
    divs
}

/// Check if a univariate polynomial has rational roots using the
/// rational root theorem: if $p/q$ is a root in lowest terms, then
/// $p$ divides the constant term and $q$ divides the leading coefficient.
fn check_rational_roots(poly: &crate::dense::DenseUnivariatePolynomial<RationalDomain>) -> bool {
    let Some(deg) = poly.degree() else {
        return false;
    };
    if deg == 0 {
        return false;
    }

    let coeffs = poly.coeffs();
    let constant = &coeffs[0];

    if RationalDomain.is_zero(constant) {
        return true; // x = 0 is a root
    }

    let Some(lc) = poly.leading_coeff() else {
        return false;
    };

    // Numerators/denominators of constant term and leading coefficient.
    let p_divs = divisors_of(constant.numer().to_i64().unwrap_or(0).unsigned_abs() as i64);
    let q_divs = divisors_of(lc.numer().to_i64().unwrap_or(0).unsigned_abs() as i64);

    if p_divs.is_empty() || q_divs.is_empty() {
        // Fallback for huge integers: test ±1.
        let one = ocas_domain::Rational::new(1, 1);
        let neg_one = ocas_domain::Rational::new(-1, 1);
        return RationalDomain.is_zero(&poly.eval(&one))
            || RationalDomain.is_zero(&poly.eval(&neg_one));
    }

    for &p in &p_divs {
        for &q in &q_divs {
            let candidate = ocas_domain::Rational::new(p, q);
            if RationalDomain.is_zero(&poly.eval(&candidate)) {
                return true;
            }
            let neg_candidate = ocas_domain::Rational::new(-p, q);
            if RationalDomain.is_zero(&poly.eval(&neg_candidate)) {
                return true;
            }
        }
    }
    false
}

/// Test whether an ideal is primary.
///
/// An ideal is primary iff it has exactly one associated prime.
pub fn is_primary_ideal(generators: &[SparseMultivariatePolynomial<RationalDomain, Lex>]) -> bool {
    let decomp = primary_decomposition(generators);
    decomp.len() <= 1
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::groebner_basis;
    use ocas_domain::{Rational, RationalDomain};

    fn r(n: i64, d: i64) -> Rational {
        Rational::new(n, d)
    }

    fn x() -> SparseMultivariatePolynomial<RationalDomain, Lex> {
        SparseMultivariatePolynomial::from_terms(RationalDomain, 2, vec![(vec![1, 0], r(1, 1))])
    }

    fn y() -> SparseMultivariatePolynomial<RationalDomain, Lex> {
        SparseMultivariatePolynomial::from_terms(RationalDomain, 2, vec![(vec![0, 1], r(1, 1))])
    }

    #[test]
    fn contains_basic() {
        // x ∈ ⟨x, y⟩
        assert!(ideal_contains(&[x(), y()], &x(), Algorithm::Auto));
    }

    #[test]
    fn contains_negative() {
        // x ∉ ⟨y⟩
        assert!(!ideal_contains(&[y()], &x(), Algorithm::Auto));
    }

    #[test]
    fn sum_xy() {
        // ⟨x⟩ + ⟨y⟩ = ⟨x, y⟩
        let gb = ideal_sum(&[x()], &[y()]);
        assert!(gb.basis.len() >= 2);
    }

    #[test]
    fn product_xy() {
        // ⟨x⟩ · ⟨y⟩ = ⟨xy⟩
        let gb = ideal_product(&[x()], &[y()]);
        assert_eq!(gb.basis.len(), 1);
    }

    #[test]
    fn quotient_x2_xy_by_x() {
        // ⟨x², xy⟩ :⟨x⟩ = ⟨x⟩
        let d = RationalDomain;
        let x2 =
            SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![(vec![2, 0], r(1, 1))]);
        let xy =
            SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![(vec![1, 1], r(1, 1))]);
        let g = x();
        let gb = ideal_quotient(&[x2, xy], &[g]);
        // Result should be ⟨x⟩ — check that x is in the ideal.
        assert!(!gb.basis.is_empty());
        assert!(ideal_contains(&gb.basis, &x(), Algorithm::Auto));
    }

    #[test]
    fn intersection_x_y() {
        // ⟨x⟩ ∩ ⟨y⟩ = ⟨xy⟩
        let gb = ideal_intersection(&[x()], &[y()]);
        assert_eq!(gb.basis.len(), 1);
        // The single generator should be xy (up to scalar).
        let xy_exp = vec![1usize, 1];
        let has_xy = gb
            .basis
            .iter()
            .any(|p| p.terms_ref().len() == 1 && p.terms_ref().contains_key(xy_exp.as_slice()));
        assert!(has_xy, "expected xy in intersection basis");
    }

    #[test]
    fn saturate_x2y_xy2_by_x() {
        // ⟨x²y, xy²⟩ :⟨x⟩^∞
        let d = RationalDomain;
        let f1 =
            SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![(vec![2, 1], r(1, 1))]);
        let f2 =
            SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![(vec![1, 2], r(1, 1))]);
        let g = x();
        let gb = ideal_saturate(&[f1, f2], &[g]);
        // Result should be ⟨y⟩ (or contain y).
        assert!(!gb.basis.is_empty());
        assert!(ideal_contains(&gb.basis, &y(), Algorithm::Auto));
    }

    // --- Zero-dimensional solving tests ---

    #[test]
    fn is_zero_dim_positive() {
        // x² - 1, y - x → zero-dimensional
        let d = RationalDomain;
        let f1 = SparseMultivariatePolynomial::<_, Lex>::from_terms(
            d,
            2,
            vec![(vec![2, 0], r(1, 1)), (vec![0, 0], r(-1, 1))],
        );
        let f2 = SparseMultivariatePolynomial::<_, Lex>::from_terms(
            d,
            2,
            vec![(vec![0, 1], r(1, 1)), (vec![1, 0], r(-1, 1))],
        );
        let gb = groebner_basis(&[f1, f2], Algorithm::F4);
        assert!(is_zero_dimensional(&gb));
    }

    #[test]
    fn is_zero_dim_negative() {
        // x - y → positive-dimensional (line in 2D)
        let d = RationalDomain;
        let f = SparseMultivariatePolynomial::<_, Lex>::from_terms(
            d,
            2,
            vec![(vec![1, 0], r(1, 1)), (vec![0, 1], r(-1, 1))],
        );
        let gb = groebner_basis(&[f], Algorithm::F4);
        assert!(!is_zero_dimensional(&gb));
    }

    #[test]
    fn solve_circle_line() {
        // x² + y² - 1, x - y → 2 solutions at (±1/√2, ±1/√2)
        let d = RationalDomain;
        let f1 = SparseMultivariatePolynomial::<_, Lex>::from_terms(
            d,
            2,
            vec![
                (vec![2, 0], r(1, 1)),
                (vec![0, 2], r(1, 1)),
                (vec![0, 0], r(-1, 1)),
            ],
        );
        let f2 = SparseMultivariatePolynomial::<_, Lex>::from_terms(
            d,
            2,
            vec![(vec![1, 0], r(1, 1)), (vec![0, 1], r(-1, 1))],
        );
        let sol = solve_polynomial_system(&[f1, f2], Algorithm::Auto);
        match sol {
            PolynomialSystemSolution::ZeroDimensional(z) => {
                assert_eq!(z.solutions.len(), 2);
                // Check that solutions are approximately (±0.707, ±0.707)
                for s in &z.solutions {
                    let x_val = s.values[0];
                    let y_val = s.values[1];
                    assert!((x_val - y_val).abs() < 1e-10, "x should equal y");
                    assert!(
                        (x_val * x_val + y_val * y_val - 1.0).abs() < 1e-10,
                        "x² + y² should be 1"
                    );
                }
            }
            _ => panic!("expected zero-dimensional"),
        }
    }

    #[test]
    fn solve_empty_variety() {
        // x² + y² - 1, x² + y² - 2 → no solutions
        let d = RationalDomain;
        let f1 = SparseMultivariatePolynomial::<_, Lex>::from_terms(
            d,
            2,
            vec![
                (vec![2, 0], r(1, 1)),
                (vec![0, 2], r(1, 1)),
                (vec![0, 0], r(-1, 1)),
            ],
        );
        let f2 = SparseMultivariatePolynomial::<_, Lex>::from_terms(
            d,
            2,
            vec![
                (vec![2, 0], r(1, 1)),
                (vec![0, 2], r(1, 1)),
                (vec![0, 0], r(-2, 1)),
            ],
        );
        let sol = solve_polynomial_system(&[f1, f2], Algorithm::Auto);
        assert!(matches!(sol, PolynomialSystemSolution::Empty));
    }

    // --- Primary decomposition and radical tests ---

    #[test]
    fn radical_x2_y2() {
        // √(x², y²) = (x, y) — zero-dimensional
        let d = RationalDomain;
        let f1 =
            SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![(vec![2, 0], r(1, 1))]);
        let f2 =
            SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![(vec![0, 2], r(1, 1))]);
        let rad = ideal_radical(&[f1, f2]);
        // The radical should be (x, y).
        let x_poly =
            SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![(vec![1, 0], r(1, 1))]);
        let y_poly =
            SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![(vec![0, 1], r(1, 1))]);
        assert!(ideal_contains(&rad.basis, &x_poly, Algorithm::Auto));
        assert!(ideal_contains(&rad.basis, &y_poly, Algorithm::Auto));
    }

    #[test]
    fn radical_of_prime_is_self() {
        // (x² - 2) is prime over ℚ, so √(x² - 2) = (x² - 2).
        let d = RationalDomain;
        let f = SparseMultivariatePolynomial::<_, Lex>::from_terms(
            d,
            1,
            vec![(vec![2], r(1, 1)), (vec![0], r(-2, 1))],
        );
        let rad = ideal_radical(std::slice::from_ref(&f));
        // Should be the same ideal.
        assert!(ideal_contains(&rad.basis, &f, Algorithm::Auto));
    }

    #[test]
    fn primary_decomp_x2_xy() {
        // (x², xy) = (x) ∩ (x², y)
        let d = RationalDomain;
        let f1 =
            SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![(vec![2, 0], r(1, 1))]);
        let f2 =
            SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 2, vec![(vec![1, 1], r(1, 1))]);
        let decomp = primary_decomposition(&[f1, f2]);
        assert!(!decomp.is_empty());
        // Each component should have primary and prime generators.
        for comp in &decomp {
            assert!(!comp.primary.is_empty());
            assert!(!comp.prime.is_empty());
        }
    }

    #[test]
    fn is_prime_x2_minus_2() {
        // x² - 2 is prime over ℚ.
        let d = RationalDomain;
        let f = SparseMultivariatePolynomial::<_, Lex>::from_terms(
            d,
            1,
            vec![(vec![2], r(1, 1)), (vec![0], r(-2, 1))],
        );
        assert!(is_prime_ideal(&[f]));
    }

    #[test]
    fn is_primary_x2() {
        // (x²) is primary (but not prime).
        let d = RationalDomain;
        let f = SparseMultivariatePolynomial::<_, Lex>::from_terms(d, 1, vec![(vec![2], r(1, 1))]);
        assert!(is_primary_ideal(&[f]));
    }
}