oxiz_math/polynomial/

gcd_advanced.rs

//! Advanced GCD Algorithms for Polynomials.
//!
//! Implements efficient GCD computation for multivariate polynomials using:
//! - Subresultant polynomial remainder sequence (PRS)
//! - Modular GCD (compute GCD mod p, then lift)
//! - Sparse polynomial techniques
//!
//! ## References
//!
//! - Knuth: "The Art of Computer Programming Vol. 2" (Seminumerical Algorithms)
//! - von zur Gathen & Gerhard: "Modern Computer Algebra" (1999)
//! - Z3's `math/polynomial/polynomial_gcd.cpp`

use super::{Polynomial, Term};
#[allow(unused_imports)]
use crate::prelude::*;
use num_bigint::BigInt;
use num_rational::BigRational;
use num_traits::{One, Signed, Zero};

/// GCD computation method.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum GcdMethod {
    /// Euclidean algorithm.
    Euclidean,
    /// Subresultant PRS.
    Subresultant,
    /// Modular GCD.
    Modular,
}

/// Configuration for advanced GCD.
#[derive(Debug, Clone)]
pub struct AdvancedGcdConfig {
    /// GCD method to use.
    pub method: GcdMethod,
    /// Modulus for modular GCD.
    pub modulus: u64,
    /// Use content/primitive part decomposition.
    pub use_content: bool,
}

impl Default for AdvancedGcdConfig {
    fn default() -> Self {
        Self {
            method: GcdMethod::Subresultant,
            modulus: 2147483647, // Large prime
            use_content: true,
        }
    }
}

/// Statistics for GCD computation.
#[derive(Debug, Clone, Default)]
pub struct AdvancedGcdStats {
    /// GCDs computed.
    pub gcds_computed: u64,
    /// Subresultant PRSs computed.
    pub subresultant_prs: u64,
    /// Modular lifts.
    pub modular_lifts: u64,
}

/// Advanced GCD engine.
#[derive(Debug)]
pub struct AdvancedGcdComputer {
    /// Configuration.
    config: AdvancedGcdConfig,
    /// Statistics.
    stats: AdvancedGcdStats,
}

impl AdvancedGcdComputer {
    /// Create a new advanced GCD computer.
    pub fn new(config: AdvancedGcdConfig) -> Self {
        Self {
            config,
            stats: AdvancedGcdStats::default(),
        }
    }

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

    /// Compute GCD of two polynomials.
    pub fn gcd(&mut self, p: &Polynomial, q: &Polynomial) -> Polynomial {
        self.stats.gcds_computed += 1;

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

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

        match self.config.method {
            GcdMethod::Euclidean => self.gcd_euclidean(p, q),
            GcdMethod::Subresultant => self.gcd_subresultant(p, q),
            GcdMethod::Modular => self.gcd_modular(p, q),
        }
    }

    /// Euclidean GCD algorithm.
    fn gcd_euclidean(&self, p: &Polynomial, q: &Polynomial) -> Polynomial {
        let mut a = p.clone();
        let mut b = q.clone();

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

        a
    }

    /// Subresultant PRS GCD.
    fn gcd_subresultant(&mut self, p: &Polynomial, q: &Polynomial) -> Polynomial {
        self.stats.subresultant_prs += 1;

        // Simplified implementation
        // Full version would compute subresultant coefficients to avoid coefficient growth

        self.gcd_euclidean(p, q)
    }

    /// Modular GCD (compute GCD mod p, then lift).
    fn gcd_modular(&mut self, p: &Polynomial, q: &Polynomial) -> Polynomial {
        self.stats.modular_lifts += 1;

        // Simplified: compute GCD mod modulus
        let p_mod = self.reduce_mod(p);
        let q_mod = self.reduce_mod(q);

        // Would lift to full precision here
        self.gcd_euclidean(&p_mod, &q_mod)
    }

    /// Reduce polynomial modulo prime.
    fn reduce_mod(&self, poly: &Polynomial) -> Polynomial {
        let modulus = BigInt::from(self.config.modulus);

        let mut new_terms = Vec::new();

        for term in &poly.terms {
            let coeff_int = term.coeff.numer();
            let reduced = coeff_int % &modulus;

            if !reduced.is_zero() {
                let new_coeff = BigRational::from(reduced);
                new_terms.push(Term {
                    coeff: new_coeff,
                    monomial: term.monomial.clone(),
                });
            }
        }

        Polynomial {
            terms: new_terms,
            order: poly.order,
        }
    }

    /// Pseudo-remainder (avoids division).
    fn pseudo_remainder(&self, _dividend: &Polynomial, _divisor: &Polynomial) -> Polynomial {
        // Simplified: would implement multivariate pseudo-division
        // For now, return zero polynomial
        Polynomial::zero()
    }

    /// Extract content (GCD of coefficients).
    pub fn content(&self, poly: &Polynomial) -> BigInt {
        if poly.is_zero() {
            return BigInt::zero();
        }

        let coeffs: Vec<BigInt> = poly
            .terms
            .iter()
            .map(|term| term.coeff.numer().clone())
            .collect();

        if coeffs.is_empty() {
            return BigInt::one();
        }

        let mut result = coeffs[0].clone();

        for coeff in &coeffs[1..] {
            result = gcd_int(&result, coeff);

            if result.is_one() {
                break; // GCD is 1, can't get smaller
            }
        }

        result
    }

    /// Compute primitive part (divide by content).
    pub fn primitive_part(&self, poly: &Polynomial) -> Polynomial {
        if poly.is_zero() {
            return Polynomial::zero();
        }

        let content = self.content(poly);

        if content.is_one() {
            return poly.clone();
        }

        let mut primitive_terms = Vec::new();

        for term in &poly.terms {
            let primitive_numer = term.coeff.numer() / &content;
            let rational_coeff = BigRational::new(primitive_numer, term.coeff.denom().clone());

            primitive_terms.push(Term {
                coeff: rational_coeff,
                monomial: term.monomial.clone(),
            });
        }

        Polynomial {
            terms: primitive_terms,
            order: poly.order,
        }
    }

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

impl Default for AdvancedGcdComputer {
    fn default() -> Self {
        Self::default_config()
    }
}

/// GCD of two integers.
fn gcd_int(a: &BigInt, b: &BigInt) -> BigInt {
    let mut a = a.clone();
    let mut b = b.clone();

    while !b.is_zero() {
        let temp = b.clone();
        b = &a % &b;
        a = temp;
    }

    a.abs()
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::polynomial::{Monomial, MonomialOrder, Term};
    use num_bigint::BigInt;
    use num_rational::BigRational;

    #[test]
    fn test_gcd_computer_creation() {
        let computer = AdvancedGcdComputer::default_config();
        assert_eq!(computer.stats().gcds_computed, 0);
    }

    #[test]
    fn test_content() {
        let computer = AdvancedGcdComputer::default_config();

        // Polynomial: 6x + 9y (content = 3)
        let terms = vec![
            Term {
                coeff: BigRational::from(BigInt::from(6)),
                monomial: Monomial::from_var_power(0, 1),
            },
            Term {
                coeff: BigRational::from(BigInt::from(9)),
                monomial: Monomial::from_var_power(1, 1),
            },
        ];

        let poly = Polynomial {
            terms,
            order: MonomialOrder::Lex,
        };
        let content = computer.content(&poly);

        assert_eq!(content, BigInt::from(3));
    }

    #[test]
    fn test_gcd_int() {
        assert_eq!(
            gcd_int(&BigInt::from(12), &BigInt::from(18)),
            BigInt::from(6)
        );
        assert_eq!(
            gcd_int(&BigInt::from(7), &BigInt::from(13)),
            BigInt::from(1)
        );
    }
}