add S-polynomials

src/poly/flat.rs

@@ -116,6 +116,34 @@ impl<V: Var> Mono<V> {
         true
     }
 
+    /// Returns the lcm of self and other (element-wise max of exponents).
+    pub fn lcm(&self, other: &Mono<V>) -> Mono<V> {
+        let mut self_it = self.term.iter().peekable();
+        let mut other_it = other.term.iter().peekable();
+        let mut result: Vec<(V, u32)> = vec![];
+
+        loop {
+            match (self_it.peek(), other_it.peek()) {
+                (None, None) => break,
+                (Some(_), None) => result.push(self_it.next().unwrap().clone()),
+                (None, Some(_)) => result.push(other_it.next().unwrap().clone()),
+                (Some((s_var, _)), Some((o_var, _))) => {
+                    if s_var < o_var {
+                        result.push(self_it.next().unwrap().clone());
+                    } else if s_var > o_var {
+                        result.push(other_it.next().unwrap().clone());
+                    } else {
+                        let (var, s_exp) = self_it.next().unwrap();
+                        let (_, o_exp) = other_it.next().unwrap();
+                        result.push((var.clone(), (*s_exp).max(*o_exp)));
+                    }
+                }
+            }
+        }
+
+        Mono { term: result }
+    }
+
     /// Divides self by other. Assumes `self.contains(other)`.
     pub fn div(self, other: &Mono<V>) -> Mono<V> {
         let mut self_it = self.term.into_iter().peekable();

src/poly/ops.rs

@@ -128,6 +128,38 @@ impl<V: Var> Sub for Poly<V> {
     }
 }
 
+fn gcd(a: i32, b: i32) -> i32 {
+    let (mut a, mut b) = (a.abs(), b.abs());
+    while b != 0 {
+        let t = b;
+        b = a % b;
+        a = t;
+    }
+    a
+}
+
+impl<V: Var> Poly<V> {
+    /// S-polynomial of `self` and `other` under lex order.
+    ///
+    /// Let `L = lcm(LM(f), LM(g))` and `d = gcd(lc(f), lc(g))`. Returns:
+    /// `(lc(g)/d) * (L/LM(f)) * f - (lc(f)/d) * (L/LM(g)) * g`
+    ///
+    /// The leading terms cancel by construction. Panics if either polynomial is zero.
+    pub fn s_poly(&self, other: &Poly<V>) -> Poly<V> {
+        let (lm_f, lc_f) = self.leading_term_lex().expect("f must be nonzero");
+        let (lm_g, lc_g) = other.leading_term_lex().expect("g must be nonzero");
+
+        let l = lm_f.lcm(&lm_g);
+        let t_f: Poly<V> = 1 * l.clone().div(&lm_f);
+        let t_g: Poly<V> = 1 * l.div(&lm_g);
+
+        let d = gcd(lc_f, lc_g);
+        let left = (lc_g / d) * (t_f * self.clone());
+        let right = (lc_f / d) * (t_g * other.clone());
+        left - right
+    }
+}
+
 impl<V: Var> Poly<V> {
     /// Pseudo-division with remainder.
     ///

src/poly/tests.rs

@@ -145,6 +145,61 @@ fn test_mono_div() {
     assert_eq!(a.div(&b), Mono { term: vec![] });
 }
 
+#[test]
+fn test_mono_lcm() {
+    // lcm(x², xy) = x²y
+    let a: Mono<StaticVar> = [("x", 2)].into();
+    let b: Mono<StaticVar> = [("x", 1), ("y", 1)].into();
+    assert_eq!(a.lcm(&b), Mono::from([("x", 2), ("y", 1)]));
+
+    // lcm(x, y) = xy  (disjoint)
+    let a: Mono<StaticVar> = [("x", 1)].into();
+    let b: Mono<StaticVar> = [("y", 1)].into();
+    assert_eq!(a.lcm(&b), Mono::from([("x", 1), ("y", 1)]));
+
+    // lcm(x², x²) = x²  (identical)
+    let a: Mono<StaticVar> = [("x", 2)].into();
+    assert_eq!(a.lcm(&a.clone()), a);
+}
+
+#[test]
+fn test_s_poly() {
+    // S(x², xy) = 0: S-poly of two monomials always vanishes
+    let f: Poly<StaticVar> = [(1, [("x", 2)])].into();
+    let g: Poly<StaticVar> = [(1, [("x", 1), ("y", 1)])].into();
+    assert!(f.s_poly(&g).is_zero());
+
+    // f = x² + y, g = xy + z  (lex: LM(f)=x², LM(g)=xy)
+    // L = x²y, t_f = y, t_g = x, d = gcd(1,1) = 1
+    // S = y*(x²+y) - x*(xy+z) = x²y + y² - x²y - xz = y² - xz
+    let f: Poly<StaticVar> = [
+        (1i32, Mono::from([("x", 2u32)])),
+        (1i32, Mono::from([("y", 1u32)])),
+    ].into_iter().collect();
+    let g: Poly<StaticVar> = [
+        (1i32, Mono::from([("x", 1u32), ("y", 1u32)])),
+        (1i32, Mono::from([("z", 1u32)])),
+    ].into_iter().collect();
+    let expected: Poly<StaticVar> = [
+        (1i32, Mono::from([("y", 2u32)])),
+        (-1i32, Mono::from([("x", 1u32), ("z", 1u32)])),
+    ].into_iter().collect();
+    assert_eq!(f.s_poly(&g), expected);
+
+    // f = 2x + y, g = 3x + z  (same LM=x, d=gcd(2,3)=1)
+    // S = (3/1)*f - (2/1)*g = 3(2x+y) - 2(3x+z) = 3y - 2z
+    let f: Poly<StaticVar> = [(2, [("x", 1)]), (1, [("y", 1)])].into();
+    let g: Poly<StaticVar> = [(3, [("x", 1)]), (1, [("z", 1)])].into();
+    let expected: Poly<StaticVar> = [(3, [("y", 1)]), (-2, [("z", 1)])].into();
+    assert_eq!(f.s_poly(&g), expected);
+
+    // f = 4x, g = 6x  (d = gcd(4,6) = 2)
+    // S = (6/2)*x*4x - (4/2)*x*6x = 3*4x - 2*6x = 12x - 12x = 0
+    let f: Poly<StaticVar> = [(4, [("x", 1)])].into();
+    let g: Poly<StaticVar> = [(6, [("x", 1)])].into();
+    assert!(f.s_poly(&g).is_zero());
+}
+
 fn make_const_poly(c: i32) -> Poly<StaticVar> {
     Poly { mono: [(Mono { term: vec![] }, c)].into_iter().collect() }
 }