feanor_math/algorithms/poly_gcd/
mod.rs

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

feanor_math/algorithms/poly_gcd/
mod.rs
