From 68bcff26fb40c9a9b12eb39700573a20fc194213 Mon Sep 17 00:00:00 2001
From: Viet Dinh <hoangviet2695@gmail.com>
Date: Thu, 5 Nov 2015 19:19:56 -0500
Subject: [PATCH] fix solveCubic

The original solution did not handle correctly when delta = 0,
resulting as nan errors. I also wrote a test case to test solving
equation x^3 = 0 after fixing.
---
 modules/core/src/mathfuncs.cpp  | 20 ++++++++++++++++++--
 modules/core/test/test_math.cpp | 11 +++++++++++
 2 files changed, 29 insertions(+), 2 deletions(-)

diff --git a/modules/core/src/mathfuncs.cpp b/modules/core/src/mathfuncs.cpp
index ed7b61e6d4..3b84e9d923 100644
--- a/modules/core/src/mathfuncs.cpp
+++ b/modules/core/src/mathfuncs.cpp
@@ -2424,7 +2424,7 @@ int cv::solveCubic( InputArray _coeffs, OutputArray _roots )
         double Qcubed = Q * Q * Q;
         double d = Qcubed - R * R;
 
-        if( d >= 0 )
+        if( d > 0 )
         {
             double theta = acos(R / sqrt(Qcubed));
             double sqrtQ = sqrt(Q);
@@ -2436,11 +2436,27 @@ int cv::solveCubic( InputArray _coeffs, OutputArray _roots )
             x2 = t0 * cos(t1 + (4.*CV_PI/3)) - t2;
             n = 3;
         }
+        else if( d == 0)
+        {
+            if(R >= 0)
+            {
+                x0 = pow(R/4, 1./3) - a1/3;
+                x1 = -pow(2*R, 1./3) - a1/3;
+            }
+            else
+            {
+                x0 = -pow(-R/4, 1./3) - a1/3;
+                x1 = pow(-2*R, 1./3) - a1/3;
+            }
+            x2 = 0;
+            n = x0 == x1 ? 1 : 2;
+            x1 = x0 == x1 ? 0 : x1;
+        }
         else
         {
             double e;
             d = sqrt(-d);
-            e = pow(d + fabs(R), 0.333333333333);
+            e = pow(d + fabs(R), 1./3);
             if( R > 0 )
                 e = -e;
             x0 = (e + Q / e) - a1 * (1./3);
diff --git a/modules/core/test/test_math.cpp b/modules/core/test/test_math.cpp
index d2c0695306..76a844406c 100644
--- a/modules/core/test/test_math.cpp
+++ b/modules/core/test/test_math.cpp
@@ -2329,6 +2329,17 @@ void Core_SolvePolyTest::run( int )
             pass = pass && div < err_eps;
         }
 
+        //test bug #5623 - solves equation x^3 = 0
+        cv::Mat coeffs_5623(4, 1, CV_64FC1);
+        cv::Mat r_5623(3, 1, CV_64FC2);
+        coeffs_5623.at<double>(0) = 1;
+        coeffs_5623.at<double>(1) = 0;
+        coeffs_5623.at<double>(2) = 0;
+        coeffs_5623.at<double>(3) = 0;
+        double prec_5623 = cv::solveCubic(coeffs_5623, r_5623);
+        pass = pass && r_5623.at<double>(0) == 0 && r_5623.at<double>(1) == 0 && r_5623.at<double>(2) == 0;
+        pass = pass && prec_5623 == 1;
+
         if (!pass)
         {
             ts->set_failed_test_info(cvtest::TS::FAIL_INVALID_OUTPUT);
-- 
GitLab