Merge pull request #28117 from satyam102006:fix/solveCubic-numerical-instability

core: fix solveCubic numerical instability via coefficient normalization (fixes #27748) #28117

 Summary
This PR fixes numerical instability in `cv::solveCubic` when the leading coefficient `a` is non-zero but extremely small relative to other coefficients (Issue #27748). 

It introduces a **normalization step** that scales all coefficients by their maximum magnitude before solving. This ensures robust detection of when the equation should degenerate to a quadratic solver, without breaking valid cubic equations that happen to have small coefficients (e.g., scaled by 1e-9).

The Problem (Issue #27748)
The previous implementation checked `if (a == 0)` to decide whether to use the cubic or quadratic formula. 
- When `a` is extremely small (e.g., 1e-17) but not exactly zero, and other coefficients are normal (e.g., 5.0), the standard cubic formula suffers from catastrophic cancellation and overflow, producing incorrect roots (e.g., 1e14).

 The Fix
1. Normalization: The solver now finds `max_coeff = max(|a|, |b|, |c|, |d|)` and scales all coefficients by `1.0 / max_coeff`.
2. Relative Threshold: It then checks `if (abs(a) < epsilon)` on the *normalized* coefficients. 

 Why this is better than previous attempts
In a previous attempt (PR #28057), a simple absolute check `abs(a) < epsilon` was proposed. That approach was rejected because it failed for scaled equations. 


Fixes #27748

modules/core/src/mathfuncs.cpp

@@ -46,6 +46,8 @@
 #include <atomic>
 #include <limits>
 #include <iostream>
+#include <algorithm>
+#include <cmath>
 #include "mathfuncs.hpp"
 
 namespace cv
@@ -1586,12 +1588,21 @@ int cv::solveCubic( InputArray _coeffs, OutputArray _roots )
         a3 = coeffs.at<double>(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<double>::epsilon())
+        return -1;
+
+    // Check if the cubic term is negligible relative to the largest coefficient
+    if( std::abs(a0) < std::numeric_limits<double>::epsilon() * max_coeff )
     {
-        if( a1 == 0 )
+        if( std::abs(a1) < std::numeric_limits<double>::epsilon() * max_coeff )
         {
-            if( a2 == 0 ) // constant
-                n = a3 == 0 ? -1 : 0;
+            if( std::abs(a2) < std::numeric_limits<double>::epsilon() * max_coeff )
+                n = std::abs(a3) < std::numeric_limits<double>::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`

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_<double>(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<double>(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