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 - unwrap or fallback to original
232            if let (Ok((q1, r1)), Ok((q2, r2))) = (
233                crate::algebra::polynomial_division::polynomial_div(self, &gcd_val, var),
234                crate::algebra::polynomial_division::polynomial_div(other, &gcd_val, var),
235            ) {
236                if r1.is_zero() && r2.is_zero() {
237                    return (gcd_val, q1, q2);
238                }
239            }
240        }
241
242        (gcd_val, self.clone(), other.clone())
243    }
244
245    fn div_polynomial(&self, divisor: &Expression, var: &Symbol) -> (Expression, Expression) {
246        crate::algebra::polynomial_division::polynomial_div(self, divisor, var)
247            .unwrap_or_else(|_| (Expression::undefined(), Expression::undefined()))
248    }
249
250    fn quo_polynomial(&self, divisor: &Expression, var: &Symbol) -> Expression {
251        self.div_polynomial(divisor, var).0
252    }
253
254    fn rem_polynomial(&self, divisor: &Expression, var: &Symbol) -> Expression {
255        self.div_polynomial(divisor, var).1
256    }
257}
258
259/// Symbolic Euclidean GCD algorithm
260///
261/// Minimal fallback for rational coefficient polynomials.
262/// Uses polynomial remainder operations.
263fn symbolic_gcd_euclidean(p1: &Expression, p2: &Expression, var: &Symbol) -> Expression {
264    let mut a = p1.clone();
265    let mut b = p2.clone();
266
267    // Euclidean algorithm with at most 10 iterations
268    for _ in 0..10 {
269        if b.is_zero() {
270            // Normalize: for coprime polynomials, ensure we return +1, not -1
271            if let Expression::Number(Number::Integer(n)) = &a {
272                if *n < 0 {
273                    return Expression::integer(n.abs());
274                }
275            }
276            return a;
277        }
278
279        let remainder =
280            crate::algebra::polynomial_division::polynomial_rem(&a, &b, var).unwrap_or(b.clone());
281        a = b;
282        b = remainder;
283    }
284
285    // Fallback: return 1 if algorithm doesn't converge
286    Expression::integer(1)
287}
288
289/// Polynomial GCD - routes to best algorithm automatically
290///
291/// Uses fast-path optimization to avoid expensive classification for simple cases.
292///
293/// # Arguments
294///
295/// * `p1` - First expression
296/// * `p2` - Second expression
297///
298/// # Returns
299///
300/// GCD of the two expressions
301///
302/// # Examples
303///
304/// ```rust
305/// use mathhook_core::algebra::gcd::polynomial_gcd;
306/// use mathhook_core::core::Expression;
307///
308/// let a = Expression::integer(12);
309/// let b = Expression::integer(18);
310/// let gcd = polynomial_gcd(&a, &b).unwrap();
311/// assert_eq!(gcd, Expression::integer(6));
312/// ```
313pub fn polynomial_gcd(p1: &Expression, p2: &Expression) -> Result<Expression, PolynomialError> {
314    use crate::expr;
315
316    // Fast path 1: Both integers
317    if let (Expression::Number(Number::Integer(n1)), Expression::Number(Number::Integer(n2))) =
318        (p1, p2)
319    {
320        return Ok(Expression::integer(n1.gcd(n2)));
321    }
322
323    // Fast path 2: Either is zero
324    if p1.is_zero() {
325        return Ok(p2.clone());
326    }
327    if p2.is_zero() {
328        return Ok(p1.clone());
329    }
330
331    // Fast path 3: Either is one
332    if p1.is_one() || p2.is_one() {
333        return Ok(expr!(1));
334    }
335
336    // Fast path 4: Both are same symbol
337    if let (Expression::Symbol(s1), Expression::Symbol(s2)) = (p1, p2) {
338        if s1 == s2 {
339            return Ok(p1.clone());
340        }
341        return Ok(expr!(1));
342    }
343
344    // Fast path 5: IntPoly for univariate integer polynomials
345    let vars_p1 = p1.find_variables();
346    if vars_p1.len() == 1 {
347        let var = &vars_p1[0];
348        if IntPoly::can_convert(p1, var) && IntPoly::can_convert(p2, var) {
349            if let (Some(poly1), Some(poly2)) = (
350                IntPoly::try_from_expression(p1, var),
351                IntPoly::try_from_expression(p2, var),
352            ) {
353                return Ok(poly1
354                    .gcd_i64(&poly2)
355                    .map_err(|e| PolynomialError::GcdComputationFailed {
356                        reason: format!("{:?}", e),
357                    })?
358                    .to_expression(var));
359            }
360        }
361    }
362
363    // Fallback to Expression::gcd()
364    Ok(p1.gcd(p2))
365}
366
367/// Univariate polynomial GCD with IntPoly fast-path
368///
369/// Uses IntPoly for integer coefficient polynomials, falls back to Expression::gcd().
370///
371/// # Arguments
372///
373/// * `p1` - First univariate polynomial
374/// * `p2` - Second univariate polynomial
375/// * `var` - The variable of the polynomials
376///
377/// # Returns
378///
379/// GCD of the two polynomials
380pub fn univariate_gcd(
381    p1: &Expression,
382    p2: &Expression,
383    var: &Symbol,
384) -> Result<Expression, PolynomialError> {
385    // IntPoly fast-path
386    if IntPoly::can_convert(p1, var) && IntPoly::can_convert(p2, var) {
387        if let (Some(poly1), Some(poly2)) = (
388            IntPoly::try_from_expression(p1, var),
389            IntPoly::try_from_expression(p2, var),
390        ) {
391            let gcd_poly =
392                poly1
393                    .gcd_i64(&poly2)
394                    .map_err(|e| PolynomialError::GcdComputationFailed {
395                        reason: format!("{:?}", e),
396                    })?;
397            return Ok(gcd_poly.to_expression(var));
398        }
399    }
400
401    // Fallback to Expression::gcd()
402    Ok(p1.gcd(p2))
403}
404
405/// Univariate polynomial GCD with cofactors
406///
407/// Returns (gcd, cofactor_p1, cofactor_p2) using IntPoly fast-path.
408///
409/// # Arguments
410///
411/// * `p1` - First univariate polynomial
412/// * `p2` - Second univariate polynomial
413/// * `var` - The variable of the polynomials
414///
415/// # Returns
416///
417/// Tuple of (gcd, cofactor_p1, cofactor_p2)
418pub fn univariate_gcd_modular(
419    p1: &Expression,
420    p2: &Expression,
421    var: &Symbol,
422) -> Result<(Expression, Expression, Expression), PolynomialError> {
423    // IntPoly fast-path with cofactors
424    if IntPoly::can_convert(p1, var) && IntPoly::can_convert(p2, var) {
425        if let (Some(poly1), Some(poly2)) = (
426            IntPoly::try_from_expression(p1, var),
427            IntPoly::try_from_expression(p2, var),
428        ) {
429            let gcd_poly =
430                poly1
431                    .gcd_i64(&poly2)
432                    .map_err(|e| PolynomialError::GcdComputationFailed {
433                        reason: format!("{:?}", e),
434                    })?;
435
436            // Compute cofactors
437            let (cof1, _) = poly1
438                .div_rem(&gcd_poly)
439                .map_err(|_| PolynomialError::DivisionByZero)?;
440            let (cof2, _) = poly2
441                .div_rem(&gcd_poly)
442                .map_err(|_| PolynomialError::DivisionByZero)?;
443
444            return Ok((
445                gcd_poly.to_expression(var),
446                cof1.to_expression(var),
447                cof2.to_expression(var),
448            ));
449        }
450    }
451
452    // Fallback
453    let gcd = p1.gcd(p2);
454    Ok((gcd.clone(), p1.clone(), p2.clone()))
455}
456
457#[cfg(test)]
458mod tests {
459    use super::*;
460    use crate::expr;
461    use crate::symbol;
462
463    #[test]
464    fn test_number_gcd() {
465        let a = Expression::integer(12);
466        let b = Expression::integer(8);
467        let result = a.gcd(&b);
468        assert_eq!(result, Expression::integer(4));
469
470        let a = Expression::integer(17);
471        let b = Expression::integer(13);
472        let result = a.gcd(&b);
473        assert_eq!(result, Expression::integer(1));
474    }
475
476    #[test]
477    fn test_gcd_with_zero() {
478        let a = Expression::integer(5);
479        let zero = Expression::integer(0);
480
481        let result = a.gcd(&zero);
482        assert_eq!(result, Expression::integer(5));
483
484        let result = zero.gcd(&a);
485        assert_eq!(result, Expression::integer(5));
486    }
487
488    #[test]
489    fn test_identical_expressions() {
490        let x = expr!(x);
491        let result = x.gcd(&x);
492        assert_eq!(result, x);
493    }
494
495    #[test]
496    fn test_gcd_performance_benchmark() {
497        use std::time::Instant;
498
499        let start = Instant::now();
500
501        for i in 1..10_000 {
502            let a = Expression::integer(i * 6);
503            let b = Expression::integer(i * 9);
504            let _result = a.gcd(&b);
505        }
506
507        let duration = start.elapsed();
508        let ops_per_sec = 10_000.0 / duration.as_secs_f64();
509
510        println!("GCD Performance: {:.2}M ops/sec", ops_per_sec / 1_000_000.0);
511
512        assert!(
513            ops_per_sec > 100_000.0,
514            "Expected >100K ops/sec, got {:.2}",
515            ops_per_sec
516        );
517    }
518
519    #[test]
520    fn test_polynomial_gcd_basic() {
521        let x = symbol!(x);
522
523        let poly1 = Expression::mul(vec![Expression::integer(6), Expression::symbol(x.clone())]);
524        let poly2 = Expression::mul(vec![Expression::integer(9), Expression::symbol(x.clone())]);
525
526        let result = poly1.gcd(&poly2);
527
528        println!("Polynomial GCD result: {}", result);
529        assert!(!result.is_zero());
530    }
531
532    #[test]
533    fn test_factor_gcd() {
534        let x = symbol!(x);
535
536        let term1 = Expression::mul(vec![Expression::integer(6), Expression::symbol(x.clone())]);
537        let term2 = Expression::mul(vec![Expression::integer(9), Expression::symbol(x.clone())]);
538        let sum = Expression::add(vec![term1, term2]);
539
540        let gcd_factor = sum.factor_gcd();
541        println!("Factored GCD: {}", gcd_factor);
542
543        assert!(!gcd_factor.is_zero());
544    }
545
546    #[test]
547    fn test_lcm_basic() {
548        let a = Expression::integer(6);
549        let b = Expression::integer(8);
550        let result = a.lcm(&b);
551
552        println!("LCM result: {}", result);
553        assert!(!result.is_zero());
554    }
555
556    #[test]
557    fn test_int_poly_fast_path() {
558        let _x = symbol!(x);
559
560        let poly1 = expr!((x ^ 5) + (2 * (x ^ 4)) + (3 * (x ^ 3)) + (4 * (x ^ 2)) + (5 * x) + 6);
561        let poly2 =
562            expr!((2 * (x ^ 5)) + (4 * (x ^ 4)) + (6 * (x ^ 3)) + (8 * (x ^ 2)) + (10 * x) + 12);
563
564        let result = poly1.gcd(&poly2);
565
566        println!("IntPoly fast-path GCD: {}", result);
567        assert!(!result.is_zero());
568    }
569
570    #[test]
571    fn test_intpoly_gcd_direct() {
572        let _x = symbol!(x);
573
574        let p1 = expr!((x ^ 2) - 1);
575        let p2 = expr!(x - 1);
576
577        let gcd = p1.gcd(&p2);
578
579        println!("Direct IntPoly GCD: {}", gcd);
580        assert!(!gcd.is_zero());
581    }
582
583    #[test]
584    fn test_cofactors_intpoly() {
585        let _x = symbol!(x);
586
587        let p1 = expr!((x ^ 2) - 1);
588        let p2 = expr!(x - 1);
589
590        let (gcd, cof1, cof2) = p1.cofactors(&p2);
591
592        println!("GCD: {}, Cofactor1: {}, Cofactor2: {}", gcd, cof1, cof2);
593        assert!(!gcd.is_zero());
594    }
595
596    #[test]
597    fn test_polynomial_gcd_function() {
598        let a = Expression::integer(12);
599        let b = Expression::integer(18);
600        let gcd = polynomial_gcd(&a, &b).unwrap();
601        assert_eq!(gcd, Expression::integer(6));
602    }
603
604    #[test]
605    fn test_univariate_gcd_function() {
606        let _x = symbol!(x);
607        let p1 = expr!((x ^ 2) - 1);
608        let p2 = expr!(x - 1);
609        let gcd = univariate_gcd(&p1, &p2, &_x).unwrap();
610        assert!(!gcd.is_zero());
611    }
612
613    #[test]
614    fn test_gcd_coprime_expressions() {
615        let _x = symbol!(x);
616        let a = expr!(x + 1);
617        let b = expr!(x + 2);
618        let result = a.gcd(&b);
619        println!("GCD result: {:?}", result);
620        assert_eq!(result, Expression::integer(1));
621    }
622}
mathhook_core/algebra/gcd.rs

mathhook_core/algebra/
gcd.rs
