Merge pull request #5651 from hoangviet1985:fix_solvePoly_3.0.0

modules/core/src/mathfuncs.cpp

@@ -2108,12 +2108,53 @@ double cv::solvePoly( InputArray _coeffs0, OutputArray _roots0, int maxIters )
         {
             p = roots[i];
             C num = coeffs[n], denom = coeffs[n];
+            int num_same_root = 1;
             for( j = 0; j < n; j++ )
             {
                 num = num*p + coeffs[n-j-1];
-                if( j != i ) denom = denom * (p - roots[j]);
+                if( j != i )
+                {
+                    if ( (p - roots[j]).re != 0 || (p - roots[j]).im != 0 )
+                        denom = denom * (p - roots[j]);
+                    else
+                        num_same_root++;
+                }
             }
             num /= denom;
+            if( num_same_root > 1)
+            {
+                double old_num_re = num.re;
+                double old_num_im = num.im;
+                int square_root_times = num_same_root % 2 == 0 ? num_same_root / 2 : num_same_root / 2 - 1;
+
+                for( j = 0; j < square_root_times; j++)
+                {
+                    num.re = old_num_re*old_num_re + old_num_im*old_num_im;
+                    num.re = sqrt(num.re);
+                    num.re += old_num_re;
+                    num.im = num.re - old_num_re;
+                    num.re /= 2;
+                    num.re = sqrt(num.re);
+
+                    num.im /= 2;
+                    num.im = sqrt(num.im);
+                    if( old_num_re < 0 ) num.im = -num.im;
+                }
+
+                if( num_same_root % 2 != 0){
+                    Mat cube_coefs(4, 1, CV_64FC1);
+                    Mat cube_roots(3, 1, CV_64FC2);
+                    cube_coefs.at<double>(3) = -(pow(old_num_re, 3));
+                    cube_coefs.at<double>(2) = -(15*pow(old_num_re, 2) + 27*pow(old_num_im, 2));
+                    cube_coefs.at<double>(1) = -48*old_num_re;
+                    cube_coefs.at<double>(0) = 64;
+                    solveCubic(cube_coefs, cube_roots);
+
+                    if(cube_roots.at<double>(0) >= 0) num.re = pow(cube_roots.at<double>(0), 1./3);
+                    else num.re = -pow(-cube_roots.at<double>(0), 1./3);
+                    num.im = sqrt(pow(num.re, 2) / 3 - old_num_re / (3*num.re));
+                }
+            }
             roots[i] = p - num;
             maxDiff = std::max(maxDiff, cv::abs(num));
         }

modules/core/test/test_math.cpp

@@ -2360,6 +2360,55 @@ void Core_SolvePolyTest::run( int )
     }
 }
 
+template<typename T>
+static void checkRoot(Mat& r, T re, T im)
+{
+    for (int i = 0; i < r.cols*r.rows; i++)
+    {
+        Vec<T, 2> v = *(Vec<T, 2>*)r.ptr(i);
+        if (fabs(re - v[0]) < 1e-6 && fabs(im - v[1]) < 1e-6)
+        {
+            v[0] = std::numeric_limits<T>::quiet_NaN();
+            v[1] = std::numeric_limits<T>::quiet_NaN();
+            return;
+        }
+    }
+    GTEST_NONFATAL_FAILURE_("Can't find root") << "(" << re << ", " << im << ")";
+}
+TEST(Core_SolvePoly, regression_5599)
+{
+    // x^4 - x^2 = 0, roots: 1, -1, 0, 0
+    cv::Mat coefs = (cv::Mat_<float>(1,5) << 0, 0, -1, 0, 1 );
+    {
+        cv::Mat r;
+        double prec;
+        prec = cv::solvePoly(coefs, r);
+        EXPECT_LE(prec, 1e-6);
+        EXPECT_EQ(4, r.total());
+        //std::cout << "Preciseness = " << prec << std::endl;
+        //std::cout << "roots:\n" << r << "\n" << std::endl;
+        ASSERT_EQ(CV_32FC2, r.type());
+        checkRoot<float>(r, 1, 0);
+        checkRoot<float>(r, -1, 0);
+        checkRoot<float>(r, 0, 0);
+        checkRoot<float>(r, 0, 0);
+    }
+    // x^2 - 2x + 1 = 0,  roots: 1, 1
+    coefs = (cv::Mat_<float>(1,3) << 1, -2, 1 );
+    {
+        cv::Mat r;
+        double prec;
+        prec = cv::solvePoly(coefs, r);
+        EXPECT_LE(prec, 1e-6);
+        EXPECT_EQ(2, r.total());
+        //std::cout << "Preciseness = " << prec << std::endl;
+        //std::cout << "roots:\n" << r << "\n" << std::endl;
+        ASSERT_EQ(CV_32FC2, r.type());
+        checkRoot<float>(r, 1, 0);
+        checkRoot<float>(r, 1, 0);
+    }
+}
+
 class Core_PhaseTest : public cvtest::BaseTest
 {
 public: