Implement div and rem for Poly

Using the long division algorithm.

1 files changed, 42 insertions, 17 deletions

diff --git a/src/poly.rs b/src/poly.rs
index 0e1cd9e..488b0ea 100644
--- a/src/poly.rs
+++ b/src/poly.rs
@@ -3,9 +3,9 @@ mod iter;
 use crate::number::Number;
 use iter::Iter;
 use std::cmp::Ordering;
-use std::ops::{Add, Sub, Mul, Rem};
+use std::ops::{Add, Sub, Mul, Div, Rem};
 
-#[derive(PartialEq, Debug)]
+#[derive(PartialEq, Debug, Clone)]
 pub struct Poly(pub Vec<Number>);
 
 impl Poly {
@@ -77,26 +77,44 @@ impl Rem for &Poly {
     type Output = Poly;
 
     fn rem(self, other: Self) -> Poly {
-        let ad = self.degree();
-        let bd = other.degree();
-        let an = self.0[ad];
-        let bn = other.0[bd];
-        match ad.cmp(&bd) {
-            Ordering::Equal => self - &(other * &Poly(vec![an / bn])),
-            Ordering::Greater => {
-                let mut qv = vec![0.0; ad + 1];
-                qv[ad - bd] = an / bn;
-                self - &(other * &Poly(qv))
-            },
-            Ordering::Less => Poly(vec![0.0])
+        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);
+            q = &q + &s;
+            r = &r - &(&s * other);
         }
+        r
+    }
+}
+
+impl Div for &Poly {
+    type Output = Poly;
+
+    fn div(self, other: Self) -> 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);
+            q = &q + &s;
+            r = &r - &(&s * other);
+        }
+        q
     }
 }
 
 pub fn gcd(a: &Poly, b: &Poly) -> Poly {
     let c = a % b;
     if c.is_zero() {
-        Poly(b.0.clone())
+        b.clone()
     } else {
         gcd(b, &c)
     }
@@ -107,6 +125,13 @@ mod tests {
     use super::*;
 
     #[test]
+    fn div_test() {
+        let a = Poly(vec![1.0, 4.0, 7.0, 6.0]);
+        let b = Poly(vec![1.0, 2.0, 3.0]);
+        assert_eq!(&a / &b, Poly(vec![1.0, 2.0]));
+    }
+
+    #[test]
     fn mul_test() {
         let a = Poly(vec![1.0, 2.0, 3.0]);
         let b = Poly(vec![1.0, 2.0]);
@@ -164,13 +189,13 @@ mod tests {
     fn rem_greater() {
         let a = Poly(vec![-6.0, -5.0, 1.0]);
         let b = Poly(vec![12.0, 12.0]);
-        assert_eq!(&a % &b, Poly(vec![-6.0, -6.0, 0.0]));
+        assert_eq!(&a % &b, Poly(vec![0.0, 0.0]));
     }
 
     #[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![-6.0, -6.0, 0.0]));
+        assert_eq!(gcd(&a, &b), Poly(vec![12.0, 12.0, 0.0]));
     }
 }