diff --git a/modules/core/src/mathfuncs.cpp b/modules/core/src/mathfuncs.cpp index 7af66c5723..90713fb27d 100644 --- a/modules/core/src/mathfuncs.cpp +++ b/modules/core/src/mathfuncs.cpp @@ -1631,9 +1631,16 @@ int cv::solveCubic( InputArray _coeffs, OutputArray _roots ) a3 *= a0; double Q = (a1 * a1 - 3 * a2) * (1./9); - double R = (2 * a1 * a1 * a1 - 9 * a1 * a2 + 27 * a3) * (1./54); + double R = (a1 * (2 * a1 * a1 - 9 * a2) + 27 * a3) * (1./54); double Qcubed = Q * Q * Q; - double d = Qcubed - R * R; + /* + Here we expand expression `Qcubed - R * R` for `d` variable + to reduce common terms `a1^6 / 729` and `-a1^4 * a2 / 81` + and thus decrease rounding error (in case of quite big coefficients). + + And then we additionally group terms to further reduce rounding error. + */ + double d = (a1 * a1 * (a2 * a2 - 4 * a1 * a3) + 2 * a2 * (9 * a1 * a3 - 2 * a2 * a2) - 27 * a3 * a3) * (1./108); if( d > 0 ) { diff --git a/modules/core/test/test_math.cpp b/modules/core/test/test_math.cpp index 6ff4d4984d..c6686cf1c4 100644 --- a/modules/core/test/test_math.cpp +++ b/modules/core/test/test_math.cpp @@ -2522,7 +2522,7 @@ TEST(Core_SolveCubicCubic, accuracy) const auto num_roots = solveCubic(coeffs, roots); EXPECT_EQ(num_roots, 1); - EXPECT_EQ(roots[0], 1.); + EXPECT_NEAR(roots[0], 1., 1e-8); } { @@ -2568,7 +2568,7 @@ TEST(Core_SolveCubicNormalizedCubic, accuracy) const auto num_roots = solveCubic(coeffs, roots); EXPECT_EQ(num_roots, 1); - EXPECT_EQ(roots[0], 1.); + EXPECT_NEAR(roots[0], 1., 1e-8); } { @@ -2597,6 +2597,27 @@ TEST(Core_SolveCubicNormalizedCubic, accuracy) } } +TEST(Core_SolveCubic, regression_27323) +{ + { + const std::vector coeffs{2e-13, 1, -2, 1}; + std::vector roots; + const auto num_roots = solveCubic(coeffs, roots); + + EXPECT_EQ(num_roots, 1); + EXPECT_EQ(roots[0], -5e12 - 2.); + } + + { + const std::vector coeffs{5e12, -1e13, 5e12}; + std::vector roots; + const auto num_roots = solveCubic(coeffs, roots); + + EXPECT_EQ(num_roots, 1); + EXPECT_EQ(roots[0], -5e12 - 2.); + } +} + TEST(Core_SolvePoly, regression_5599) { // x^4 - x^2 = 0, roots: 1, -1, 0, 0