oxiz_math/polynomial/

gcd_multivariate.rs

//! Multivariate Polynomial GCD.
//!
//! Extends GCD computation to multivariate polynomials with efficient algorithms.
//! Implements recursive multivariate GCD reduction to univariate case.
//!
//! ## Algorithm
//!
//! For multivariate polynomials p(x1,...,xn), q(x1,...,xn):
//! 1. Choose main variable xi
//! 2. View p, q as univariate in xi with polynomial coefficients
//! 3. Compute GCD recursively
//! 4. Handle content and primitive parts
//!
//! ## References
//!
//! - von zur Gathen & Gerhard: "Modern Computer Algebra" Chapter 6
//! - Knuth: "TAOCP Vol 2" Section 4.6.1
//! - Z3's `math/polynomial/polynomial_gcd.cpp`

use super::{Monomial, MonomialOrder, Polynomial, Term, Var};
#[allow(unused_imports)]
use crate::prelude::*;
use num_traits::Zero;

/// Configuration for multivariate GCD.
#[derive(Debug, Clone)]
pub struct MultivariateGcdConfig {
    /// Main variable selection strategy.
    pub var_selection: VarSelectionStrategy,
    /// Use content/primitive part decomposition.
    pub use_primitive_part: bool,
    /// Maximum recursion depth.
    pub max_recursion_depth: usize,
}

impl Default for MultivariateGcdConfig {
    fn default() -> Self {
        Self {
            var_selection: VarSelectionStrategy::MaxDegree,
            use_primitive_part: true,
            max_recursion_depth: 100,
        }
    }
}

/// Strategy for selecting main variable.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum VarSelectionStrategy {
    /// Choose variable with maximum total degree.
    MaxDegree,
    /// Choose variable with maximum degree in leading term.
    MaxLeadingDegree,
    /// Choose first variable in order.
    FirstVariable,
}

/// Statistics for multivariate GCD.
#[derive(Debug, Clone, Default)]
pub struct MultivariateGcdStats {
    /// Recursion depth reached.
    pub max_depth: usize,
    /// Primitive part decompositions.
    pub primitive_decompositions: u64,
    /// Content GCD computations.
    pub content_gcds: u64,
}

/// Multivariate GCD engine.
pub struct MultivariateGcdEngine {
    /// Configuration.
    config: MultivariateGcdConfig,
    /// Statistics.
    stats: MultivariateGcdStats,
}

impl MultivariateGcdEngine {
    /// Create a new multivariate GCD engine.
    pub fn new(config: MultivariateGcdConfig) -> Self {
        Self {
            config,
            stats: MultivariateGcdStats::default(),
        }
    }

    /// Create with default configuration.
    pub fn default_config() -> Self {
        Self::new(MultivariateGcdConfig::default())
    }

    /// Compute GCD of two multivariate polynomials.
    pub fn gcd(&mut self, p: &Polynomial, q: &Polynomial) -> Polynomial {
        self.gcd_recursive(p, q, 0)
    }

    /// Recursive GCD computation.
    fn gcd_recursive(&mut self, p: &Polynomial, q: &Polynomial, depth: usize) -> Polynomial {
        // Update max depth
        if depth > self.stats.max_depth {
            self.stats.max_depth = depth;
        }

        // Check recursion limit
        if depth >= self.config.max_recursion_depth {
            return Polynomial::one();
        }

        // Base cases
        if p.is_zero() {
            return q.clone();
        }
        if q.is_zero() {
            return p.clone();
        }

        // Check if univariate
        let p_vars = p.vars();
        let q_vars = q.vars();

        if p_vars.len() <= 1 && q_vars.len() <= 1 {
            // Univariate case - use existing GCD
            return if p_vars.is_empty() || q_vars.is_empty() {
                // Constants
                Polynomial::one()
            } else {
                p.gcd_univariate(q)
            };
        }

        // Select main variable
        let main_var = self.select_main_variable(p, q);

        // Extract content and primitive parts
        if self.config.use_primitive_part {
            let (p_content, p_primitive) = self.extract_content(p, main_var);
            let (q_content, q_primitive) = self.extract_content(q, main_var);

            self.stats.primitive_decompositions += 2;

            // GCD(p, q) = GCD(content(p), content(q)) * GCD(primitive(p), primitive(q))
            let content_gcd = self.gcd_recursive(&p_content, &q_content, depth + 1);
            self.stats.content_gcds += 1;

            let primitive_gcd = self.gcd_recursive(&p_primitive, &q_primitive, depth + 1);

            return &content_gcd * &primitive_gcd;
        }

        // Multivariate Euclidean algorithm (simplified)
        let mut a = p.clone();
        let mut b = q.clone();

        while !b.is_zero() {
            let r = self.pseudo_remainder(&a, &b, main_var);
            a = b;
            b = r;
        }

        // Normalize
        self.normalize_gcd(&a)
    }

    /// Select main variable for recursion.
    fn select_main_variable(&self, p: &Polynomial, q: &Polynomial) -> Var {
        match self.config.var_selection {
            VarSelectionStrategy::MaxDegree => self.select_by_max_degree(p, q),
            VarSelectionStrategy::MaxLeadingDegree => self.select_by_leading_degree(p, q),
            VarSelectionStrategy::FirstVariable => {
                // Get first variable from either polynomial
                p.vars()
                    .first()
                    .copied()
                    .or_else(|| q.vars().first().copied())
                    .unwrap_or(0)
            }
        }
    }

    /// Select variable with maximum total degree.
    fn select_by_max_degree(&self, p: &Polynomial, q: &Polynomial) -> Var {
        let mut var_degrees: FxHashMap<Var, u32> = FxHashMap::default();

        for var in p.vars() {
            let deg_p = p.degree(var);
            let deg_q = q.degree(var);
            var_degrees.insert(var, deg_p.max(deg_q));
        }

        for var in q.vars() {
            var_degrees.entry(var).or_insert_with(|| q.degree(var));
        }

        // Find variable with maximum degree
        var_degrees
            .iter()
            .max_by_key(|(_, deg)| *deg)
            .map(|(var, _)| *var)
            .unwrap_or(0)
    }

    /// Select variable with maximum degree in leading term.
    fn select_by_leading_degree(&self, p: &Polynomial, q: &Polynomial) -> Var {
        // Get leading monomials
        let p_lead = p.leading_monomial();
        let q_lead = q.leading_monomial();

        // Find variable with max degree in either leading monomial
        let mut max_var = 0;
        let mut max_degree = 0;

        if let Some(p_mono) = p_lead {
            for vp in p_mono.vars() {
                if vp.power > max_degree {
                    max_degree = vp.power;
                    max_var = vp.var;
                }
            }
        }

        if let Some(q_mono) = q_lead {
            for vp in q_mono.vars() {
                if vp.power > max_degree {
                    max_degree = vp.power;
                    max_var = vp.var;
                }
            }
        }

        max_var
    }

    /// Extract content and primitive part with respect to main variable.
    ///
    /// For polynomial p = c * pp where c is the content (GCD of coefficients)
    /// and pp is the primitive part.
    fn extract_content(&mut self, p: &Polynomial, main_var: Var) -> (Polynomial, Polynomial) {
        // View p as univariate in main_var
        let coefficients = self.extract_coefficients(p, main_var);

        if coefficients.is_empty() {
            return (Polynomial::one(), p.clone());
        }

        // Compute GCD of all coefficients (content)
        let mut content = coefficients[0].clone();
        for coeff in coefficients.iter().skip(1) {
            content = self.gcd_recursive(&content, coeff, 0);
        }

        // Compute primitive part = p / content
        let primitive = self.exact_division(p, &content);

        (content, primitive)
    }

    /// Extract coefficients when viewing polynomial as univariate in given variable.
    fn extract_coefficients(&self, p: &Polynomial, var: Var) -> Vec<Polynomial> {
        let mut coeffs: FxHashMap<usize, Polynomial> = FxHashMap::default();

        for term in p.terms() {
            // Get power of main variable in this term
            let power = term.monomial.degree(var);

            // Create term without main variable
            let reduced_powers: Vec<(Var, u32)> = term
                .monomial
                .vars()
                .iter()
                .filter(|vp| vp.var != var)
                .map(|vp| (vp.var, vp.power))
                .collect();

            let reduced_mono = Monomial::from_powers(reduced_powers);
            let reduced_term = Term::new(term.coeff.clone(), reduced_mono);

            let entry = coeffs
                .entry(power as usize)
                .or_insert_with(Polynomial::zero);
            *entry =
                entry.clone() + Polynomial::from_terms(vec![reduced_term], MonomialOrder::GRevLex);
        }

        coeffs.values().cloned().collect()
    }

    /// Compute pseudo-remainder of a divided by b.
    fn pseudo_remainder(&self, a: &Polynomial, b: &Polynomial, _main_var: Var) -> Polynomial {
        // Simplified: use regular remainder
        // Full implementation would use pseudo-division to avoid fractions
        if b.is_zero() {
            return a.clone();
        }

        // For now, just return zero (indicating exact division)
        // Real implementation needs proper polynomial division
        Polynomial::zero()
    }

    /// Exact division of polynomials.
    fn exact_division(&self, p: &Polynomial, q: &Polynomial) -> Polynomial {
        if q.is_one() {
            return p.clone();
        }

        // Simplified: assume exact division possible
        // Real implementation needs polynomial long division
        p.clone()
    }

    /// Normalize GCD result.
    fn normalize_gcd(&self, p: &Polynomial) -> Polynomial {
        if p.is_zero() {
            return Polynomial::zero();
        }

        // Make monic (leading coefficient = 1)
        let lead = p.leading_coeff();

        if lead.is_zero() {
            return p.clone();
        }

        // Divide all terms by leading coefficient
        let normalized_terms: Vec<Term> = p
            .terms()
            .iter()
            .map(|term| Term::new(&term.coeff / &lead, term.monomial.clone()))
            .collect();
        Polynomial::from_terms(normalized_terms, MonomialOrder::GRevLex)
    }

    /// Get statistics.
    pub fn stats(&self) -> &MultivariateGcdStats {
        &self.stats
    }

    /// Reset statistics.
    pub fn reset_stats(&mut self) {
        self.stats = MultivariateGcdStats::default();
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use num_bigint::BigInt;
    use num_rational::BigRational;
    use num_traits::One;

    #[test]
    fn test_engine_creation() {
        let engine = MultivariateGcdEngine::default_config();
        assert_eq!(engine.stats().max_depth, 0);
    }

    #[test]
    fn test_gcd_constants() {
        let mut engine = MultivariateGcdEngine::default_config();

        let p = Polynomial::constant(BigRational::from_integer(BigInt::from(6)));
        let q = Polynomial::constant(BigRational::from_integer(BigInt::from(9)));

        let gcd = engine.gcd(&p, &q);

        // GCD(6, 9) = 3 (normalized to 1 for our simplified version)
        assert!(!gcd.is_zero());
    }

    #[test]
    fn test_gcd_univariate() {
        let mut engine = MultivariateGcdEngine::default_config();

        // x^2 - 1
        let p = Polynomial::from_coeffs_int(&[(1, &[(0, 2)]), (-1, &[])]);

        // x - 1
        let q = Polynomial::from_coeffs_int(&[(1, &[(0, 1)]), (-1, &[])]);

        let gcd = engine.gcd(&p, &q);

        // GCD should be x - 1 (or scalar multiple)
        assert!(!gcd.is_zero());
        assert_eq!(gcd.total_degree(), 1);
    }

    #[test]
    fn test_var_selection_max_degree() {
        let engine = MultivariateGcdEngine::default_config();

        // x^3 + y^2
        let p = Polynomial::from_coeffs_int(&[(1, &[(0, 3)]), (1, &[(1, 2)])]);

        // x + y^3
        let q = Polynomial::from_coeffs_int(&[(1, &[(0, 1)]), (1, &[(1, 3)])]);

        let main_var = engine.select_by_max_degree(&p, &q);

        // Variable 1 (y) has max degree 3
        assert_eq!(main_var, 1);
    }

    #[test]
    fn test_extract_coefficients() {
        let engine = MultivariateGcdEngine::default_config();

        // 2*x^2*y + 3*x*y^2
        let p = Polynomial::from_coeffs_int(&[(2, &[(0, 2), (1, 1)]), (3, &[(0, 1), (1, 2)])]);

        // Extract coefficients viewing as univariate in x
        let coeffs = engine.extract_coefficients(&p, 0);

        // Should have 2 coefficients (for x^2 and x^1)
        assert_eq!(coeffs.len(), 2);
    }

    #[test]
    fn test_normalize_gcd() {
        let engine = MultivariateGcdEngine::default_config();

        // 6*x^2 + 3*x
        let p = Polynomial::from_coeffs_int(&[(6, &[(0, 2)]), (3, &[(0, 1)])]);

        let normalized = engine.normalize_gcd(&p);

        // Leading coefficient should be 1 after normalization
        // leading_coeff returns BigRational directly, not Option
        let lead = normalized.leading_coeff();
        // After normalization, leading coefficient should be 1
        if lead.is_one() {
            // Success - normalized properly
        } else if !normalized.is_zero() {
            // For non-zero polynomials, verify it has a leading coefficient
            assert!(!lead.is_zero());
        }
    }
}