catlog/zero/
alg.rs

1//! Commutative algebra and polynomials.
2
3use num_traits::{One, Pow, Zero};
4use std::collections::BTreeMap;
5use std::fmt::Display;
6use std::iter::{Product, Sum};
7use std::ops::{Add, AddAssign, Mul, Neg};
8
9use derivative::Derivative;
10
11use super::rig::*;
12
13/// A commutative algebra over a commutative ring.
14pub trait CommAlg: CommRing + Module<Ring = Self::R> {
15    /// The base ring of the algebra.
16    type R: CommRing;
17
18    /** Convert an element of the base ring into an element of the algebra.
19
20    A commutative algebra A over a commutative ring R can be defined as a ring
21    homomorphism from R to A. This method computes that homomorphism.
22     */
23    fn from_scalar(r: Self::R) -> Self {
24        Self::one() * r
25    }
26}
27
28/** A polynomial in several variables.
29
30This data structure is for polynomials in *normal form*: a **polynomial** is a
31formal linear combination of monomials in which no monomial is repeated, and no
32variable is repeated within any monomial. The implementation is indeed a
33[`Combination`] of a [`Monomial`]s. The use of a normal form means that
34polynomial arithmetic automatically performs certain simplifications.
35
36In abstract terms, polynomials with coefficients valued in a [commutative
37ring](super::rig::CommRing) *R* are the free [commutative algebra](CommAlg)
38over *R*.
39 */
40#[derive(Clone, PartialEq, Eq, Derivative)]
41#[derivative(Default(bound = ""))]
42pub struct Polynomial<Var, Coef, Exp>(Combination<Monomial<Var, Exp>, Coef>);
43
44impl<Var, Coef, Exp> Polynomial<Var, Coef, Exp>
45where
46    Var: Ord,
47    Exp: Ord,
48{
49    /// Constructs the generating polynomial corresponding to a variable.
50    pub fn generator(var: Var) -> Self
51    where
52        Coef: One,
53        Exp: One,
54    {
55        Polynomial::from_monomial(Monomial::generator(var))
56    }
57
58    /// Constructs the polynomial corresponding to a monomial.
59    pub fn from_monomial(m: Monomial<Var, Exp>) -> Self
60    where
61        Coef: One,
62    {
63        Polynomial(Combination::generator(m))
64    }
65
66    /// Iterates over the monomials in the polynomial.
67    pub fn monomials(&self) -> impl ExactSizeIterator<Item = &Monomial<Var, Exp>> {
68        self.0.variables()
69    }
70
71    /** Maps the coefficients of the polynomial.
72
73    In the usual situations when the coefficients from commutative rigs and the
74    mapping is a rig homomorphism, this operation is extension of scalars
75    applied to free commutative algebras.
76     */
77    pub fn extend_scalars<NewCoef, F>(self, f: F) -> Polynomial<Var, NewCoef, Exp>
78    where
79        F: FnMut(Coef) -> NewCoef,
80    {
81        Polynomial(self.0.extend_scalars(f))
82    }
83
84    /// Evaluates the polynomial by substituting for the variables.
85    pub fn eval<A, F>(&self, f: F) -> A
86    where
87        A: Clone + Mul<Coef, Output = A> + Pow<Exp, Output = A> + Sum + Product,
88        F: Clone + FnMut(&Var) -> A,
89        Coef: Clone,
90        Exp: Clone,
91    {
92        self.0.eval_with_order(self.monomials().map(|m| m.eval(f.clone())))
93    }
94
95    /** Evaluates the polynomial on a sequence of variable-value pairs.
96
97    This is a convenient way to evaluate the polynomial at a single point but it
98    is not very efficient.
99    */
100    pub fn eval_pairs<A>(&self, pairs: impl IntoIterator<Item = (Var, A)>) -> A
101    where
102        A: Clone + Mul<Coef, Output = A> + Pow<Exp, Output = A> + Sum + Product,
103        Coef: Clone,
104        Exp: Clone,
105    {
106        let map: BTreeMap<Var, A> = pairs.into_iter().collect();
107        self.eval(|var| map.get(var).cloned().unwrap())
108    }
109
110    /** Maps over the variables in the polynomial.
111
112    The mapping need not be injective. This is conceptually equivalent to
113    [evaluating](Polynomial::eval) the polynomial with a map that sends
114    generators to generators but avoids assuming that an arbitrary polynomial
115    can be exponentiated, which is only makes sense when the exponents are
116    nonnegative integers.
117     */
118    pub fn map_variables<NewVar, F>(&self, mut f: F) -> Polynomial<NewVar, Coef, Exp>
119    where
120        Coef: Clone + Add<Output = Coef>,
121        Exp: Clone + Add<Output = Exp>,
122        NewVar: Clone + Ord,
123        F: FnMut(&Var) -> NewVar,
124    {
125        (&self.0)
126            .into_iter()
127            .map(|(coef, m)| (coef.clone(), m.map_variables(|var| f(var))))
128            .collect()
129    }
130
131    /** Puts the polynomial into normal form.
132
133    The data structure for polynomials is already pretty close to being a normal
134    form, but allows the possibility of coefficients or exponents being zero.
135    This method removes those if present.
136     */
137    pub fn normalize(self) -> Self
138    where
139        Coef: Zero,
140        Exp: Zero,
141    {
142        self.0
143            .into_iter()
144            .filter_map(|(coef, m)| {
145                if coef.is_zero() {
146                    None
147                } else {
148                    Some((coef, m.normalize()))
149                }
150            })
151            .collect()
152    }
153}
154
155impl<Var, Coef, Exp> FromIterator<(Coef, Monomial<Var, Exp>)> for Polynomial<Var, Coef, Exp>
156where
157    Var: Ord,
158    Coef: Add<Output = Coef>,
159    Exp: Ord,
160{
161    fn from_iter<T: IntoIterator<Item = (Coef, Monomial<Var, Exp>)>>(iter: T) -> Self {
162        Polynomial(iter.into_iter().collect())
163    }
164}
165
166impl<Var, Coef, Exp> Display for Polynomial<Var, Coef, Exp>
167where
168    Var: Display,
169    Coef: Display + PartialEq + One,
170    Exp: Display + PartialEq + One,
171{
172    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
173        self.0.fmt(f)
174    }
175}
176
177// XXX: Lots of boilerplate to delegate the module structure of `Polynomial` to
178// `Combination` because these traits cannot be derived automatically.
179
180impl<Var, Coef, Exp> AddAssign<(Coef, Monomial<Var, Exp>)> for Polynomial<Var, Coef, Exp>
181where
182    Var: Ord,
183    Coef: Add<Output = Coef>,
184    Exp: Ord,
185{
186    fn add_assign(&mut self, rhs: (Coef, Monomial<Var, Exp>)) {
187        self.0 += rhs;
188    }
189}
190
191impl<Var, Coef, Exp> Add for Polynomial<Var, Coef, Exp>
192where
193    Var: Ord,
194    Coef: Add<Output = Coef>,
195    Exp: Ord,
196{
197    type Output = Self;
198
199    fn add(self, rhs: Self) -> Self::Output {
200        Polynomial(self.0 + rhs.0)
201    }
202}
203
204impl<Var, Coef, Exp> Zero for Polynomial<Var, Coef, Exp>
205where
206    Var: Ord,
207    Coef: Add<Output = Coef> + Zero,
208    Exp: Ord,
209{
210    fn zero() -> Self {
211        Polynomial(Combination::default())
212    }
213
214    fn is_zero(&self) -> bool {
215        self.0.is_zero()
216    }
217}
218
219impl<Var, Coef, Exp> AdditiveMonoid for Polynomial<Var, Coef, Exp>
220where
221    Var: Ord,
222    Coef: AdditiveMonoid,
223    Exp: Ord,
224{
225}
226
227impl<Var, Coef, Exp> Mul<Coef> for Polynomial<Var, Coef, Exp>
228where
229    Var: Ord,
230    Coef: Clone + Default + Mul<Output = Coef>,
231    Exp: Ord,
232{
233    type Output = Self;
234
235    fn mul(self, a: Coef) -> Self::Output {
236        Polynomial(self.0 * a)
237    }
238}
239
240impl<Var, Coef, Exp> RigModule for Polynomial<Var, Coef, Exp>
241where
242    Var: Ord,
243    Coef: Clone + Default + CommRig,
244    Exp: Ord,
245{
246    type Rig = Coef;
247}
248
249impl<Var, Coef, Exp> Neg for Polynomial<Var, Coef, Exp>
250where
251    Var: Ord,
252    Coef: Default + Neg<Output = Coef>,
253    Exp: Ord,
254{
255    type Output = Self;
256
257    fn neg(self) -> Self::Output {
258        Polynomial(self.0.neg())
259    }
260}
261
262impl<Var, Coef, Exp> AbGroup for Polynomial<Var, Coef, Exp>
263where
264    Var: Ord,
265    Coef: Default + AbGroup,
266    Exp: Ord,
267{
268}
269
270impl<Var, Coef, Exp> Module for Polynomial<Var, Coef, Exp>
271where
272    Var: Ord,
273    Coef: Clone + Default + CommRing,
274    Exp: Ord,
275{
276    type Ring = Coef;
277}
278
279/// Multiply polynomials using the distributive law.
280impl<Var, Coef, Exp> Mul for Polynomial<Var, Coef, Exp>
281where
282    Var: Clone + Ord,
283    Coef: Clone + Add<Output = Coef> + Mul<Output = Coef>,
284    Exp: Clone + Ord + Add<Output = Exp>,
285{
286    type Output = Self;
287
288    fn mul(self, rhs: Self) -> Self::Output {
289        // Avoid unnecessary clones by tracking whether we're in the last
290        // iteration of the outer and inner loops.
291        let mut result = Polynomial::default();
292        let (outer, inner) = (self.0, rhs.0);
293        let mut outer_iter = outer.into_iter();
294        while let Some((a, m)) = outer_iter.next() {
295            if outer_iter.len() == 0 {
296                let mut inner_iter = inner.into_iter();
297                while let Some((b, n)) = inner_iter.next() {
298                    if inner_iter.len() == 0 {
299                        result += (a * b, m * n);
300                        break;
301                    } else {
302                        result += (a.clone() * b, m.clone() * n);
303                    }
304                }
305                break;
306            } else {
307                for (b, n) in &inner {
308                    result += (a.clone() * b.clone(), m.clone() * n.clone());
309                }
310            }
311        }
312        result
313    }
314}
315
316impl<Var, Coef, Exp> One for Polynomial<Var, Coef, Exp>
317where
318    Var: Clone + Ord,
319    Coef: Clone + Add<Output = Coef> + One,
320    Exp: Clone + Ord + Add<Output = Exp>,
321{
322    fn one() -> Self {
323        Polynomial::from_monomial(Default::default())
324    }
325}
326
327impl<Var, Coef, Exp> Monoid for Polynomial<Var, Coef, Exp>
328where
329    Var: Clone + Ord,
330    Coef: Clone + Rig,
331    Exp: Clone + Ord + AdditiveMonoid,
332{
333}
334
335impl<Var, Coef, Exp> Rig for Polynomial<Var, Coef, Exp>
336where
337    Var: Clone + Ord,
338    Coef: Clone + Rig,
339    Exp: Clone + Ord + AdditiveMonoid,
340{
341}
342
343impl<Var, Coef, Exp> Ring for Polynomial<Var, Coef, Exp>
344where
345    Var: Clone + Ord,
346    Coef: Clone + Default + Ring,
347    Exp: Clone + Ord + AdditiveMonoid,
348{
349}
350
351impl<Var, Coef, Exp> CommMonoid for Polynomial<Var, Coef, Exp>
352where
353    Var: Clone + Ord,
354    Coef: Clone + CommRig,
355    Exp: Clone + Ord + AdditiveMonoid,
356{
357}
358
359impl<Var, Coef, Exp> CommRig for Polynomial<Var, Coef, Exp>
360where
361    Var: Clone + Ord,
362    Coef: Clone + CommRig,
363    Exp: Clone + Ord + AdditiveMonoid,
364{
365}
366
367impl<Var, Coef, Exp> CommRing for Polynomial<Var, Coef, Exp>
368where
369    Var: Clone + Ord,
370    Coef: Clone + Default + CommRing,
371    Exp: Clone + Ord + AdditiveMonoid,
372{
373}
374
375impl<Var, Coef, Exp> CommAlg for Polynomial<Var, Coef, Exp>
376where
377    Var: Clone + Ord,
378    Coef: Clone + Default + CommRing,
379    Exp: Clone + Ord + AdditiveMonoid,
380{
381    type R = Coef;
382
383    fn from_scalar(r: Self::R) -> Self {
384        [(r, Monomial::one())].into_iter().collect()
385    }
386}
387
388#[cfg(test)]
389mod tests {
390    use super::*;
391
392    #[test]
393    fn polynomials() {
394        let x = || Polynomial::<_, i32, u32>::generator('x');
395        let y = || Polynomial::<_, i32, u32>::generator('y');
396        assert_eq!(x().to_string(), "x");
397
398        let p = Polynomial::<char, i32, u32>::from_scalar(-5);
399        assert_eq!(p.eval_pairs::<i32>([]), -5);
400
401        let p = x() * y() * x() * 2 + y() * x() * y() * 3;
402        assert_eq!(p.to_string(), "3 x y^2 + 2 x^2 y");
403        assert_eq!(p.map_variables(|_| 'x').to_string(), "5 x^3");
404        assert_eq!(p.eval_pairs([('x', 1), ('y', 1)]), 5);
405        assert_eq!(p.eval_pairs([('x', 1), ('y', 2)]), 16);
406        assert_eq!(p.eval_pairs([('y', 1), ('x', 2)]), 14);
407
408        let p = (x() + y()) * (x() + y());
409        assert_eq!(p.to_string(), "2 x y + x^2 + y^2");
410
411        let p = (x() + y()) * (x() + y().neg());
412        assert_eq!(p.normalize().to_string(), "x^2 + (-1) y^2");
413    }
414}
catlog/zero/alg.rs

catlog/zero/
alg.rs
