feanor_math/algorithms/
poly_div.rs

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

feanor_math/algorithms/
poly_div.rs
