feanor_math/rings/poly/
sparse_poly.rs

1use crate::algorithms;
2use crate::algorithms::convolution::ConvolutionAlgorithm;
3use crate::algorithms::poly_gcd::PolyTFracGCDRing;
4use crate::divisibility::*;
5use crate::integer::IntegerRing;
6use crate::integer::IntegerRingStore;
7use crate::pid::*;
8use crate::field::Field;
9use crate::seq::VectorViewMut;
10use crate::ring::*;
11use crate::rings::poly::*;
12use crate::seq::*;
13use crate::seq::sparse::*;
14
15use std::alloc::Allocator;
16use std::cmp::max;
17use std::rc::Rc;
18
19///
20/// The univariate polynomial ring `R[X]`. Polynomials are stored as sparse vectors of
21/// coefficients, thus giving improved performance in the case that most coefficients are
22/// zero.
23/// 
24/// # Example
25/// ```
26/// # use feanor_math::ring::*;
27/// # use feanor_math::homomorphism::*;
28/// # use feanor_math::rings::poly::*;
29/// # use feanor_math::rings::poly::sparse_poly::*;
30/// # use feanor_math::primitive_int::*;
31/// 
32/// let ZZ = StaticRing::<i32>::RING;
33/// let P = SparsePolyRing::new(ZZ, "X");
34/// let x10_plus_1 = P.add(P.pow(P.indeterminate(), 10), P.int_hom().map(1));
35/// let power = P.pow(x10_plus_1, 10);
36/// assert_eq!(0, *P.coefficient_at(&power, 1));
37/// ```
38/// This ring has a [`CanIsoFromTo`] to [`dense_poly::DensePolyRingBase`].
39/// ```
40/// # use feanor_math::assert_el_eq;
41/// # use feanor_math::homomorphism::*;
42/// # use feanor_math::ring::*;
43/// # use feanor_math::rings::poly::*;
44/// # use feanor_math::rings::poly::dense_poly::*;
45/// # use feanor_math::rings::poly::sparse_poly::*;
46/// # use feanor_math::primitive_int::*;
47/// 
48/// let ZZ = StaticRing::<i32>::RING;
49/// let P = SparsePolyRing::new(ZZ, "X");
50/// let P2 = DensePolyRing::new(ZZ, "X");
51/// let high_power_of_x = P.pow(P.indeterminate(), 10);
52/// assert_el_eq!(P2, P2.pow(P2.indeterminate(), 10), &P.can_iso(&P2).unwrap().map(high_power_of_x));
53/// ```
54/// 
55pub struct SparsePolyRingBase<R: RingStore> {
56    base_ring: Rc<R>,
57    unknown_name: &'static str,
58    zero: El<R>
59}
60
61impl<R: RingStore + Clone> Clone for SparsePolyRingBase<R> {
62    
63    fn clone(&self) -> Self {
64        SparsePolyRingBase {
65            base_ring: self.base_ring.clone(), 
66            unknown_name: self.unknown_name, 
67            zero: self.base_ring.zero()
68        }
69    }
70}
71
72///
73/// The univariate polynomial ring `R[X]`, with polynomials being stored as sparse vectors of coefficients.
74/// For details, see [`SparsePolyRingBase`].
75/// 
76#[allow(type_alias_bounds)]
77pub type SparsePolyRing<R: RingStore> = RingValue<SparsePolyRingBase<R>>;
78
79impl<R: RingStore> SparsePolyRing<R> {
80
81    pub fn new(base_ring: R, unknown_name: &'static str) -> Self {
82        let zero = base_ring.zero();
83        Self::from(SparsePolyRingBase { 
84            base_ring: Rc::new(base_ring), 
85            unknown_name: unknown_name, 
86            zero: zero
87        })
88    }
89}
90
91impl<R: RingStore> SparsePolyRingBase<R> {
92
93    fn degree_truncate(&self, el: &mut SparseMapVector<Rc<R>>) {
94        for i in (0..el.len()).rev() {
95            if !self.base_ring.is_zero(&el.at(i)) {
96                el.set_len(i + 1);
97                return;
98            }
99        }
100        el.set_len(0);
101    }
102
103    fn poly_div<F>(&self, lhs: &mut <Self as RingBase>::Element, rhs: &<Self as RingBase>::Element, mut left_div_lc: F) -> Option<<Self as RingBase>::Element>
104        where F: FnMut(El<R>) -> Option<El<R>>
105    {
106        let lhs_val = std::mem::replace(lhs, self.zero());
107        let (quo, rem) = algorithms::poly_div::poly_div_rem(
108            RingRef::new(self), 
109            lhs_val, 
110            rhs, 
111            |x| left_div_lc(self.base_ring().clone_el(x)).ok_or(())
112        ).ok()?;
113        *lhs = rem;
114        return Some(quo);
115    }
116}
117
118pub struct SparsePolyRingEl<R: RingStore> {
119    data: SparseMapVector<Rc<R>>
120}
121
122impl<R: RingStore> RingBase for SparsePolyRingBase<R> {
123    
124    type Element = SparsePolyRingEl<R>;
125
126    fn clone_el(&self, val: &Self::Element) -> Self::Element {
127        SparsePolyRingEl {
128            data: val.data.clone()
129        }
130    }
131
132    fn add_assign_ref(&self, lhs: &mut Self::Element, rhs: &Self::Element) {
133        lhs.data.set_len(max(lhs.data.len(), rhs.data.len()));
134        for (i, c) in rhs.data.nontrivial_entries() {
135            self.base_ring.add_assign_ref(lhs.data.at_mut(i), c);
136        }
137        self.degree_truncate(&mut lhs.data);
138    }
139
140    fn add_assign(&self, lhs: &mut Self::Element, rhs: Self::Element) {
141        self.add_assign_ref(lhs, &rhs);
142    }
143
144    fn sub_assign_ref(&self, lhs: &mut Self::Element, rhs: &Self::Element) {
145        lhs.data.set_len(max(lhs.data.len(), rhs.data.len()));
146        for (i, c) in rhs.data.nontrivial_entries() {
147            self.base_ring.sub_assign_ref(lhs.data.at_mut(i), c);
148        }
149        self.degree_truncate(&mut lhs.data);
150    }
151
152    fn negate_inplace(&self, lhs: &mut Self::Element) {
153        lhs.data.scan(|_, c| self.base_ring.negate_inplace(c));
154    }
155
156    fn mul_assign(&self, lhs: &mut Self::Element, rhs: Self::Element) {
157        self.mul_assign_ref(lhs, &rhs);
158    }
159
160    fn mul_assign_ref(&self, lhs: &mut Self::Element, rhs: &Self::Element) {
161        *lhs = self.mul_ref(lhs, rhs);
162    }
163
164    fn zero(&self) -> Self::Element {
165        SparsePolyRingEl {
166            data: SparseMapVector::new(0, self.base_ring.clone())
167        }
168    }
169    
170    fn from_int(&self, value: i32) -> Self::Element {
171        self.from(self.base_ring().get_ring().from_int(value))
172    }
173
174    fn eq_el(&self, lhs: &Self::Element, rhs: &Self::Element) -> bool {
175        if lhs.data.len() != rhs.data.len() {
176            return false;
177        }
178        for (i, c) in lhs.data.nontrivial_entries() {
179            if !self.base_ring.eq_el(rhs.data.at(i), c) {
180                return false;
181            }
182        }
183        for (i, c) in rhs.data.nontrivial_entries() {
184            if !self.base_ring.eq_el(lhs.data.at(i), c) {
185                return false;
186            }
187        }
188        return true;
189    }
190
191    fn is_commutative(&self) -> bool {
192        self.base_ring.is_commutative()
193    }
194
195    fn is_noetherian(&self) -> bool {
196        // by Hilbert's basis theorem
197        self.base_ring.is_noetherian()
198    }
199
200    fn dbg_within<'a>(&self, value: &Self::Element, out: &mut std::fmt::Formatter<'a>, env: EnvBindingStrength) -> std::fmt::Result {
201        super::generic_impls::dbg_poly(self, value, out, self.unknown_name, env)
202    }
203
204    fn dbg<'a>(&self, value: &Self::Element, out: &mut std::fmt::Formatter<'a>) -> std::fmt::Result {
205        self.dbg_within(value, out, EnvBindingStrength::Weakest)
206    }
207
208    fn square(&self, value: &mut Self::Element) {
209        *value = self.mul_ref(&value, &value);
210    }
211
212    fn mul_ref(&self, lhs: &Self::Element, rhs: &Self::Element) -> Self::Element {
213        if lhs.data.len() == 0 || rhs.data.len() == 0 {
214            return self.zero();
215        }
216        let mut result = SparseMapVector::new(lhs.data.len() + rhs.data.len() - 1, self.base_ring.clone());
217        for (i, c1) in lhs.data.nontrivial_entries() {
218            for (j, c2) in rhs.data.nontrivial_entries() {
219                self.base_ring.add_assign(result.at_mut(i + j), self.base_ring.mul_ref(c1, c2));
220            }
221        }
222        // if the ring is not zero-divisor free
223        self.degree_truncate(&mut result);
224        return SparsePolyRingEl {
225            data: result
226        };
227    }
228
229    fn mul_assign_int(&self, lhs: &mut Self::Element, rhs: i32) {
230        if rhs == 0 {
231            *lhs = self.zero();
232        } else {
233            lhs.data.scan(|_, c| self.base_ring.int_hom().mul_assign_map(c, rhs));
234        }
235    }
236    
237    fn characteristic<I: IntegerRingStore + Copy>(&self, ZZ: I) -> Option<El<I>>
238        where I::Type: IntegerRing
239    {
240        self.base_ring().characteristic(ZZ)
241    }
242    
243    fn is_approximate(&self) -> bool {
244        self.base_ring().get_ring().is_approximate()
245    }
246}
247
248impl<R> PartialEq for SparsePolyRingBase<R> 
249    where R: RingStore
250{
251    fn eq(&self, other: &Self) -> bool {
252        self.base_ring.get_ring() == other.base_ring.get_ring()
253    }
254}
255
256impl<R: RingStore> RingExtension for SparsePolyRingBase<R> {
257    
258    type BaseRing = R;
259
260    fn base_ring<'a>(&'a self) -> &'a Self::BaseRing {
261        &self.base_ring
262    }
263
264    fn from(&self, x: El<Self::BaseRing>) -> Self::Element {
265        let mut result = self.zero();
266        if !self.base_ring().is_zero(&x) {
267            result.data.set_len(1);
268            *result.data.at_mut(0) = x;
269        }
270        return result;
271    }
272}
273
274pub trait ImplGenericCanIsoFromToMarker: PolyRing {}
275
276impl<R, A, C> ImplGenericCanIsoFromToMarker for dense_poly::DensePolyRingBase<R, A, C> 
277    where R: RingStore, A: Allocator + Clone, C: ConvolutionAlgorithm<R::Type>
278{}
279
280impl<R, P> CanHomFrom<P> for SparsePolyRingBase<R> 
281    where R: RingStore, R::Type: CanHomFrom<<P::BaseRing as RingStore>::Type>, P: ImplGenericCanIsoFromToMarker
282{
283    type Homomorphism = super::generic_impls::Homomorphism<P, Self>;
284
285    fn has_canonical_hom(&self, from: &P) -> Option<Self::Homomorphism> {
286        super::generic_impls::has_canonical_hom(from, self)
287    }
288
289    fn map_in(&self, from: &P, el: P::Element, hom: &Self::Homomorphism) -> Self::Element {
290        super::generic_impls::map_in(from, self, el, hom)
291    }
292}
293
294impl<R1, R2> CanHomFrom<SparsePolyRingBase<R1> > for SparsePolyRingBase<R2> 
295    where R1: RingStore, R2: RingStore, R2::Type: CanHomFrom<R1::Type>
296{
297    type Homomorphism = <R2::Type as CanHomFrom<R1::Type>>::Homomorphism;
298
299    fn has_canonical_hom(&self, from: &SparsePolyRingBase<R1>) -> Option<Self::Homomorphism> {
300        self.base_ring().get_ring().has_canonical_hom(from.base_ring().get_ring())
301    }
302
303    fn map_in_ref(&self, from: &SparsePolyRingBase<R1> , el: &SparsePolyRingEl<R1>, hom: &Self::Homomorphism) -> Self::Element {
304        let mut result = SparseMapVector::new(el.data.len(), self.base_ring.clone());
305        for (j, c) in el.data.nontrivial_entries() {
306            *result.at_mut(j) = self.base_ring().get_ring().map_in_ref(from.base_ring().get_ring(), c, hom);
307        }
308        return SparsePolyRingEl {
309            data: result
310        };
311    }
312
313    fn map_in(&self, from: &SparsePolyRingBase<R1>, el: <SparsePolyRingBase<R1> as RingBase>::Element, hom: &Self::Homomorphism) -> Self::Element {
314        self.map_in_ref(from, &el, hom)
315    }
316}
317
318impl<R, P> CanIsoFromTo<P> for SparsePolyRingBase<R> 
319    where R: RingStore, R::Type: CanIsoFromTo<<P::BaseRing as RingStore>::Type>, P: ImplGenericCanIsoFromToMarker
320{
321    type Isomorphism = super::generic_impls::Isomorphism<P, Self>;
322
323    fn has_canonical_iso(&self, from: &P) -> Option<Self::Isomorphism> {
324        self.base_ring().get_ring().has_canonical_iso(from.base_ring().get_ring())
325    }
326
327    fn map_out(&self, from: &P, el: Self::Element, iso: &Self::Isomorphism) -> P::Element {
328        super::generic_impls::map_out(from, self, el, iso)
329    }
330}
331
332impl<R1, R2> CanIsoFromTo<SparsePolyRingBase<R1>> for SparsePolyRingBase<R2> 
333    where R1: RingStore, R2: RingStore, R2::Type: CanIsoFromTo<R1::Type>
334{
335    type Isomorphism = <R2::Type as CanIsoFromTo<R1::Type>>::Isomorphism;
336
337    fn has_canonical_iso(&self, from: &SparsePolyRingBase<R1>) -> Option<Self::Isomorphism> {
338        self.base_ring().get_ring().has_canonical_iso(from.base_ring().get_ring())
339    }
340
341    fn map_out(&self, from: &SparsePolyRingBase<R1>, el: Self::Element, iso: &Self::Isomorphism) -> SparsePolyRingEl<R1> {
342        let mut result = SparseMapVector::new(el.data.len(), from.base_ring.clone());
343        for (j, c) in el.data.nontrivial_entries() {
344            *result.at_mut(j) = self.base_ring().get_ring().map_out(from.base_ring().get_ring(), self.base_ring().clone_el(c), iso);
345        }
346        return SparsePolyRingEl {
347            data: result
348        };
349    }
350}
351
352pub struct TermIterator<'a, R>
353    where R: RingStore
354{
355    iter: sparse::SparseMapVectorIter<'a, Rc<R>>
356}
357
358impl<'a, R> Iterator for TermIterator<'a, R>
359    where R: RingStore
360{
361    type Item = (&'a El<R>, usize);
362
363    fn next(&mut self) -> Option<Self::Item> {
364        if let Some((i, c)) = self.iter.next() {
365            Some((c, i))
366        } else {
367            None
368        }
369    }
370}
371
372impl<R> PolyRing for SparsePolyRingBase<R> 
373    where R: RingStore
374{
375    type TermsIterator<'a> = TermIterator<'a, R>
376        where Self: 'a;
377
378    fn indeterminate(&self) -> Self::Element {
379        let mut result = self.zero();
380        result.data.set_len(2);
381        *result.data.at_mut(1) = self.base_ring.one();
382        return result;
383    }
384
385    fn terms<'a>(&'a self, f: &'a Self::Element) -> TermIterator<'a, R> {
386        TermIterator {
387            iter: f.data.nontrivial_entries()
388        }
389    }
390
391    fn add_assign_from_terms<I>(&self, lhs: &mut Self::Element, rhs: I)
392        where I: IntoIterator<Item = (El<Self::BaseRing>, usize)>
393    {
394        for (c, i) in rhs {
395            lhs.data.set_len(max(lhs.data.len(), i + 1));
396            self.base_ring().add_assign(lhs.data.at_mut(i), c);
397        }
398        // if a previously set entry is then set to zero or adds up to zero, this might be not truncated
399        self.degree_truncate(&mut lhs.data);
400    }
401
402    fn coefficient_at<'a>(&'a self, f: &'a Self::Element, i: usize) -> &'a El<Self::BaseRing> {
403        if i < f.data.len() {
404            return f.data.at(i);
405        } else {
406            return &self.zero;
407        }
408    }
409
410    fn degree(&self, f: &Self::Element) -> Option<usize> {
411        f.data.len().checked_sub(1)
412    }
413
414    fn div_rem_monic(&self, mut lhs: Self::Element, rhs: &Self::Element) -> (Self::Element, Self::Element) {
415        assert!(self.base_ring().is_one(self.coefficient_at(rhs, self.degree(rhs).unwrap())));
416        let quo = self.poly_div(&mut lhs, rhs, |x| Some(x)).unwrap();
417        return (quo, lhs);
418    }
419}
420
421impl<R> Domain for SparsePolyRingBase<R> 
422    where R: RingStore, R::Type: Domain
423{}
424
425impl<R,> DivisibilityRing for SparsePolyRingBase<R> 
426    where R: RingStore, 
427        R::Type: DivisibilityRing + Domain
428{
429    fn checked_left_div(&self, lhs: &Self::Element, rhs: &Self::Element) -> Option<Self::Element> {
430        if let Some(d) = self.degree(rhs) {
431            let lc = rhs.data.at(d);
432            let mut lhs_copy = self.clone_el(&lhs);
433            let quo = self.poly_div(&mut lhs_copy, rhs, |x| self.base_ring().checked_left_div(&x, lc))?;
434            if self.is_zero(&lhs_copy) {
435                Some(quo)
436            } else {
437                None
438            }
439        } else if self.is_zero(lhs) {
440            Some(self.zero())
441        } else {
442            None
443        }
444    }
445}
446
447impl<R> PrincipalIdealRing for SparsePolyRingBase<R>
448    where R: RingStore, R::Type: Field + PolyTFracGCDRing
449{
450    fn checked_div_min(&self, lhs: &Self::Element, rhs: &Self::Element) -> Option<Self::Element> {
451        // base ring is a field, so everything is fine
452        if self.is_zero(rhs) && self.is_zero(lhs) {
453            return Some(self.one());
454        } else if self.is_zero(rhs) {
455            return None;
456        }
457        let (quo, rem) = self.euclidean_div_rem(self.clone_el(lhs), rhs);
458        if self.is_zero(&rem) {
459            return Some(quo);
460        } else {
461            return None;
462        }
463    }
464
465    fn extended_ideal_gen(&self, lhs: &Self::Element, rhs: &Self::Element) -> (Self::Element, Self::Element, Self::Element) {
466        algorithms::eea::eea(self.clone_el(lhs), self.clone_el(rhs), RingRef::new(self))
467    }
468    
469    fn ideal_gen(&self, lhs: &Self::Element, rhs: &Self::Element) -> Self::Element {
470        <_ as PolyTFracGCDRing>::gcd(RingRef::new(self), lhs, rhs)
471    }
472}
473
474impl<R> EuclideanRing for SparsePolyRingBase<R> 
475    where R: RingStore, R::Type: Field + PolyTFracGCDRing
476{
477    fn euclidean_div_rem(&self, mut lhs: Self::Element, rhs: &Self::Element) -> (Self::Element, Self::Element) {
478        let lc_inv = self.base_ring.invert(rhs.data.at(self.degree(rhs).unwrap())).unwrap();
479        let quo = self.poly_div(&mut lhs, rhs, |x| Some(self.base_ring().mul_ref_snd(x, &lc_inv))).unwrap();
480        return (quo, lhs);
481    }
482
483    fn euclidean_deg(&self, val: &Self::Element) -> Option<usize> {
484        return Some(self.degree(val).map(|x| x + 1).unwrap_or(0));
485    }
486}
487
488#[cfg(test)]
489use crate::rings::zn::*;
490#[cfg(test)]
491use crate::rings::zn::zn_static::{Zn, Fp};
492#[cfg(test)]
493use crate::rings::finite::FiniteRingStore;
494#[cfg(test)]
495use super::dense_poly::DensePolyRing;
496#[cfg(test)]
497use crate::primitive_int::StaticRing;
498
499#[cfg(test)]
500fn edge_case_elements<P: PolyRingStore>(poly_ring: P) -> impl Iterator<Item = El<P>>
501    where P::Type: PolyRing
502{
503    let base_ring = poly_ring.base_ring();
504    vec![ 
505        poly_ring.from_terms([].into_iter()),
506        poly_ring.from_terms([(base_ring.int_hom().map(1), 0)].into_iter()),
507        poly_ring.from_terms([(base_ring.int_hom().map(1), 1)].into_iter()),
508        poly_ring.from_terms([(base_ring.int_hom().map(1), 0), (base_ring.int_hom().map(1), 1)].into_iter()),
509        poly_ring.from_terms([(base_ring.int_hom().map(-1), 0)].into_iter()),
510        poly_ring.from_terms([(base_ring.int_hom().map(-1), 1)].into_iter()),
511        poly_ring.from_terms([(base_ring.int_hom().map(-1), 0), (base_ring.int_hom().map(1), 1)].into_iter()),
512        poly_ring.from_terms([(base_ring.int_hom().map(1), 0), (base_ring.int_hom().map(-1), 1)].into_iter()),
513        poly_ring.from_terms([(base_ring.int_hom().map(-1), 0), (base_ring.int_hom().map(1), 2)].into_iter()),
514        poly_ring.from_terms([(base_ring.int_hom().map(1), 0), (base_ring.int_hom().map(-1), 2)].into_iter()),
515        poly_ring.from_terms([(base_ring.int_hom().map(1), 0), (base_ring.int_hom().map(-1), 2), (base_ring.int_hom().map(0), 2)].into_iter())
516    ].into_iter()
517}
518
519#[test]
520fn test_ring_axioms() {
521    let poly_ring = SparsePolyRing::new(Zn::<7>::RING, "X");
522    crate::ring::generic_tests::test_ring_axioms(&poly_ring, edge_case_elements(&poly_ring));
523}
524
525#[test]
526fn test_poly_ring_axioms() {
527    let poly_ring = SparsePolyRing::new(Zn::<7>::RING, "X");
528    super::generic_tests::test_poly_ring_axioms(poly_ring, Zn::<7>::RING.elements());
529}
530
531#[test]
532fn test_canonical_iso_axioms_different_base_ring() {
533    let poly_ring1 = SparsePolyRing::new(zn_big::Zn::new(StaticRing::<i128>::RING, 7), "X");
534    let poly_ring2 = SparsePolyRing::new(zn_64::Zn::new(7), "X");
535    crate::ring::generic_tests::test_hom_axioms(&poly_ring1, &poly_ring2, edge_case_elements(&poly_ring1));
536    crate::ring::generic_tests::test_iso_axioms(&poly_ring1, &poly_ring2, edge_case_elements(&poly_ring1));
537}
538
539#[test]
540fn test_canonical_iso_dense_poly_ring() {
541    let poly_ring1 = SparsePolyRing::new(zn_64::Zn::new(7), "X");
542    let poly_ring2 = DensePolyRing::new(zn_64::Zn::new(7), "X");
543    crate::ring::generic_tests::test_hom_axioms(&poly_ring2, &poly_ring1, edge_case_elements(&poly_ring2));
544    crate::ring::generic_tests::test_iso_axioms(&poly_ring2, &poly_ring1, edge_case_elements(&poly_ring2));
545}
546
547#[test]
548fn test_divisibility_ring_axioms() {
549    let poly_ring = SparsePolyRing::new(Fp::<7>::RING, "X");
550    crate::divisibility::generic_tests::test_divisibility_axioms(&poly_ring, edge_case_elements(&poly_ring));
551}
552
553#[test]
554fn test_euclidean_ring_axioms() {
555    let poly_ring = SparsePolyRing::new(Fp::<7>::RING, "X");
556    crate::pid::generic_tests::test_euclidean_ring_axioms(&poly_ring, edge_case_elements(&poly_ring));
557}
558
559#[test]
560fn test_principal_ideal_ring_axioms() {
561    let poly_ring = SparsePolyRing::new(Fp::<7>::RING, "X");
562    crate::pid::generic_tests::test_principal_ideal_ring_axioms(&poly_ring, edge_case_elements(&poly_ring));
563}
feanor_math/rings/poly/sparse_poly.rs

feanor_math/rings/poly/
sparse_poly.rs
