polynomial_ring/lib.rs
1#![cfg_attr(feature = "__internal_inject_debug", recursion_limit = "16")]
2#![doc = include_str!("../README.md")]
3
4mod sealed {
5 #[cfg(feature = "__internal_inject_debug")]
6 pub trait SizedExt: std::marker::Sized + std::fmt::Debug + std::fmt::Display {}
7 #[cfg(feature = "__internal_inject_debug")]
8 impl<T> SizedExt for T where T: std::marker::Sized + std::fmt::Debug + std::fmt::Display {}
9 #[cfg(not(feature = "__internal_inject_debug"))]
10 pub use std::marker::Sized;
11 #[cfg(feature = "__internal_inject_debug")]
12 pub use SizedExt as Sized;
13}
14use num_traits::{One, Zero};
15use ring_algorithm::{gcd, RingNormalize};
16use sealed::Sized;
17use std::fmt;
18use std::ops::{AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Rem, Sub, SubAssign};
19mod ops;
20#[cfg(feature = "serde")]
21use serde::{Deserialize, Serialize};
22
23/** Polynomial ring $`R[x]`$
24
25```
26use num::Rational64;
27use polynomial_ring::Polynomial;
28let f = Polynomial::new(vec![3, 1, 4, 1, 5].into_iter().map(|x| Rational64::from_integer(x)).collect());
29let g = Polynomial::new(vec![2, 7, 1].into_iter().map(|x| Rational64::from_integer(x)).collect());
30let mut r = f.clone();
31let q = r.division(&g);
32assert_eq!(f, q * g + r);
33```
34*/
35#[derive(Clone, Debug, PartialEq, Eq, Default)]
36#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
37pub struct Polynomial<T> {
38 coeff: Vec<T>,
39}
40
41impl<T: crate::Sized> Polynomial<T> {
42 fn len(&self) -> usize {
43 self.coeff.len()
44 }
45
46 /** The degree of the polynomial.
47 This is the highest power, or [None] for the zero polynomial.
48
49 ```
50 use polynomial_ring::Polynomial;
51 let p = Polynomial::new(vec![3, 2, 1]); // 3+2x+x^2
52 assert_eq!(p.deg(), Some(2));
53 let q = Polynomial::new(vec![0]); // 0
54 assert_eq!(q.deg(), None);
55 ```
56 */
57 pub fn deg(&self) -> Option<usize> {
58 if self.coeff.is_empty() {
59 None
60 } else {
61 Some(self.len() - 1)
62 }
63 }
64
65 /** The leading coefficent.
66 This is the coeffient of the highest power, or [None] for the zero polynomial.
67
68 ```
69 use polynomial_ring::Polynomial;
70 let p = Polynomial::new(vec![3, 2, 1]); // 3+2x+x^2
71 assert_eq!(p.lc(), Some(&1));
72 let q = Polynomial::new(vec![0]); // 0
73 assert_eq!(q.lc(), None);
74 ```
75 */
76 pub fn lc(&self) -> Option<&T> {
77 self.deg().map(|d| &self.coeff[d])
78 }
79
80 /** Get the coefficents.
81
82 ```
83 use polynomial_ring::Polynomial;
84 let p = Polynomial::new(vec![3, 2, 1]); // 3+2x+x^2
85 assert_eq!(p.coeffs(), &[3, 2, 1]);
86 let q = Polynomial::new(vec![0]); // 0
87 assert_eq!(q.coeffs(), &[]);
88 ```
89 */
90 pub fn coeffs(&self) -> &[T] {
91 &self.coeff
92 }
93}
94
95// additive monoid
96impl<M: crate::Sized + Zero> Polynomial<M> {
97 fn trim_zero(&mut self) {
98 let len = self
99 .coeff
100 .iter()
101 .rposition(|x| !x.is_zero())
102 .map(|pos| pos + 1)
103 .unwrap_or(0);
104 self.coeff.truncate(len);
105 }
106
107 /** Construct a polynomial from a [Vec] of coefficients.
108
109 ```
110 use polynomial_ring::Polynomial;
111 let p = Polynomial::new(vec![3, 2, 1]);
112 assert_eq!(p.to_string(), "x^2+2*x+3");
113 ```
114 */
115 pub fn new(coeff: Vec<M>) -> Self {
116 let mut poly = Self { coeff };
117 poly.trim_zero();
118 poly
119 }
120
121 fn extend(&mut self, len: usize) {
122 if self.len() < len {
123 self.coeff.resize_with(len, M::zero);
124 }
125 }
126
127 /** Construct a monomial $`cx^d`$ from a coefficient and a degree ($`c`$=coefficent, $`d`$=degree).
128
129 ```
130 use polynomial_ring::Polynomial;
131 let p = Polynomial::from_monomial(3, 2);
132 let q = Polynomial::new(vec![0, 0, 3]);
133 assert_eq!(p, q);
134 ```
135 */
136 pub fn from_monomial(coefficent: M, degree: usize) -> Self {
137 let coeff = if coefficent.is_zero() {
138 Vec::new()
139 } else {
140 let mut coeff = Vec::with_capacity(degree + 1);
141 for _ in 0..degree {
142 coeff.push(M::zero());
143 }
144 coeff.push(coefficent);
145 coeff
146 };
147 Self { coeff }
148 }
149}
150
151#[macro_export]
152/** Make a [Polynomial] from a list of coefficients (like `vec!`).
153
154```
155use polynomial_ring::{Polynomial, polynomial};
156let p = Polynomial::new(vec![3, 2, 1]);
157let q = polynomial![3, 2, 1];
158assert_eq!(p, q);
159```
160*/
161macro_rules! polynomial {
162 ($($x:expr),*) => {
163 Polynomial::new(vec![$($x), *])
164 }
165}
166
167// unitary ring
168impl<R> fmt::Display for Polynomial<R>
169where
170 R: PartialEq + fmt::Display + Zero + One,
171{
172 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
173 let vec = &self.coeff;
174 if vec.is_empty() {
175 return write!(f, "0");
176 }
177 let mut is_first = true;
178 for (i, c) in vec.iter().enumerate().rev() {
179 if c.is_zero() {
180 continue;
181 }
182 if is_first {
183 is_first = false;
184 } else {
185 write!(f, "+")?;
186 }
187 if c.is_one() {
188 match i {
189 0 => write!(f, "1")?,
190 1 => write!(f, "x")?,
191 _ => write!(f, "x^{i}")?,
192 }
193 } else {
194 match i {
195 0 => write!(f, "{c}")?,
196 1 => write!(f, "{c}*x")?,
197 _ => write!(f, "{c}*x^{i}")?,
198 }
199 }
200 }
201 Ok(())
202 }
203}
204
205impl<R: Sized> Polynomial<R> {
206 /** Evaluate a polynomial at some point, using Horner's method.
207
208 ```
209 use polynomial_ring::Polynomial;
210 let p = Polynomial::new(vec![3, 2, 1]); // 3+2x+x^2
211 assert_eq!(p.eval(&1), 6);
212 assert_eq!(p.eval(&2), 11);
213 ```
214 */
215 pub fn eval<'a>(&self, x: &'a R) -> R
216 where
217 R: Sized + Zero + for<'x> AddAssign<&'x R> + MulAssign<&'a R>,
218 {
219 let mut sum = R::zero();
220 for coeff in self.coeff.iter().rev() {
221 sum *= x;
222 sum += coeff;
223 }
224 sum
225 }
226
227 /** Calculate the derivative.
228
229 ```
230 use polynomial_ring::{Polynomial, polynomial};
231 let p = polynomial![1i32, 2, 3, 2, 1]; // 1+2x+3x^2+2x^3+x^4
232 assert_eq!(p.derivative(), polynomial![2, 6, 6, 4]); // 2+6x+6x^2+4x^3
233 ```
234 */
235 #[must_use]
236 pub fn derivative(&self) -> Self
237 where
238 R: Sized + Zero + One + for<'x> AddAssign<&'x R> + for<'x> From<<&'x R as Mul>::Output>,
239 for<'x> &'x R: Mul,
240 {
241 let n = self.coeff.len();
242 let n = if n > 0 { n - 1 } else { 0 };
243 let mut coeff = Vec::with_capacity(n);
244 let mut i = R::one();
245 for c in self.coeff.iter().skip(1) {
246 coeff.push(R::from(&i * c));
247 i += &R::one();
248 }
249 Polynomial::new(coeff)
250 }
251
252 /** Pseudo division.
253
254 Let $`R`$ be an [integral domain](https://en.wikipedia.org/wiki/Integral_domain).
255 Let $`f, g \in R[x]`$, where $`g \neq 0`$.
256 This function calculates $`s \in R`$, $`q, r \in R[x]`$ s.t. $`sf=qg+r`$,
257 where $`r=0`$ or $`\deg(r)<\deg(g)`$.
258 ```
259 use polynomial_ring::{polynomial, Polynomial};
260
261 let f = polynomial![1i32, 3, 1]; // 1+3x+x^2 ∈ Z[x]
262 let g = polynomial![5, 2]; // 5+2x ∈ Z[x]
263 let mut r = f.clone();
264 let (s, q) = r.pseudo_division(&g);
265 assert_eq!(f * s, q * g + r);
266 ```
267 */
268 pub fn pseudo_division(&mut self, other: &Self) -> (R, Self)
269 where
270 R: Sized
271 + Clone
272 + Zero
273 + One
274 + for<'x> AddAssign<&'x R>
275 + for<'x> MulAssign<&'x R>
276 + for<'x> From<<&'x R as Sub>::Output>
277 + for<'x> From<<&'x R as Mul>::Output>,
278 for<'x> &'x R: Sub + Mul,
279 {
280 let g_deg = other.deg().expect("Division by zero");
281 let f_deg = self.deg();
282 if f_deg < other.deg() {
283 return (R::one(), Self::zero());
284 }
285 let f_deg = f_deg.unwrap();
286 debug_assert!(f_deg >= g_deg);
287 let k = f_deg - g_deg + 1;
288 let lc = other.lc().unwrap();
289 let mut coeff = vec![R::zero(); k];
290 let mut scale = R::one();
291 while self.deg() >= other.deg() {
292 let d = self.deg().unwrap() - g_deg;
293 let c = self.lc().unwrap().clone();
294 for i in 0..other.len() - 1 {
295 self.coeff[i + d] =
296 R::from(&R::from(lc * &self.coeff[i + d]) - &R::from(&c * &other.coeff[i]));
297 }
298 for i in 0..d {
299 self.coeff[i] *= lc;
300 }
301 self.coeff.pop(); // new deg < prev deg
302 self.trim_zero();
303 for c_i in coeff.iter_mut().skip(d + 1) {
304 *c_i *= lc;
305 }
306 coeff[d] = c;
307 scale *= lc;
308 }
309 (scale, Self { coeff })
310 }
311
312 /** Calculate the [resultant](https://en.wikipedia.org/wiki/Resultant)
313
314 ```
315 use polynomial_ring::{polynomial, Polynomial};
316
317 let f = polynomial![0i32, 2, 0, 1]; // 2x+x^3 ∈ Z[x]
318 let g = polynomial![2, 3, 5]; // 2+3x+5x^2 ∈ Z[x]
319 let r = f.resultant(g);
320 assert_eq!(r, 164);
321 ```
322 */
323 pub fn resultant(mut self, other: Self) -> R
324 where
325 R: Sized
326 + Clone
327 + Zero
328 + One
329 + for<'x> AddAssign<&'x R>
330 + for<'x> MulAssign<&'x R>
331 + Neg<Output = R>
332 + for<'x> From<<&'x R as Sub>::Output>
333 + for<'x> From<<&'x R as Mul>::Output>
334 + for<'x> From<<&'x R as Div>::Output>,
335 for<'x> &'x R: Sub + Mul + Div,
336 {
337 let f_deg = self.deg();
338 let g_deg = other.deg();
339 match (f_deg, g_deg) {
340 (Some(0), Some(0)) => R::one(),
341 (Some(0), None) => R::one(),
342 (None, Some(0)) => R::one(),
343 (None, None) => R::zero(),
344 (None, Some(_)) => R::zero(),
345 (Some(_), None) => R::zero(),
346 (Some(0), Some(m)) => {
347 ring_algorithm::power::<R, u64>(self.lc().unwrap().clone(), m as u64)
348 }
349 (Some(n), Some(0)) => {
350 ring_algorithm::power::<R, u64>(other.lc().unwrap().clone(), n as u64)
351 }
352 (Some(n), Some(m)) => {
353 debug_assert!(n >= 1);
354 debug_assert!(m >= 1);
355 let (scale, _) = self.pseudo_division(&other);
356 if let Some(l) = self.deg() {
357 let sign = if n * m % 2 == 0 { R::one() } else { -R::one() };
358 let mul = ring_algorithm::power::<R, u64>(
359 other.lc().unwrap().clone(),
360 (n - l) as u64,
361 );
362 let div = ring_algorithm::power::<R, u64>(scale, m as u64);
363 R::from(&(other.resultant(self) * sign * mul) / &div)
364 } else {
365 // g | f, gcd(f, g) = g
366 R::zero()
367 }
368 }
369 }
370 }
371
372 /** Calculate the primitive part of the input polynomial,
373 that is divide the polynomial by the GCD of its coefficents.
374 ```
375 use polynomial_ring::{polynomial, Polynomial};
376 use num_traits::{One, Zero};
377 let mut f = polynomial![2i32, 4, -2, 6]; // 2+4x+2x^2+6x^3 ∈ Z[x]
378 f.primitive_part_mut();
379 assert_eq!(f, polynomial![1, 2, -1, 3]);// 1+2x+x^2+3x^3 ∈ Z[x]
380 let mut g = polynomial![polynomial![1i32, 1], polynomial![1, 2, 1], polynomial![3, 4, 1], polynomial![-1, -1]]; // (1+x)+(1+2x+x^2)y+(3+4x+x^2)y^2+(-1-x)y^3 ∈ Z[x][y]
381 g.primitive_part_mut();
382 assert_eq!(g, polynomial![polynomial![1], polynomial![1, 1], polynomial![3, 1], polynomial![-1]]); // 1+(1+x)y+(3+x)y^2-y^3 ∈ Z[x][y]
383 ```
384 */
385 pub fn primitive_part_mut(&mut self)
386 where
387 R: Sized
388 + Clone
389 + Zero
390 + for<'x> DivAssign<&'x R>
391 + RingNormalize
392 + for<'x> From<<&'x R as Rem>::Output>,
393 for<'x> &'x R: Rem,
394 {
395 if self.deg().is_none() {
396 return;
397 }
398 let mut g = self.coeff[0].clone();
399 for c in &self.coeff[1..] {
400 g = gcd::<R>(g, c.clone());
401 }
402 g.normalize_mut();
403 *self /= &g;
404 }
405}
406
407// field
408impl<K: Sized> Polynomial<K> {
409 /** Make the polynomial monic, that is divide it by its leading coefficient.
410
411 ```
412 use num::Rational64;
413 use polynomial_ring::Polynomial;
414 let mut p = Polynomial::new(vec![1, 2, 3].into_iter().map(|x| Rational64::from_integer(x)).collect());
415 p.monic();
416 let q = Polynomial::new(vec![(1, 3), (2, 3), (1, 1)].into_iter().map(|(n, d)| Rational64::new(n, d)).collect());
417 assert_eq!(p, q);
418 ```
419 */
420 pub fn monic(&mut self)
421 where
422 K: Clone + for<'x> DivAssign<&'x K>,
423 {
424 if let Some(lc) = self.lc() {
425 let lc = lc.clone();
426 *self /= &lc;
427 }
428 }
429
430 /** Polynomial division.
431
432 ```
433 use num::Rational64;
434 use polynomial_ring::Polynomial;
435 let f = Polynomial::new(vec![3, 1, 4, 1, 5].into_iter().map(|x| Rational64::from_integer(x)).collect());
436 let g = Polynomial::new(vec![2, 7, 1].into_iter().map(|x| Rational64::from_integer(x)).collect());
437 let mut r = f.clone();
438 let q = r.division(&g);
439 assert_eq!(f, q * g + r);
440 ```
441 */
442 pub fn division(&mut self, other: &Self) -> Self
443 where
444 K: Sized
445 + Clone
446 + Zero
447 + for<'x> AddAssign<&'x K>
448 + for<'x> SubAssign<&'x K>
449 + for<'x> MulAssign<&'x K>
450 + for<'x> DivAssign<&'x K>,
451 {
452 let g_deg = other.deg().expect("Division by zero");
453 if self.deg() < other.deg() {
454 return Self::zero();
455 }
456 let lc = other.lc().unwrap();
457 let mut coeff = vec![K::zero(); self.len() - other.len() + 1];
458 while self.deg() >= other.deg() {
459 let d = self.deg().unwrap() - g_deg;
460 let mut c = self.lc().unwrap().clone();
461 c /= lc;
462 for i in 0..other.len() - 1 {
463 let mut c = c.clone();
464 c *= &other.coeff[i];
465 self.coeff[i + d] -= &c;
466 }
467 self.coeff.pop(); // new deg < prev deg
468 self.trim_zero();
469 coeff[d] = c;
470 }
471 Self { coeff }
472 }
473
474 /** Calculate the [square-free decomposition](https://en.wikipedia.org/wiki/Square-free_polynomial).
475
476 ```
477 use polynomial_ring::{Polynomial, polynomial};
478 use num::Rational64;
479 let f = polynomial![Rational64::from(1), Rational64::from(1)];
480 let g = polynomial![Rational64::from(1), Rational64::from(1), Rational64::from(1)];
481 let p = &f * &f * &f * &g * &g; // (x+1)^3(x^2+x+1)^2
482 assert_eq!(p.square_free(), &f * &g); // (x+1)(x^2+x+1)
483 ```
484 */
485 #[must_use]
486 pub fn square_free(&self) -> Self
487 where
488 K: Sized
489 + Clone
490 + Zero
491 + One
492 + for<'x> AddAssign<&'x K>
493 + for<'x> SubAssign<&'x K>
494 + for<'x> MulAssign<&'x K>
495 + for<'x> DivAssign<&'x K>
496 + for<'x> From<<&'x K as Mul>::Output>
497 + for<'x> From<<&'x K as Div>::Output>,
498 for<'x> &'x K: Mul + Div,
499 {
500 let d = self.derivative().into_normalize();
501 let f = ring_algorithm::gcd::<Self>(self.clone(), d).into_normalize();
502 (self / &f).into_normalize()
503 }
504}