mathhook_core/core/polynomial/poly/
division.rs

1use super::Poly;
2use crate::core::polynomial::traits::EuclideanDomain;
3use crate::error::MathError;
4
5/// Integer GCD using Euclidean algorithm
6#[inline(always)]
7fn gcd_i64(mut a: i64, mut b: i64) -> i64 {
8    while b != 0 {
9        let t = b;
10        b = a % b;
11        a = t;
12    }
13    a.abs()
14}
15
16impl<T: EuclideanDomain> Poly<T> {
17    /// Polynomial division with remainder
18    ///
19    /// Returns (quotient, remainder) such that self = quotient * divisor + remainder
20    /// and degree(remainder) < degree(divisor)
21    ///
22    /// # Errors
23    /// Returns `MathError::DivisionByZero` if divisor is zero
24    #[inline(always)]
25    pub fn div_rem(&self, divisor: &Self) -> Result<(Self, Self), MathError> {
26        if divisor.is_zero() {
27            return Err(MathError::DivisionByZero);
28        }
29
30        if self.is_zero() {
31            return Ok((Self::zero(), Self::zero()));
32        }
33
34        let self_deg = match self.degree() {
35            Some(d) => d,
36            None => return Ok((Self::zero(), Self::zero())),
37        };
38
39        let divisor_deg = match divisor.degree() {
40            Some(d) => d,
41            None => return Err(MathError::DivisionByZero),
42        };
43
44        if self_deg < divisor_deg {
45            return Ok((Self::zero(), self.clone()));
46        }
47
48        let divisor_lc = divisor.leading_coeff();
49        let mut remainder = self.coeffs.clone();
50        let mut quotient = Vec::with_capacity(self_deg - divisor_deg + 1);
51        for _ in 0..=self_deg - divisor_deg {
52            quotient.push(T::zero());
53        }
54
55        for i in (0..=self_deg - divisor_deg).rev() {
56            let rem_idx = i + divisor_deg;
57            let rem_coeff = remainder[rem_idx].clone();
58
59            let (q_coeff, r) = rem_coeff.div_rem(&divisor_lc);
60            if !r.is_zero() {
61                break;
62            }
63
64            quotient[i] = q_coeff.clone();
65
66            for (j, d_coeff) in divisor.coeffs.iter().enumerate() {
67                let sub_val = q_coeff.wrapping_mul(d_coeff);
68                remainder[i + j] = remainder[i + j].wrapping_sub(&sub_val);
69            }
70        }
71
72        Ok((Self::from_coeffs(quotient), Self::from_coeffs(remainder)))
73    }
74
75    /// Polynomial pseudo-division
76    ///
77    /// Returns (quotient, remainder, multiplier) such that:
78    /// multiplier * self = quotient * divisor + remainder
79    ///
80    /// This avoids fractions by multiplying by powers of the leading coefficient.
81    ///
82    /// # Errors
83    /// Returns `MathError::DivisionByZero` if divisor is zero
84    #[inline(always)]
85    pub fn pseudo_div_rem(&self, divisor: &Self) -> Result<(Self, Self, T), MathError> {
86        if divisor.is_zero() {
87            return Err(MathError::DivisionByZero);
88        }
89
90        if self.is_zero() {
91            return Ok((Self::zero(), Self::zero(), T::one()));
92        }
93
94        let self_deg = match self.degree() {
95            Some(d) => d,
96            None => return Ok((Self::zero(), Self::zero(), T::one())),
97        };
98
99        let divisor_deg = match divisor.degree() {
100            Some(d) => d,
101            None => return Err(MathError::DivisionByZero),
102        };
103
104        if self_deg < divisor_deg {
105            return Ok((Self::zero(), self.clone(), T::one()));
106        }
107
108        let divisor_lc = divisor.leading_coeff();
109        let mut remainder = self.coeffs.clone();
110        let mut quotient = Vec::with_capacity(self_deg - divisor_deg + 1);
111        for _ in 0..=self_deg - divisor_deg {
112            quotient.push(T::zero());
113        }
114        let mut multiplier = T::one();
115
116        for i in (0..=self_deg - divisor_deg).rev() {
117            let rem_idx = i + divisor_deg;
118
119            for c in remainder.iter_mut() {
120                *c = c.wrapping_mul(&divisor_lc);
121            }
122            for c in quotient.iter_mut() {
123                *c = c.wrapping_mul(&divisor_lc);
124            }
125            multiplier = multiplier.wrapping_mul(&divisor_lc);
126
127            let (q_coeff, _) = remainder[rem_idx].div_rem(&divisor_lc);
128            quotient[i] = q_coeff.clone();
129
130            for (j, d_coeff) in divisor.coeffs.iter().enumerate() {
131                let sub_val = q_coeff.wrapping_mul(d_coeff);
132                remainder[i + j] = remainder[i + j].wrapping_sub(&sub_val);
133            }
134        }
135
136        Ok((
137            Self::from_coeffs(quotient),
138            Self::from_coeffs(remainder),
139            multiplier,
140        ))
141    }
142
143    /// GCD using Euclidean algorithm
144    ///
145    /// Returns the greatest common divisor of two polynomials.
146    /// The result is made primitive (content = 1).
147    ///
148    /// # Errors
149    /// Returns `MathError::DivisionByZero` if division by zero occurs during computation
150    #[inline(always)]
151    pub fn gcd(&self, other: &Self) -> Result<Self, MathError> {
152        if self.is_zero() {
153            return other.primitive_part();
154        }
155        if other.is_zero() {
156            return self.primitive_part();
157        }
158
159        let mut a = self.primitive_part()?;
160        let mut b = other.primitive_part()?;
161
162        while !b.is_zero() {
163            let (_, rem, _) = a.pseudo_div_rem(&b)?;
164            a = b;
165            b = if rem.is_zero() {
166                rem
167            } else {
168                rem.primitive_part()?
169            };
170        }
171
172        Ok(a)
173    }
174
175    /// Content: GCD of all coefficients
176    #[inline(always)]
177    pub fn content(&self) -> T {
178        if self.coeffs.is_empty() {
179            return T::zero();
180        }
181
182        let mut g = self.coeffs[0].abs();
183        for c in &self.coeffs[1..] {
184            g = g.gcd(&c.abs());
185            if g.is_one() {
186                break;
187            }
188        }
189        g
190    }
191
192    /// Primitive part: polynomial / content
193    ///
194    /// # Errors
195    /// Returns `MathError::DivisionByZero` if content is zero
196    #[inline(always)]
197    pub fn primitive_part(&self) -> Result<Self, MathError> {
198        let c = self.content();
199        if c.is_zero() {
200            return Err(MathError::DivisionByZero);
201        }
202        if c.is_one() {
203            return Ok(self.clone());
204        }
205
206        let coeffs: Vec<T> = self
207            .coeffs
208            .iter()
209            .map(|x| {
210                let (quot, _) = x.div_rem(&c);
211                quot
212            })
213            .collect();
214        Ok(Self { coeffs })
215    }
216}
217
218impl super::IntPoly {
219    /// Check if polynomial is monic (leading coefficient = 1)
220    #[inline(always)]
221    pub fn is_monic(&self) -> bool {
222        self.leading_coeff() == 1
223    }
224
225    /// Make polynomial monic by dividing by leading coefficient
226    /// Returns None if leading coefficient doesn't divide all coefficients
227    #[inline(always)]
228    pub fn try_make_monic(&self) -> Option<Self> {
229        let lc = self.leading_coeff();
230        if lc == 0 {
231            return None;
232        }
233        if lc == 1 {
234            return Some(self.clone());
235        }
236
237        for &c in &self.coeffs {
238            if c % lc != 0 {
239                return None;
240            }
241        }
242
243        let coeffs: Vec<i64> = self.coeffs.iter().map(|&c| c / lc).collect();
244        Some(Self { coeffs })
245    }
246
247    /// Normalize IntPoly (make leading coefficient positive)
248    #[inline(always)]
249    pub fn normalize(&self) -> Self {
250        if self.leading_coeff() < 0 {
251            self.negate()
252        } else {
253            self.clone()
254        }
255    }
256
257    /// Content for IntPoly (optimized)
258    #[inline(always)]
259    pub fn content_i64(&self) -> i64 {
260        if self.coeffs.is_empty() {
261            return 0;
262        }
263
264        let mut g = self.coeffs[0].abs();
265        for &c in &self.coeffs[1..] {
266            g = gcd_i64(g, c.abs());
267            if g == 1 {
268                break;
269            }
270        }
271        g
272    }
273
274    /// Primitive part for IntPoly (optimized)
275    ///
276    /// # Errors
277    /// Returns `MathError::DivisionByZero` if content is zero
278    #[inline(always)]
279    pub fn primitive_part_i64(&self) -> Result<Self, MathError> {
280        let c = self.content_i64();
281        if c == 0 {
282            return Err(MathError::DivisionByZero);
283        }
284        if c == 1 {
285            return Ok(self.clone());
286        }
287
288        let coeffs: Vec<i64> = self.coeffs.iter().map(|&x| x / c).collect();
289        Ok(Self { coeffs })
290    }
291
292    /// GCD for IntPoly (optimized)
293    ///
294    /// # Errors
295    /// Returns `MathError::DivisionByZero` if division by zero occurs during computation
296    #[inline(always)]
297    pub fn gcd_i64(&self, other: &Self) -> Result<Self, MathError> {
298        if self.is_zero() {
299            return Ok(other.primitive_part_i64()?.normalize());
300        }
301        if other.is_zero() {
302            return Ok(self.primitive_part_i64()?.normalize());
303        }
304
305        let mut a = self.primitive_part_i64()?;
306        let mut b = other.primitive_part_i64()?;
307
308        while !b.is_zero() {
309            let (_, rem, _) = a.pseudo_div_rem(&b)?;
310            a = b;
311            b = if rem.is_zero() {
312                rem
313            } else {
314                rem.primitive_part_i64()?
315            };
316        }
317
318        Ok(a.normalize())
319    }
320}
321
322#[cfg(test)]
323mod tests {
324    use crate::core::polynomial::poly::IntPoly;
325
326    #[test]
327    fn test_division_exact() {
328        let dividend = IntPoly::from_coeffs(vec![-1, 0, 1]);
329        let divisor = IntPoly::from_coeffs(vec![-1, 1]);
330        let (quot, rem) = dividend.div_rem(&divisor).unwrap();
331
332        assert_eq!(quot.coefficients(), &[1, 1]);
333        assert!(rem.is_zero());
334    }
335
336    #[test]
337    fn test_gcd() {
338        let p1 = IntPoly::from_coeffs(vec![0, 0, 6]);
339        let p2 = IntPoly::from_coeffs(vec![0, 9]);
340        let g = p1.gcd_i64(&p2).unwrap();
341
342        assert_eq!(g.degree(), Some(1));
343        assert_eq!(g.coeff(0), 0);
344        assert!(g.coeff(1) > 0);
345    }
346
347    #[test]
348    fn test_gcd_coprime() {
349        let p1 = IntPoly::from_coeffs(vec![1, 1]);
350        let p2 = IntPoly::from_coeffs(vec![2, 1]);
351        let g = p1.gcd_i64(&p2).unwrap();
352
353        assert!(g.is_constant());
354        assert_eq!(g.coefficients(), &[1]);
355    }
356
357    #[test]
358    fn test_content() {
359        let p = IntPoly::from_coeffs(vec![6, 12, 18]);
360        assert_eq!(p.content_i64(), 6);
361    }
362
363    #[test]
364    fn test_primitive_part() {
365        let p = IntPoly::from_coeffs(vec![6, 12, 18]);
366        let pp = p.primitive_part_i64().unwrap();
367        assert_eq!(pp.coefficients(), &[1, 2, 3]);
368    }
369
370    #[test]
371    fn test_division_by_zero() {
372        let dividend = IntPoly::from_coeffs(vec![1, 2, 3]);
373        let divisor = IntPoly::from_coeffs(vec![]);
374        assert!(dividend.div_rem(&divisor).is_err());
375    }
376}
mathhook_core/core/polynomial/poly/division.rs

mathhook_core/core/polynomial/poly/
division.rs
