Fix div not stopping at remainder 0

1 files changed, 17 insertions, 12 deletions

diff --git a/src/poly.rs b/src/poly.rs
index 6ffaf7d..60eb707 100644
--- a/src/poly.rs
+++ b/src/poly.rs
@@ -8,6 +8,12 @@ use std::ops::{Add, Sub, Mul, Div};
 pub struct Poly(pub Vec<Number>);
 
 impl Poly {
+    fn mono(degree: usize, coefficient: Number) -> Poly {
+        let mut p = vec![0.0; degree + 1];
+        p[degree] = coefficient;
+        Poly(p)
+    }
+
     fn degree(&self) -> usize {
         let mut i = self.0.len() - 1;
         while self.0[i] == 0.0 && i > 0 {
@@ -16,6 +22,10 @@ impl Poly {
         i
     }
 
+    fn lc(&self) -> Number {
+        self.0[self.degree()]
+    }
+
     fn iter(&self) -> Iter {
         Iter::new(self.0.clone(), self.degree())
     }
@@ -75,17 +85,13 @@ impl Mul for &Poly {
 impl Div for &Poly {
     type Output = (Poly, Poly);
 
-    fn div(self, other: Self) -> (Poly, Poly) {
+    fn div(self, divisor: Self) -> (Poly, Poly) {
         let mut q = Poly(vec![0.0]);
         let mut r = self.clone();
-        let d = other.degree();
-        let c = other.0[d];
-        while r.degree() >= d {
-            let mut p = vec![0.0; r.degree() - d + 1];
-            p[r.degree() - d] = r.0[r.degree()] / c;
-            let s = Poly(p);
+        while !r.is_zero() && r.degree() >= divisor.degree() {
+            let s = Poly::mono(r.degree() - divisor.degree(), r.lc() / divisor.lc());
             q = &q + &s;
-            r = &r - &(&s * other);
+            r = &r - &(&s * divisor);
         }
         (q, r)
     }
@@ -93,7 +99,6 @@ impl Div for &Poly {
 
 pub fn gcd(a: &Poly, b: &Poly) -> Poly {
     let (_, c) = a / b;
-    println!("a: {:?}, b: {:?}, c: {:?}", a, b, c);
     if c.is_zero() {
         b.clone()
     } else {
@@ -175,8 +180,8 @@ mod tests {
 
     #[test]
     fn gcd_test() {
-        let a = Poly(vec![6.0, 7.0, 1.0]);
-        let b = Poly(vec![-6.0, -5.0, 1.0]);
-        assert_eq!(gcd(&a, &b), Poly(vec![12.0, 12.0, 0.0]));
+        let a = Poly(vec![2.0, -8.0, 8.0]);
+        let b = Poly(vec![-1.0, 4.0, -1.0]);
+        assert_eq!(gcd(&a, &b), Poly(vec![-0.0625, 0.0]));
     }
 }