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

catlog/zero/
alg.rs
