Implement polynomial gcd

This was done to find the common roots between two polynomials.

2 files changed, 126 insertions, 2 deletions

diff --git a/src/lib.rs b/src/lib.rs
index 31f230e..c70a40f 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -1,3 +1,5 @@
+use std::cmp::Ordering;
+
 pub type Number = f32;
 
 pub type Poly = Vec<Number>;
@@ -29,13 +31,13 @@ pub fn lp(l: Box<Lerp>) -> Poly {
         Lerp::Node(a, b) => {
             let a = lp(a);
             let b = lp(b);
-            let c = poly_sub(&b, &a);
+            let c = sub(&b, &a);
             skewed_sum(a, c)
         }
     }
 }
 
-fn poly_sub(a: &Poly, b: &Poly) -> Poly {
+fn sub(a: &Poly, b: &Poly) -> Poly {
     let mut r = a.clone();
     for i in 0..r.len() {
         r[i] -= b[i];
@@ -51,3 +53,121 @@ fn skewed_sum(a: Poly, b: Poly) -> Poly {
     r.push(b[b.len() - 1]);
     r
 }
+
+fn is_zero(p: &Poly) -> bool {
+    for i in 0..p.len() {
+        if p[i] != 0.0 {
+            return false
+        }
+    }
+    true
+}
+
+fn degree(p: &Poly) -> usize {
+    let mut i = p.len() - 1;
+    while p[i] == 0.0 && i > 0 {
+        i -= 1;
+    }
+    i
+}
+
+fn prod(p: &Poly, n: Number) -> Poly {
+    let mut r = p.clone();
+    for i in 0..r.len() {
+        r[i] *= n;
+    }
+    r
+}
+
+fn shift(p: &Poly, amount: usize) -> Poly {
+    let mut r = vec![0.0; p.len()];
+    for i in 0..p.len() {
+        if i + amount < r.len() {
+            r[i + amount] = p[i];
+        }
+    }
+    r
+}
+
+fn rem(a: &Poly, b: &Poly) -> Poly {
+    let an = a[degree(a)];
+    let bn = b[degree(b)];
+    match degree(a).cmp(&degree(b)) {
+        Ordering::Equal => sub(a, &prod(b, an / bn)),
+        Ordering::Greater => sub(a, &prod(&shift(b, (degree(a) - degree(b)) as usize), an / bn)),
+        Ordering::Less => vec![0.0; b.len()]
+    }
+}
+
+pub fn gcd(a: &Poly, b: &Poly) -> Poly {
+    let c = rem(a, b);
+    if is_zero(&c) {
+        b.clone()
+    } else {
+        gcd(b, &c)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn poly_is_zero() {
+        let p = vec![0.0; 5];
+        assert_eq!(is_zero(&p), true);
+    }
+
+    #[test]
+    fn poly_is_not_zero() {
+        let mut p = vec![0.0; 5];
+        p[2] = 8.0;
+        assert_eq!(is_zero(&p), false);
+    }
+
+    #[test]
+    fn degree_is_five() {
+        let mut p = vec![0.0; 6];
+        p[5] = 2.0;
+        assert_eq!(degree(&p), 5);
+    }
+
+    #[test]
+    fn degree_is_zero() {
+        let p = vec![0.0; 6];
+        assert_eq!(degree(&p), 0);
+    }
+
+    #[test]
+    fn prod_test() {
+        let p = vec![1.0, 2.0, 3.0, 4.0, 5.0];
+        assert_eq!(prod(&p, 2.0), vec![2.0, 4.0, 6.0, 8.0, 10.0]);
+    }
+
+    #[test]
+    fn shift_test() {
+        let p = vec![1.0, 2.0, 3.0, 0.0, 0.0];
+        assert_eq!(shift(&p, 2), vec![0.0, 0.0, 1.0, 2.0, 3.0]);
+    }
+
+    #[test]
+    fn rem_equal() {
+        let a = vec![6.0, 7.0, 1.0];
+        let b = vec![-6.0, -5.0, 1.0];
+        assert_eq!(rem(&a, &b), vec![12.0, 12.0, 0.0]);
+    }
+
+    #[test]
+    fn rem_greater() {
+        let a = vec![-6.0, -5.0, 1.0];
+        let b = vec![12.0, 12.0, 0.0];
+        assert_eq!(rem(&a, &b), vec![-6.0, -6.0, 0.0]);
+    }
+
+    #[test]
+    fn gcd_test() {
+        let a = vec![6.0, 7.0, 1.0];
+        let b = vec![-6.0, -5.0, 1.0];
+        assert_eq!(gcd(&a, &b), vec![-6.0, -6.0, 0.0]);
+    }
+}
diff --git a/src/main.rs b/src/main.rs
index 13f13da..8831976 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -23,6 +23,7 @@ fn main() {
     let b = Lerp::new(vec![400.0, 140.0, 500.0, 300.0]);
     let pa = poly::lp(a);
     let pb = poly::lp(b);
+    let p = poly::gcd(&pa, &pb);
 
     let sdl_context = sdl2::init().unwrap();
     let video_subsystem = sdl_context.video().unwrap();
@@ -51,10 +52,13 @@ fn main() {
         for t in 0..800 {
             let x = eval_poly(&pa, t as Number / 800.0);
             let y = eval_poly(&pb, t as Number / 800.0);
+            let z = eval_poly(&p,  t as Number / 800.0);
             canvas.set_draw_color(Color::RGB(180, 20, 20));
             canvas.fill_rect(Rect::new(t, x as i32, 5, 5)).unwrap();
             canvas.set_draw_color(Color::RGB(20, 180, 20));
             canvas.fill_rect(Rect::new(t, y as i32, 5, 5)).unwrap();
+            canvas.set_draw_color(Color::RGB(20, 20, 180));
+            canvas.fill_rect(Rect::new(t, z as i32, 5, 5)).unwrap();
         }
 
         canvas.present();