From 93dc8a70a70a20fc3970a891c9f0c213b3f4637d Mon Sep 17 00:00:00 2001
From: Guillaume Bittoun <guillaume.bittoun@gmail.com>
Date: Wed, 13 May 2015 07:19:28 +0000
Subject: [PATCH] polynomial: 3rd order solving v1.0

---
 src/lol/math/polynomial.h | 30 ++++++++----------------------
 src/t/math/polynomial.cpp | 12 +++++++++---
 2 files changed, 17 insertions(+), 25 deletions(-)

diff --git a/src/lol/math/polynomial.h b/src/lol/math/polynomial.h
index 5695ecb1..9d130c99 100644
--- a/src/lol/math/polynomial.h
+++ b/src/lol/math/polynomial.h
@@ -20,8 +20,6 @@
 
 #include <functional>
 
-// #define ENABLE_3SOLVE
-// #include <iostream>
 
 namespace lol
 {
@@ -111,14 +109,9 @@ struct polynomial
 
     array<T> roots() const
     {
-        /* For now we can only solve polynomials of degrees 0, 1 or 2. */
-#ifdef ENABLE_3SOLVE
+        /* For now we can only solve polynomials of degrees 0, 1, 2 or 3. */
         ASSERT(degree() >= 0 && degree() <= 3,
                "roots() called on polynomial of degree %d", degree());
-#else
-        ASSERT(degree() >= 0 && degree() <= 2,
-               "roots() called on polynomial of degree %d", degree());
-#endif
 
         if (degree() == 0)
         {
@@ -156,7 +149,6 @@ struct polynomial
                  return array<T> { -k };
             }
         }
-#ifdef ENABLE_3SOLVE
         else if (degree() == 3)
         {
             /* p(x) = ax³ + bx² + cx + d */
@@ -226,21 +218,15 @@ struct polynomial
                     pow(v3_norm, 1.f / 3.f) * cos(v_angle);
             }
 
-            if (delta < 0)
+            if (delta < 0) // 3 real solutions
                 return array<T> {solutions[0] - k, solutions[1] - k, solutions[2] - k};
-            else
-            {
-                if (u3_norm > 0)
-                {
-                    return array<T> {solutions[0] - k, solutions[1] - k};
-                }
-                else
-                {
-                    return array<T> {solutions[0] - k};
-                }
-            }
+            else if (delta > 0) // 1 real solution
+                return array<T> {solutions[0] - k};
+            else if (u3_norm > 0) // 2 real solutions
+                return array<T> {solutions[0] - k, solutions[1] - k};
+            else // one triple solution
+                return array<T> {solutions[0] - k};
         }
-#endif
 
         /* It is an error to reach this point. */
         ASSERT(false);
diff --git a/src/t/math/polynomial.cpp b/src/t/math/polynomial.cpp
index fd02772d..c4884d32 100644
--- a/src/t/math/polynomial.cpp
+++ b/src/t/math/polynomial.cpp
@@ -272,13 +272,13 @@ lolunit_declare_fixture(PolynomialTest)
         lolunit_assert_equal(roots2[1], 7.f);
     }
 
-#ifdef ENABLE_3SOLVE
     lolunit_declare_test(RootsDegree3TripleSolution)
     {
         polynomial<float> p { 1.f, 3.f, 3.f, 1.f };
         auto roots1 = p.roots();
 
-        lolunit_assert_equal(roots1.count(), 3);
+        lolunit_assert_equal(roots1.count(), 1);
+        lolunit_assert_doubles_equal(roots1[0], -1, 0);
     }
 
     lolunit_declare_test(RootsDegree3DoubleSolution)
@@ -286,7 +286,11 @@ lolunit_declare_fixture(PolynomialTest)
         polynomial<float> p { 2.f, 5.f, 4.f, 1.f };
         auto roots1 = p.roots();
 
+        // Should have 2 solutions only, but precision leads to 3 solutions
         lolunit_assert_equal(roots1.count(), 3);
+        lolunit_assert_doubles_equal(roots1[0], -1, 1e-3);
+        lolunit_assert_doubles_equal(roots1[1], -2, 1e-6);
+        lolunit_assert_doubles_equal(roots1[2], -1, 1e-3);
     }
 
     lolunit_declare_test(RootsDegree3SingleSolutions)
@@ -295,8 +299,10 @@ lolunit_declare_fixture(PolynomialTest)
         auto roots1 = p.roots();
 
         lolunit_assert_equal(roots1.count(), 3);
+        lolunit_assert_doubles_equal(roots1[0], -1, 1e-8);
+        lolunit_assert_doubles_equal(roots1[1], -3, 1e-8);
+        lolunit_assert_doubles_equal(roots1[2], -2, 1e-8);
     }
-#endif
 
     lolunit_declare_test(Chebyshev)
     {