feanor_math/algorithms/
poly_div.rs

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

feanor_math/algorithms/
poly_div.rs
