diff --git a/src/lol/math/polynomial.h b/src/lol/math/polynomial.h index 986b3130..d5a306b5 100644 --- a/src/lol/math/polynomial.h +++ b/src/lol/math/polynomial.h @@ -225,20 +225,22 @@ struct polynomial solutions[i] = u_norm * cos(u_angle) + v_norm * cos(v_angle); } - if (delta < 0) // 3 real solutions - return array { solutions[0] - k, - solutions[1] - k, - solutions[2] - k }; - - if (delta > 0) // 1 real solution + if (a == d && b == c && b == T(3)*a) // triple solution return array { solutions[0] - k }; - if (u_norm > 0) // 2 real solutions - return array { solutions[0] - k, - solutions[1] - k }; + // if root of the derivative is also root of the current polynomial, we have a double root. + for (auto root : polynomial{c, T(2)*b, T(3)*a}.roots()) + if (eval(root) == T(0)) + return array { solutions[0] - k, + solutions[1] - k }; - // one triple solution - return array { solutions[0] - k }; + // we have 3 or 1 root depending on delta sign + if (delta > 0) + return array { solutions[0] - k }; + else // if (delta < 0) 3 real solutions + return array { solutions[0] - k, + solutions[1] - k, + solutions[2] - k }; } /* It is an error to reach this point. */ diff --git a/src/t/math/polynomial.cpp b/src/t/math/polynomial.cpp index d431d6f8..86892cf6 100644 --- a/src/t/math/polynomial.cpp +++ b/src/t/math/polynomial.cpp @@ -287,10 +287,9 @@ lolunit_declare_fixture(polynomial_test) auto roots1 = p.roots(); // Should have 2 solutions only, but precision leads to 3 solutions - lolunit_assert_equal(roots1.count(), 3); + lolunit_assert_equal(roots1.count(), 2); 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(degree_3_three_roots)