feanor_math/algorithms/poly_gcd/
finite.rs

1use std::mem::swap;
2use std::cmp::max;
3
4use crate::algorithms::matmul::strassen::{dispatch_strassen_impl, strassen_mem_size};
5use crate::algorithms::int_factor::is_prime_power;
6use crate::algorithms::poly_div::{poly_div_rem_finite_reduced, PolyDivRemReducedError};
7use crate::computation::*;
8use crate::matrix::*;
9use crate::primitive_int::*;
10use crate::ring::*;
11use crate::rings::poly::*;
12use crate::field::*;
13use crate::pid::*;
14use crate::divisibility::*;
15use crate::homomorphism::*;
16use crate::integer::*;
17use crate::rings::finite::*;
18
19///
20/// Returns a list of `(fi, ki)` such that the `fi` are monic, square-free and pairwise coprime, and
21/// `f = a prod_i fi^ki` for a unit `a` of the base field.
22/// 
23#[stability::unstable(feature = "enable")]
24pub fn poly_power_decomposition_finite_field<P, Controller>(poly_ring: P, poly: &El<P>, controller: Controller) -> Vec<(El<P>, usize)>
25    where P: RingStore + Copy,
26        P::Type: PolyRing + EuclideanRing,
27        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: FiniteRing + Field,
28        Controller: ComputationController
29{
30    assert!(!poly_ring.is_zero(&poly));
31    let squarefree_part = poly_squarefree_part_finite_field(poly_ring, poly, controller.clone());
32    if poly_ring.degree(&squarefree_part).unwrap() == poly_ring.degree(&poly).unwrap() {
33        return vec![(squarefree_part, 1)];
34    } else {
35        let square_part = poly_ring.checked_div(&poly, &squarefree_part).unwrap();
36        let square_part_decomposition = poly_power_decomposition_finite_field(poly_ring, &square_part, controller);
37        let mut result = square_part_decomposition;
38        let mut degree = 0;
39        for (g, k) in &mut result {
40            *k += 1;
41            degree += poly_ring.degree(g).unwrap() * *k;
42        }
43        if degree != poly_ring.degree(&poly).unwrap() {
44            let remaining_part = poly_ring.checked_div(&poly, &poly_ring.prod(result.iter().map(|(g, e)| poly_ring.pow(poly_ring.clone_el(g), *e)))).unwrap();
45            result.insert(0, (poly_ring.normalize(remaining_part), 1));
46        }
47        return result;
48    }
49}
50
51///
52/// Computes the square-free part of a polynomial `f`, i.e. the greatest (w.r.t.
53/// divisibility) polynomial `g | f` that is square-free.
54/// 
55/// The returned polynomial is always monic, and with this restriction, it
56/// is unique.
57/// 
58#[stability::unstable(feature = "enable")]
59pub fn poly_squarefree_part_finite_field<P, Controller>(poly_ring: P, poly: &El<P>, controller: Controller) -> El<P>
60    where P: RingStore,
61        P::Type: PolyRing + PrincipalIdealRing,
62        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: FiniteRing + Field,
63        Controller: ComputationController
64{
65    assert!(!poly_ring.is_zero(&poly));
66    if poly_ring.degree(poly).unwrap() == 0 {
67        return poly_ring.one();
68    }
69    let derivate = derive_poly(&poly_ring, poly);
70    if poly_ring.is_zero(&derivate) {
71        let q = poly_ring.base_ring().size(&BigIntRing::RING).unwrap();
72        let (p, e) = is_prime_power(BigIntRing::RING, &q).unwrap();
73        let p = int_cast(p, StaticRing::<i64>::RING, BigIntRing::RING) as usize;
74        assert!(p > 0);
75        let undo_frobenius = Frobenius::new(poly_ring.base_ring(), e - 1);
76        let base_poly = poly_ring.from_terms(poly_ring.terms(poly).map(|(c, i)| {
77            debug_assert!(i % p == 0);
78            (undo_frobenius.map_ref(c), i / p)
79        }));
80        return poly_squarefree_part_finite_field(poly_ring, &base_poly, controller);
81    } else {
82        let square_part = poly_ring.ideal_gen(poly, &derivate);
83        let result = poly_ring.checked_div(poly, &square_part).unwrap();
84        return poly_ring.normalize(result);
85    }
86}
87
88const FAST_POLY_EEA_THRESHOLD: usize = 32;
89
90///
91/// Computes linearly independent vectors `(s, t)` and `(s', t')` and `a, a'` 
92/// such that `a = s * lhs + t * rhs` and `a' = s' * lhs + t' * rhs` are both of
93/// degree at most `target_deg`, or alternatively one of them is zero (in which
94/// case the degree cannot be reduced further).
95/// 
96/// The degrees of `s, t, s', t'` are bounded as
97/// ```text
98///   deg(s) < deg(rhs) - deg(a)
99///   deg(t) < deg(lhs) - deg(a)
100///   deg(s') < deg(rhs) - deg(a')
101///   deg(t') < deg(lhs) - deg(a')
102/// ```
103///
104fn partial_eea<P>(ring: P, lhs: El<P>, rhs: El<P>, target_deg: usize) -> ([El<P>; 4], [El<P>; 2])
105    where P: RingStore + Copy,
106        P::Type: PolyRing + EuclideanRing,
107        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: Field
108{
109    if ring.is_zero(&lhs) || ring.is_zero(&rhs) {
110        return ([ring.one(), ring.zero(), ring.zero(), ring.one()], [lhs, rhs]);
111    }
112    let (mut a, mut b) = (ring.clone_el(&lhs), ring.clone_el(&rhs));
113    let (mut sa, mut ta) = (ring.one(), ring.zero());
114    let (mut sb, mut tb) = (ring.zero(), ring.one());
115    
116    if ring.degree(&a).unwrap() < ring.degree(&b).unwrap() {
117        swap(&mut a, &mut b);
118        swap(&mut sa, &mut sb);
119        swap(&mut ta, &mut tb);
120    }
121
122    while ring.degree(&a).unwrap() > target_deg && !ring.is_zero(&b) {
123        debug_assert!(ring.eq_el(&a, &ring.add(ring.mul_ref(&sa, &lhs), ring.mul_ref(&ta, &rhs))));
124        debug_assert!(ring.eq_el(&b, &ring.add(ring.mul_ref(&sb, &lhs), ring.mul_ref(&tb, &rhs))));
125
126        let (quo, rem) = ring.euclidean_div_rem(a, &b);
127        ta = ring.sub(ta, ring.mul_ref(&quo, &tb));
128        sa = ring.sub(sa, ring.mul_ref_snd(quo, &sb));
129        a = rem;
130
131        swap(&mut a, &mut b);
132        swap(&mut sa, &mut sb);
133        swap(&mut ta, &mut tb);
134        
135        debug_assert_eq!(ring.degree(&sb).unwrap(), ring.degree(&rhs).unwrap() - ring.degree(&a).unwrap());
136        debug_assert_eq!(ring.degree(&tb).unwrap(), ring.degree(&lhs).unwrap() - ring.degree(&a).unwrap());
137    }
138    return ([sa, ta, sb, tb], [a, b]);
139}
140
141///
142/// Computes a Bezout identity for polynomials, using a fast divide-and-conquer
143/// polynomial gcd algorithm. Unless you are implementing [`crate::pid::PrincipalIdealRing`]
144/// for a custom type, you should use [`crate::pid::PrincipalIdealRing::extended_ideal_gen()`]
145/// to get a Bezout identity instead.
146/// 
147/// A Bezout identity is exactly as specified by [`crate::pid::PrincipalIdealRing::extended_ideal_gen()`],
148/// i.e. `s, t, d` such that `d` is the gcd of `lhs` and `rhs`, and `d = lhs * s + rhs * t`.
149/// Note that this algorithm does not try to avoid coefficient growth, and thus is only fast
150/// over finite fields. Furthermore, it will fall back to a slightly less efficient variant of
151/// the standard Euclidean algorithm on small inputs.
152/// 
153#[stability::unstable(feature = "enable")]
154pub fn fast_poly_eea<P, Controller>(poly_ring: P, lhs: El<P>, rhs: El<P>, controller: Controller) -> (El<P>, El<P>, El<P>)
155    where P: RingStore + Copy,
156        P::Type: PolyRing + EuclideanRing,
157        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: Field,
158        Controller: ComputationController
159{
160    fn fast_poly_eea_impl<P, Controller>(poly_ring: P, lhs: El<P>, rhs: El<P>, target_deg: usize, controller: Controller, memory: &mut [El<P>]) -> ([El<P>; 4], [El<P>; 2])
161        where P: RingStore + Copy,
162            P::Type: PolyRing + EuclideanRing,
163            <<P::Type as RingExtension>::BaseRing as RingStore>::Type: Field,
164            Controller: ComputationController
165    {
166        if poly_ring.is_zero(&lhs) || poly_ring.is_zero(&rhs) {
167            return ([poly_ring.one(), poly_ring.zero(), poly_ring.zero(), poly_ring.one()], [lhs, rhs]);
168        }
169        let ldeg = poly_ring.degree(&lhs).unwrap();
170        let rdeg = poly_ring.degree(&rhs).unwrap();
171        if ldeg < target_deg + FAST_POLY_EEA_THRESHOLD || rdeg < target_deg + FAST_POLY_EEA_THRESHOLD {
172            log_progress!(controller, ".");
173            return partial_eea(poly_ring, lhs, rhs, target_deg);
174        } else if ldeg >= 2 * rdeg {
175            let (mut q, r) = poly_ring.euclidean_div_rem(lhs, &rhs);
176            poly_ring.negate_inplace(&mut q);
177            let (transform, rest) = fast_poly_eea_impl(poly_ring, r, rhs, target_deg, controller, memory);
178            let mut transform: (_, _, _, _) = transform.into();
179            transform.1 = poly_ring.fma(&q, &transform.0, transform.1);
180            transform.3 = poly_ring.fma(&q, &transform.2, transform.3);
181            return (transform.into(), rest);
182        } else if rdeg >= 2 * ldeg {
183            let (transform, rest) = fast_poly_eea_impl(poly_ring, rhs, lhs, target_deg, controller, memory);
184            let transform: (_, _, _, _) = transform.into();
185            return ([transform.1, transform.0, transform.3, transform.2], rest);
186        }
187        let split_deg = max(ldeg, rdeg) / 3;
188        assert!(2 * split_deg + 1 < max(ldeg, rdeg));
189        let part_target_deg = max(split_deg, target_deg.saturating_sub(split_deg));
190
191        let lhs_upper = poly_ring.from_terms(poly_ring.terms(&lhs).filter(|(_, i)| *i >= split_deg).map(|(c, i)| (poly_ring.base_ring().clone_el(c), i - split_deg)));
192        let mut lhs_lower = lhs;
193        poly_ring.truncate_monomials(&mut lhs_lower, split_deg);
194        let rhs_upper = poly_ring.from_terms(poly_ring.terms(&rhs).filter(|(_, i)| *i >= split_deg).map(|(c, i)| (poly_ring.base_ring().clone_el(c), i - split_deg)));
195        let mut rhs_lower = rhs;
196        poly_ring.truncate_monomials(&mut rhs_lower, split_deg);
197
198        log_progress!(controller, "({},{})", max(poly_ring.degree(&lhs_upper).unwrap(), poly_ring.degree(&rhs_upper).unwrap()), part_target_deg);
199        let (fst_transform, [mut lhs_rest, mut rhs_rest]) = fast_poly_eea_impl(poly_ring, lhs_upper, rhs_upper, part_target_deg, controller.clone(), memory);
200        poly_ring.mul_assign_monomial(&mut lhs_rest, split_deg);
201        poly_ring.mul_assign_monomial(&mut rhs_rest, split_deg);
202
203        lhs_rest = poly_ring.fma(&fst_transform[0], &lhs_lower, lhs_rest);
204        lhs_rest = poly_ring.fma(&fst_transform[1], &rhs_lower, lhs_rest);
205        rhs_rest = poly_ring.fma(&fst_transform[2], &lhs_lower, rhs_rest);
206        rhs_rest = poly_ring.fma(&fst_transform[3], &rhs_lower, rhs_rest);
207
208        log_progress!(controller, "({},{})", max(poly_ring.degree(&lhs_rest).unwrap_or(0), poly_ring.degree(&rhs_rest).unwrap_or(0)), target_deg);
209        let (snd_transform, rest) = fast_poly_eea_impl(poly_ring, lhs_rest, rhs_rest, target_deg, controller, memory);
210
211        // multiply snd_transform * fst_transform
212        let mut result = [poly_ring.zero(), poly_ring.zero(), poly_ring.zero(), poly_ring.zero()];
213        dispatch_strassen_impl::<_, _, _, _, false, _, _, _>(
214            1, 
215            0, 
216            TransposableSubmatrix::from(Submatrix::from_1d(&snd_transform, 2, 2)),
217            TransposableSubmatrix::from(Submatrix::from_1d(&fst_transform, 2, 2)),
218            TransposableSubmatrixMut::from(SubmatrixMut::from_1d(&mut result, 2, 2)), 
219            poly_ring, 
220            memory
221        );
222
223        return (result, rest);
224    }
225
226    if poly_ring.is_zero(&lhs) {
227        return (poly_ring.zero(), poly_ring.one(), rhs);
228    } else if poly_ring.is_zero(&rhs) {
229        return (poly_ring.one(), poly_ring.zero(), lhs);
230    }
231
232    controller.run_computation(format_args!("fast_poly_eea(ldeg={}, rdeg={})", poly_ring.degree(&lhs).unwrap(), poly_ring.degree(&rhs).unwrap()), |controller| {
233        let ([s1, t1, s2, t2], [a1, a2]) = fast_poly_eea_impl(poly_ring, lhs, rhs, 0, controller,  &mut (0..strassen_mem_size(false, 2, 0)).map(|_| poly_ring.zero()).collect::<Vec<_>>());
234        if poly_ring.is_zero(&a1) {
235            return (s2, t2, a2);
236        } else {
237            assert!(poly_ring.is_zero(&a2));
238            return (s1, t1, a1);
239        }
240    })
241}
242
243///
244/// Computes the gcd of two polynomials over a finite and reduced ring.
245/// 
246/// This is well-defined, since a finite reduced ring is always a product of
247/// finite fields.
248/// 
249/// If the ring is not reduced, this function may fail and return `Err(nil)`, where
250/// `nil` is a nilpotent element of the ring. However, as long as the gcd of `lhs` and
251/// `rhs` exists, this function may alternatively return it, even in cases where the ring is 
252/// not reduced (note however that over a non-reduced ring, the gcd does not always exist).
253/// 
254#[stability::unstable(feature = "enable")]
255pub fn poly_gcd_finite_reduced<P>(poly_ring: P, mut lhs: El<P>, mut rhs: El<P>) -> Result<El<P>, El<<P::Type as RingExtension>::BaseRing>>
256    where P: RingStore + Copy,
257        P::Type: PolyRing,
258        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: FiniteRing + PrincipalIdealRing
259{
260    while !poly_ring.is_zero(&rhs) {
261        match poly_div_rem_finite_reduced(poly_ring, poly_ring.clone_el(&lhs), &rhs) {
262            Ok((_q, r)) => {
263                lhs = r;
264                std::mem::swap(&mut lhs, &mut rhs);
265            },
266            Err(PolyDivRemReducedError::NotReduced(nilpotent)) => return Err(nilpotent),
267            Err(PolyDivRemReducedError::NotDivisibleByContent(content_rhs)) => {
268                // we find a decomposition `R ~ R/c x R/Ann(c)` for the content `c` of `rhs`.
269                // clearly the gcd must be `lhs` modulo `c` (since `rhs = 0 mod c`); furthermore,
270                // modulo `Ann(c)` the content `c` is a unit, so `gcd(lhs, rhs) = gcd(c * lhs, rhs) mod Ann(c)`
271                let content_ann = poly_ring.base_ring().annihilator(&content_rhs);
272                if !poly_ring.base_ring().is_unit(&poly_ring.base_ring().ideal_gen(&content_rhs, &content_ann)) {
273                    return Err(poly_ring.base_ring().annihilator(&poly_ring.base_ring().ideal_gen(&content_rhs, &content_ann)));
274                }
275                let mod_content_gcd = poly_gcd_finite_reduced(poly_ring, poly_ring.inclusion().mul_ref_map(&lhs, &content_rhs), rhs)?;
276                debug_assert!(poly_ring.terms(&mod_content_gcd).all(|(c, _)| poly_ring.base_ring().divides(c, &content_rhs)));
277                return Ok(poly_ring.add(
278                    poly_ring.inclusion().mul_ref_map(&lhs, &content_ann),
279                    mod_content_gcd
280                ));
281            }
282        }
283    }
284    return Ok(lhs);
285}
286
287#[cfg(test)]
288use crate::rings::zn::zn_64;
289#[cfg(test)]
290use crate::rings::zn::zn_rns::*;
291#[cfg(test)]
292use crate::rings::poly::dense_poly::DensePolyRing;
293#[cfg(test)]
294use crate::rings::zn::*;
295#[cfg(test)]
296use crate::seq::VectorView;
297#[cfg(test)]
298use crate::algorithms::poly_div::poly_checked_div_finite_reduced;
299
300#[test]
301fn test_poly_gcd_finite_reduced() {
302    let base_ring = Zn::new([5, 7, 11, 13].into_iter().map(|p| zn_64::Zn::new(p).as_field().ok().unwrap()).collect(), StaticRing::<i64>::RING);
303    let poly_ring = DensePolyRing::new(&base_ring, "X");
304    let component_poly_rings: [_; 4] = std::array::from_fn(|i| DensePolyRing::new(base_ring.at(i), "X"));
305
306    let [f0, g0, expected0] = component_poly_rings[0].with_wrapped_indeterminate(|X| [
307        (X.pow_ref(2) + 2) * (X.pow_ref(3) + X + 1),
308        (X + 1) * (X + 2) * (X.pow_ref(3) + X + 1),
309        X.pow_ref(3) + X + 1
310    ]);
311
312    let f1 = component_poly_rings[1].zero();
313    let [g1, expected1] = component_poly_rings[1].with_wrapped_indeterminate(|X| [
314        X.pow_ref(3) + X + 1,
315        X.pow_ref(3) + X + 1
316    ]);
317
318    let f2 = component_poly_rings[2].int_hom().map(1);
319    let g2 = component_poly_rings[2].zero();
320    let expected2 = component_poly_rings[2].one();
321    
322    let f3 = component_poly_rings[3].zero();
323    let g3 = component_poly_rings[3].zero();
324    let expected3 = component_poly_rings[3].zero();
325
326    fn reconstruct<'a, 'b, R>(polys: [El<DensePolyRing<&'a R>>; 4], poly_rings: &[DensePolyRing<&'a R>; 4], poly_ring: &DensePolyRing<&'b Zn<R, StaticRing<i64>>>) -> El<DensePolyRing<&'b Zn<R, StaticRing<i64>>>>
327        where R: RingStore,
328            R::Type: ZnRing + CanHomFrom<StaticRingBase<i64>>
329    {
330        poly_ring.from_terms((0..10).map(|i| (poly_ring.base_ring().from_congruence(polys.iter().zip(poly_rings.iter()).map(|(f, P)| P.base_ring().clone_el(P.coefficient_at(f, i)))), i)))
331    }
332    
333    let f = reconstruct([f0, f1, f2, f3], &component_poly_rings, &poly_ring);
334    let g = reconstruct([g0, g1, g2, g3], &component_poly_rings, &poly_ring);
335    let expected = reconstruct([expected0, expected1, expected2, expected3], &component_poly_rings, &poly_ring);
336    let actual = poly_gcd_finite_reduced(&poly_ring, poly_ring.clone_el(&f), poly_ring.clone_el(&g)).ok().unwrap();
337
338    assert!(poly_checked_div_finite_reduced(&poly_ring, poly_ring.clone_el(&actual), poly_ring.clone_el(&expected)).ok().unwrap().is_some());
339    assert!(poly_checked_div_finite_reduced(&poly_ring, poly_ring.clone_el(&expected), poly_ring.clone_el(&actual)).ok().unwrap().is_some());
340}
341
342#[test]
343fn test_partial_eea() {
344    let field = zn_64::Zn::new(65537).as_field().ok().unwrap();
345    let poly_ring = DensePolyRing::new(field, "X");
346    let [f, g] = poly_ring.with_wrapped_indeterminate(|X| [
347        X.pow_ref(9) - X.pow_ref(7) + 3 * X.pow_ref(2) - 1,
348        X.pow_ref(10) + X.pow_ref(6) + 1,
349    ]);
350
351    for k in (1..10).rev() {
352        let ([s1, t1, s2, t2], [a, b]) = partial_eea(&poly_ring, poly_ring.clone_el(&f), poly_ring.clone_el(&g), k);
353        assert_el_eq!(&poly_ring, poly_ring.add(poly_ring.mul_ref(&s1, &f), poly_ring.mul_ref(&t1, &g)), &a);
354        assert_el_eq!(&poly_ring, poly_ring.add(poly_ring.mul_ref(&s2, &f), poly_ring.mul_ref(&t2, &g)), &b);
355        assert_eq!(k, poly_ring.degree(&a).unwrap());
356        assert_eq!(k - 1, poly_ring.degree(&b).unwrap());
357        assert!(poly_ring.degree(&s1).is_none() || poly_ring.degree(&s1).unwrap() < 10 - poly_ring.degree(&a).unwrap());
358        assert!(poly_ring.degree(&t1).is_none() || poly_ring.degree(&t1).unwrap() < 9 - poly_ring.degree(&a).unwrap());
359        assert!(poly_ring.degree(&s2).is_none() || poly_ring.degree(&s2).unwrap() < 10 - poly_ring.degree(&b).unwrap());
360        assert!(poly_ring.degree(&t2).is_none() || poly_ring.degree(&t2).unwrap() < 9 - poly_ring.degree(&b).unwrap());
361    }
362}
363
364#[test]
365fn test_fast_poly_eea() {
366    let field = zn_64::Zn::new(65537).as_field().ok().unwrap();
367    let poly_ring = DensePolyRing::new(field, "X");
368    let [f, g] = poly_ring.with_wrapped_indeterminate(|X| [
369        X.pow_ref(90) - X.pow_ref(70) + 3 * X.pow_ref(20) - 1,
370        X.pow_ref(100) + X.pow_ref(60) + 1,
371    ]);
372
373    let (s, t, d) = fast_poly_eea(&poly_ring, poly_ring.clone_el(&f), poly_ring.clone_el(&g), LOG_PROGRESS);
374    assert!(poly_ring.is_unit(&d));
375    assert_el_eq!(&poly_ring, &d, poly_ring.add(poly_ring.mul_ref(&s, &f), poly_ring.mul_ref(&t, &g)));
376
377    let [f, g] = poly_ring.with_wrapped_indeterminate(|X| [
378        X.pow_ref(9) - X.pow_ref(7) + 3 * X.pow_ref(2) - 1,
379        X.pow_ref(100) + X.pow_ref(60) + 1,
380    ]);
381
382    let (s, t, d) = fast_poly_eea(&poly_ring, poly_ring.clone_el(&f), poly_ring.clone_el(&g), LOG_PROGRESS);
383    assert!(poly_ring.is_unit(&d));
384    assert_el_eq!(&poly_ring, &d, poly_ring.add(poly_ring.mul_ref(&s, &f), poly_ring.mul_ref(&t, &g)));
385}
feanor_math/algorithms/poly_gcd/finite.rs

feanor_math/algorithms/poly_gcd/
finite.rs
