math: polynomial division.

src/lol/math/polynomial.h

@@ -91,7 +91,7 @@ struct polynomial
      * as the value instead of a scalar. */
     template<typename U> U eval(U x) const
     {
-        U ret((*this)[degree()]);
+        U ret(leading());
         for (int i = degree() - 1; i >= 0; --i)
             ret = ret * x + U(m_coefficients[i]);
         return ret;
@@ -166,6 +166,12 @@ struct polynomial
         return m_coefficients[n];
     }
 
+    /* Return the leading coefficient */
+    inline T leading() const
+    {
+        return (*this)[degree()];
+    }
+
     /* Unary plus */
     polynomial<T> operator +() const
     {
@@ -264,6 +270,30 @@ struct polynomial
         return ret;
     }
 
+    /* Divide a polynomial by another one. There is no /= variant because
+     * the return value contains both the quotient and the remainder. */
+    tuple<polynomial<T>, polynomial<T>> operator /(polynomial<T> p) const
+    {
+        ASSERT(p.degree() >= 0);
+
+        tuple<polynomial<T>, polynomial<T>> ret;
+        polynomial<T> &quotient = ret.m1;
+        polynomial<T> &remainder = ret.m2;
+
+        remainder = *this / p.leading();
+        p /= p.leading();
+
+        for (int n = remainder.degree() - p.degree(); n >= 0; --n)
+        {
+            quotient.set(n, remainder.leading());
+            remainder.m_coefficients.Pop();
+            for (int i = 0; i < p.degree(); ++i)
+                remainder.m_coefficients[n + i] -= remainder.leading() * p[i];
+        }
+
+        return ret;
+    }
+
 private:
     /* Enforce the non-zero leading coefficient rule. */
     void reduce_degree()

src/t/math/polynomial.cpp

@@ -183,6 +183,22 @@ lolunit_declare_fixture(PolynomialTest)
         lolunit_assert_equal(r[3], 15.f);
     }
 
+    lolunit_declare_test(Division)
+    {
+        /* p(x) = -4 - 2x² + x³ */
+        /* q(x) = -3 + x */
+        polynomial<float> p { -4.f, 0.f, -2.f, 1.f };
+        polynomial<float> q { -3.f, 1.f };
+
+        auto r = p / q;
+        lolunit_assert_equal(r.m1.degree(), 2);
+        lolunit_assert_doubles_equal(r.m1[0], 3.f, 1e-5f);
+        lolunit_assert_doubles_equal(r.m1[1], 1.f, 1e-5f);
+        lolunit_assert_doubles_equal(r.m1[2], 1.f, 1e-5f);
+        lolunit_assert_equal(r.m2.degree(), 0);
+        lolunit_assert_doubles_equal(r.m2[0], 5.f, 1e-5f);
+    }
+
     lolunit_declare_test(Composition1)
     {
         /* p(x) = 1 + x² */