Merge pull request #27347 from MaximSmolskiy:improve-solveCubic-accuracy

MaximSmolskiy · web-flow · commit 023d14ecc412 · 2025-05-24T09:56:48.000+03:00
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
diff --git 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
@@ -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<double> coeffs{2e-13, 1, -2, 1};
+        std::vector<double> roots;
+        const auto num_roots = solveCubic(coeffs, roots);
+
+        EXPECT_EQ(num_roots, 1);
+        EXPECT_EQ(roots[0], -5e12 - 2.);
+    }
+
+    {
+        const std::vector<double> coeffs{5e12, -1e13, 5e12};
+        std::vector<double> 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

Original file line number	Diff line number	Diff line change
`@@ -2522,7 +2522,7 @@ TEST(Core_SolveCubicCubic, accuracy)`
`2522`	`2522`	`const auto num_roots = solveCubic(coeffs, roots);`
`2523`	`2523`
`2524`	`2524`	`EXPECT_EQ(num_roots, 1);`
`2525`		`- EXPECT_EQ(roots[0], 1.);`
	`2525`	`+ EXPECT_NEAR(roots[0], 1., 1e-8);`
`2526`	`2526`	`}`
`2527`	`2527`
`2528`	`2528`	`{`
`@@ -2568,7 +2568,7 @@ TEST(Core_SolveCubicNormalizedCubic, accuracy)`
`2568`	`2568`	`const auto num_roots = solveCubic(coeffs, roots);`
`2569`	`2569`
`2570`	`2570`	`EXPECT_EQ(num_roots, 1);`
`2571`		`- EXPECT_EQ(roots[0], 1.);`
	`2571`	`+ EXPECT_NEAR(roots[0], 1., 1e-8);`
`2572`	`2572`	`}`
`2573`	`2573`
`2574`	`2574`	`{`
`@@ -2597,6 +2597,27 @@ TEST(Core_SolveCubicNormalizedCubic, accuracy)`
`2597`	`2597`	`}`
`2598`	`2598`	`}`
`2599`	`2599`
	`2600`	`+TEST(Core_SolveCubic, regression_27323)`
	`2601`	`+{`
	`2602`	`+ {`
	`2603`	`+ const std::vector<double> coeffs{2e-13, 1, -2, 1};`
	`2604`	`+ std::vector<double> roots;`
	`2605`	`+ const auto num_roots = solveCubic(coeffs, roots);`
	`2606`	`+`
	`2607`	`+ EXPECT_EQ(num_roots, 1);`
	`2608`	`+ EXPECT_EQ(roots[0], -5e12 - 2.);`
	`2609`	`+ }`
	`2610`	`+`
	`2611`	`+ {`
	`2612`	`+ const std::vector<double> coeffs{5e12, -1e13, 5e12};`
	`2613`	`+ std::vector<double> roots;`
	`2614`	`+ const auto num_roots = solveCubic(coeffs, roots);`
	`2615`	`+`
	`2616`	`+ EXPECT_EQ(num_roots, 1);`
	`2617`	`+ EXPECT_EQ(roots[0], -5e12 - 2.);`
	`2618`	`+ }`
	`2619`	`+}`
	`2620`	`+`
`2600`	`2621`	`TEST(Core_SolvePoly, regression_5599)`
`2601`	`2622`	`{`
`2602`	`2623`	`// x^4 - x^2 = 0, roots: 1, -1, 0, 0`