From 1c93dabbad95861d30a6c63feae25976b31742f3 Mon Sep 17 00:00:00 2001
From: Guillaume Bittoun <guillaume.bittoun@gmail.com>
Date: Mon, 11 May 2015 21:11:28 +0000
Subject: [PATCH] polynomial: more 3rd order fixes

---
 src/lol/math/polynomial.h | 15 ++++++++++-----
 src/t/math/polynomial.cpp |  4 ++--
 2 files changed, 12 insertions(+), 7 deletions(-)

diff --git a/src/lol/math/polynomial.h b/src/lol/math/polynomial.h
index 19692388..e5ccc9c7 100644
--- a/src/lol/math/polynomial.h
+++ b/src/lol/math/polynomial.h
@@ -20,7 +20,10 @@
 
 #include <functional>
 
-#if 1
+// #define ENABLE_3SOLVE
+// #include <iostream>
+
+#ifdef ENABLE_3SOLVE
 #include <complex>
 #endif
 
@@ -169,7 +172,7 @@ struct polynomial
             /* Using x = (X - k) so that p2(X) = p(X - k) = aX³ + 0×X² + mX + n */
             T const k = b / (3 * a);
             T const m = 3 * a * k * k - 2 * b * k + c;
-            T const n = -a * k * k * k + b * k * k + c * k + d;
+            T const n = -a * k * k * k + b * k * k - c * k + d;
 
             std::cout << "k,m,n: " << k << ", " << m << ", " << n << std::endl;
 
@@ -184,8 +187,8 @@ struct polynomial
              * u³ - v³ = √((n/a)² + 4m³/27)
              * u³ + v³ = -n/a
              *
-             * u³ = -n/2a + √((n/a)² + 4m³/27)/2
-             * v³ = -n/2a - √((n/a)² + 4m³/27)/2
+             * u³ = (-n/a + √((n/a)² + 4m³/27))/2
+             * v³ = (-n/a - √((n/a)² + 4m³/27))/2
              */
             T const delta = (n * n) / (a * a) + 4 * m * m * m / 27;
 
@@ -229,9 +232,11 @@ struct polynomial
                     pow(v3_norm, 1.f / 3.f) * std::complex<T>(cos(v_angle), sin(v_angle));
             }
 
+            std::cout << "delta: " << delta << std::endl;
+
             std::cout << "complex_solutions: " << complex_solutions[0] << ", " << complex_solutions[1] << ", " << complex_solutions[2] << std::endl;
 
-            return array<T> {complex_solutions[0].real(), complex_solutions[1].real(), complex_solutions[2].real()};
+            return array<T> {complex_solutions[0].real() - k, complex_solutions[1].real() - k, complex_solutions[2].real() - k};
         }
 #endif
 
diff --git a/src/t/math/polynomial.cpp b/src/t/math/polynomial.cpp
index ee9512f0..5a2b7ae7 100644
--- a/src/t/math/polynomial.cpp
+++ b/src/t/math/polynomial.cpp
@@ -275,7 +275,7 @@ lolunit_declare_fixture(PolynomialTest)
 #ifdef ENABLE_3SOLVE // Development in progress
     lolunit_declare_test(RootsDegree3)
     {
-        polynomial<float> p { 2.f, 0.f, 0.f, 1.f };
+        polynomial<float> p { 1.f, 0.f, 0.f, 1.f };
         auto roots1 = p.roots();
 
         lolunit_assert_equal(roots1.count(), 3);
@@ -285,7 +285,7 @@ lolunit_declare_fixture(PolynomialTest)
 
     lolunit_declare_test(RootsDegree3_2)
     {
-        polynomial<float> p { -1.f, 0.f, 0.f, 1.f };
+        polynomial<float> p { 0.f, -2.f, 0.f, 1.f };
         auto roots1 = p.roots();
 
         lolunit_assert_equal(roots1.count(), 3);