From 1eb5f9f361e0d90bc9b78d568a4db3c318173859 Mon Sep 17 00:00:00 2001
From: Sam Hocevar <sam@hocevar.net>
Date: Mon, 8 Dec 2014 14:11:05 +0000
Subject: [PATCH] math: add a roots() method to find polynomial roots for
 degrees 1 and 2.

---
 src/lol/math/polynomial.h | 46 +++++++++++++++++++++++++++++++++++++++
 src/t/math/polynomial.cpp | 37 +++++++++++++++++++++++++++++++
 2 files changed, 83 insertions(+)

diff --git a/src/lol/math/polynomial.h b/src/lol/math/polynomial.h
index 09e988bd..b3463612 100644
--- a/src/lol/math/polynomial.h
+++ b/src/lol/math/polynomial.h
@@ -103,6 +103,52 @@ struct polynomial
         return ret;
     }
 
+    array<T> roots() const
+    {
+        ASSERT(degree() >= 0, "roots() called on the null polynomial");
+
+        if (degree() == 0)
+        {
+             /* p(x) = a > 0 */
+             return array<T> {};
+        }
+        else if (degree() == 1)
+        {
+             /* p(x) = ax + b */
+             T const &a = m_coefficients[1];
+             T const &b = m_coefficients[0];
+             return array<T> { -b / a };
+        }
+        else if (degree() == 2)
+        {
+             /* p(x) = ax² + bx + c */
+             T const &a = m_coefficients[2];
+             T const &b = m_coefficients[1];
+             T const &c = m_coefficients[0];
+
+             T const k = b / (a + a);
+             T const delta = k * k - c / a;
+
+             if (delta < T(0))
+             {
+                 return array<T> {};
+             }
+             else if (delta > T(0))
+             {
+                 T const sqrt_delta = sqrt(delta);
+                 return array<T> { -k - sqrt_delta, -k + sqrt_delta };
+             }
+             else
+             {
+                 return array<T> { -k };
+             }
+        }
+
+        ASSERT(false, "roots() called on polynomial of degree %d > 2",
+                      degree());
+        return array<T> {};
+    }
+
     inline T operator[](ptrdiff_t n) const
     {
         if (n < 0 || n > degree())
diff --git a/src/t/math/polynomial.cpp b/src/t/math/polynomial.cpp
index 0757aab9..9fcb85de 100644
--- a/src/t/math/polynomial.cpp
+++ b/src/t/math/polynomial.cpp
@@ -214,6 +214,43 @@ lolunit_declare_fixture(PolynomialTest)
         lolunit_assert_equal(r[2], 1.f);
     }
 
+    lolunit_declare_test(RootsDegree0)
+    {
+        /* p(x) = 42 */
+        polynomial<float> p { 42.f };
+        auto roots = p.roots();
+
+        lolunit_assert_equal(roots.Count(), 0);
+    }
+
+    lolunit_declare_test(RootsDegree1)
+    {
+        /* p(x) = -6 + 2x */
+        polynomial<float> p { -6.f, 2.f };
+        auto roots = p.roots();
+
+        lolunit_assert_equal(roots.Count(), 1);
+        lolunit_assert_equal(roots[0], 3.f);
+    }
+
+    lolunit_declare_test(RootsDegree2)
+    {
+        /* p(x) = 81 - 18x + x² */
+        polynomial<float> p { 81.f, -18.f, 1.f };
+        auto roots1 = p.roots();
+
+        lolunit_assert_equal(roots1.Count(), 1);
+        lolunit_assert_equal(roots1[0], 9.f);
+
+        /* p(x) = 42 - 20x + 2x² */
+        polynomial<float> q { 42.f, -20.f, 2.f };
+        auto roots2 = q.roots();
+
+        lolunit_assert_equal(roots2.Count(), 2);
+        lolunit_assert_equal(roots2[0], 3.f);
+        lolunit_assert_equal(roots2[1], 7.f);
+    }
+
     lolunit_declare_test(Chebyshev)
     {
         polynomial<float> t0 = polynomial<float>::chebyshev(0);