diff --git a/modules/core/src/mathfuncs.cpp b/modules/core/src/mathfuncs.cpp index 90713fb27d..7108f2e81d 100644 --- a/modules/core/src/mathfuncs.cpp +++ b/modules/core/src/mathfuncs.cpp @@ -46,6 +46,8 @@ #include #include #include +#include +#include #include "mathfuncs.hpp" namespace cv @@ -1586,12 +1588,21 @@ int cv::solveCubic( InputArray _coeffs, OutputArray _roots ) a3 = coeffs.at(i+3); } - if( a0 == 0 ) + // Fix for Issue #27748: Normalize coefficients to avoid overflow/underflow + // and correctly detect negligible cubic terms. + double max_coeff = std::max({std::abs(a0), std::abs(a1), std::abs(a2), std::abs(a3)}); + + // If max_coeff is effectively zero, the equation is 0=0 + if (max_coeff < std::numeric_limits::epsilon()) + return -1; + + // Check if the cubic term is negligible relative to the largest coefficient + if( std::abs(a0) < std::numeric_limits::epsilon() * max_coeff ) { - if( a1 == 0 ) + if( std::abs(a1) < std::numeric_limits::epsilon() * max_coeff ) { - if( a2 == 0 ) // constant - n = a3 == 0 ? -1 : 0; + if( std::abs(a2) < std::numeric_limits::epsilon() * max_coeff ) + n = std::abs(a3) < std::numeric_limits::epsilon() * max_coeff ? -1 : 0; else { // linear equation @@ -1608,7 +1619,7 @@ int cv::solveCubic( InputArray _coeffs, OutputArray _roots ) d = std::sqrt(d); double q1 = (-a2 + d) * 0.5; double q2 = (a2 + d) * -0.5; - if( fabs(q1) > fabs(q2) ) + if( std::abs(q1) > std::abs(q2) ) { x0 = q1 / a1; x1 = a3 / q1; @@ -1625,6 +1636,13 @@ int cv::solveCubic( InputArray _coeffs, OutputArray _roots ) else { // cubic equation + // Normalize coefficients in-place to prevent overflow/instability + double scale = 1.0 / max_coeff; + a0 *= scale; + a1 *= scale; + a2 *= scale; + a3 *= scale; + a0 = 1./a0; a1 *= a0; a2 *= a0; @@ -1633,6 +1651,7 @@ int cv::solveCubic( InputArray _coeffs, OutputArray _roots ) double Q = (a1 * a1 - 3 * a2) * (1./9); double R = (a1 * (2 * a1 * a1 - 9 * a2) + 27 * a3) * (1./54); double Qcubed = Q * Q * Q; + /* Here we expand expression `Qcubed - R * R` for `d` variable to reduce common terms `a1^6 / 729` and `-a1^4 * a2 / 81` diff --git a/modules/core/test/test_math.cpp b/modules/core/test/test_math.cpp index c6686cf1c4..c7c5dceb31 100644 --- a/modules/core/test/test_math.cpp +++ b/modules/core/test/test_math.cpp @@ -2618,6 +2618,33 @@ TEST(Core_SolveCubic, regression_27323) } } +TEST(Core_SolveCubic, regression_27748) +{ + // a is extremely small relative to others (approx 1.8e-19 ratio), + // causing instability in standard cubic formula. + double a = 1.56041e-17; + double b = 84.4504; + double c = -96.795; + double d = 13.6826; + + Mat coeffs = (Mat_(1, 4) << a, b, c, d); + Mat roots; + + int n = solveCubic(coeffs, roots); + + // Expecting quadratic behavior (2 roots) + EXPECT_GE(n, 2); + + // Verify roots satisfy the FULL cubic equation + for(int i = 0; i < n; i++) + { + double x = roots.at(i); + double val = a*x*x*x + b*x*x + c*x + d; + // Check residual is small + EXPECT_LE(fabs(val), 1e-6) << "Root " << x << " does not satisfy the equation"; + } +} + TEST(Core_SolvePoly, regression_5599) { // x^4 - x^2 = 0, roots: 1, -1, 0, 0