feanor_math/algorithms/poly_gcd/
mod.rs

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

feanor_math/algorithms/poly_gcd/
mod.rs
