From 3916a44c589e04170b7cd20772cfb784553709c1 Mon Sep 17 00:00:00 2001
From: Sam Hocevar <sam@hocevar.net>
Date: Mon, 8 Dec 2014 14:11:01 +0000
Subject: [PATCH] math: add derive() method to polynomial and allow to call
 eval() with polynomials as arguments so as to compose them together.

---
 src/lol/math/polynomial.h | 25 +++++++++++++++---
 src/t/math/polynomial.cpp | 54 +++++++++++++++++++++++++++++++++++++++
 2 files changed, 76 insertions(+), 3 deletions(-)

diff --git a/src/lol/math/polynomial.h b/src/lol/math/polynomial.h
index 48f182a4..09e988bd 100644
--- a/src/lol/math/polynomial.h
+++ b/src/lol/math/polynomial.h
@@ -26,6 +26,13 @@ struct polynomial
     /* The zero polynomial */
     explicit inline polynomial() {}
 
+    /* A constant polynomial */
+    explicit inline polynomial(T const &a)
+    {
+        m_coefficients.Push(a);
+        reduce_degree();
+    }
+
     /* Create a polynomial from a list of coefficients */
     explicit polynomial(std::initializer_list<T> const &init)
     {
@@ -76,11 +83,23 @@ struct polynomial
         reduce_degree();
     }
 
-    T eval(T x) const
+    /* Evaluate polynomial at a given value. This method can also
+     * be used to compose polynomials, i.e. using another polynomial
+     * as the value instead of a scalar. */
+    template<typename U> U eval(U x) const
     {
-        T ret = this->operator[](degree());
+        U ret((*this)[degree()]);
         for (int i = degree() - 1; i >= 0; --i)
-            ret = ret * x + m_coefficients[i];
+            ret = ret * x + U(m_coefficients[i]);
+        return ret;
+    }
+
+    polynomial<T> derive() const
+    {
+        /* No need to reduce the degree after deriving. */
+        polynomial<T> ret;
+        for (int i = 1; i <= degree(); ++i)
+            ret.m_coefficients.Push(m_coefficients[i] * T(i));
         return ret;
     }
 
diff --git a/src/t/math/polynomial.cpp b/src/t/math/polynomial.cpp
index 72de936b..0757aab9 100644
--- a/src/t/math/polynomial.cpp
+++ b/src/t/math/polynomial.cpp
@@ -50,6 +50,29 @@ lolunit_declare_fixture(PolynomialTest)
         lolunit_assert_equal(s.degree(), -1);
     }
 
+    lolunit_declare_test(Derive)
+    {
+        polynomial<float> p {};
+        p = p.derive();
+        lolunit_assert_equal(p.degree(), -1);
+
+        polynomial<float> q { 1.f };
+        q = q.derive();
+        lolunit_assert_equal(q.degree(), -1);
+
+        polynomial<float> r { 1.f, 2.f };
+        r = r.derive();
+        lolunit_assert_equal(r.degree(), 0);
+        lolunit_assert_equal(r[0], 2.f);
+
+        polynomial<float> s { 1.f, 2.f, 3.f, 4.f };
+        s = s.derive();
+        lolunit_assert_equal(s.degree(), 2);
+        lolunit_assert_equal(s[0], 2.f);
+        lolunit_assert_equal(s[1], 6.f);
+        lolunit_assert_equal(s[2], 12.f);
+    }
+
     lolunit_declare_test(Eval)
     {
         /* Special null polynomial */
@@ -160,6 +183,37 @@ lolunit_declare_fixture(PolynomialTest)
         lolunit_assert_equal(r[3], 15.f);
     }
 
+    lolunit_declare_test(Composition1)
+    {
+        /* p(x) = 1 + x² */
+        polynomial<float> p({ 1, 0, 1 });
+
+        /* q(x) = (p o p)(x) = 2 + 2x² + x⁴ */
+        polynomial<float> q = p.eval(p);
+        lolunit_assert_equal(q.degree(), 4);
+        lolunit_assert_equal(q[0], 2.f);
+        lolunit_assert_equal(q[1], 0.f);
+        lolunit_assert_equal(q[2], 2.f);
+        lolunit_assert_equal(q[3], 0.f);
+        lolunit_assert_equal(q[4], 1.f);
+    }
+
+    lolunit_declare_test(Composition2)
+    {
+        /* p(x) = 1 + x */
+        polynomial<float> p({ 1, 1 });
+
+        /* q(x) = 1 + x + x² */
+        polynomial<float> q({ 1, 1, 1 });
+
+        /* r(x) = (q o p)(x) = 3 + 3x + x² */
+        polynomial<float> r = q.eval(p);
+        lolunit_assert_equal(r.degree(), 2);
+        lolunit_assert_equal(r[0], 3.f);
+        lolunit_assert_equal(r[1], 3.f);
+        lolunit_assert_equal(r[2], 1.f);
+    }
+
     lolunit_declare_test(Chebyshev)
     {
         polynomial<float> t0 = polynomial<float>::chebyshev(0);