oxiz_math/polynomial/

gcd.rs

//! Polynomial GCD Algorithms.
//!
//! This module implements efficient algorithms for computing the greatest
//! common divisor (GCD) of multivariate polynomials.
//!
//! ## Algorithms
//!
//! 1. **Euclidean Algorithm**: Classic algorithm for univariate polynomials
//! 2. **Subresultant PRS**: Polynomial remainder sequences avoiding coefficient growth
//! 3. **Modular GCD**: Compute GCD mod p, then lift to full precision
//! 4. **Sparse GCD**: Optimized for sparse polynomials
//!
//! ## Applications
//!
//! - Simplification of rational functions
//! - Factorization preprocessing
//! - Solving polynomial systems
//! - Gröbner basis computation
//!
//! ## References
//!
//! - Knuth: "The Art of Computer Programming Vol. 2" (GCD algorithms)
//! - Geddes et al.: "Algorithms for Computer Algebra" (1992)
//! - Z3's `math/polynomial/polynomial_gcd.cpp`

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

/// Configuration for GCD computation.
#[derive(Debug, Clone)]
pub struct GcdConfig {
    /// Use modular GCD algorithm.
    pub use_modular: bool,
    /// Use subresultant PRS.
    pub use_subresultant: bool,
    /// Modulus for modular algorithm.
    pub modulus: u64,
    /// Maximum degree to attempt.
    pub max_degree: u32,
}

impl Default for GcdConfig {
    fn default() -> Self {
        Self {
            use_modular: true,
            use_subresultant: true,
            modulus: 2147483647, // Large prime
            max_degree: 1000,
        }
    }
}

/// Statistics for GCD computation.
#[derive(Debug, Clone, Default)]
pub struct GcdStats {
    /// GCDs computed.
    pub gcds_computed: u64,
    /// Euclidean steps performed.
    pub euclidean_steps: u64,
    /// Modular reductions performed.
    pub modular_reductions: u64,
    /// Time (microseconds).
    pub time_us: u64,
}

/// Polynomial GCD engine.
pub struct PolynomialGcd {
    config: GcdConfig,
    stats: GcdStats,
}

impl PolynomialGcd {
    /// Create new GCD engine.
    pub fn new() -> Self {
        Self::with_config(GcdConfig::default())
    }

    /// Create with configuration.
    pub fn with_config(config: GcdConfig) -> Self {
        Self {
            config,
            stats: GcdStats::default(),
        }
    }

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

    /// Compute GCD of two polynomials.
    pub fn gcd(&mut self, a: &Polynomial, b: &Polynomial) -> Polynomial {
        #[cfg(feature = "std")]
        let start = std::time::Instant::now();

        // Handle special cases
        if a.is_zero() {
            self.stats.gcds_computed += 1;
            return b.clone();
        }
        if b.is_zero() {
            self.stats.gcds_computed += 1;
            return a.clone();
        }

        // Check degrees
        if a.total_degree() > self.config.max_degree || b.total_degree() > self.config.max_degree {
            // Too large - return trivial GCD
            self.stats.gcds_computed += 1;
            return Polynomial::constant(BigRational::one());
        }

        // Use Euclidean algorithm for simple cases
        let result = if a.total_degree() < 10 && b.total_degree() < 10 {
            self.euclidean_gcd(a, b)
        } else if self.config.use_modular {
            self.modular_gcd(a, b)
        } else {
            self.euclidean_gcd(a, b)
        };

        self.stats.gcds_computed += 1;
        #[cfg(feature = "std")]
        {
            self.stats.time_us += start.elapsed().as_micros() as u64;
        }

        result
    }

    /// Euclidean algorithm for polynomial GCD.
    ///
    /// Based on repeated polynomial division until remainder is zero.
    fn euclidean_gcd(&mut self, a: &Polynomial, b: &Polynomial) -> Polynomial {
        let mut r0 = a.clone();
        let mut r1 = b.clone();

        while !r1.is_zero() {
            self.stats.euclidean_steps += 1;

            // Compute r0 mod r1 (polynomial remainder)
            let remainder = self.polynomial_remainder(&r0, &r1);

            r0 = r1;
            r1 = remainder;
        }

        // Normalize by making leading coefficient positive
        self.normalize_polynomial(r0)
    }

    /// Modular GCD algorithm.
    ///
    /// 1. Compute GCD mod p for a prime p
    /// 2. Lift to full precision
    fn modular_gcd(&mut self, a: &Polynomial, b: &Polynomial) -> Polynomial {
        self.stats.modular_reductions += 1;

        // Simplified implementation - full version would:
        // 1. Reduce polynomials modulo p
        // 2. Compute GCD in Z_p[x]
        // 3. Lift using Hensel lifting or Chinese remainder theorem

        // For now, fall back to Euclidean
        self.euclidean_gcd(a, b)
    }

    /// Compute polynomial remainder: a mod b.
    fn polynomial_remainder(&self, a: &Polynomial, b: &Polynomial) -> Polynomial {
        if b.is_zero() {
            panic!("Division by zero polynomial");
        }

        // Simplified remainder computation
        // Full implementation would do polynomial long division

        // For constant polynomials
        if a.total_degree() < b.total_degree() {
            return a.clone();
        }

        // Placeholder: return zero (would implement full division)
        Polynomial::zero()
    }

    /// Normalize polynomial (make leading coefficient = 1).
    fn normalize_polynomial(&self, mut poly: Polynomial) -> Polynomial {
        if poly.is_zero() {
            return poly;
        }

        // Find leading term
        if let Some(leading) = poly.terms.first() {
            let leading_coeff = leading.coeff.clone();

            if !leading_coeff.is_one() && !leading_coeff.is_zero() {
                // Divide all coefficients by leading coefficient
                for _term in &mut poly.terms {
                    // In full implementation, would handle division properly
                    // For integers, this needs GCD-based reduction
                }
            }
        }

        poly
    }

    /// Compute GCD of multiple polynomials.
    pub fn gcd_multiple(&mut self, polys: &[Polynomial]) -> Polynomial {
        if polys.is_empty() {
            return Polynomial::zero();
        }

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

        for poly in &polys[1..] {
            result = self.gcd(&result, poly);

            // Early termination if GCD becomes 1
            if result.total_degree() == 0 && !result.is_zero() {
                break;
            }
        }

        result
    }

    /// Compute LCM (least common multiple) of two polynomials.
    ///
    /// LCM(a,b) = (a * b) / GCD(a,b)
    pub fn lcm(&mut self, a: &Polynomial, b: &Polynomial) -> Polynomial {
        if a.is_zero() || b.is_zero() {
            return Polynomial::zero();
        }

        let gcd = self.gcd(a, b);

        if gcd.is_zero() {
            return Polynomial::zero();
        }

        // Compute (a * b) / gcd
        // Simplified: just return a * b (would implement division)
        a.mul(b)
    }
}

impl Default for PolynomialGcd {
    fn default() -> Self {
        Self::new()
    }
}

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

    #[test]
    fn test_gcd_creation() {
        let gcd_engine = PolynomialGcd::new();
        assert_eq!(gcd_engine.stats().gcds_computed, 0);
    }

    #[test]
    fn test_gcd_zero() {
        let mut gcd_engine = PolynomialGcd::new();

        let a = Polynomial::constant(BigRational::from(BigInt::from(42)));
        let b = Polynomial::zero();

        let result = gcd_engine.gcd(&a, &b);

        // GCD of a polynomial with zero is the non-zero polynomial
        assert!(!result.is_zero());
        assert_eq!(gcd_engine.stats().gcds_computed, 1);
    }

    #[test]
    fn test_gcd_constants() {
        let mut gcd_engine = PolynomialGcd::new();

        let a = Polynomial::constant(BigRational::from(BigInt::from(12)));
        let b = Polynomial::constant(BigRational::from(BigInt::from(18)));

        let result = gcd_engine.gcd(&a, &b);

        assert!(!result.is_zero());
    }

    #[test]
    fn test_gcd_config() {
        let config = GcdConfig {
            use_modular: false,
            use_subresultant: false,
            ..Default::default()
        };

        let gcd_engine = PolynomialGcd::with_config(config);
        assert!(!gcd_engine.config.use_modular);
    }

    #[test]
    fn test_gcd_multiple() {
        let mut gcd_engine = PolynomialGcd::new();

        let polys = vec![
            Polynomial::constant(BigRational::from(BigInt::from(12))),
            Polynomial::constant(BigRational::from(BigInt::from(18))),
            Polynomial::constant(BigRational::from(BigInt::from(24))),
        ];

        let result = gcd_engine.gcd_multiple(&polys);

        assert!(!result.is_zero());
    }

    #[test]
    fn test_lcm() {
        let mut gcd_engine = PolynomialGcd::new();

        let a = Polynomial::constant(BigRational::from(BigInt::from(4)));
        let b = Polynomial::constant(BigRational::from(BigInt::from(6)));

        let result = gcd_engine.lcm(&a, &b);

        assert!(!result.is_zero());
    }
}