1use std::cmp::Ordering;
4
5pub mod coef;
6pub mod interval;
7pub mod mono;
8
9mod arithmetic;
10mod display;
11mod macros;
12
13pub use arithmetic::monomial_div;
14pub use macros::poly;
15
16use coef::Coef;
17pub use mono::{Mono, Pow, Var};
18
19#[derive(PartialEq, Eq, Hash, Clone)]
30pub struct Poly {
31 terms: Vec<(Coef, Mono)>,
32}
33
34impl Poly {
35 pub fn zero() -> Self {
36 Self { terms: vec![] }
37 }
38
39 pub fn is_zero(&self) -> bool {
40 self.terms.is_empty()
41 }
42
43 pub fn constant<T: Into<Coef>>(c: T) -> Self {
44 let c = c.into();
45 if c == 0.into() {
46 Poly::zero()
47 } else {
48 Poly {
49 terms: vec![(c, Mono::unit())],
50 }
51 }
52 }
53
54 pub fn var<T: Into<Var>>(v: T, p: Pow) -> Self {
55 if p == 0 {
56 Self::constant(1)
57 } else {
58 Self::term(1, Mono::new(vec![(v.into(), p)]))
59 }
60 }
61
62 pub fn term<T: Into<Coef>>(c: T, m: Mono) -> Self {
63 let c = c.into();
64 if c == 0.into() {
65 Poly::zero()
66 } else {
67 Poly {
68 terms: vec![(c, m)],
69 }
70 }
71 }
72
73 pub fn terms(&self) -> Vec<(Coef, Mono)> {
74 self.terms.clone()
75 }
76
77 pub fn from_terms(terms: impl IntoIterator<Item = (Coef, Mono)>) -> Self {
78 let mut raw: Vec<_> = terms.into_iter().collect();
79 raw.sort_unstable_by(|(_, m1), (_, m2)| m2.cmp(m1));
80 let mut terms = Vec::with_capacity(raw.len());
81 for (c, m) in raw {
82 match terms.last_mut() {
83 Some((lc, lm)) if *lm == m => {
84 *lc = *lc + c;
85 if *lc == 0.into() {
86 terms.pop();
87 }
88 }
89 _ => {
90 if c != 0.into() {
91 terms.push((c, m));
92 }
93 }
94 }
95 }
96 Self { terms }.debug_checked()
97 }
98
99 pub fn total_degree(&self) -> Pow {
100 self.terms
101 .first()
102 .map(|(_, m)| m.total_degree())
103 .unwrap_or(0)
104 }
105
106 pub fn degree_in<T: Into<Var>>(&self, v: T) -> Pow {
107 let v = v.into();
108 let mut deg = 0;
109 for (_, mono) in &self.terms {
110 let mut term_deg = 0;
111 for (var, pow) in mono.exps() {
112 term_deg += pow;
113 match var.cmp(&v) {
114 Ordering::Greater => continue,
115 Ordering::Equal => {
116 deg = deg.max(*pow);
117 break;
118 }
119 Ordering::Less => break,
120 }
121 }
122 if term_deg < deg {
123 break;
124 }
125 }
126
127 deg
128 }
129
130 pub fn leading_term(&self) -> (Coef, Mono) {
131 match self.terms.first() {
132 Some(m) => m.clone(),
133 None => (Coef::new(0), Mono::unit()),
134 }
135 }
136
137 pub fn eval(&self, mut point: impl FnMut(Var) -> i128) -> Coef {
138 let mut acc = Coef::new(0);
139 for (c, m) in &self.terms {
140 let mut t = *c;
141 for &(v, e) in m.exps() {
142 t = t * Coef::new(point(v).checked_pow(e).unwrap());
143 }
144 acc = acc + t;
145 }
146 acc
147 }
148
149 pub fn map_vars(&self, mut f: impl FnMut(Var) -> Self) -> Self {
150 let mut acc = Poly::zero();
151 for (c, m) in &self.terms {
152 let mut t = Poly::constant(*c);
153 for &(v, e) in m.exps() {
154 t = t.mul(&f(v).pow(e));
155 }
156 acc = acc.add(&t);
157 }
158 acc
159 }
160
161 pub fn substitute<T: Into<Var>>(&self, v: T, q: &Self) -> Self {
162 let v = v.into();
163 self.map_vars(|w| if w == v { q.clone() } else { Self::var(w, 1) })
164 }
165
166 pub fn coef_gcd(&self) -> Coef {
167 self.terms
168 .iter()
169 .fold(Coef::from(0), |g, &(c, _)| g.gcd(&c))
170 .abs()
171 }
172
173 pub fn coef_lcm(&self) -> Coef {
174 self.terms
175 .iter()
176 .fold(Coef::from(1), |g, &(c, _)| g.lcm(&c))
177 .abs()
178 }
179
180 pub fn divide_coefs<T: Into<Coef>>(&self, d: T) -> Option<Self> {
181 let d = d.into();
182 if d.is_zero() {
183 return None;
184 }
185 let mut terms = Vec::with_capacity(self.terms.len());
186 for (c, m) in &self.terms {
187 let (quot, rem) = c.divrem(d);
188 if rem.is_zero() {
189 terms.push((quot, m.clone()));
190 } else {
191 return None;
192 }
193 }
194 Some(Poly { terms }.debug_checked())
195 }
196
197 pub fn as_constant(&self) -> Option<Coef> {
198 match self.terms.as_slice() {
199 [] => Some(Coef::new(0)),
200 [(c, m)] if m.is_unit() => Some(*c),
201 _ => None,
202 }
203 }
204
205 pub fn vars(&self) -> Vec<Var> {
206 let mut vs = vec![];
207
208 for term in &self.terms {
209 for (var, _) in term.1.exps() {
210 if !vs.contains(var) {
211 vs.push(*var);
212 }
213 }
214 }
215
216 vs.sort_unstable();
217 vs
218 }
219
220 pub fn always_nonneg(&self, nonneg: &dyn Fn(Var) -> bool) -> bool {
221 self.terms()
222 .iter()
223 .all(|(c, m)| c.is_positive() && m.always_nonneg(nonneg))
224 }
225
226 pub fn assert_canonical(&self) {
227 for w in self.terms.windows(2) {
228 assert!(
229 w[0].1 > w[1].1,
230 "terms not strictly descending: {:?} then {:?}",
231 w[0],
232 w[1]
233 );
234 }
235 for (c, m) in &self.terms {
236 assert!(*c != 0.into(), "zero coefficient on {m:?}");
237 m.assert_canonical();
238 }
239 }
240
241 fn debug_checked(self) -> Self {
242 #[cfg(debug_assertions)]
243 self.assert_canonical();
244 self
245 }
246}
247
248#[cfg(test)]
249mod tests {
250 use super::{Coef, Poly, monomial_div, poly};
251
252 #[test]
253 fn constants() {
254 let x = poly!(3 * x ^ 2 + 5 * z - 2 + 1);
255 let y = poly!(4 * y * z);
256 assert_eq!(
257 x.mul(&y),
258 poly!(12 * x ^ 2 * y * z + 20 * y * z ^ 2 - 4 * y * z)
259 );
260 assert_eq!(poly!(), poly!(x - x));
261 }
262
263 #[test]
264 fn arithmetic_sanity() {
265 let x = poly!(a ^ 2 + 2 * b + c);
266 let y = poly!(2 * a ^ 2 - c ^ 3 + d);
267 assert!(x.add(&y) == poly!(-1 * c ^ 3 + 3 * a ^ 2 + 2 * b + c + d));
268
269 let x = poly!(a ^ 4 - b ^ 4);
270 let y = poly!(a ^ 2 + b ^ 2);
271 assert_eq!(x.try_divide(&y), Some(poly!(a ^ 2 - b ^ 2)));
272
273 let x = poly!(a ^ 2 - 2 * a * b + b ^ 2);
274 let y = poly!(a - b);
275 assert_eq!(x.try_divide(&y), Some(poly!(a - b)));
276
277 let x = poly!(-4 * b);
278 let y = poly!(a);
279 assert_eq!(x.try_divide(&y), None);
280 }
281
282 #[test]
283 fn arithmetic_fuzz() {
284 use rand::prelude::*;
285
286 use super::Poly;
287
288 let mut rng = SmallRng::seed_from_u64(1);
289
290 fn create_random_poly(rng: &mut SmallRng, term_max: i32) -> Poly {
291 let mut p = Poly::zero();
292
293 for _ in 0..rng.gen_range(0..term_max + 1) {
294 let coef = rng.gen_range(-7..7);
295 let xpow = rng.gen_range(0..2);
296 let ypow = rng.gen_range(0..2);
297 let zpow = rng.gen_range(0..4);
298 let coef = Poly::constant(coef);
299 let xpow = Poly::var('x', xpow);
300 let ypow = Poly::var('y', ypow);
301 let zpow = Poly::var('z', zpow);
302
303 p = p.add(&coef.mul(&xpow).mul(&ypow).mul(&zpow));
304 }
305
306 p
307 }
308
309 for _ in 0..10_000 {
310 let dividend = create_random_poly(&mut rng, 10);
311 let n_divs = rng.gen_range(0..4);
312 let mut divisors: Vec<_> = std::iter::repeat_with(|| create_random_poly(&mut rng, 6))
313 .take(n_divs)
314 .collect();
315
316 let (quotients, rem) = dividend.compound_divide(&mut divisors);
317
318 println!("-------------------------");
319 println!("calculated {:?} / {:?}", dividend, divisors);
320 println!("got {:?} rem {:?}", quotients, rem);
321 println!("-------------------------");
322
323 let calculated_dividend = quotients
324 .into_iter()
325 .zip(divisors.clone())
326 .fold(Poly::zero(), |acc, (x, y)| acc.add(&mut x.mul(&y)))
327 .add(&mut rem.clone());
328
329 assert_eq!(calculated_dividend, dividend);
330 }
331 }
332
333 pub fn s_poly(p: &Poly, q: &Poly) -> Poly {
334 let p_lt = Poly {
335 terms: vec![p.leading_term()],
336 };
337 let q_lt = Poly {
338 terms: vec![q.leading_term()],
339 };
340
341 let lcm_ltp_ltq = Poly {
342 terms: vec![(
343 p_lt.terms[0].0 * q_lt.terms[0].0,
344 p_lt.terms[0].1.lcm(&q_lt.terms[0].1),
345 )],
346 };
347
348 let coef_p = lcm_ltp_ltq.try_divide(&p_lt).unwrap();
349 let coef_q = lcm_ltp_ltq.try_divide(&q_lt).unwrap();
350
351 let a = coef_p.mul(&p);
352 let b = coef_q.mul(&q);
353 let res = a.sub(&b);
354
355 if res.is_zero() {
356 res
357 } else {
358 res.divide_coefs(res.coef_gcd()).unwrap()
359 }
360 }
361
362 pub fn groebner_basis(sys: &[Poly]) -> Vec<Poly> {
363 let mut sys = sys.to_vec();
364
365 let mut combs = {
368 let mut combs = vec![];
369 for i in 0..sys.len() {
370 for j in 0..sys.len() {
371 if i != j {
372 combs.push((sys[i].clone(), sys[j].clone()));
373 }
374 }
375 }
376
377 combs
378 };
379
380 while let Some((a, b)) = combs.pop() {
381 let s = s_poly(&a, &b);
382 let (_, rem) = s.compound_divide(&sys);
383
384 if !rem.is_zero() {
385 for member in &sys {
386 combs.push((member.clone(), rem.clone()));
387 }
388 sys.push(rem);
389 }
390 }
391
392 let mut keep = vec![];
395
396 for i in 0..sys.len() {
397 let mut divides_any = false;
398
399 for j in 0..sys.len() {
400 if i != j {
401 let i_lt = (Coef::new(1), sys[i].leading_term().1);
402 let j_lt = (Coef::new(1), sys[j].leading_term().1);
403 if let Some((_, m)) = monomial_div(&i_lt, &j_lt) {
404 if m.is_unit() {
405 divides_any = i > j;
406 } else {
407 divides_any = true;
408 }
409
410 if divides_any {
411 break;
412 }
413 }
414 }
415 }
416
417 if !divides_any {
418 keep.push(sys[i].clone());
419 }
420 }
421
422 keep.sort_by(|p, q| p.leading_term().1.cmp(&q.leading_term().1).reverse());
423 println!("sys: {:?}", keep);
424
425 let mut keep2 = vec![];
426
427 for (i, k) in keep.iter().enumerate() {
428 let all_except = keep
429 .iter()
430 .enumerate()
431 .filter_map(|(j, p)| if j != i { Some(p.clone()) } else { None })
432 .collect::<Vec<_>>();
433
434 let all_except_lcm = all_except
435 .iter()
436 .fold(Coef::new(1), |c, p| c.lcm(&p.coef_lcm()));
437
438 let (_, rem) = k.mul_scalar(all_except_lcm).compound_divide(&all_except);
439 if !rem.is_zero() {
440 keep2.push(rem.divide_coefs(rem.coef_gcd()).unwrap());
441 }
442 }
443
444 for p in keep2.iter_mut() {
445 if !p.leading_term().0.is_positive() {
446 *p = p.mul_scalar(-1);
447 }
448 }
449
450 keep2
451 }
452
453 #[test]
454 fn groebner_basis_validation() {
455 let sys = [
456 poly!(x + y ^ 2 + z),
457 poly!(x - y + 3 * z + 5),
458 poly!(x - 2 * y + 3),
459 ];
460
461 let gb = groebner_basis(&sys);
462
463 let correct = [
464 poly!(9 * z ^ 2 + 7 * z - 3),
465 poly!(x + 6 * z + 7),
466 poly!(y + 3 * z + 2),
467 ];
468
469 println!("{:?}", gb);
470 println!("{:?}", correct);
471
472 for (g, c) in gb.iter().zip(correct.iter()) {
473 assert_eq!(g, c, "{:?} = {:?}", g, c);
474 }
475
476 let sys = [
477 poly!(x ^ 2 * y + 1),
478 poly!(2 * x + y * z - 1),
479 poly!(x - y ^ 2 * z ^ 2 + 1),
480 ];
481
482 let gb = groebner_basis(&sys);
483
484 let correct = [poly!(4 * x - 5), poly!(25 * y + 16), poly!(32 * z - 75)];
485
486 println!("{:?}", gb);
487 println!("{:?}", correct);
488
489 for (g, c) in gb.iter().zip(correct.iter()) {
490 assert_eq!(g, c, "{:?} = {:?}", g, c);
491 }
492 }
493}