feanor_math/algorithms/
poly_div.rs

1use crate::divisibility::*;
2use crate::computation::*;
3use crate::pid::*;
4use crate::ring::*;
5use crate::homomorphism::*;
6use crate::rings::finite::FiniteRing;
7use crate::rings::poly::*;
8
9use std::cmp::max;
10
11///
12/// Computes the polynomial division of `lhs` by `rhs`, i.e. `lhs = q * rhs + r` with
13/// `deg(r) < deg(rhs)`.
14/// 
15/// Note that this function does not compute the proper polynomial division if the leading
16/// coefficient of `rhs` is a zero-divisor in the ring. See [`poly_div_rem_finite_reduced()`]
17/// for details.
18/// 
19/// This requires a function `left_div_lc` that computes the division of an element of the 
20/// base ring by the leading coefficient of `rhs`. If the base ring is a field, this can
21/// just be standard division. In other cases, this depends on the exact situation you are
22/// in - e.g. `rhs` might be monic or in in a specific context, it might be guaranteed that the 
23/// division always works. If this is not the case, look also at [`poly_div_rem_domain()`], which
24/// implicitly performs the polynomial division over the field of fractions.
25/// 
26#[stability::unstable(feature = "enable")]
27pub fn poly_div_rem<P, F, E>(poly_ring: P, mut lhs: El<P>, rhs: &El<P>, mut left_div_lc: F) -> Result<(El<P>, El<P>), E>
28    where P: RingStore,
29        P::Type: PolyRing,
30        F: FnMut(&El<<P::Type as RingExtension>::BaseRing>) -> Result<El<<P::Type as RingExtension>::BaseRing>, E>
31{
32    assert!(poly_ring.degree(rhs).is_some());
33
34    let rhs_deg = poly_ring.degree(rhs).unwrap();
35    if poly_ring.degree(&lhs).is_none() {
36        return Ok((poly_ring.zero(), lhs));
37    }
38    let lhs_deg = poly_ring.degree(&lhs).unwrap();
39    if lhs_deg < rhs_deg {
40        return Ok((poly_ring.zero(), lhs));
41    }
42    let result = poly_ring.try_from_terms(
43        (0..(lhs_deg + 1 - rhs_deg)).rev().map(|i| {
44            let quo = left_div_lc(poly_ring.coefficient_at(&lhs, i +  rhs_deg))?;
45            let neg_quo = poly_ring.base_ring().negate(quo);
46            if !poly_ring.base_ring().is_zero(&neg_quo) {
47                poly_ring.get_ring().add_assign_from_terms(
48                    &mut lhs, 
49                    poly_ring.terms(rhs).map(|(c, j)| 
50                        (poly_ring.base_ring().mul_ref(&neg_quo, c), i + j)
51                    )
52                );
53            }
54            Ok((poly_ring.base_ring().negate(neg_quo), i))
55        })
56    )?;
57    return Ok((result, lhs));
58}
59
60const FAST_POLY_DIV_THRESHOLD: usize = 32;
61
62///
63/// Computes the polynomial division of `lhs` by `rhs`, i.e. `lhs = q * rhs + r` with
64/// `deg(r) < deg(rhs)`, i.e. is functionally equivalent to [`poly_div_rem()`].
65/// 
66/// As opposed to [`poly_div_rem()`], this function uses a fast polynomial division algorithm,
67/// which is faster for large inputs.
68/// 
69pub fn fast_poly_div_rem<P, F, E, Controller>(poly_ring: P, f: El<P>, g: &El<P>, mut left_div_lc: F, controller: Controller)-> Result<(El<P>, El<P>), E>
70    where P: RingStore + Copy,
71        P::Type: PolyRing,
72        F: FnMut(&El<<P::Type as RingExtension>::BaseRing>) -> Result<El<<P::Type as RingExtension>::BaseRing>, E>,
73        Controller: ComputationController
74{
75    fn fast_poly_div_impl<P, F, E, Controller>(poly_ring: P, f: El<P>, g: &El<P>, left_div_lc: &mut F, controller: Controller)-> Result<(El<P>, El<P>), E>
76        where P: RingStore + Copy,
77            P::Type: PolyRing,
78            F: FnMut(&El<<P::Type as RingExtension>::BaseRing>) -> Result<El<<P::Type as RingExtension>::BaseRing>, E>,
79            Controller: ComputationController
80    {
81        let deg_g = poly_ring.degree(g).unwrap();
82        if poly_ring.degree(&f).is_none() || poly_ring.degree(&f).unwrap() < deg_g {
83            return Ok((poly_ring.zero(), f));
84        } 
85        let deg_f = poly_ring.degree(&f).unwrap();
86        if deg_g < FAST_POLY_DIV_THRESHOLD || (deg_f - deg_g) < FAST_POLY_DIV_THRESHOLD {
87            return poly_div_rem(poly_ring, f, g, left_div_lc);
88        }
89
90        let (split_degree_f, split_degree_g) = if deg_f >= 3 * deg_g {
91            (deg_f / 3, 0)
92        } else if 2 * (deg_f / 3) < deg_g {
93            (deg_g / 2, deg_g / 2)
94        } else {
95            (deg_f / 3, deg_g - deg_f / 3)
96        };
97        assert!(split_degree_f >= split_degree_g);
98        assert!(split_degree_f <= deg_f);
99        assert!(split_degree_g <= deg_g);
100
101        let f_upper = poly_ring.from_terms(poly_ring.terms(&f).filter(|(_, i)| *i >= split_degree_f).map(|(c, i)| (poly_ring.base_ring().clone_el(c), i - split_degree_f)));
102        let mut f_lower = f;
103        poly_ring.truncate_monomials(&mut f_lower, split_degree_f);
104        let g_upper = poly_ring.from_terms(poly_ring.terms(&g).filter(|(_, i)| *i >= split_degree_g).map(|(c, i)| (poly_ring.base_ring().clone_el(c), i - split_degree_g)));
105        let mut g_lower = poly_ring.clone_el(g);
106        poly_ring.truncate_monomials(&mut g_lower, split_degree_g);
107
108        let (q_upper, r) = fast_poly_div_impl(poly_ring, poly_ring.clone_el(&f_upper), &g_upper, &mut *left_div_lc, controller.clone())?;
109        debug_assert!(poly_ring.degree(&q_upper).is_none() || poly_ring.degree(&q_upper).unwrap() <= deg_f + split_degree_g - split_degree_f - deg_g);
110        debug_assert!(poly_ring.degree(&r).is_none() || poly_ring.degree(&r).unwrap() <= deg_g - split_degree_g - 1);
111
112        poly_ring.get_ring().add_assign_from_terms(&mut f_lower, poly_ring.terms(&r).map(|(c, i)| (poly_ring.base_ring().clone_el(c), i + split_degree_f)));
113        debug_assert!(poly_ring.degree(&f_lower).is_none() || poly_ring.degree(&f_lower).unwrap() <= deg_g + split_degree_f - split_degree_g);
114        poly_ring.mul_assign_ref(&mut g_lower, &q_upper);
115        poly_ring.get_ring().add_assign_from_terms(&mut f_lower, poly_ring.terms(&g_lower).map(|(c, i)| (poly_ring.base_ring().negate(poly_ring.base_ring().clone_el(c)), i + split_degree_f - split_degree_g)));
116        debug_assert!(poly_ring.degree(&f_lower).is_none() || poly_ring.degree(&f_lower).unwrap() <= max(deg_f + split_degree_g - deg_g, deg_g + split_degree_f - split_degree_g));
117
118        let (mut q_lower, r) = fast_poly_div_impl(poly_ring, poly_ring.clone_el(&f_lower), g, &mut *left_div_lc, controller)?;
119
120        poly_ring.get_ring().add_assign_from_terms(&mut q_lower, poly_ring.terms(&q_upper).map(|(c, i)| (poly_ring.base_ring().clone_el(c), i + split_degree_f - split_degree_g)));
121        return Ok((q_lower, r));
122    }
123
124    assert!(!poly_ring.is_zero(g));
125    if poly_ring.is_zero(&f) {
126        return Ok((poly_ring.zero(), f));
127    }
128    controller.run_computation(format_args!("fast_poly_div(ldeg={}, rdeg={})", poly_ring.degree(&f).unwrap(), poly_ring.degree(g).unwrap()), |controller| 
129        fast_poly_div_impl(poly_ring, f, g, &mut left_div_lc, controller)
130    )
131}
132
133///
134/// Computes `(q, r, a)` such that `a * lhs = q * rhs + r` and `deg(r) < deg(rhs)`.
135/// The chosen factor `a` is in the base ring and is the smallest possible w.r.t.
136/// divisibility.
137/// 
138#[stability::unstable(feature = "enable")]
139pub fn poly_div_rem_domain<P>(ring: P, mut lhs: El<P>, rhs: &El<P>) -> (El<P>, El<P>, El<<P::Type as RingExtension>::BaseRing>)
140    where P: RingStore,
141        P::Type: PolyRing,
142        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: Domain + PrincipalIdealRing
143{
144    assert!(!ring.is_zero(rhs));
145    let d = ring.degree(rhs).unwrap();
146    let base_ring = ring.base_ring();
147    let rhs_lc = ring.lc(rhs).unwrap();
148
149    let mut current_scale = base_ring.one();
150    let mut terms = Vec::new();
151    while let Some(lhs_deg) = ring.degree(&lhs) {
152        if lhs_deg < d {
153            break;
154        }
155        let lhs_lc = base_ring.clone_el(ring.lc(&lhs).unwrap());
156        let gcd = base_ring.ideal_gen(&lhs_lc, &rhs_lc);
157        let additional_scale = base_ring.checked_div(&rhs_lc, &gcd).unwrap();
158
159        base_ring.mul_assign_ref(&mut current_scale, &additional_scale);
160        terms.iter_mut().for_each(|(c, _)| base_ring.mul_assign_ref(c, &additional_scale));
161        ring.inclusion().mul_assign_map(&mut lhs, additional_scale);
162
163        let factor = base_ring.checked_div(ring.lc(&lhs).unwrap(), rhs_lc).unwrap();
164        ring.get_ring().add_assign_from_terms(&mut lhs,
165            ring.terms(rhs).map(|(c, i)| (base_ring.negate(base_ring.mul_ref(c, &factor)), i + lhs_deg - d))
166        );
167        terms.push((factor, lhs_deg - d));
168    }
169    return (ring.from_terms(terms.into_iter()), lhs, current_scale);
170}
171
172///
173/// Possible errors that might be returned by [`poly_div_rem_finite_reduced()`].
174/// 
175#[stability::unstable(feature = "enable")]
176pub enum PolyDivRemReducedError<R>
177    where R: ?Sized + RingBase
178{
179    NotReduced(R::Element),
180    NotDivisibleByContent(R::Element)
181}
182
183///
184/// Given polynomials `f, g` over a finite and reduced ring `R`, tries to compute the
185/// polynomial division of `f` by `g`, i.e. values `q, r in R[X]` with `f = q g + r`
186/// and `deg(r) < deg(g)`.
187/// 
188/// As opposed to the case when `R` is a field, it is possible that this fails if `cont(g)`
189/// does not divide `cont(f)`. Hence, the possible results of this function are
190///  - success, i.e. `q, r` such that `f = q g + r` and `deg(r) < deg(g)`
191///  - [`PolyDivRemReducedError::NotDivisibleByContent`] with the content `c in R`, if `c` does not divide `cont(f)`
192///  - [`PolyDivRemReducedError::NotReduced`] with a nilpotent element `x != 0` if `x^2 = 0`, which means that 
193///    the ring is not reduced
194/// 
195/// Since a finite reduced ring is always a product of finite fields, one could (in theory)
196/// compute the polynomial division in each field and reconstruct the result. However, this
197/// function finds the result without computing the ring factors, which can be very difficult
198/// in some situations (e.g. it might require factoring the modulus `n`). Note however that
199/// this function does not necessarily give the same result as the div-in-field-and-reconstruct
200/// approach. See "lack of uniqueness" later.
201/// 
202/// Note that sometimes, even if `R` is not reduced, or `cont(g)` does not divide `cont(f)`, there
203/// might still exist suitable `q, r`. In these cases, it is unspecified whether the function aborts
204/// or returns suitable `q, r`.
205/// 
206/// # Why reduced?
207/// 
208/// Polynomial division as above makes sense over finite reduced rings (since they are always
209/// a product of finite fields), but cannot be properly defined over unreduced rings anymore.
210/// The reason is that the divisors of some polynomial `f` don't need to have degree `<= deg(f)`,
211/// e.g.
212/// ```text
213///   1 = (5 X^2 + 5 X + 1) (-5 X^2 - 5 X + 1) mod 5^2
214/// ```
215/// 
216/// # Lack of Uniqueness
217/// 
218/// The result of this function does not have to be unique.
219/// Consider for example
220/// ```rust
221/// # use feanor_math::assert_el_eq;
222/// # use feanor_math::ring::*;
223/// # use feanor_math::rings::zn::zn_64::*;
224/// # use feanor_math::rings::poly::dense_poly::*;
225/// # use feanor_math::rings::poly::*;
226/// # use feanor_math::algorithms::poly_div::*;
227/// let Z6 = Zn::new(6);
228/// let Z6X = DensePolyRing::new(Z6, "X");
229/// let [g, q1, q2, r1, r2] = Z6X.with_wrapped_indeterminate(|X| [
230///     X.pow_ref(2) * 2 + 1,
231///     X.clone(),
232///     4 * X,
233///     X + 1,
234///     4 * X + 1
235/// ]);
236/// assert_el_eq!(&Z6X, Z6X.add_ref(&Z6X.mul_ref(&g, &q1), &r1), Z6X.add_ref(&Z6X.mul_ref(&g, &q2), &r2));
237/// ```
238/// In particular, `g | f` does not not imply that the remainder of the division of `f` by `g` is `0`.
239/// Clearly `g` must divide the remainder `r`, but if the ring is not a domain, `g` might divide polynomials
240/// of smaller degree. Hence, use [`poly_checked_div_finite_reduced()`] to check for divisibility.
241/// 
242#[stability::unstable(feature = "enable")]
243pub fn poly_div_rem_finite_reduced<P>(ring: P, mut lhs: El<P>, rhs: &El<P>) -> Result<(El<P>, El<P>), PolyDivRemReducedError<<<P::Type as RingExtension>::BaseRing as RingStore>::Type>>
244    where P: RingStore,
245        P::Type: PolyRing,
246        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: FiniteRing + PrincipalIdealRing
247{
248    if ring.is_zero(rhs) && ring.is_zero(&lhs) {
249        return Ok((ring.zero(), ring.zero()));
250    } else if ring.is_zero(rhs) {
251        return Err(PolyDivRemReducedError::NotDivisibleByContent(ring.base_ring().zero()));
252    }
253    let rhs_deg = ring.degree(rhs).unwrap();
254    let mut result = ring.zero();
255    let zero = ring.base_ring().zero();
256    while ring.degree(&lhs).is_some() && ring.degree(&lhs).unwrap() >= rhs_deg {
257        let lhs_deg = ring.degree(&lhs).unwrap();
258        let lcf = ring.lc(&lhs).unwrap();
259        let mut h = ring.zero();
260        let mut annihilator = ring.base_ring().one();
261        let mut i: i64 = rhs_deg.try_into().unwrap();
262        let mut d = ring.base_ring().zero();
263        while ring.base_ring().checked_div(lcf, &d).is_none() {
264            debug_assert!(ring.base_ring().eq_el(ring.lc(&ring.mul_ref(&h, rhs)).unwrap_or(&zero), &d));
265            if i == -1 {
266                return Err(PolyDivRemReducedError::NotDivisibleByContent(d));
267            }
268            let (s, t, new_d) = ring.base_ring().extended_ideal_gen(&d, &ring.base_ring().mul_ref(&annihilator, ring.coefficient_at(&rhs, i as usize)));
269            ring.inclusion().mul_assign_map(&mut h, s);
270            ring.add_assign(&mut h, ring.from_terms([(ring.base_ring().mul_ref(&annihilator, &t), lhs_deg - i as usize)]));
271            annihilator = ring.base_ring().annihilator(&new_d);
272            d = new_d;
273            i = i - 1;
274            if !ring.base_ring().is_unit(&ring.base_ring().ideal_gen(&annihilator, &d)) {
275                let nilpotent = ring.base_ring().annihilator(&ring.base_ring().ideal_gen(&annihilator, &d));
276                debug_assert!(!ring.base_ring().is_zero(&nilpotent));
277                debug_assert!(ring.base_ring().is_zero(&ring.base_ring().mul_ref(&nilpotent, &nilpotent)));
278                return Err(PolyDivRemReducedError::NotReduced(nilpotent));
279            }
280            debug_assert!(ring.base_ring().eq_el(ring.lc(&ring.mul_ref(&h, rhs)).unwrap_or(&zero), &d));
281        }
282        let scale = ring.base_ring().checked_div(lcf, &d).unwrap();
283        ring.inclusion().mul_assign_map(&mut h, scale);
284        ring.sub_assign(&mut lhs, ring.mul_ref(&h, rhs));
285        ring.add_assign(&mut result, h);
286        debug_assert!(ring.degree(&lhs).map(|d| d + 1).unwrap_or(0) <= lhs_deg);
287    }
288    return Ok((result, lhs));
289}
290
291///
292/// Checks whether `rhs | lhs` and returns a quotient if one exists, assuming that
293/// the base ring is reduced. If it is not, the function may fail with a nilpotent
294/// element of the base ring.
295/// 
296#[stability::unstable(feature = "enable")]
297pub fn poly_checked_div_finite_reduced<P>(ring: P, mut lhs: El<P>, mut rhs: El<P>) -> Result<Option<El<P>>, El<<P::Type as RingExtension>::BaseRing>>
298    where P: RingStore + Copy,
299        P::Type: PolyRing,
300        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: FiniteRing + PrincipalIdealRing
301{
302    let mut result = ring.zero();
303    while !ring.is_zero(&lhs) {
304        match poly_div_rem_finite_reduced(ring, lhs, &rhs) {
305            Ok((q, r)) => {
306                ring.add_assign(&mut result, q);
307                lhs = r;
308            },
309            Err(PolyDivRemReducedError::NotReduced(nilpotent)) => {
310                return Err(nilpotent);
311            },
312            Err(PolyDivRemReducedError::NotDivisibleByContent(_)) => {
313                return Ok(None);
314            }
315        }
316        let annihilate_lc = ring.base_ring().annihilator(ring.lc(&rhs).unwrap());
317        ring.inclusion().mul_assign_map(&mut rhs, annihilate_lc);
318    }
319    return Ok(Some(result));
320}
321
322#[cfg(test)]
323use crate::rings::zn::zn_64::*;
324#[cfg(test)]
325use crate::rings::zn::*;
326#[cfg(test)]
327use crate::integer::*;
328#[cfg(test)]
329use crate::primitive_int::*;
330#[cfg(test)]
331use dense_poly::DensePolyRing;
332
333#[test]
334fn test_poly_div_rem_finite_reduced() {
335    let base_ring = Zn::new(5 * 7 * 11);
336    let ring = DensePolyRing::new(base_ring, "X");
337
338    let [f, g, _q, _r] = ring.with_wrapped_indeterminate(|X| [
339        X.pow_ref(2),
340        X.pow_ref(2) * 5 + X * 7 + 11,
341        X * (-77) + 108,
342        91 * X - 33
343    ]);
344    let (q, r) = poly_div_rem_finite_reduced(&ring, ring.clone_el(&f), &g).ok().unwrap();
345    assert_eq!(1, ring.degree(&r).unwrap());
346    assert_el_eq!(&ring, &f, ring.add(ring.mul(q, g), r));
347
348    let [f, g] = ring.with_wrapped_indeterminate(|X| [
349        5 * X.pow_ref(2),
350        X * 5 * 11 + X * 7 * 11,
351    ]);
352    if let Err(PolyDivRemReducedError::NotDivisibleByContent(content)) = poly_div_rem_finite_reduced(&ring, f, &g) {
353        assert!(base_ring.checked_div(&content, &base_ring.int_hom().map(11)).is_some());
354        assert!(base_ring.checked_div(&base_ring.int_hom().map(11), &content).is_some());
355    } else {
356        assert!(false);
357    }
358
359    let base_ring = Zn::new(5 * 5 * 11);
360    let ring = DensePolyRing::new(base_ring, "X");
361
362    // note that `11 = (1 - 5 X - 5 X^2)(1 + 55 X + 55 X^2) mod 5^2 * 11`
363    let [g] = ring.with_wrapped_indeterminate(|X| [
364        1 - 5 * X - 5 * X.pow_ref(2)
365    ]);
366    let f = ring.from_terms([(base_ring.int_hom().map(11), 2)]);
367    if let Err(PolyDivRemReducedError::NotReduced(nilpotent)) = poly_div_rem_finite_reduced(&ring, ring.clone_el(&f), &g) {
368        assert!(!base_ring.is_zero(&nilpotent));
369        assert!(base_ring.is_zero(&base_ring.pow(nilpotent, 2)));
370    } else {
371        assert!(false);
372    }
373}
374
375#[test]
376fn test_poly_div_rem_finite_reduced_nonmonic() {
377    let base_ring = zn_big::Zn::new(BigIntRing::RING, int_cast(8589934594, BigIntRing::RING, StaticRing::<i64>::RING));
378    let poly_ring = DensePolyRing::new(&base_ring, "X");
379    let [f, g] = poly_ring.with_wrapped_indeterminate(|X| [
380        1431655767 + -1431655765 * X,
381        -2 + X + X.pow_ref(2)
382    ]);
383    let (q, r) = poly_div_rem_finite_reduced(&poly_ring, poly_ring.clone_el(&g), &f).ok().unwrap();
384    assert_eq!(0, poly_ring.degree(&r).unwrap_or(0));
385    assert_el_eq!(&poly_ring, &g, poly_ring.add(poly_ring.mul_ref(&q, &f), r));
386}
387
388#[test]
389fn test_fast_poly_div() {
390    let ZZ = BigIntRing::RING;
391    let ZZX = DensePolyRing::new(ZZ, "X");
392    let [f, g] = ZZX.with_wrapped_indeterminate(|X| [X.pow_ref(80) - 1, X.pow_ref(40) - 2 * X.pow_ref(33) + X.pow_ref(21) - X + 10]);
393    assert_el_eq!(&ZZX, ZZX.div_rem_monic(ZZX.clone_el(&f), &g).0, fast_poly_div_rem(&ZZX, ZZX.clone_el(&f), &g, |c| Ok(ZZ.clone_el(c)), LOG_PROGRESS).unwrap_or_else(no_error).0);
394}
feanor_math/algorithms/poly_div.rs

feanor_math/algorithms/
poly_div.rs
