Merge pull request #5605 from hoangviet1985:fix_bug_5599

modules/core/src/mathfuncs.cpp

@@ -2512,7 +2512,6 @@ double cv::solvePoly( InputArray _coeffs0, OutputArray _roots0, int maxIters )
     }
 
     C p(1, 0), r(1, 1);
-
     for( i = 0; i < n; i++ )
     {
         roots[i] = p;
@@ -2527,12 +2526,54 @@ 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 = max(maxDiff, abs(num));
         }

modules/core/test/test_math.cpp

@@ -2359,6 +2359,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_CheckRange_Empty : public cvtest::BaseTest
 {
 public:

samples/cpp/polynominal_equations.cpp

@@ -0,0 +1,38 @@
+ /*
+ * solves equations x^4 - x^2 = 0, and x^2 - 2x + 1 = 0
+ */
+
+#include "opencv2/core/core.hpp"
+#include <iostream>
+
+
+int main(void)
+{
+        cv::Mat coefs = (cv::Mat_<float>(1,5) << 0, 0, -1, 0, 1 );
+        std::cout << "x^4 - x^2 = 0\n\n" << "Coefficients: " << coefs << "\n" << std::endl;
+        cv::Mat r;
+
+        double prec;
+        prec = cv::solvePoly(coefs, r);
+        std::cout << "Preciseness = " << prec << std::endl;
+        std::cout << "roots after 1000 iterations:\n" << r << "\n" << std::endl;
+
+        prec = cv::solvePoly(coefs, r, 9999);
+        std::cout << "Preciseness = " << prec << std::endl;
+        std::cout << "roots after 9999 iterations:\n" << r << "\n" << std::endl;
+
+        std::cout << "\n---------------------------------------\n" << std::endl;
+
+        coefs = (cv::Mat_<float>(1,3) << 1, -2, 1 );
+        std::cout << "x^2 - 2x + 1 = 0\n\n" << "Coefficients: " << coefs << "\n" << std::endl;
+
+        prec = cv::solvePoly(coefs, r);
+        std::cout << "Preciseness = " << prec << std::endl;
+        std::cout << "roots after 1000 iterations:\n" << r << "\n" << std::endl;
+
+        prec = cv::solvePoly(coefs, r, 9999);
+        std::cout << "Preciseness = " << prec << std::endl;
+        std::cout << "roots after 9999 iterations:\n" << r << "\n" << std::endl;
+
+    return 0;
+}