From 023d14ecc412c1c475c5756f6e0f5d7c0ba46648 Mon Sep 17 00:00:00 2001 From: Maxim Smolskiy Date: Sat, 24 May 2025 09:56:48 +0300 Subject: [PATCH] Merge pull request #27347 from MaximSmolskiy:improve-solveCubic-accuracy Improve solveCubic accuracy #27347 ### Pull Request Readiness Checklist Fix #27323 ``` 2e-13 * x^3 + x^2 - 2 * x + 1 = 0 -> x^3 + 5e12 * x^2 - 1e13 * x + 5e12 = 0 ``` The problem that coefficients have quite big magnitudes and current calculations are subject to round-off error ``` Q = (a1 * a1 - 3 * a2) * (1./9) R = (2 * a1 * a1 * a1 - 9 * a1 * a2 + 27 * a3) * (1./54) Qcubed = Q * Q * Q = a1^6/729 - (a1^4 a2)/81 + (a1^2 a2^2)/27 - a2^3/27 R * R = R^2 = a1^6/729 - (a1^4 a2)/81 + (a1^2 a2^2)/36 + (a1^3 a3)/27 - (a1 a2 a3)/6 + a3^2/4 d = Qcubed - R * R ``` Let `a1`, `a2`, `a3` have quite big same magnitudes, then we see that `Qcubed` and `R * R` have same terms `a1^6/729` and `-(a1^4 a2)/81` (which will be reduced in `d`), but they level out the other terms (these terms have `6`th and `5`th degree and other terms - less or equal than `4`th degree). So, if these terms will participate in the calculation, this will lead to a huge round-off error. But if we expand the expression, then round-off error should be less ``` d = Qcubed - R * R = 1/108 (a1^2 a2^2 - 4 a2^3 - 4 a1^3 a3 + 18 a1 a2 a3 - 27 a3^2) ``` See details at https://github.com/opencv/opencv/wiki/How_to_contribute#making-a-good-pull-request - [x] I agree to contribute to the project under Apache 2 License. - [x] To the best of my knowledge, the proposed patch is not based on a code under GPL or another license that is incompatible with OpenCV - [x] The PR is proposed to the proper branch - [x] There is a reference to the original bug report and related work - [x] There is accuracy test, performance test and test data in opencv_extra repository, if applicable Patch to opencv_extra has the same branch name. - [x] The feature is well documented and sample code can be built with the project CMake --- modules/core/src/mathfuncs.cpp | 11 +++++++++-- modules/core/test/test_math.cpp | 25 +++++++++++++++++++++++-- 2 files changed, 32 insertions(+), 4 deletions(-) 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