From c55f25f8216e983896649b2ac7e360a36dd592a3 Mon Sep 17 00:00:00 2001
From: Guillaume Bittoun <guillaume.bittoun@gmail.com>
Date: Sat, 9 May 2015 21:14:24 +0000
Subject: [PATCH] =?UTF-8?q?drafting=20polynomial=203rd=20order=20solving.?=
 =?UTF-8?q?=20To=20be=20continued=E2=80=A6?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/lol/math/polynomial.h | 79 +++++++++++++++++++++++++++++++++++++++
 src/t/math/polynomial.cpp | 12 ++++++
 2 files changed, 91 insertions(+)

diff --git a/src/lol/math/polynomial.h b/src/lol/math/polynomial.h
index 4f0ec570..5efdb87f 100644
--- a/src/lol/math/polynomial.h
+++ b/src/lol/math/polynomial.h
@@ -20,6 +20,11 @@
 
 #include <functional>
 
+#if 1
+#include <complex>
+#endif
+
+
 namespace lol
 {
 
@@ -109,8 +114,13 @@ struct polynomial
     array<T> roots() const
     {
         /* For now we can only solve polynomials of degrees 0, 1 or 2. */
+#ifdef ENABLE_3SOLVE
+        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)
         {
@@ -148,6 +158,75 @@ struct polynomial
                  return array<T> { -k };
             }
         }
+#ifdef ENABLE_3SOLVE // Development in progress
+        else if (degree() == 3)
+        {
+            /* p(x) = ax³ + bx² + cx + d */
+            T const &a = m_coefficients[3];
+            T const &b = m_coefficients[2];
+            T const &c = m_coefficients[1];
+            T const &d = m_coefficients[0];
+
+            /* 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;
+
+            std::cout << "k,m,n: " << k << ", " << m << ", " << n << std::endl;
+
+            /* Assuming X = u + v and 3uv = -m, then
+             * p2(u + v) = a(u + v) + n
+             *
+             * We then need to solve the following system:
+             *    u³v³ = -m³/27
+             * u³ + v³ = -n/a
+             *
+             * which gives :
+             * 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
+             */
+            T const delta = (m * m) / (a * a) + 4 * m * m * m / 27;
+
+            std::complex<T> u3, v3;
+
+            if (delta < 0)
+            {
+                u3 = std::complex<T>(-n/(2.0f*a),  sqrt(abs(delta)));
+                v3 = std::complex<T>(-n/(2.0f*a), -sqrt(abs(delta)));
+            }
+            else
+            {
+                u3 = std::complex<T>(-n/(2.0f*a) + sqrt(delta), 0);
+                v3 = std::complex<T>(-n/(2.0f*a) - sqrt(delta), 0);
+            }
+
+            std::cout << "delta,u3,v3: " << delta << ", " << u3 << "," << v3 << std::endl;
+
+            T const u3_angle = atan2(u3.imag(), u3.real());
+            T const v3_angle = atan2(v3.imag(), v3.real());
+
+            std::cout << "u3_angle,v3_angle: " << u3_angle << "," << v3_angle << std::endl;
+
+            std::complex<T> complex_solutions[3];
+
+            for (int i = 0 ; i < 3 ; ++i)
+            {
+                T u_angle = u3_angle / 3 + i * 2 * M_PI / 3;
+                T v_angle = v3_angle / 3 - i * 2 * M_PI / 3;
+
+                std::cout << i << " => u_angle,v_angle: " << u_angle << "," << v_angle << std::endl;
+
+                complex_solutions[i] =
+                    pow(abs(u3), 1.0f/3.0f) * std::complex<T>(cos(u_angle), sin(u_angle)) +
+                    pow(abs(v3), 1.0f/3.0f) * std::complex<T>(cos(v_angle), sin(v_angle));
+            }
+
+            return array<T> {complex_solutions[0].real(), complex_solutions[1].real(), complex_solutions[2].real()};
+        }
+#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 ffbccd03..622a55fe 100644
--- a/src/t/math/polynomial.cpp
+++ b/src/t/math/polynomial.cpp
@@ -272,6 +272,18 @@ lolunit_declare_fixture(PolynomialTest)
         lolunit_assert_equal(roots2[1], 7.f);
     }
 
+#ifdef ENABLE_3SOLVE // Development in progress
+    lolunit_declare_test(RootsDegree3)
+    {
+        polynomial<float> p { 2.f, 0.f, 0.f, 1.f };
+        auto roots1 = p.roots();
+
+        lolunit_assert_equal(roots1.count(), 3);
+
+        std::cout << roots1[0] << ", " << roots1[1] << ", " << roots1[2] << std::endl;
+    }
+#endif
+
     lolunit_declare_test(Chebyshev)
     {
         polynomial<float> t0 = polynomial<float>::chebyshev(0);