mathhook_core/algebra/
gcd.rs

1//! Greatest Common Divisor operations for polynomials and expressions
2
3use crate::core::polynomial::IntPoly;
4use crate::core::polynomial::PolynomialError;
5use crate::core::{Expression, Number, Symbol};
6use num_integer::Integer;
7
8/// Trait for GCD operations on expressions
9pub trait PolynomialGcd {
10    fn gcd(&self, other: &Self) -> Self;
11    fn lcm(&self, other: &Self) -> Self;
12    fn factor_gcd(&self) -> Self;
13    fn cofactors(&self, other: &Self) -> (Expression, Expression, Expression);
14
15    /// Divides this polynomial by another, returning (quotient, remainder).
16    ///
17    /// Performs polynomial long division with respect to the specified variable.
18    /// Returns a tuple (quotient, remainder) satisfying the division identity:
19    /// `dividend = divisor * quotient + remainder` where `degree(remainder) < degree(divisor)`.
20    ///
21    /// # Arguments
22    ///
23    /// * `divisor` - The polynomial to divide by (must be non-zero)
24    /// * `var` - The variable to treat as the polynomial variable
25    ///
26    /// # Examples
27    ///
28    /// ```rust
29    /// use mathhook_core::{symbol, expr, Expression, algebra::gcd::PolynomialGcd};
30    ///
31    /// let x = symbol!(x);
32    /// // Divide (x^2 + 3x + 2) by (x + 1)
33    /// // Expected: (x + 1)(x + 2) = x^2 + 3x + 2
34    /// let dividend = expr!((x^2) + (3*x) + 2);
35    /// let divisor = expr!(x + 1);
36    /// let (quotient, remainder) = dividend.div_polynomial(&divisor, &x);
37    /// assert_eq!(remainder, Expression::integer(0));
38    /// ```
39    ///
40    /// # Returns
41    ///
42    /// Returns `(quotient, remainder)` tuple where both are expressions
43    fn div_polynomial(&self, divisor: &Expression, var: &Symbol) -> (Expression, Expression);
44
45    /// Returns the quotient of polynomial division.
46    ///
47    /// Computes only the quotient part of polynomial division, discarding the remainder.
48    /// Equivalent to `div_polynomial(divisor, var).0`.
49    ///
50    /// # Arguments
51    ///
52    /// * `divisor` - The polynomial to divide by
53    /// * `var` - The variable to treat as the polynomial variable
54    ///
55    /// # Examples
56    ///
57    /// ```rust
58    /// use mathhook_core::{symbol, expr, algebra::gcd::PolynomialGcd};
59    ///
60    /// let x = symbol!(x);
61    /// let dividend = expr!((x^2) - 1);
62    /// let divisor = expr!(x - 1);
63    /// let quotient = dividend.quo_polynomial(&divisor, &x);
64    /// // quotient = x + 1
65    /// ```
66    ///
67    /// # Returns
68    ///
69    /// Returns the quotient expression
70    fn quo_polynomial(&self, divisor: &Expression, var: &Symbol) -> Expression;
71
72    /// Returns the remainder of polynomial division.
73    ///
74    /// Computes only the remainder part of polynomial division, discarding the quotient.
75    /// Equivalent to `div_polynomial(divisor, var).1`.
76    ///
77    /// # Arguments
78    ///
79    /// * `divisor` - The polynomial to divide by
80    /// * `var` - The variable to treat as the polynomial variable
81    ///
82    /// # Examples
83    ///
84    /// ```rust
85    /// use mathhook_core::{symbol, expr, Expression, algebra::gcd::PolynomialGcd};
86    ///
87    /// let x = symbol!(x);
88    /// let dividend = expr!((x^2) + 1);
89    /// let divisor = expr!(x - 1);
90    /// let remainder = dividend.rem_polynomial(&divisor, &x);
91    /// assert_eq!(remainder, Expression::integer(2));
92    /// ```
93    ///
94    /// # Returns
95    ///
96    /// Returns the remainder expression
97    fn rem_polynomial(&self, divisor: &Expression, var: &Symbol) -> Expression;
98}
99
100impl PolynomialGcd for Expression {
101    /// GCD using IntPoly fast-path (primary path)
102    ///
103    /// Converts to IntPoly at API boundary, performs pure numeric GCD,
104    /// converts result back. NO Expression tree walking in fast path.
105    #[inline(always)]
106    fn gcd(&self, other: &Self) -> Self {
107        // Numeric GCD (most common case)
108        if let (Expression::Number(Number::Integer(a)), Expression::Number(Number::Integer(b))) =
109            (self, other)
110        {
111            return Expression::integer(a.gcd(b));
112        }
113
114        if self == other {
115            return self.clone();
116        }
117
118        if self.is_zero() {
119            return other.clone();
120        }
121        if other.is_zero() {
122            return self.clone();
123        }
124
125        // IntPoly fast-path: PRIMARY PATH for polynomial GCD
126        let vars = self.find_variables();
127        if vars.len() == 1 {
128            let var = &vars[0];
129            if IntPoly::can_convert(self, var) && IntPoly::can_convert(other, var) {
130                if let (Some(poly1), Some(poly2)) = (
131                    IntPoly::try_from_expression(self, var),
132                    IntPoly::try_from_expression(other, var),
133                ) {
134                    if let Ok(gcd_poly) = poly1.gcd_i64(&poly2) {
135                        let result = gcd_poly.to_expression(var);
136                        // Ensure GCD of coprime polynomials returns positive 1, not -1
137                        if let Expression::Number(Number::Integer(n)) = &result {
138                            if *n < 0 {
139                                return Expression::integer(n.abs());
140                            }
141                        }
142                        return result;
143                    }
144                }
145            }
146        }
147
148        // Minimal symbolic fallback for rational coefficient polynomials
149        if vars.len() == 1 {
150            let var = &vars[0];
151            let result = symbolic_gcd_euclidean(self, other, var);
152            // Final normalization: ensure we return +1 not -1 for coprime
153            if let Expression::Number(Number::Integer(n)) = &result {
154                if *n < 0 {
155                    return Expression::integer(n.abs());
156                }
157            }
158            return result;
159        }
160
161        // For multivariate or non-polynomial cases, return 1 (coprime)
162        Expression::integer(1)
163    }
164
165    /// Least Common Multiple
166    #[inline(always)]
167    fn lcm(&self, other: &Self) -> Self {
168        let gcd_val = self.gcd(other);
169
170        if gcd_val.is_zero() {
171            return Expression::integer(0);
172        }
173
174        let product = Expression::mul(vec![self.clone(), other.clone()]);
175        Expression::div(product, gcd_val)
176    }
177
178    /// Factor out GCD from expression
179    #[inline(always)]
180    fn factor_gcd(&self) -> Self {
181        match self {
182            Expression::Add(terms) => {
183                if terms.len() < 2 {
184                    return self.clone();
185                }
186
187                let mut common_gcd = terms[0].clone();
188                for term in &terms[1..] {
189                    common_gcd = common_gcd.gcd(term);
190                    if common_gcd.is_one() {
191                        return self.clone();
192                    }
193                }
194
195                common_gcd
196            }
197            Expression::Mul(_factors) => self.clone(),
198            _ => self.clone(),
199        }
200    }
201
202    /// Compute GCD and cofactors: returns (gcd, a/gcd, b/gcd)
203    fn cofactors(&self, other: &Self) -> (Expression, Expression, Expression) {
204        let gcd_val = self.gcd(other);
205
206        if gcd_val.is_zero() || gcd_val.is_one() {
207            return (gcd_val, self.clone(), other.clone());
208        }
209
210        // IntPoly fast-path for cofactors
211        let vars = self.find_variables();
212        if vars.len() == 1 {
213            let var = &vars[0];
214            if IntPoly::can_convert(self, var)
215                && IntPoly::can_convert(other, var)
216                && IntPoly::can_convert(&gcd_val, var)
217            {
218                if let (Some(poly1), Some(poly2), Some(gcd_poly)) = (
219                    IntPoly::try_from_expression(self, var),
220                    IntPoly::try_from_expression(other, var),
221                    IntPoly::try_from_expression(&gcd_val, var),
222                ) {
223                    if let (Ok((cof1, _)), Ok((cof2, _))) =
224                        (poly1.div_rem(&gcd_poly), poly2.div_rem(&gcd_poly))
225                    {
226                        return (gcd_val, cof1.to_expression(var), cof2.to_expression(var));
227                    }
228                }
229            }
230
231            // Symbolic fallback for cofactors
232            let (q1, r1) = crate::algebra::polynomial_division::polynomial_div(self, &gcd_val, var);
233            let (q2, r2) =
234                crate::algebra::polynomial_division::polynomial_div(other, &gcd_val, var);
235            if r1.is_zero() && r2.is_zero() {
236                return (gcd_val, q1, q2);
237            }
238        }
239
240        (gcd_val, self.clone(), other.clone())
241    }
242
243    fn div_polynomial(&self, divisor: &Expression, var: &Symbol) -> (Expression, Expression) {
244        crate::algebra::polynomial_division::polynomial_div(self, divisor, var)
245    }
246
247    fn quo_polynomial(&self, divisor: &Expression, var: &Symbol) -> Expression {
248        self.div_polynomial(divisor, var).0
249    }
250
251    fn rem_polynomial(&self, divisor: &Expression, var: &Symbol) -> Expression {
252        self.div_polynomial(divisor, var).1
253    }
254}
255
256/// Symbolic Euclidean GCD algorithm
257///
258/// Minimal fallback for rational coefficient polynomials.
259/// Uses polynomial remainder operations.
260fn symbolic_gcd_euclidean(p1: &Expression, p2: &Expression, var: &Symbol) -> Expression {
261    let mut a = p1.clone();
262    let mut b = p2.clone();
263
264    // Euclidean algorithm with at most 10 iterations
265    for _ in 0..10 {
266        if b.is_zero() {
267            // Normalize: for coprime polynomials, ensure we return +1, not -1
268            if let Expression::Number(Number::Integer(n)) = &a {
269                if *n < 0 {
270                    return Expression::integer(n.abs());
271                }
272            }
273            return a;
274        }
275
276        let remainder = crate::algebra::polynomial_division::polynomial_rem(&a, &b, var);
277        a = b;
278        b = remainder;
279    }
280
281    // Fallback: return 1 if algorithm doesn't converge
282    Expression::integer(1)
283}
284
285/// Polynomial GCD - routes to best algorithm automatically
286///
287/// Uses fast-path optimization to avoid expensive classification for simple cases.
288///
289/// # Arguments
290///
291/// * `p1` - First expression
292/// * `p2` - Second expression
293///
294/// # Returns
295///
296/// GCD of the two expressions
297///
298/// # Examples
299///
300/// ```rust
301/// use mathhook_core::algebra::gcd::polynomial_gcd;
302/// use mathhook_core::core::Expression;
303///
304/// let a = Expression::integer(12);
305/// let b = Expression::integer(18);
306/// let gcd = polynomial_gcd(&a, &b).unwrap();
307/// assert_eq!(gcd, Expression::integer(6));
308/// ```
309pub fn polynomial_gcd(p1: &Expression, p2: &Expression) -> Result<Expression, PolynomialError> {
310    use crate::expr;
311
312    // Fast path 1: Both integers
313    if let (Expression::Number(Number::Integer(n1)), Expression::Number(Number::Integer(n2))) =
314        (p1, p2)
315    {
316        return Ok(Expression::integer(n1.gcd(n2)));
317    }
318
319    // Fast path 2: Either is zero
320    if p1.is_zero() {
321        return Ok(p2.clone());
322    }
323    if p2.is_zero() {
324        return Ok(p1.clone());
325    }
326
327    // Fast path 3: Either is one
328    if p1.is_one() || p2.is_one() {
329        return Ok(expr!(1));
330    }
331
332    // Fast path 4: Both are same symbol
333    if let (Expression::Symbol(s1), Expression::Symbol(s2)) = (p1, p2) {
334        if s1 == s2 {
335            return Ok(p1.clone());
336        }
337        return Ok(expr!(1));
338    }
339
340    // Fast path 5: IntPoly for univariate integer polynomials
341    let vars_p1 = p1.find_variables();
342    if vars_p1.len() == 1 {
343        let var = &vars_p1[0];
344        if IntPoly::can_convert(p1, var) && IntPoly::can_convert(p2, var) {
345            if let (Some(poly1), Some(poly2)) = (
346                IntPoly::try_from_expression(p1, var),
347                IntPoly::try_from_expression(p2, var),
348            ) {
349                return Ok(poly1
350                    .gcd_i64(&poly2)
351                    .map_err(|e| PolynomialError::GcdComputationFailed {
352                        reason: format!("{:?}", e),
353                    })?
354                    .to_expression(var));
355            }
356        }
357    }
358
359    // Fallback to Expression::gcd()
360    Ok(p1.gcd(p2))
361}
362
363/// Univariate polynomial GCD with IntPoly fast-path
364///
365/// Uses IntPoly for integer coefficient polynomials, falls back to Expression::gcd().
366///
367/// # Arguments
368///
369/// * `p1` - First univariate polynomial
370/// * `p2` - Second univariate polynomial
371/// * `var` - The variable of the polynomials
372///
373/// # Returns
374///
375/// GCD of the two polynomials
376pub fn univariate_gcd(
377    p1: &Expression,
378    p2: &Expression,
379    var: &Symbol,
380) -> Result<Expression, PolynomialError> {
381    // IntPoly fast-path
382    if IntPoly::can_convert(p1, var) && IntPoly::can_convert(p2, var) {
383        if let (Some(poly1), Some(poly2)) = (
384            IntPoly::try_from_expression(p1, var),
385            IntPoly::try_from_expression(p2, var),
386        ) {
387            let gcd_poly =
388                poly1
389                    .gcd_i64(&poly2)
390                    .map_err(|e| PolynomialError::GcdComputationFailed {
391                        reason: format!("{:?}", e),
392                    })?;
393            return Ok(gcd_poly.to_expression(var));
394        }
395    }
396
397    // Fallback to Expression::gcd()
398    Ok(p1.gcd(p2))
399}
400
401/// Univariate polynomial GCD with cofactors
402///
403/// Returns (gcd, cofactor_p1, cofactor_p2) using IntPoly fast-path.
404///
405/// # Arguments
406///
407/// * `p1` - First univariate polynomial
408/// * `p2` - Second univariate polynomial
409/// * `var` - The variable of the polynomials
410///
411/// # Returns
412///
413/// Tuple of (gcd, cofactor_p1, cofactor_p2)
414pub fn univariate_gcd_modular(
415    p1: &Expression,
416    p2: &Expression,
417    var: &Symbol,
418) -> Result<(Expression, Expression, Expression), PolynomialError> {
419    // IntPoly fast-path with cofactors
420    if IntPoly::can_convert(p1, var) && IntPoly::can_convert(p2, var) {
421        if let (Some(poly1), Some(poly2)) = (
422            IntPoly::try_from_expression(p1, var),
423            IntPoly::try_from_expression(p2, var),
424        ) {
425            let gcd_poly =
426                poly1
427                    .gcd_i64(&poly2)
428                    .map_err(|e| PolynomialError::GcdComputationFailed {
429                        reason: format!("{:?}", e),
430                    })?;
431
432            // Compute cofactors
433            let (cof1, _) = poly1
434                .div_rem(&gcd_poly)
435                .map_err(|_| PolynomialError::DivisionByZero)?;
436            let (cof2, _) = poly2
437                .div_rem(&gcd_poly)
438                .map_err(|_| PolynomialError::DivisionByZero)?;
439
440            return Ok((
441                gcd_poly.to_expression(var),
442                cof1.to_expression(var),
443                cof2.to_expression(var),
444            ));
445        }
446    }
447
448    // Fallback
449    let gcd = p1.gcd(p2);
450    Ok((gcd.clone(), p1.clone(), p2.clone()))
451}
452
453#[cfg(test)]
454mod tests {
455    use super::*;
456    use crate::expr;
457    use crate::symbol;
458
459    #[test]
460    fn test_number_gcd() {
461        let a = Expression::integer(12);
462        let b = Expression::integer(8);
463        let result = a.gcd(&b);
464        assert_eq!(result, Expression::integer(4));
465
466        let a = Expression::integer(17);
467        let b = Expression::integer(13);
468        let result = a.gcd(&b);
469        assert_eq!(result, Expression::integer(1));
470    }
471
472    #[test]
473    fn test_gcd_with_zero() {
474        let a = Expression::integer(5);
475        let zero = Expression::integer(0);
476
477        let result = a.gcd(&zero);
478        assert_eq!(result, Expression::integer(5));
479
480        let result = zero.gcd(&a);
481        assert_eq!(result, Expression::integer(5));
482    }
483
484    #[test]
485    fn test_identical_expressions() {
486        let x = expr!(x);
487        let result = x.gcd(&x);
488        assert_eq!(result, x);
489    }
490
491    #[test]
492    fn test_gcd_performance_benchmark() {
493        use std::time::Instant;
494
495        let start = Instant::now();
496
497        for i in 1..10_000 {
498            let a = Expression::integer(i * 6);
499            let b = Expression::integer(i * 9);
500            let _result = a.gcd(&b);
501        }
502
503        let duration = start.elapsed();
504        let ops_per_sec = 10_000.0 / duration.as_secs_f64();
505
506        println!("GCD Performance: {:.2}M ops/sec", ops_per_sec / 1_000_000.0);
507
508        assert!(
509            ops_per_sec > 100_000.0,
510            "Expected >100K ops/sec, got {:.2}",
511            ops_per_sec
512        );
513    }
514
515    #[test]
516    fn test_polynomial_gcd_basic() {
517        let x = symbol!(x);
518
519        let poly1 = Expression::mul(vec![Expression::integer(6), Expression::symbol(x.clone())]);
520        let poly2 = Expression::mul(vec![Expression::integer(9), Expression::symbol(x.clone())]);
521
522        let result = poly1.gcd(&poly2);
523
524        println!("Polynomial GCD result: {}", result);
525        assert!(!result.is_zero());
526    }
527
528    #[test]
529    fn test_factor_gcd() {
530        let x = symbol!(x);
531
532        let term1 = Expression::mul(vec![Expression::integer(6), Expression::symbol(x.clone())]);
533        let term2 = Expression::mul(vec![Expression::integer(9), Expression::symbol(x.clone())]);
534        let sum = Expression::add(vec![term1, term2]);
535
536        let gcd_factor = sum.factor_gcd();
537        println!("Factored GCD: {}", gcd_factor);
538
539        assert!(!gcd_factor.is_zero());
540    }
541
542    #[test]
543    fn test_lcm_basic() {
544        let a = Expression::integer(6);
545        let b = Expression::integer(8);
546        let result = a.lcm(&b);
547
548        println!("LCM result: {}", result);
549        assert!(!result.is_zero());
550    }
551
552    #[test]
553    fn test_int_poly_fast_path() {
554        let _x = symbol!(x);
555
556        let poly1 = expr!((x ^ 5) + (2 * (x ^ 4)) + (3 * (x ^ 3)) + (4 * (x ^ 2)) + (5 * x) + 6);
557        let poly2 =
558            expr!((2 * (x ^ 5)) + (4 * (x ^ 4)) + (6 * (x ^ 3)) + (8 * (x ^ 2)) + (10 * x) + 12);
559
560        let result = poly1.gcd(&poly2);
561
562        println!("IntPoly fast-path GCD: {}", result);
563        assert!(!result.is_zero());
564    }
565
566    #[test]
567    fn test_intpoly_gcd_direct() {
568        let _x = symbol!(x);
569
570        let p1 = expr!((x ^ 2) - 1);
571        let p2 = expr!(x - 1);
572
573        let gcd = p1.gcd(&p2);
574
575        println!("Direct IntPoly GCD: {}", gcd);
576        assert!(!gcd.is_zero());
577    }
578
579    #[test]
580    fn test_cofactors_intpoly() {
581        let _x = symbol!(x);
582
583        let p1 = expr!((x ^ 2) - 1);
584        let p2 = expr!(x - 1);
585
586        let (gcd, cof1, cof2) = p1.cofactors(&p2);
587
588        println!("GCD: {}, Cofactor1: {}, Cofactor2: {}", gcd, cof1, cof2);
589        assert!(!gcd.is_zero());
590    }
591
592    #[test]
593    fn test_polynomial_gcd_function() {
594        let a = Expression::integer(12);
595        let b = Expression::integer(18);
596        let gcd = polynomial_gcd(&a, &b).unwrap();
597        assert_eq!(gcd, Expression::integer(6));
598    }
599
600    #[test]
601    fn test_univariate_gcd_function() {
602        let _x = symbol!(x);
603        let p1 = expr!((x ^ 2) - 1);
604        let p2 = expr!(x - 1);
605        let gcd = univariate_gcd(&p1, &p2, &_x).unwrap();
606        assert!(!gcd.is_zero());
607    }
608
609    #[test]
610    fn test_gcd_coprime_expressions() {
611        let _x = symbol!(x);
612        let a = expr!(x + 1);
613        let b = expr!(x + 2);
614        let result = a.gcd(&b);
615        println!("GCD result: {:?}", result);
616        assert_eq!(result, Expression::integer(1));
617    }
618}
mathhook_core/algebra/gcd.rs

mathhook_core/algebra/
gcd.rs
