feanor_math/rings/multivariate/
mod.rs

1use crate::homomorphism::Homomorphism;
2use crate::ring::*;
3use crate::seq::VectorFn;
4use crate::wrapper::RingElementWrapper;
5
6///
7/// Contains an implementation [`multivariate_impl::MultivariatePolyRingImpl`] of a multivariate polynomial
8/// ring.
9/// 
10pub mod multivariate_impl;
11
12use std::any::Any;
13use std::cmp::{max, Ordering};
14
15///
16/// Type of coefficients of multivariate polynomials of the given ring.
17/// 
18pub type PolyCoeff<P> = El<<<P as RingStore>::Type as RingExtension>::BaseRing>;
19
20///
21/// Type of monomials of multivariate polynomials of the given ring.
22/// 
23pub type PolyMonomial<P> = <<P as RingStore>::Type as MultivariatePolyRing>::Monomial;
24
25///
26/// Trait for multivariate polynomial rings.
27/// 
28pub trait MultivariatePolyRing: RingExtension {
29
30    type Monomial;
31    type TermIter<'a>: Iterator<Item = (&'a El<Self::BaseRing>, &'a Self::Monomial)>
32        where Self: 'a;
33
34    ///
35    /// Returns the number of variables of this polynomial ring, i.e. the transcendence degree
36    /// of the base ring.
37    /// 
38    fn indeterminate_count(&self) -> usize;
39
40    ///
41    /// Returns the monomial `Xi`, where `Xi` is the `i`-th generator of this ring.
42    /// 
43    fn indeterminate(&self, i: usize) -> Self::Monomial {
44        assert!(i < self.indeterminate_count());
45        self.create_monomial((0..self.indeterminate_count()).map(|j| if i == j { 1 } else { 0 }))
46    }
47
48    ///
49    /// Creates a monomial with the given exponents.
50    /// 
51    /// Note that when building a polynomial, the most convenient method is usually
52    /// to use [`MultivariatePolyRingStore::with_wrapped_indeterminates()`].
53    /// 
54    /// # Example
55    /// ```rust
56    /// # use feanor_math::ring::*;
57    /// # use feanor_math::primitive_int::*;
58    /// # use feanor_math::rings::multivariate::*;
59    /// # use feanor_math::rings::multivariate::multivariate_impl::*;
60    /// let poly_ring = MultivariatePolyRingImpl::new(StaticRing::<i64>::RING, 3);
61    /// let x_as_monomial = poly_ring.create_monomial([1, 0, 0]);
62    /// let x_as_poly = poly_ring.create_term(1, x_as_monomial);
63    /// assert_eq!("X0", format!("{}", poly_ring.format(&x_as_poly)));
64    /// ```
65    /// 
66    fn create_monomial<I>(&self, exponents: I) -> Self::Monomial
67        where I: IntoIterator<Item = usize>,
68            I::IntoIter: ExactSizeIterator;
69
70    ///
71    /// Multiplies the given polynomial with the given monomial.
72    /// 
73    fn mul_assign_monomial(&self, f: &mut Self::Element, monomial: Self::Monomial);
74
75    ///
76    /// Returns the coefficient corresponding to the given monomial in the given polynomial.
77    /// If the polynomial does not contain a term with that monomial, zero is returned.
78    /// 
79    fn coefficient_at<'a>(&'a self, f: &'a Self::Element, m: &Self::Monomial) -> &'a El<Self::BaseRing>;
80
81    ///
82    /// Returns the power of the `var_index`-th variable in the given monomial.
83    /// In other words, this maps `X1^i1 ... Xm^im` to `i(var_index)`.
84    /// 
85    fn exponent_at(&self, m: &Self::Monomial, var_index: usize) -> usize;
86
87    ///
88    /// Writes the powers of each variable in the given monomial into the given 
89    /// output slice. 
90    /// 
91    /// This is equivalent to performing `out[i] = self.exponent_at(m, i)` for
92    /// every `i` in `0..self.indeterminate_count()`.
93    /// 
94    fn expand_monomial_to(&self, m: &Self::Monomial, out: &mut [usize]) {
95        assert_eq!(out.len(), self.indeterminate_count());
96        for i in 0..self.indeterminate_count() {
97            out[i] = self.exponent_at(m, i)
98        }
99    }
100
101    ///
102    /// Returns an iterator over all nonzero terms of the given polynomial.
103    /// 
104    fn terms<'a>(&'a self, f: &'a Self::Element) -> Self::TermIter<'a>;
105
106    ///
107    /// Creates a new single-term polynomial.
108    /// 
109    fn create_term(&self, coeff: El<Self::BaseRing>, monomial: Self::Monomial) -> Self::Element {
110        let mut result = self.from(coeff);
111        self.mul_assign_monomial(&mut result, monomial);
112        return result;
113    }
114
115    ///
116    /// Returns the **L**eading **T**erm of `f`, i.e. the term whose monomial is largest w.r.t. the given order.
117    /// 
118    fn LT<'a, O: MonomialOrder>(&'a self, f: &'a Self::Element, order: O) -> Option<(&'a El<Self::BaseRing>, &'a Self::Monomial)> {
119        self.terms(f).max_by(|l, r| order.compare(RingRef::new(self), &l.1, &r.1))
120    }
121
122    ///
123    /// Returns the term of `f` whose monomial is largest (w.r.t. the given order) among all monomials smaller than `lt_than`.
124    /// 
125    fn largest_term_lt<'a, O: MonomialOrder>(&'a self, f: &'a Self::Element, order: O, lt_than: &Self::Monomial) -> Option<(&'a El<Self::BaseRing>, &'a Self::Monomial)> {
126        self.terms(f).filter(|(_, m)| order.compare(RingRef::new(self), m, lt_than) == Ordering::Less).max_by(|l, r| order.compare(RingRef::new(self), &l.1, &r.1))
127    }
128
129    fn add_assign_from_terms<I>(&self, lhs: &mut Self::Element, rhs: I)
130        where I: IntoIterator<Item = (El<Self::BaseRing>, Self::Monomial)>
131    {
132        let self_ring = RingRef::new(self);
133        self.add_assign(lhs, self_ring.sum(
134            rhs.into_iter().map(|(c, m)| self.create_term(c, m))
135        ));
136    }
137
138    ///
139    /// Applies the given homomorphism `R -> S` to each coefficient of the given polynomial
140    /// in `R[X1, ..., Xm]` to produce a monomial in `S[X1, ..., Xm]`.
141    /// 
142    fn map_terms<P, H>(&self, from: &P, el: &P::Element, hom: H) -> Self::Element
143        where P: ?Sized + MultivariatePolyRing,
144            H: Homomorphism<<P::BaseRing as RingStore>::Type, <Self::BaseRing as RingStore>::Type>
145    {
146        assert!(self.base_ring().get_ring() == hom.codomain().get_ring());
147        assert!(from.base_ring().get_ring() == hom.domain().get_ring());
148        assert_eq!(self.indeterminate_count(), from.indeterminate_count());
149        let mut exponents_storage = (0..self.indeterminate_count()).map(|_| 0).collect::<Vec<_>>();
150        return RingRef::new(self).from_terms(from.terms(el).map(|(c, m)| {
151            from.expand_monomial_to(m, &mut exponents_storage);
152            (hom.map_ref(c), self.create_monomial(exponents_storage.iter().copied()))
153        }));
154    }
155
156    fn clone_monomial(&self, mon: &Self::Monomial) -> Self::Monomial {
157        self.create_monomial((0..self.indeterminate_count()).map(|i| self.exponent_at(mon, i)))
158    }
159
160    ///
161    /// Returns a list of all variables appearing in the given polynomial.
162    /// Associated with each variable is the highest degree in which it
163    /// appears in some term.
164    /// 
165    /// # Example
166    /// ```rust
167    /// # use feanor_math::primitive_int::*;
168    /// # use feanor_math::ring::*;
169    /// # use feanor_math::rings::multivariate::*;
170    /// # use feanor_math::rings::multivariate::multivariate_impl::*;
171    /// let poly_ring = MultivariatePolyRingImpl::new(StaticRing::<i64>::RING, 2);
172    /// let [f, g] = poly_ring.with_wrapped_indeterminates(|[X, Y]| [1 + X + X.pow_ref(2) * Y, X.pow_ref(3)]);
173    /// assert_eq!(vec![(0, 2), (1, 1)], poly_ring.appearing_indeterminates(&f));
174    /// assert_eq!(vec![(0, 3)], poly_ring.appearing_indeterminates(&g));
175    /// ```
176    /// 
177    fn appearing_indeterminates(&self, f: &Self::Element) -> Vec<(usize, usize)> {
178        let mut result = (0..self.indeterminate_count()).map(|_| 0).collect::<Vec<_>>();
179        for (_, m) in self.terms(f) {
180            for i in 0..self.indeterminate_count() {
181                result[i] = max(result[i], self.exponent_at(m, i));
182            }
183        }
184        return result.into_iter().enumerate().filter(|(_, e)| *e > 0).collect();
185    }
186
187    ///
188    /// Multiplies two monomials.
189    /// 
190    fn monomial_mul(&self, lhs: Self::Monomial, rhs: &Self::Monomial) -> Self::Monomial {
191        self.create_monomial((0..self.indeterminate_count()).map(|i| self.exponent_at(&lhs, i) + self.exponent_at(rhs, i)))
192    }
193
194    ///
195    /// Returns the degree of a monomial, i.e. the sum of the exponents of all variables.
196    /// 
197    fn monomial_deg(&self, mon: &Self::Monomial) -> usize {
198        (0..self.indeterminate_count()).map(|i| self.exponent_at(mon, i)).sum()
199    }
200
201    ///
202    /// Returns the least common multiple of two monomials.
203    /// 
204    fn monomial_lcm(&self, lhs: Self::Monomial, rhs: &Self::Monomial) -> Self::Monomial {
205        self.create_monomial((0..self.indeterminate_count()).map(|i| max(self.exponent_at(&lhs, i), self.exponent_at(rhs, i))))
206    }
207
208    ///
209    /// Computes the quotient of two monomials.
210    /// 
211    /// If `lhs` does not divide `rhs`, this returns `Result::Err` with the monomial
212    /// `lhs / gcd(rhs, lhs)`.
213    /// 
214    fn monomial_div(&self, lhs: Self::Monomial, rhs: &Self::Monomial) -> Result<Self::Monomial, Self::Monomial> {
215        let mut failed = false;
216        let result = self.create_monomial((0..self.indeterminate_count()).map(|i| {
217            if let Some(res) = self.exponent_at(&lhs, i).checked_sub(self.exponent_at(rhs, i)) {
218                res
219            } else {
220                failed = true;
221                0
222            }
223        }));
224        if failed {
225            Err(result)
226        } else {
227            Ok(result)
228        }
229    }
230
231    ///
232    /// Evaluates the given polynomial at the given values.
233    /// 
234    /// # Example
235    /// ```rust
236    /// # use feanor_math::ring::*;
237    /// # use feanor_math::rings::multivariate::*;
238    /// # use feanor_math::rings::multivariate::multivariate_impl::*;
239    /// # use feanor_math::primitive_int::*; 
240    /// # use feanor_math::seq::*; 
241    /// let poly_ring = MultivariatePolyRingImpl::new(StaticRing::<i64>::RING, 2);
242    /// let [f] = poly_ring.with_wrapped_indeterminates(|[X, Y]| [1 + X + X.pow_ref(2) * Y]);
243    /// assert_eq!(1 + 5 + 5 * 5 * 8, poly_ring.evaluate(&f, [5, 8].clone_ring_els(StaticRing::<i64>::RING), &poly_ring.base_ring().identity()));
244    /// ```
245    /// 
246    fn evaluate<R, V, H>(&self, f: &Self::Element, value: V, hom: H) -> R::Element
247        where R: ?Sized + RingBase,
248            H: Homomorphism<<Self::BaseRing as RingStore>::Type, R>,
249            V: VectorFn<R::Element>
250    {
251        assert_eq!(self.indeterminate_count(), value.len());
252        assert!(hom.domain().get_ring() == self.base_ring().get_ring());
253        hom.codomain().sum(self.terms(f).map(|(c, m)| hom.mul_map(
254            hom.codomain().prod((0..self.indeterminate_count()).map(|i| hom.codomain().pow(value.at(i), self.exponent_at(m, i)))),
255            hom.domain().clone_el(c)
256        )))
257    }
258
259    ///
260    /// Replaces the given indeterminate in the given polynomial by the value `val`.
261    /// 
262    /// Conceptually, this is similar to [`MultivariatePolyRing::evaluate()`], but less general,
263    /// which can allow a faster implementation sometimes. In particular, this only replaces a single
264    /// indeterminate, and does not change the ring.
265    /// 
266    fn specialize(&self, f: &Self::Element, var: usize, val: &Self::Element) -> Self::Element {
267        assert!(var < self.indeterminate_count());
268        let mut parts = Vec::new();
269        for (c, m) in self.terms(f) {
270            while self.exponent_at(m, var) >= parts.len() {
271                parts.push(Vec::new());
272            }
273            let new_m = self.create_monomial((0..self.indeterminate_count()).map(|i| if i == var { 0 } else { self.exponent_at(m, i) }));
274            parts[self.exponent_at(m, var)].push((self.base_ring().clone_el(c), new_m));
275        }
276        if let Some(first) = parts.pop() {
277            let mut current = self.zero();
278            self.add_assign_from_terms(&mut current, first);
279            while let Some(new) = parts.pop() {
280                let mut next = self.zero();
281                self.add_assign_from_terms(&mut next, new);
282                self.mul_assign_ref(&mut current, val);
283                self.add_assign(&mut current, next);
284            }
285            return current;
286        } else {
287            return self.zero();
288        }
289    }
290}
291
292///
293/// [`RingStore`] for [`MultivariatePolyRing`]
294/// 
295pub trait MultivariatePolyRingStore: RingStore
296    where Self::Type: MultivariatePolyRing
297{
298    delegate!{ MultivariatePolyRing, fn indeterminate_count(&self) -> usize }
299    delegate!{ MultivariatePolyRing, fn indeterminate(&self, i: usize) -> <Self::Type as MultivariatePolyRing>::Monomial }
300    delegate!{ MultivariatePolyRing, fn create_term(&self, coeff: PolyCoeff<Self>, monomial: PolyMonomial<Self>) -> El<Self> }
301    delegate!{ MultivariatePolyRing, fn exponent_at(&self, m: &PolyMonomial<Self>, var_index: usize) -> usize }
302    delegate!{ MultivariatePolyRing, fn expand_monomial_to(&self, m: &PolyMonomial<Self>, out: &mut [usize]) -> () }
303    delegate!{ MultivariatePolyRing, fn monomial_mul(&self, lhs: PolyMonomial<Self>, rhs: &PolyMonomial<Self>) -> PolyMonomial<Self> }
304    delegate!{ MultivariatePolyRing, fn monomial_lcm(&self, lhs: PolyMonomial<Self>, rhs: &PolyMonomial<Self>) -> PolyMonomial<Self> }
305    delegate!{ MultivariatePolyRing, fn monomial_div(&self, lhs: PolyMonomial<Self>, rhs: &PolyMonomial<Self>) -> Result<PolyMonomial<Self>, PolyMonomial<Self>> }
306    delegate!{ MultivariatePolyRing, fn monomial_deg(&self, val: &PolyMonomial<Self>) -> usize }
307    delegate!{ MultivariatePolyRing, fn mul_assign_monomial(&self, f: &mut El<Self>, monomial: PolyMonomial<Self>) -> () }
308    delegate!{ MultivariatePolyRing, fn appearing_indeterminates(&self, f: &El<Self>) -> Vec<(usize, usize)> }
309    delegate!{ MultivariatePolyRing, fn specialize(&self, f: &El<Self>, var: usize, val: &El<Self>) -> El<Self> }
310
311    fn expand_monomial(&self, m: &PolyMonomial<Self>) -> Vec<usize> {
312        let mut result = (0..self.indeterminate_count()).map(|_| 0).collect::<Vec<_>>();
313        self.expand_monomial_to(m, &mut result);
314        return result;
315    }
316
317    ///
318    /// Returns the term of `f` whose monomial is largest (w.r.t. the given order) among all monomials smaller than `lt_than`.
319    /// 
320    fn largest_term_lt<'a, O: MonomialOrder>(&'a self, f: &'a El<Self>, order: O, lt_than: &PolyMonomial<Self>) -> Option<(&'a PolyCoeff<Self>, &'a PolyMonomial<Self>)> {
321        self.get_ring().largest_term_lt(f, order, lt_than)
322    }
323    
324    ///
325    /// Returns the **L**eading **T**erm of `f`, i.e. the term whose monomial is largest w.r.t. the given order.
326    /// 
327    fn LT<'a, O: MonomialOrder>(&'a self, f: &'a El<Self>, order: O) -> Option<(&'a PolyCoeff<Self>, &'a PolyMonomial<Self>)> {
328        self.get_ring().LT(f, order)
329    }
330
331    ///
332    /// Creates a new monomial with given exponents.
333    /// 
334    /// For details, see [`MultivariatePolyRing::create_monomial()`].
335    /// 
336    fn create_monomial<I>(&self, exponents: I) -> PolyMonomial<Self>
337        where I: IntoIterator<Item = usize>,
338            I::IntoIter: ExactSizeIterator
339    {
340        self.get_ring().create_monomial(exponents)
341    }
342
343    fn clone_monomial(&self, mon: &PolyMonomial<Self>) -> PolyMonomial<Self> {
344        self.get_ring().clone_monomial(mon)
345    }
346
347    ///
348    /// Returns the coefficient of a polynomial corresponding to a monomial.
349    /// 
350    /// For details, see [`MultivariatePolyRing::coefficient_at()`].
351    /// 
352    fn coefficient_at<'a>(&'a self, f: &'a El<Self>, m: &PolyMonomial<Self>) -> &'a PolyCoeff<Self> {
353        self.get_ring().coefficient_at(f, m)
354    }
355
356    ///
357    /// Returns an iterator over all nonzero terms of the polynomial.
358    /// 
359    /// For details, see [`MultivariatePolyRing::terms()`].
360    /// 
361    fn terms<'a>(&'a self, f: &'a El<Self>) -> <Self::Type as MultivariatePolyRing>::TermIter<'a> {
362        self.get_ring().terms(f)
363    }
364
365    ///
366    /// Creates a new polynomial by summing up all the given terms.
367    /// 
368    fn from_terms<I>(&self, terms: I) -> El<Self>
369        where I: IntoIterator<Item = (PolyCoeff<Self>, PolyMonomial<Self>)>
370    {
371        let mut result = self.zero();
372        self.get_ring().add_assign_from_terms(&mut result, terms);
373        return result;
374    }
375
376    ///
377    /// Evaluates the polynomial at the given values.
378    /// 
379    /// For details, see [`MultivariatePolyRing::evaluate()`].
380    /// 
381    fn evaluate<R, V, H>(&self, f: &El<Self>, value: V, hom: H) -> R::Element
382        where R: ?Sized + RingBase,
383            H: Homomorphism<<<Self::Type as RingExtension>::BaseRing as RingStore>::Type, R>,
384            V: VectorFn<R::Element>
385    {
386        self.get_ring().evaluate(f, value, hom)
387    }
388    
389    ///
390    /// Returns the homomorphism `R[X1, ..., Xm] -> S[X1, ..., Xm]` that is induced by
391    /// applying the given homomorphism `R -> S` coefficient-wise.
392    /// 
393    fn into_lifted_hom<P, H>(self, from: P, hom: H) -> CoefficientHom<P, Self, H>
394        where P: RingStore,
395            P::Type: MultivariatePolyRing,
396            H: Homomorphism<<<P::Type as RingExtension>::BaseRing as RingStore>::Type, <<Self::Type as RingExtension>::BaseRing as RingStore>::Type>
397    {
398        CoefficientHom {
399            from: from,
400            to: self,
401            hom: hom
402        }
403    }
404
405    ///
406    /// Returns the homomorphism `R[X1, ..., Xm] -> S[X1, ..., Xm]` that is induced by
407    /// applying the given homomorphism `R -> S` coefficient-wise.
408    /// 
409    /// If the ownership of this ring should be transferred to the homomorphism, consider
410    /// using [`MultivariatePolyRingStore::into_lifted_hom()`].
411    /// 
412    fn lifted_hom<'a, P, H>(&'a self, from: P, hom: H) -> CoefficientHom<P, &'a Self, H>
413        where P: RingStore,
414            P::Type: MultivariatePolyRing,
415            H: Homomorphism<<<P::Type as RingExtension>::BaseRing as RingStore>::Type, <<Self::Type as RingExtension>::BaseRing as RingStore>::Type>
416    {
417        self.into_lifted_hom(from, hom)
418    }
419
420    ///
421    /// Invokes the function with a wrapped version of the indeterminates of this poly ring.
422    /// Use for convenient creation of polynomials.
423    /// 
424    /// Note however that [`MultivariatePolyRingStore::from_terms()`] might be more performant.
425    /// 
426    /// # Example
427    /// ```rust
428    /// use feanor_math::assert_el_eq;
429    /// use feanor_math::homomorphism::*;
430    /// use feanor_math::ring::*;
431    /// use feanor_math::rings::multivariate::*;
432    /// use feanor_math::rings::zn::zn_64::*;
433    /// use feanor_math::rings::multivariate::multivariate_impl::*;
434    /// let base_ring = Zn::new(7);
435    /// let poly_ring = MultivariatePolyRingImpl::new(base_ring, 3);
436    /// let f_version1 = poly_ring.from_terms([(base_ring.int_hom().map(3), poly_ring.create_monomial([0, 0, 0])), (base_ring.int_hom().map(2), poly_ring.create_monomial([0, 1, 1])), (base_ring.one(), poly_ring.create_monomial([2, 0, 0]))].into_iter());
437    /// let f_version2 = poly_ring.with_wrapped_indeterminates_dyn(|[x, y, z]| [3 + 2 * y * z + x.pow_ref(2)]).into_iter().next().unwrap();
438    /// let [f_version3] = poly_ring.with_wrapped_indeterminates(|[x, y, z]| [3 + 2 * y * z + x.pow_ref(2)]);
439    /// assert_el_eq!(poly_ring, f_version1, f_version2);
440    /// ```
441    /// 
442    #[stability::unstable(feature = "enable")]
443    fn with_wrapped_indeterminates_dyn<'a, F, T, const N: usize>(&'a self, f: F) -> Vec<El<Self>>
444        where F: FnOnce([&RingElementWrapper<&'a Self>; N]) -> T,
445            T: IntoIterator<Item = RingElementWrapper<&'a Self>>
446    {
447        assert_eq!(self.indeterminate_count(), N);
448        let wrapped_indets: [_; N] = std::array::from_fn(|i| RingElementWrapper::new(self, self.create_term(self.base_ring().one(), self.indeterminate(i))));
449        f(std::array::from_fn(|i| &wrapped_indets[i])).into_iter().map(|f| f.unwrap()).collect()
450    }
451
452    ///
453    /// Same as [`MultivariatePolyRingStore::with_wrapped_indeterminates_dyn()`], but returns result
454    /// as an array to allow pattern matching.
455    /// 
456    /// # Example
457    /// ```rust
458    /// use feanor_math::assert_el_eq;
459    /// use feanor_math::homomorphism::*;
460    /// use feanor_math::ring::*;
461    /// use feanor_math::rings::multivariate::*;
462    /// use feanor_math::rings::zn::zn_64::*;
463    /// use feanor_math::rings::multivariate::multivariate_impl::*;
464    /// let base_ring = Zn::new(7);
465    /// let poly_ring = MultivariatePolyRingImpl::new(base_ring, 3);
466    /// let f_version1 = poly_ring.from_terms([(base_ring.int_hom().map(3), poly_ring.create_monomial([0, 0, 0])), (base_ring.int_hom().map(2), poly_ring.create_monomial([0, 1, 1])), (base_ring.one(), poly_ring.create_monomial([2, 0, 0]))].into_iter());
467    /// let f_version2 = poly_ring.with_wrapped_indeterminates_dyn(|[x, y, z]| [3 + 2 * y * z + x.pow_ref(2)]).into_iter().next().unwrap();
468    /// let [f_version3] = poly_ring.with_wrapped_indeterminates(|[x, y, z]| [3 + 2 * y * z + x.pow_ref(2)]);
469    /// assert_el_eq!(poly_ring, f_version1, f_version2);
470    /// ```
471    /// 
472    #[stability::unstable(feature = "enable")]
473    fn with_wrapped_indeterminates<'a, F, const N: usize, const M: usize>(&'a self, f: F) -> [El<Self>; M]
474        where F: FnOnce([&RingElementWrapper<&'a Self>; N]) -> [RingElementWrapper<&'a Self>; M]
475    {
476        assert_eq!(self.indeterminate_count(), N);
477        let wrapped_indets: [_; N] = std::array::from_fn(|i| RingElementWrapper::new(self, self.create_term(self.base_ring().one(), self.indeterminate(i))));
478        let mut result_it = f(std::array::from_fn(|i| &wrapped_indets[i])).into_iter().map(|f| f.unwrap());
479        let result = std::array::from_fn(|_| result_it.next().unwrap());
480        debug_assert!(result_it.next().is_none());
481        return result;
482    }
483}
484
485impl<P> MultivariatePolyRingStore for P
486    where P: RingStore,
487        P::Type: MultivariatePolyRing {}
488
489///
490/// Trait for monomial orders, i.e. orderings on the monomials `X1^i1 ... Xm^im` of a multivariate
491/// polynomial ring that are compatible with multiplication. These are also sometimes called term
492/// orders.
493/// 
494/// To be precise, a total order on monomials is a monomial order, if for all `u, v, w` with `u <= v`
495/// also `uw <= vw`.
496/// 
497pub trait MonomialOrder: Clone {
498
499    fn compare<P>(&self, ring: P, lhs: &PolyMonomial<P>, rhs: &PolyMonomial<P>) -> Ordering
500        where P: RingStore,
501            P::Type: MultivariatePolyRing;
502
503    ///
504    /// Checks whether two monomials are equal.
505    /// 
506    /// This may be faster than [`MonomialOrder::compare()`], but clearly must have
507    /// a compatible behavior.
508    /// 
509    fn eq_mon<P>(&self, ring: P, lhs: &PolyMonomial<P>, rhs: &PolyMonomial<P>) -> bool
510        where P: RingStore,
511            P::Type: MultivariatePolyRing
512    {
513        self.compare(ring, lhs, rhs) == Ordering::Equal
514    }
515
516    ///
517    /// Whether this order is the same as the given other order, i.e. [`MonomialOrder::compare()`]
518    /// gives the same output on all inputs.
519    /// 
520    /// Many monomial orders are likely to be implemented as zero-sized types with only a single
521    /// instance. In this case, the default implementation is sufficient.
522    /// 
523    fn is_same<O>(&self, rhs: &O) -> bool
524        where O: MonomialOrder
525    {
526        assert!(self.as_any().is_some());
527        assert!(std::mem::size_of::<Self>() == 0);
528        if let Some(rhs_as_any) = rhs.as_any() {
529            self.as_any().unwrap().type_id() == rhs_as_any.type_id()
530        } else {
531            false
532        }
533    }
534
535    ///
536    /// Upcasts this reference to `&dyn Any`, which is sometimes required to compare monomial order
537    /// objects of different types.
538    /// 
539    fn as_any(&self) -> Option<&dyn Any>;
540}
541
542///
543/// Trait for [`MonomialOrder`]s that are graded, i.e. for `v, w` with `deg(v) < deg(w)`
544/// they always satisfy `v < w`.
545/// 
546pub trait GradedMonomialOrder: MonomialOrder {}
547
548///
549/// The graded reverse lexicographic order. This is the most important monomial order, since
550/// computing a Groebner basis w.r.t. this order is usually more efficient than for other orders.
551/// Also sometimes referred to as "grevlex".
552/// 
553/// Concretely, this is the ordering of monomials we get by first comparing monomial degrees, and
554/// in case of a tie reverse the outcome of a lexicographic comparison, using a reversed order
555/// of variables.
556/// 
557/// # Example
558/// ```rust
559/// # use feanor_math::ring::*;
560/// # use feanor_math::rings::multivariate::*;
561/// # use feanor_math::rings::multivariate::multivariate_impl::*;
562/// # use feanor_math::primitive_int::*; 
563/// # use std::cmp::Ordering;
564/// let poly_ring = MultivariatePolyRingImpl::new(StaticRing::<i64>::RING, 3);
565/// let monomials_descending = [
566///     [2, 0, 0], // x1^2
567///     [1, 1, 0], // x1 x2
568///     [0, 2, 0], // x2^2
569///     [1, 0, 1], // x1 x3
570///     [0, 1, 1], // x2 x3
571///     [0, 0, 2], // x3^2
572/// ].into_iter().map(|m| poly_ring.create_monomial(m)).collect::<Vec<_>>();
573/// for i in 1..6 {
574///     assert!(DegRevLex.compare(&poly_ring, &monomials_descending[i - 1], &monomials_descending[i]) == Ordering::Greater);
575/// }
576/// ```
577/// 
578#[derive(Clone, Copy, PartialEq, Eq, Debug)]
579pub struct DegRevLex;
580
581impl MonomialOrder for DegRevLex {
582
583    fn as_any(&self) -> Option<&dyn Any> {
584        Some(self as &dyn Any)
585    }
586
587    fn compare<P>(&self, ring: P, lhs: &PolyMonomial<P>, rhs: &PolyMonomial<P>) -> Ordering
588        where P: RingStore,
589            P::Type:MultivariatePolyRing
590    {
591        let lhs_deg = ring.monomial_deg(lhs);
592        let rhs_deg = ring.monomial_deg(rhs);
593        if lhs_deg < rhs_deg {
594            return Ordering::Less;
595        } else if lhs_deg > rhs_deg {
596            return Ordering::Greater;
597        } else {
598            for i in (0..ring.indeterminate_count()).rev() {
599                if ring.exponent_at(lhs, i) > ring.exponent_at(rhs, i) {
600                    return Ordering::Less
601                } else if ring.exponent_at(lhs, i) < ring.exponent_at(rhs, i) {
602                    return Ordering::Greater;
603                }
604            }
605            return Ordering::Equal;
606        }
607    }
608}
609
610impl GradedMonomialOrder for DegRevLex {}
611
612///
613/// Lexicographic ordering of monomials.
614/// 
615#[derive(Clone, Copy, PartialEq, Eq, Debug)]
616pub struct Lex;
617
618impl MonomialOrder for Lex {
619
620    fn as_any(&self) -> Option<&dyn Any> {
621        Some(self as &dyn Any)
622    }
623
624    fn compare<P>(&self, ring: P, lhs: &PolyMonomial<P>, rhs: &PolyMonomial<P>) -> Ordering
625        where P:RingStore,
626            P::Type:MultivariatePolyRing
627    {
628        for i in 0..ring.indeterminate_count() {
629            match ring.exponent_at(lhs, i).cmp(&ring.exponent_at(rhs, i)) {
630                Ordering::Less => { return Ordering::Less; },
631                Ordering::Greater => { return Ordering::Greater; },
632                Ordering::Equal => {}
633            }
634        }
635        return Ordering::Equal;
636    }
637}
638
639pub struct CoefficientHom<PFrom, PTo, H>
640    where PFrom: RingStore,
641        PTo: RingStore,
642        PFrom::Type: MultivariatePolyRing,
643        PTo::Type: MultivariatePolyRing,
644        H: Homomorphism<<<PFrom::Type as RingExtension>::BaseRing as RingStore>::Type, <<PTo::Type as RingExtension>::BaseRing as RingStore>::Type>
645{
646    from: PFrom,
647    to: PTo,
648    hom: H
649}
650
651impl<PFrom, PTo, H> Homomorphism<PFrom::Type, PTo::Type> for CoefficientHom<PFrom, PTo, H>
652    where PFrom: RingStore,
653        PTo: RingStore,
654        PFrom::Type: MultivariatePolyRing,
655        PTo::Type: MultivariatePolyRing,
656        H: Homomorphism<<<PFrom::Type as RingExtension>::BaseRing as RingStore>::Type, <<PTo::Type as RingExtension>::BaseRing as RingStore>::Type>
657{
658    type DomainStore = PFrom;
659    type CodomainStore = PTo;
660
661    fn codomain<'a>(&'a self) -> &'a Self::CodomainStore {
662        &self.to
663    }
664
665    fn domain<'a>(&'a self) -> &'a Self::DomainStore {
666        &self.from
667    }
668
669    fn map(&self, x: <PFrom::Type as RingBase>::Element) -> <PTo::Type as RingBase>::Element {
670        self.map_ref(&x)
671    }
672
673    fn map_ref(&self, x: &<PFrom::Type as RingBase>::Element) -> <PTo::Type as RingBase>::Element {
674        self.to.get_ring().map_terms(self.from.get_ring(), x, &self.hom)
675    }
676}
677
678pub mod generic_impls {
679    use std::fmt::{Formatter, Result};
680    use super::*;
681
682    #[stability::unstable(feature = "enable")]
683    pub fn print<P>(ring: P, poly: &El<P>, out: &mut Formatter, env: EnvBindingStrength) -> Result
684        where P: RingStore,
685            P::Type: MultivariatePolyRing
686    {
687        if ring.is_zero(poly) {
688            ring.base_ring().get_ring().dbg_within(&ring.base_ring().zero(), out, env)?;
689        } else {
690            if env >= EnvBindingStrength::Product {
691                write!(out, "(")?;
692            }
693            
694            let mut print_term = |c: &PolyCoeff<P>, m: &PolyMonomial<P>, with_plus: bool| {
695                if with_plus {
696                    write!(out, " + ")?;
697                }
698                let is_constant_term = ring.monomial_deg(m) == 0;
699                if !ring.base_ring().is_one(c) || is_constant_term {
700                    ring.base_ring().get_ring().dbg_within(c, out, if is_constant_term { EnvBindingStrength::Sum } else { EnvBindingStrength::Product })?;
701                    if !is_constant_term {
702                        write!(out, " * ")?;
703                    }
704                }
705                let mut needs_space = false;
706                for i in 0..ring.indeterminate_count() {
707                    if ring.exponent_at(m, i) > 0 {
708                        if needs_space {
709                            write!(out, " * ")?;
710                        }
711                        write!(out, "X{}", i)?;
712                        needs_space = true;
713                    }
714                    if ring.exponent_at(m, i) > 1 {
715                        write!(out, "^{}", ring.exponent_at(m, i))?;
716                    }
717                }
718                return Ok::<(), std::fmt::Error>(());
719            };
720            
721            for (i, (c, m)) in ring.terms(poly).enumerate() {
722                print_term(c, m, i != 0)?;
723            }
724            if env >= EnvBindingStrength::Product {
725                write!(out, ")")?;
726            }
727        }
728
729        return Ok(());
730    }
731}
732
733#[cfg(any(test, feature = "generic_tests"))]
734pub mod generic_tests {
735    use super::*;
736    use crate::seq::*;
737
738    #[stability::unstable(feature = "enable")]
739    pub fn test_poly_ring_axioms<P: RingStore, I: Iterator<Item = PolyCoeff<P>>>(ring: P, interesting_base_ring_elements: I)
740        where P::Type: MultivariatePolyRing
741    {
742        let elements = interesting_base_ring_elements.collect::<Vec<_>>();
743        let n = ring.indeterminate_count();
744        let base_ring = ring.base_ring();
745
746        // test multiplication of variables
747        for i in 0..n {
748            for j in 0..n {
749                let xi = ring.create_term(base_ring.one(), ring.indeterminate(i));
750                let xj = ring.create_term(base_ring.one(), ring.indeterminate(j));
751                let xixj = ring.create_term(base_ring.one(), ring.create_monomial((0..n).map(|k| if k == i && k == j { 2 } else if k == j || k == i { 1 } else { 0 })));
752                assert_el_eq!(ring, xixj, ring.mul(xi, xj));
753            }
754        }
755
756        // test monomial operations
757        for i in 0..n {
758            for j in 0..n {
759                let xi = ring.indeterminate(i);
760                let xj = ring.indeterminate(j);
761                let xixj_lcm = ring.create_monomial((0..n).map(|k| if k == j || k == i { 1 } else { 0 }));
762                assert_el_eq!(ring, ring.create_term(base_ring.one(), xixj_lcm), ring.create_term(base_ring.one(), ring.monomial_lcm(xi, &xj)));
763            }
764        }
765
766        // all monomials should be different
767        for i in 0..n {
768            for j in 0..n {
769                let xi = ring.create_term(base_ring.one(), ring.create_monomial((0..n).map(|k| if k == i { 1 } else { 0 })));
770                let xj = ring.create_term(base_ring.one(), ring.create_monomial((0..n).map(|k| if k == j { 1 } else { 0 })));
771                assert!((i == j) == ring.eq_el(&xi, &xj));
772            }
773        }
774
775        // monomials shouldn't be zero divisors
776        for i in 0..n {
777            for a in &elements {
778                let xi = ring.create_term(base_ring.one(), ring.indeterminate(i));
779                assert!(base_ring.is_zero(a) == ring.is_zero(&ring.inclusion().mul_ref_map(&xi, a)));
780            }
781        }
782
783        // test add_assign_from_terms
784        for i in 0..n {
785            let xi = ring.indeterminate(i);
786            let mut a = ring.create_term(base_ring.int_hom().map(3), ring.create_monomial((0..n).map(|_| 0)));
787            let terms_with_multiples = [
788                (base_ring.one(), ring.clone_monomial(&xi)),
789                (base_ring.one(), ring.clone_monomial(&xi)),
790                (base_ring.one(), ring.create_monomial((0..n).map(|_| 0))),
791                (base_ring.one(), ring.create_monomial((0..n).map(|_| 0))),
792                (base_ring.one(), ring.clone_monomial(&xi)),
793                (base_ring.one(), ring.create_monomial((0..n).map(|_| 0))),
794                (base_ring.one(), ring.clone_monomial(&xi)),
795                (base_ring.one(), ring.create_monomial((0..n).map(|_| 0))),
796            ];
797            ring.get_ring().add_assign_from_terms(&mut a, terms_with_multiples);
798            assert_el_eq!(&ring, ring.from_terms([
799                (base_ring.int_hom().map(7), ring.create_monomial((0..n).map(|_| 0))),
800                (base_ring.int_hom().map(4), xi),
801            ]), a);
802        }
803
804        if n >= 2 {
805            let one = ring.create_monomial((0..n).map(|_| 0));
806            let x0 = ring.indeterminate(0);
807            let x1 = ring.indeterminate(1);
808            let x0_v = ring.create_term(base_ring.one(), ring.clone_monomial(&x0));
809            let x1_v = ring.create_term(base_ring.one(), ring.clone_monomial(&x1));
810            let x0_2 = ring.create_monomial((0..n).map(|k| if k == 0 { 2 } else { 0 }));
811            let x0_3 = ring.create_monomial((0..n).map(|k| if k == 0 { 3 } else { 0 }));
812            let x0_4 = ring.create_monomial((0..n).map(|k| if k == 0 { 4 } else { 0 }));
813            let x0x1 = ring.create_monomial((0..n).map(|k| if k == 0 || k == 1 { 1 } else { 0 }));
814            let x0_3x1 = ring.create_monomial((0..n).map(|k| if k == 0 { 3 } else if k == 1 { 1 } else { 0 }));
815            let x1_2 = ring.create_monomial((0..n).map(|k| if k == 1 { 2 } else { 0 }));
816
817            // test product
818            for a in &elements {
819                for b in &elements {
820                    for c in &elements {
821                        let f = ring.add(ring.inclusion().mul_ref_map(&x0_v, a), ring.inclusion().mul_ref_map(&x1_v, b));
822                        let g = ring.add(ring.inclusion().mul_ref_map(&x0_v, c), ring.clone_el(&x1_v));
823                        let h = ring.from_terms([
824                            (base_ring.mul_ref(a, c), ring.clone_monomial(&x0_2)),
825                            (base_ring.add_ref_snd(base_ring.mul_ref(b, c), a), ring.clone_monomial(&x0x1)),
826                            (base_ring.clone_el(b), ring.clone_monomial(&x1_2)),
827                        ].into_iter());
828                        assert_el_eq!(ring, h, ring.mul(f, g));
829                    }
830                }
831            }
832
833            // test sum
834            for a in &elements {
835                for b in &elements {
836                    for c in &elements {
837                        let f = ring.from_terms([
838                            (base_ring.clone_el(a), ring.clone_monomial(&one)),
839                            (base_ring.clone_el(c), ring.clone_monomial(&x0_2)),
840                            (base_ring.one(), ring.clone_monomial(&x0x1))
841                        ]);
842                        let g = ring.from_terms([
843                            (base_ring.clone_el(b), ring.clone_monomial(&x0)),
844                            (base_ring.one(), ring.clone_monomial(&x0_2)),
845                        ]);
846                        let h = ring.from_terms([
847                            (base_ring.clone_el(a), ring.clone_monomial(&one)),
848                            (base_ring.clone_el(b), ring.clone_monomial(&x0)),
849                            (base_ring.add_ref_fst(c, base_ring.one()), ring.clone_monomial(&x0_2)),
850                            (base_ring.one(), ring.clone_monomial(&x0x1))
851                        ]);
852                        assert_el_eq!(ring, h, ring.add(f, g));
853                    }
854                }
855            }
856
857            // test mul_assign_monomial
858            for a in &elements {
859                for b in &elements {
860                    for c in &elements {
861                        let mut f = ring.from_terms([
862                            (base_ring.clone_el(a), ring.clone_monomial(&one)),
863                            (base_ring.clone_el(b), ring.clone_monomial(&x0)),
864                            (base_ring.clone_el(c), ring.clone_monomial(&x0_2)),
865                            (base_ring.one(), ring.clone_monomial(&x0x1))
866                        ]);
867                        let h = ring.from_terms([
868                            (base_ring.clone_el(a), ring.clone_monomial(&x0_2)),
869                            (base_ring.clone_el(b), ring.clone_monomial(&x0_3)),
870                            (base_ring.clone_el(c), ring.clone_monomial(&x0_4)),
871                            (base_ring.one(), ring.clone_monomial(&x0_3x1))
872                        ]);
873                        let m = ring.clone_monomial(&x0_2);
874                        ring.mul_assign_monomial(&mut f, m);
875                        assert_el_eq!(ring, h, f);
876                    }
877                }
878            }
879
880            // test evaluate
881            for a in &elements {
882                for b in &elements {
883                    let f = ring.from_terms([
884                        (base_ring.int_hom().map(3), ring.clone_monomial(&one)),
885                        (base_ring.int_hom().map(10), ring.clone_monomial(&x0)),
886                        (base_ring.neg_one(), ring.clone_monomial(&x0_2)),
887                        (base_ring.one(), ring.clone_monomial(&x0x1))
888                    ]);
889                    let expected = <_ as RingStore>::sum(base_ring, [
890                        base_ring.int_hom().map(3),
891                        base_ring.int_hom().mul_ref_map(a, &10),
892                        base_ring.negate(base_ring.pow(base_ring.clone_el(a), 2)),
893                        base_ring.mul_ref(a, b)
894                    ]);
895                    let values = [base_ring.clone_el(a), base_ring.clone_el(b)].into_iter().chain((0..(ring.indeterminate_count() - 2)).map(|_| base_ring.zero())).collect::<Vec<_>>();
896                    assert_el_eq!(&base_ring, &expected, &ring.evaluate(&f, values.into_clone_ring_els(base_ring), &base_ring.identity()));
897                }
898            }
899        }
900    }
901}