feanor_math/algorithms/poly_gcd/
mod.rs

1use finite::{fast_poly_eea, poly_power_decomposition_finite_field};
2use gcd::poly_gcd_local;
3use squarefree_part::poly_power_decomposition_local;
4
5use crate::computation::*;
6use crate::delegate::DelegateRing;
7use crate::divisibility::*;
8use crate::homomorphism::*;
9use crate::pid::*;
10use crate::reduce_lift::poly_factor_gcd::*;
11use crate::ring::*;
12use crate::rings::field::*;
13use crate::rings::finite::*;
14use crate::rings::poly::dense_poly::*;
15use crate::rings::poly::*;
16use crate::specialization::FiniteRingOperation;
17
18/// Contains an implementation of factoring polynomials over finite fields.
19pub mod finite;
20/// Contains algorithms for computing the gcd of polynomials.
21pub mod gcd;
22/// Contains an implementation of Hensel lifting, to lift a factorization modulo
23/// a maximal ideal to a factorization modulo a power of this ideal.
24pub mod hensel;
25/// Contains algorithms for computing power decompositions and the square-free
26/// part of polynomials.
27pub mod squarefree_part;
28
29const INCREASE_EXPONENT_PER_ATTEMPT_CONSTANT: f64 = 1.5;
30
31/// Trait for domain `R` for which there is an efficient way of computing the gcd
32/// of univariate polynomials over `TFrac(R)`, where `TFrac(R)` is the total ring
33/// of fractions.
34///
35/// However, computations in `TFrac(R)` are avoided by most implementations due to
36/// performance reasons, and both inputs and outputs are polynomials over `R`. Despite
37/// this, the gcd is the gcd over `TFrac(R)` and not `R` (the gcd over `R` is often not
38/// even defined, since `R` does not have to be UFD).
39///
40/// # Example
41/// ```rust
42/// # use feanor_math::assert_el_eq;
43/// # use feanor_math::ring::*;
44/// # use feanor_math::algorithms::poly_gcd::*;
45/// # use feanor_math::rings::poly::*;
46/// # use feanor_math::rings::poly::dense_poly::*;
47/// # use feanor_math::primitive_int::*;
48/// let ZZX = DensePolyRing::new(StaticRing::<i64>::RING, "X");
49/// let [f, g, expected] =
50///     ZZX.with_wrapped_indeterminate(|X| [X.pow_ref(2) - 2 * X + 1, X.pow_ref(2) - 1, X - 1]);
51/// assert_el_eq!(&ZZX, expected, <_ as PolyTFracGCDRing>::gcd(&ZZX, &f, &g));
52/// ```
53///
54/// # Implementation notes
55///
56/// Efficient implementations for polynomial gcds are often quite complicated, since the standard
57/// euclidean algorithm is only efficient over finite fields, where no coefficient explosion
58/// happens. The general idea for other rings/fields is to reduce it to the finite case, by
59/// considering the situation modulo a finite-index ideal. The requirements for this approach are
60/// defined by the trait [`PolyGCDLocallyDomain`], and there is a blanket impl `R: PolyTFracGCDRing
61/// where R: PolyGCDLocallyDomain`.
62///
63/// Note that this blanket impl used [`crate::specialization::FiniteRingSpecializable`] to use the
64/// standard algorithm whenever the corresponding ring is actually finite. In other words, despite
65/// the fact that the blanket implementation for `PolyGCDLocallyDomain`s also applies to finite
66/// fields, the local implementation is not actually used in these cases.
67pub trait PolyTFracGCDRing {
68    /// Computes the square-free part of a polynomial `f`, which is the largest-degree squarefree
69    /// polynomial `d` such that `d | a f` for some non-zero-divisor `a` of this ring.
70    ///
71    /// This value is unique up to multiplication by units. If the base ring is a field,
72    /// we impose the additional constraint that it be monic, which makes it unique.
73    ///
74    /// # Example
75    /// ```rust
76    /// # use feanor_math::assert_el_eq;
77    /// # use feanor_math::ring::*;
78    /// # use feanor_math::algorithms::poly_gcd::*;
79    /// # use feanor_math::rings::poly::*;
80    /// # use feanor_math::rings::poly::dense_poly::*;
81    /// # use feanor_math::primitive_int::*;
82    /// let ZZX = DensePolyRing::new(StaticRing::<i64>::RING, "X");
83    /// let [f] = ZZX.with_wrapped_indeterminate(|X| [1 - X.pow_ref(2)]);
84    /// assert_el_eq!(
85    ///     &ZZX,
86    ///     &f,
87    ///     <_ as PolyTFracGCDRing>::squarefree_part(&ZZX, &ZZX.mul_ref(&f, &f))
88    /// );
89    /// ```
90    fn squarefree_part<P>(poly_ring: P, poly: &El<P>) -> El<P>
91    where
92        P: RingStore + Copy,
93        P::Type: PolyRing + DivisibilityRing,
94        <P::Type as RingExtension>::BaseRing: RingStore<Type = Self>,
95    {
96        poly_ring.prod(Self::power_decomposition(poly_ring, poly).into_iter().map(|(f, _)| f))
97    }
98
99    /// Compute square-free polynomials `f1, f2, ...` such that `a f = f1 f2^2 f3^3 ...`
100    /// for some non-zero-divisor `a` of this ring. They are returned as tuples `(fi, i)`
101    /// where `deg(fi) > 0`.
102    ///
103    /// These values are unique up to multiplication by units. If the base ring is a field,
104    /// we impose the additional constraint that all `fi` be monic, which makes them unique.
105    fn power_decomposition<P>(poly_ring: P, poly: &El<P>) -> Vec<(El<P>, usize)>
106    where
107        P: RingStore + Copy,
108        P::Type: PolyRing + DivisibilityRing,
109        <P::Type as RingExtension>::BaseRing: RingStore<Type = Self>;
110
111    /// As [`PolyTFracGCDRing::power_decomposition()`], this writes a polynomial as a product
112    /// of powers of square-free polynomials. However, it additionally accepts a
113    /// [`ComputationController`] to customize the performed computation.
114    fn power_decomposition_with_controller<P, Controller>(
115        poly_ring: P,
116        poly: &El<P>,
117        _: Controller,
118    ) -> Vec<(El<P>, usize)>
119    where
120        P: RingStore + Copy,
121        P::Type: PolyRing + DivisibilityRing,
122        <P::Type as RingExtension>::BaseRing: RingStore<Type = Self>,
123        Controller: ComputationController,
124    {
125        Self::power_decomposition(poly_ring, poly)
126    }
127
128    /// Computes the greatest common divisor of two polynomials `f, g` over the fraction field,
129    /// which is the largest-degree polynomial `d` such that `d | a f, a g` for some
130    /// non-zero-divisor `a` of this ring.
131    ///
132    /// This value is unique up to multiplication by units. If the base ring is a field,
133    /// we impose the additional constraint that it be monic, which makes it unique.
134    ///
135    /// # Example
136    /// ```rust
137    /// # use feanor_math::assert_el_eq;
138    /// # use feanor_math::ring::*;
139    /// # use feanor_math::algorithms::poly_gcd::*;
140    /// # use feanor_math::rings::poly::*;
141    /// # use feanor_math::rings::poly::dense_poly::*;
142    /// # use feanor_math::primitive_int::*;
143    /// let ZZX = DensePolyRing::new(StaticRing::<i64>::RING, "X");
144    /// let [f, g, expected] = ZZX.with_wrapped_indeterminate(|X| [X.pow_ref(2) - 2 * X + 1, 2 * X.pow_ref(2) - 2, X - 1]);
145    /// // note that `expected` is not the gcd over `ZZ[X]` (which would be `2 X - 2`), but `X - 1`, i.e. the (monic) gcd over `QQ[X]`
146    /// assert_el_eq!(&ZZX, expected, <_ as PolyTFracGCDRing>::gcd(&ZZX, &f, &g));
147    ///
148    /// // of course, the result does not have to be monic
149    /// let [f, g, expected] = ZZX.with_wrapped_indeterminate(|X| [4 * X.pow_ref(2) - 1, 4 * X.pow_ref(2) - 4 * X + 1, - 2 * X + 1]);
150    /// assert_el_eq!(&ZZX, expected, <_ as PolyTFracGCDRing>::gcd(&ZZX, &f, &g));
151    /// ```
152    fn gcd<P>(poly_ring: P, lhs: &El<P>, rhs: &El<P>) -> El<P>
153    where
154        P: RingStore + Copy,
155        P::Type: PolyRing + DivisibilityRing,
156        <P::Type as RingExtension>::BaseRing: RingStore<Type = Self>;
157
158    /// As [`PolyTFracGCDRing::gcd()`], this computes the gcd of two polynomials.
159    /// However, it additionally accepts a [`ComputationController`] to customize
160    /// the performed computation.
161    fn gcd_with_controller<P, Controller>(poly_ring: P, lhs: &El<P>, rhs: &El<P>, _: Controller) -> El<P>
162    where
163        P: RingStore + Copy,
164        P::Type: PolyRing + DivisibilityRing,
165        <P::Type as RingExtension>::BaseRing: RingStore<Type = Self>,
166        Controller: ComputationController,
167    {
168        Self::gcd(poly_ring, lhs, rhs)
169    }
170}
171
172/// Computes the map
173/// ```text
174///   R[X] -> R[X],  f(X) -> a^(deg(f) - 1) f(X / a)
175/// ```
176/// that can be used to make polynomials over a domain monic (when setting `a = lc(f)`).
177#[stability::unstable(feature = "enable")]
178pub fn evaluate_aX<P>(poly_ring: P, f: &El<P>, a: &El<<P::Type as RingExtension>::BaseRing>) -> El<P>
179where
180    P: RingStore,
181    P::Type: PolyRing,
182    <<P::Type as RingExtension>::BaseRing as RingStore>::Type: DivisibilityRing + Domain,
183{
184    if poly_ring.is_zero(f) {
185        return poly_ring.zero();
186    }
187    let ring = poly_ring.base_ring();
188    let d = poly_ring.degree(f).unwrap();
189    let result = poly_ring.from_terms(poly_ring.terms(f).map(|(c, i)| {
190        if i == d {
191            (ring.checked_div(c, a).unwrap(), d)
192        } else {
193            (ring.mul_ref_fst(c, ring.pow(ring.clone_el(a), d - i - 1)), i)
194        }
195    }));
196    return result;
197}
198
199/// Computes the inverse to [`evaluate_aX()`].
200#[stability::unstable(feature = "enable")]
201pub fn unevaluate_aX<P>(poly_ring: P, g: &El<P>, a: &El<<P::Type as RingExtension>::BaseRing>) -> El<P>
202where
203    P: RingStore,
204    P::Type: PolyRing,
205    <<P::Type as RingExtension>::BaseRing as RingStore>::Type: DivisibilityRing + Domain,
206{
207    if poly_ring.is_zero(g) {
208        return poly_ring.zero();
209    }
210    let ring = poly_ring.base_ring();
211    let d = poly_ring.degree(g).unwrap();
212    let result = poly_ring.from_terms(poly_ring.terms(g).map(|(c, i)| {
213        if i == d {
214            (ring.mul_ref(c, a), d)
215        } else {
216            (ring.checked_div(c, &ring.pow(ring.clone_el(a), d - i - 1)).unwrap(), i)
217        }
218    }));
219    return result;
220}
221
222/// Given a polynomial `f` over a PID, returns `(f/cont(f), cont(f))`, where `cont(f)`
223/// is the content of `f`, i.e. the gcd of all coefficients of `f`.
224#[stability::unstable(feature = "enable")]
225pub fn make_primitive<P>(poly_ring: P, f: &El<P>) -> (El<P>, El<<P::Type as RingExtension>::BaseRing>)
226where
227    P: RingStore,
228    P::Type: PolyRing,
229    <<P::Type as RingExtension>::BaseRing as RingStore>::Type: PrincipalIdealRing + Domain,
230{
231    if poly_ring.is_zero(f) {
232        return (poly_ring.zero(), poly_ring.base_ring().one());
233    }
234    let ring = poly_ring.base_ring();
235    let content = poly_ring
236        .terms(f)
237        .map(|(c, _)| c)
238        .fold(ring.zero(), |a, b| ring.ideal_gen(&a, b));
239    let result = poly_ring.from_terms(
240        poly_ring
241            .terms(f)
242            .map(|(c, i)| (ring.checked_div(c, &content).unwrap(), i)),
243    );
244    return (result, content);
245}
246
247/// Checks whether there exists a polynomial `g` such that `g^k = f`, and if yes,
248/// returns `g`.
249///
250/// # Example
251/// ```rust
252/// # use feanor_math::assert_el_eq;
253/// # use feanor_math::ring::*;
254/// # use feanor_math::rings::poly::*;
255/// # use feanor_math::rings::poly::dense_poly::*;
256/// # use feanor_math::primitive_int::*;
257/// # use feanor_math::algorithms::poly_gcd::*;
258/// let poly_ring = DensePolyRing::new(StaticRing::<i64>::RING, "X");
259/// let [f, f_sqrt] = poly_ring.with_wrapped_indeterminate(|X| [X.pow_ref(2) + 2 * X + 1, X + 1]);
260/// assert_el_eq!(&poly_ring, f_sqrt, poly_root(&poly_ring, &f, 2).unwrap());
261/// ```
262#[stability::unstable(feature = "enable")]
263pub fn poly_root<P>(poly_ring: P, f: &El<P>, k: usize) -> Option<El<P>>
264where
265    P: RingStore,
266    P::Type: PolyRing,
267    <<P::Type as RingExtension>::BaseRing as RingStore>::Type: DivisibilityRing + Domain,
268{
269    assert!(poly_ring.degree(f).unwrap().is_multiple_of(k));
270    let d = poly_ring.degree(f).unwrap() / k;
271    let ring = poly_ring.base_ring();
272    let k_in_ring = ring.int_hom().map(k.try_into().unwrap());
273
274    let mut result_reversed = Vec::new();
275    result_reversed.push(ring.one());
276    for i in 1..=d {
277        let g = poly_ring.pow(
278            poly_ring.from_terms((0..i).map(|j| (ring.clone_el(&result_reversed[j]), j))),
279            k,
280        );
281        let partition_sum = poly_ring.coefficient_at(&g, i);
282        let next_coeff = ring.checked_div(
283            &ring.sub_ref(poly_ring.coefficient_at(f, k * d - i), partition_sum),
284            &k_in_ring,
285        )?;
286        result_reversed.push(next_coeff);
287    }
288
289    let result = poly_ring.from_terms(result_reversed.into_iter().enumerate().map(|(i, c)| (c, d - i)));
290    if poly_ring.eq_el(f, &poly_ring.pow(poly_ring.clone_el(&result), k)) {
291        return Some(result);
292    } else {
293        return None;
294    }
295}
296
297impl<R> PolyTFracGCDRing for R
298where
299    R: ?Sized + PolyGCDLocallyDomain + SelfIso,
300{
301    default fn power_decomposition<P>(poly_ring: P, poly: &El<P>) -> Vec<(El<P>, usize)>
302    where
303        P: RingStore + Copy,
304        P::Type: PolyRing + DivisibilityRing,
305        <P::Type as RingExtension>::BaseRing: RingStore<Type = Self>,
306    {
307        Self::power_decomposition_with_controller(poly_ring, poly, DontObserve)
308    }
309
310    default fn power_decomposition_with_controller<P, Controller>(
311        poly_ring: P,
312        poly: &El<P>,
313        controller: Controller,
314    ) -> Vec<(El<P>, usize)>
315    where
316        P: RingStore + Copy,
317        P::Type: PolyRing + DivisibilityRing,
318        <P::Type as RingExtension>::BaseRing: RingStore<Type = Self>,
319        Controller: ComputationController,
320    {
321        struct PowerDecompositionOperation<'a, P, Controller>(P, &'a El<P>, Controller)
322        where
323            P: RingStore + Copy,
324            P::Type: PolyRing,
325            <<P::Type as RingExtension>::BaseRing as RingStore>::Type: DivisibilityRing + SelfIso,
326            Controller: ComputationController;
327
328        impl<'a, P, Controller> FiniteRingOperation<<<P::Type as RingExtension>::BaseRing as RingStore>::Type>
329            for PowerDecompositionOperation<'a, P, Controller>
330        where
331            P: RingStore + Copy,
332            P::Type: PolyRing + DivisibilityRing,
333            <<P::Type as RingExtension>::BaseRing as RingStore>::Type:
334                PolyGCDLocallyDomain + DivisibilityRing + SelfIso,
335            Controller: ComputationController,
336        {
337            type Output = Vec<(El<P>, usize)>;
338
339            fn execute(self) -> Vec<(El<P>, usize)>
340            where
341                <<P::Type as RingExtension>::BaseRing as RingStore>::Type: FiniteRing,
342            {
343                static_assert_impls!(<<P::Type as RingExtension>::BaseRing as RingStore>::Type: SelfIso);
344
345                let new_poly_ring = DensePolyRing::new(
346                    AsField::from(AsFieldBase::promise_is_perfect_field(self.0.base_ring())),
347                    "X",
348                );
349                let new_poly = new_poly_ring.from_terms(self.0.terms(self.1).map(|(c, i)| {
350                    (
351                        new_poly_ring
352                            .base_ring()
353                            .get_ring()
354                            .rev_delegate(self.0.base_ring().clone_el(c)),
355                        i,
356                    )
357                }));
358                poly_power_decomposition_finite_field(&new_poly_ring, &new_poly, self.2)
359                    .into_iter()
360                    .map(|(f, k)| {
361                        (
362                            self.0.from_terms(new_poly_ring.terms(&f).map(|(c, i)| {
363                                (
364                                    new_poly_ring
365                                        .base_ring()
366                                        .get_ring()
367                                        .unwrap_element(new_poly_ring.base_ring().clone_el(c)),
368                                    i,
369                                )
370                            })),
371                            k,
372                        )
373                    })
374                    .collect()
375            }
376
377            fn fallback(self) -> Self::Output {
378                let poly_ring = self.0;
379                poly_power_decomposition_local(poly_ring, poly_ring.clone_el(self.1), self.2)
380            }
381        }
382
383        R::specialize(PowerDecompositionOperation(poly_ring, poly, controller))
384    }
385
386    default fn gcd<P>(poly_ring: P, lhs: &El<P>, rhs: &El<P>) -> El<P>
387    where
388        P: RingStore + Copy,
389        P::Type: PolyRing + DivisibilityRing,
390        <P::Type as RingExtension>::BaseRing: RingStore<Type = Self>,
391    {
392        Self::gcd_with_controller(poly_ring, lhs, rhs, DontObserve)
393    }
394
395    fn gcd_with_controller<P, Controller>(poly_ring: P, lhs: &El<P>, rhs: &El<P>, controller: Controller) -> El<P>
396    where
397        P: RingStore + Copy,
398        P::Type: PolyRing + DivisibilityRing,
399        <P::Type as RingExtension>::BaseRing: RingStore<Type = Self>,
400        Controller: ComputationController,
401    {
402        struct PolyGCDOperation<'a, P, Controller>(P, &'a El<P>, &'a El<P>, Controller)
403        where
404            P: RingStore + Copy,
405            P::Type: PolyRing,
406            <<P::Type as RingExtension>::BaseRing as RingStore>::Type: DivisibilityRing + SelfIso,
407            Controller: ComputationController;
408
409        impl<'a, P, Controller> FiniteRingOperation<<<P::Type as RingExtension>::BaseRing as RingStore>::Type>
410            for PolyGCDOperation<'a, P, Controller>
411        where
412            P: RingStore + Copy,
413            P::Type: PolyRing + DivisibilityRing,
414            <<P::Type as RingExtension>::BaseRing as RingStore>::Type:
415                PolyGCDLocallyDomain + DivisibilityRing + SelfIso,
416            Controller: ComputationController,
417        {
418            type Output = El<P>;
419
420            fn execute(self) -> El<P>
421            where
422                <<P::Type as RingExtension>::BaseRing as RingStore>::Type: FiniteRing,
423            {
424                static_assert_impls!(<<P::Type as RingExtension>::BaseRing as RingStore>::Type: SelfIso);
425
426                let new_poly_ring = DensePolyRing::new(
427                    AsField::from(AsFieldBase::promise_is_perfect_field(self.0.base_ring())),
428                    "X",
429                );
430                let new_lhs = new_poly_ring.from_terms(self.0.terms(self.1).map(|(c, i)| {
431                    (
432                        new_poly_ring
433                            .base_ring()
434                            .get_ring()
435                            .rev_delegate(self.0.base_ring().clone_el(c)),
436                        i,
437                    )
438                }));
439                let new_rhs = new_poly_ring.from_terms(self.0.terms(self.2).map(|(c, i)| {
440                    (
441                        new_poly_ring
442                            .base_ring()
443                            .get_ring()
444                            .rev_delegate(self.0.base_ring().clone_el(c)),
445                        i,
446                    )
447                }));
448                let result = new_poly_ring.normalize(fast_poly_eea(&new_poly_ring, new_lhs, new_rhs, self.3).2);
449                return self.0.from_terms(new_poly_ring.terms(&result).map(|(c, i)| {
450                    (
451                        new_poly_ring
452                            .base_ring()
453                            .get_ring()
454                            .unwrap_element(new_poly_ring.base_ring().clone_el(c)),
455                        i,
456                    )
457                }));
458            }
459
460            fn fallback(self) -> Self::Output {
461                let poly_ring = self.0;
462                poly_gcd_local(
463                    poly_ring,
464                    poly_ring.clone_el(self.1),
465                    poly_ring.clone_el(self.2),
466                    self.3,
467                )
468            }
469        }
470
471        R::specialize(PolyGCDOperation(poly_ring, lhs, rhs, controller))
472    }
473}
474
475#[cfg(test)]
476use crate::integer::*;
477#[cfg(test)]
478use crate::rings::extension::galois_field::GaloisField;
479#[cfg(test)]
480use crate::rings::zn::ZnRingStore;
481#[cfg(test)]
482use crate::rings::zn::zn_64;
483
484#[test]
485fn test_poly_root() {
486    let ring = BigIntRing::RING;
487    let poly_ring = DensePolyRing::new(ring, "X");
488    let [f] = poly_ring.with_wrapped_indeterminate(|X| {
489        [X.pow_ref(7) + X.pow_ref(6) + X.pow_ref(5) + X.pow_ref(4) + X.pow_ref(3) + X.pow_ref(2) + X + 1]
490    });
491    for k in 1..5 {
492        assert_el_eq!(
493            &poly_ring,
494            &f,
495            poly_root(&poly_ring, &poly_ring.pow(poly_ring.clone_el(&f), k), k).unwrap()
496        );
497    }
498
499    let [f] = poly_ring.with_wrapped_indeterminate(|X| {
500        [X.pow_ref(5) + 2 * X.pow_ref(4) + 3 * X.pow_ref(3) + 4 * X.pow_ref(2) + 5 * X + 6]
501    });
502    for k in 1..5 {
503        assert_el_eq!(
504            &poly_ring,
505            &f,
506            poly_root(&poly_ring, &poly_ring.pow(poly_ring.clone_el(&f), k), k).unwrap()
507        );
508    }
509}
510
511#[test]
512fn test_poly_gcd_galois_field() {
513    let field = GaloisField::new(5, 3);
514    let poly_ring = DensePolyRing::new(&field, "X");
515    let [f, g, f_g_gcd] = poly_ring.with_wrapped_indeterminate(|X| {
516        [
517            (X.pow_ref(2) + 2) * (X.pow_ref(5) + 1),
518            (X.pow_ref(2) + 2) * (X + 1) * (X + 2),
519            (X.pow_ref(2) + 2) * (X + 1),
520        ]
521    });
522    assert_el_eq!(&poly_ring, &f_g_gcd, <_ as PolyTFracGCDRing>::gcd(&poly_ring, &f, &g));
523}
524
525#[test]
526fn test_poly_gcd_prime_field() {
527    let field = zn_64::Zn::new(5).as_field().ok().unwrap();
528    let poly_ring = DensePolyRing::new(&field, "X");
529    let [f, g, f_g_gcd] = poly_ring.with_wrapped_indeterminate(|X| {
530        [
531            (X.pow_ref(2) + 2) * (X.pow_ref(5) + 1),
532            (X.pow_ref(2) + 2) * (X + 1) * (X + 2),
533            (X.pow_ref(2) + 2) * (X + 1),
534        ]
535    });
536    assert_el_eq!(&poly_ring, &f_g_gcd, <_ as PolyTFracGCDRing>::gcd(&poly_ring, &f, &g));
537}
feanor_math/algorithms/poly_gcd/mod.rs

feanor_math/algorithms/poly_gcd/
mod.rs
