mathhook_core/algebra/
polynomial_division.rs

1//! Polynomial long division operations
2//!
3//! Implements polynomial division algorithms for univariate polynomials,
4//! supporting the Euclidean GCD algorithm and general polynomial arithmetic.
5//!
6//! # Algorithm
7//!
8//! Uses IntPoly fast-path for univariate integer polynomials (primary path).
9//! Falls back to symbolic polynomial division for rational coefficient cases.
10//!
11//! # Example
12//!
13//! ```rust
14//! use mathhook_core::{expr, symbol};
15//! use mathhook_core::algebra::polynomial_division::polynomial_div;
16//!
17//! let x = symbol!(x);
18//! let (quotient, remainder) = polynomial_div(&expr!((x^2) - 1), &expr!(x - 1), &x)?;
19//! # Ok::<(), mathhook_core::error::MathError>(())
20//! ```
21
22use crate::core::polynomial::IntPoly;
23use crate::core::{Expression, Number, Symbol};
24use crate::error::MathError;
25use crate::simplify::Simplify;
26
27/// Polynomial long division
28///
29/// Returns (quotient, remainder) such that:
30/// `dividend = divisor * quotient + remainder`
31/// and `degree(remainder) < degree(divisor)`
32///
33/// Uses IntPoly fast-path for univariate integer polynomials (primary path).
34///
35/// # Arguments
36///
37/// * `dividend` - Polynomial to divide
38/// * `divisor` - Polynomial to divide by (must be non-zero)
39/// * `var` - Variable to treat as polynomial variable
40///
41/// # Errors
42///
43/// Returns `MathError::DivisionByZero` if divisor is zero.
44/// Returns `MathError::NotImplemented` for complex symbolic polynomial division.
45///
46/// # Examples
47///
48/// ```rust
49/// use mathhook_core::{expr, symbol};
50/// use mathhook_core::algebra::polynomial_division::polynomial_div;
51///
52/// let x = symbol!(x);
53/// // (x^2 + 3x + 2) / (x + 1) = (x + 2) with remainder 0
54/// let dividend = expr!((x^2) + (3*x) + 2);
55/// let divisor = expr!(x + 1);
56/// let (quot, rem) = polynomial_div(&dividend, &divisor, &x)?;
57/// # Ok::<(), mathhook_core::error::MathError>(())
58/// ```
59///
60/// # Returns
61///
62/// Returns `(quotient, remainder)` tuple where both are expressions
63pub fn polynomial_div(
64    dividend: &Expression,
65    divisor: &Expression,
66    var: &Symbol,
67) -> Result<(Expression, Expression), MathError> {
68    if divisor.is_zero() {
69        return Err(MathError::DivisionByZero);
70    }
71
72    if dividend.is_zero() {
73        return Ok((Expression::integer(0), Expression::integer(0)));
74    }
75
76    if dividend == divisor {
77        return Ok((Expression::integer(1), Expression::integer(0)));
78    }
79
80    // Fast path: divisor is constant
81    if let Some(divisor_const) = extract_constant(divisor) {
82        if divisor_const.is_zero() {
83            return Err(MathError::DivisionByZero);
84        }
85        let quotient = Expression::mul(vec![
86            dividend.clone(),
87            Expression::pow(divisor.clone(), Expression::integer(-1)),
88        ])
89        .simplify();
90        return Ok((quotient, Expression::integer(0)));
91    }
92
93    // IntPoly fast-path - PRIMARY PATH
94    let vars = dividend.find_variables();
95    if vars.len() == 1 {
96        let dividend_var = &vars[0];
97        if dividend_var == var {
98            let divisor_vars = divisor.find_variables();
99            if divisor_vars.len() == 1
100                && &divisor_vars[0] == var
101                && IntPoly::can_convert(dividend, var)
102                && IntPoly::can_convert(divisor, var)
103            {
104                if let (Some(p1), Some(p2)) = (
105                    IntPoly::try_from_expression(dividend, var),
106                    IntPoly::try_from_expression(divisor, var),
107                ) {
108                    if let Ok((q, r)) = p1.div_rem(&p2) {
109                        return Ok((q.to_expression(var), r.to_expression(var)));
110                    }
111                }
112            }
113        }
114    }
115
116    symbolic_polynomial_div(dividend, divisor, var)
117}
118
119/// Extract constant value from expression if it's a constant
120fn extract_constant(expr: &Expression) -> Option<Expression> {
121    match expr {
122        Expression::Number(_) => Some(expr.clone()),
123        _ => None,
124    }
125}
126
127/// Symbolic polynomial division using Expression operations
128///
129/// This is a MINIMAL fallback for cases that cannot use IntPoly.
130/// Used for rational coefficient polynomials and special cases.
131fn symbolic_polynomial_div(
132    dividend: &Expression,
133    divisor: &Expression,
134    var: &Symbol,
135) -> Result<(Expression, Expression), MathError> {
136    let dividend_degree = polynomial_degree_in_var(dividend, var);
137    let divisor_degree = polynomial_degree_in_var(divisor, var);
138
139    if dividend_degree < divisor_degree {
140        return Ok((Expression::integer(0), dividend.clone()));
141    }
142
143    if dividend_degree == divisor_degree {
144        let dividend_lc = polynomial_leading_coefficient(dividend, var);
145        let divisor_lc = polynomial_leading_coefficient(divisor, var);
146
147        let quotient_term = Expression::mul(vec![
148            dividend_lc,
149            Expression::pow(divisor_lc, Expression::integer(-1)),
150        ])
151        .simplify();
152
153        let product = Expression::mul(vec![quotient_term.clone(), divisor.clone()]).simplify();
154        let remainder = Expression::add(vec![
155            dividend.clone(),
156            Expression::mul(vec![Expression::integer(-1), product]),
157        ])
158        .simplify();
159
160        return Ok((quotient_term, remainder));
161    }
162
163    Err(MathError::NotImplemented {
164        feature: format!(
165            "complex symbolic polynomial division (degree {} ÷ degree {})",
166            dividend_degree, divisor_degree
167        ),
168    })
169}
170
171/// Get polynomial degree with respect to a specific variable
172fn polynomial_degree_in_var(expr: &Expression, var: &Symbol) -> i64 {
173    match expr {
174        Expression::Symbol(s) if s == var => 1,
175        Expression::Number(_) => 0,
176        Expression::Pow(base, exp) => {
177            if let (Expression::Symbol(s), Expression::Number(Number::Integer(e))) =
178                (base.as_ref(), exp.as_ref())
179            {
180                if s == var {
181                    return *e;
182                }
183            }
184            0
185        }
186        Expression::Add(terms) => {
187            let mut max_degree = 0i64;
188            for term in terms.iter() {
189                let deg = polynomial_degree_in_var(term, var);
190                max_degree = max_degree.max(deg);
191            }
192            max_degree
193        }
194        Expression::Mul(factors) => {
195            let mut total_degree = 0i64;
196            for factor in factors.iter() {
197                total_degree += polynomial_degree_in_var(factor, var);
198            }
199            total_degree
200        }
201        _ => 0,
202    }
203}
204
205/// Get leading coefficient of polynomial
206fn polynomial_leading_coefficient(expr: &Expression, var: &Symbol) -> Expression {
207    let degree = polynomial_degree_in_var(expr, var);
208
209    match expr {
210        Expression::Number(_n) => expr.clone(),
211        Expression::Symbol(s) if s == var => Expression::integer(1),
212        Expression::Pow(base, exp) => {
213            if let (Expression::Symbol(s), Expression::Number(Number::Integer(_))) =
214                (base.as_ref(), exp.as_ref())
215            {
216                if s == var {
217                    return Expression::integer(1);
218                }
219            }
220            expr.clone()
221        }
222        Expression::Mul(factors) => {
223            let mut coeff = Expression::integer(1);
224            for factor in factors.iter() {
225                if polynomial_degree_in_var(factor, var) == 0 {
226                    coeff = Expression::mul(vec![coeff, factor.clone()]);
227                }
228            }
229            coeff
230        }
231        Expression::Add(terms) => {
232            for term in terms.iter() {
233                if polynomial_degree_in_var(term, var) == degree {
234                    return polynomial_leading_coefficient(term, var);
235                }
236            }
237            Expression::integer(0)
238        }
239        _ => Expression::integer(1),
240    }
241}
242
243/// Polynomial quotient (division without remainder)
244///
245/// Returns only the quotient part of polynomial division
246///
247/// # Arguments
248///
249/// * `dividend` - Polynomial to divide
250/// * `divisor` - Polynomial to divide by
251/// * `var` - Variable to treat as polynomial variable
252///
253/// # Errors
254///
255/// Returns `MathError::DivisionByZero` if divisor is zero.
256///
257/// # Examples
258///
259/// ```rust
260/// use mathhook_core::{expr, symbol};
261/// use mathhook_core::algebra::polynomial_division::polynomial_quo;
262///
263/// let x = symbol!(x);
264/// let dividend = expr!((x^2) + (3*x) + 2);
265/// let divisor = expr!(x + 1);
266/// let quot = polynomial_quo(&dividend, &divisor, &x)?;
267/// # Ok::<(), mathhook_core::error::MathError>(())
268/// ```
269pub fn polynomial_quo(
270    dividend: &Expression,
271    divisor: &Expression,
272    var: &Symbol,
273) -> Result<Expression, MathError> {
274    polynomial_div(dividend, divisor, var).map(|(q, _)| q)
275}
276
277/// Polynomial remainder
278///
279/// Returns only the remainder part of polynomial division
280///
281/// # Arguments
282///
283/// * `dividend` - Polynomial to divide
284/// * `divisor` - Polynomial to divide by
285/// * `var` - Variable to treat as polynomial variable
286///
287/// # Errors
288///
289/// Returns `MathError::DivisionByZero` if divisor is zero.
290///
291/// # Examples
292///
293/// ```rust
294/// use mathhook_core::{expr, symbol};
295/// use mathhook_core::algebra::polynomial_division::polynomial_rem;
296///
297/// let x = symbol!(x);
298/// let dividend = expr!((x^2) + 1);
299/// let divisor = expr!(x - 1);
300/// let rem = polynomial_rem(&dividend, &divisor, &x)?;
301/// # Ok::<(), mathhook_core::error::MathError>(())
302/// ```
303pub fn polynomial_rem(
304    dividend: &Expression,
305    divisor: &Expression,
306    var: &Symbol,
307) -> Result<Expression, MathError> {
308    // IntPoly fast-path - dedicated remainder computation
309    let vars = dividend.find_variables();
310    if vars.len() == 1 {
311        let dividend_var = &vars[0];
312        if dividend_var == var {
313            let divisor_vars = divisor.find_variables();
314            if divisor_vars.len() == 1
315                && &divisor_vars[0] == var
316                && IntPoly::can_convert(dividend, var)
317                && IntPoly::can_convert(divisor, var)
318            {
319                if let (Some(p1), Some(p2)) = (
320                    IntPoly::try_from_expression(dividend, var),
321                    IntPoly::try_from_expression(divisor, var),
322                ) {
323                    if let Ok((_, r)) = p1.div_rem(&p2) {
324                        return Ok(r.to_expression(var));
325                    }
326                }
327            }
328        }
329    }
330
331    polynomial_div(dividend, divisor, var).map(|(_, r)| r)
332}
333
334#[cfg(test)]
335mod tests {
336    use super::*;
337    use crate::{expr, symbol};
338
339    #[test]
340    fn test_polynomial_div_exact() {
341        let x = symbol!(x);
342
343        let dividend = expr!((x ^ 2) - 1);
344        let divisor = expr!(x - 1);
345        let (_quot, rem) = polynomial_div(&dividend, &divisor, &x).unwrap();
346
347        assert!(rem.is_zero(), "Expected zero remainder");
348    }
349
350    #[test]
351    fn test_polynomial_div_with_remainder() {
352        let x = symbol!(x);
353
354        let dividend = expr!((x ^ 2) + 1);
355        let divisor = expr!(x - 1);
356        let (_quot, rem) = polynomial_div(&dividend, &divisor, &x).unwrap();
357
358        assert!(!rem.is_zero(), "Expected non-zero remainder");
359    }
360
361    #[test]
362    fn test_polynomial_div_by_constant() {
363        let x = symbol!(x);
364
365        let dividend = Expression::add(vec![
366            Expression::pow(Expression::symbol(x.clone()), Expression::integer(2)),
367            Expression::mul(vec![Expression::integer(2), Expression::symbol(x.clone())]),
368            Expression::integer(1),
369        ]);
370        let divisor = Expression::integer(2);
371        let (_quot, rem) = polynomial_div(&dividend, &divisor, &x).unwrap();
372
373        assert!(rem.is_zero(), "Expected zero remainder");
374    }
375
376    #[test]
377    fn test_polynomial_div_identical() {
378        let x = symbol!(x);
379
380        let dividend = expr!(x + 1);
381        let divisor = expr!(x + 1);
382        let (quot, rem) = polynomial_div(&dividend, &divisor, &x).unwrap();
383
384        assert_eq!(quot, Expression::integer(1));
385        assert!(rem.is_zero());
386    }
387
388    #[test]
389    fn test_polynomial_quo() {
390        let x = symbol!(x);
391
392        let dividend = expr!((x ^ 2) - 1);
393        let divisor = expr!(x - 1);
394        let quot = polynomial_quo(&dividend, &divisor, &x).unwrap();
395
396        assert!(!quot.is_zero());
397    }
398
399    #[test]
400    fn test_polynomial_rem() {
401        let x = symbol!(x);
402
403        let dividend = expr!((x ^ 2) + 1);
404        let divisor = expr!(x - 1);
405        let rem = polynomial_rem(&dividend, &divisor, &x).unwrap();
406
407        assert!(!rem.is_zero());
408    }
409
410    #[test]
411    fn test_intpoly_fastpath() {
412        let x = symbol!(x);
413
414        let dividend = expr!((x ^ 3) + (2 * (x ^ 2)) + (3 * x) + 4);
415        let divisor = expr!((x ^ 2) + 1);
416        let (quot, rem) = polynomial_div(&dividend, &divisor, &x).unwrap();
417
418        println!("Quotient: {}, Remainder: {}", quot, rem);
419        assert_ne!(quot, Expression::undefined());
420    }
421
422    #[test]
423    fn test_polynomial_div_by_zero() {
424        let x = symbol!(x);
425        let dividend = expr!(x ^ 2);
426        let divisor = Expression::integer(0);
427
428        let result = polynomial_div(&dividend, &divisor, &x);
429        assert!(matches!(result, Err(MathError::DivisionByZero)));
430    }
431}
mathhook_core/algebra/polynomial_division.rs

mathhook_core/algebra/
polynomial_division.rs
