diff --git a/modules/core/misc/java/test/CoreTest.java b/modules/core/misc/java/test/CoreTest.java
index 62d72b706f..8d6b731616 100644
--- a/modules/core/misc/java/test/CoreTest.java
+++ b/modules/core/misc/java/test/CoreTest.java
@@ -1820,12 +1820,18 @@ public class CoreTest extends OpenCVTestCase {
 
         assertGE(1e-6, Math.abs(Core.solvePoly(coeffs, roots)));
 
-        truth = new Mat(3, 1, CvType.CV_32FC2) {
+        List<Mat> rootsReIm = new ArrayList<Mat>();
+        Core.split(roots, rootsReIm);
+        Core.sort(rootsReIm.get(0), dst, Core.SORT_EVERY_COLUMN);
+
+        Mat truthRe = new Mat(3, 1, CvType.CV_32F) {
             {
-                put(0, 0, 1, 0, 2, 0, 3, 0);
+                put(0, 0, 1, 2, 3);
             }
         };
-        assertMatEqual(truth, roots, EPS);
+        Mat truthIm = Mat.zeros(3, 1, CvType.CV_32F);
+        assertMatEqual(truthRe, dst, EPS);
+        assertMatEqual(truthIm, rootsReIm.get(1), EPS);
     }
 
     public void testSolvePolyMatMatInt() {
@@ -1836,14 +1842,20 @@ public class CoreTest extends OpenCVTestCase {
         };
         Mat roots = new Mat();
 
-        assertEquals(10.198039027185569, Core.solvePoly(coeffs, roots, 1));
+        assertGE(1e-6, Core.solvePoly(coeffs, roots, 10));
 
-        truth = new Mat(3, 1, CvType.CV_32FC2) {
+        List<Mat> rootsReIm = new ArrayList<Mat>();
+        Core.split(roots, rootsReIm);
+        Core.sort(rootsReIm.get(0), dst, Core.SORT_EVERY_COLUMN);
+
+        Mat truthRe = new Mat(3, 1, CvType.CV_32F) {
             {
-                put(0, 0, 1, 0, -1, 2, -2, 12);
+                put(0, 0, 1, 2, 3);
             }
         };
-        assertMatEqual(truth, roots, EPS);
+        Mat truthIm = Mat.zeros(3, 1, CvType.CV_32F);
+        assertMatEqual(truthRe, dst, EPS);
+        assertMatEqual(truthIm, rootsReIm.get(1), EPS);
     }
 
     public void testSort() {
diff --git a/modules/core/src/mathfuncs.cpp b/modules/core/src/mathfuncs.cpp
index abcfb54215..c3bd50e0d4 100644
--- a/modules/core/src/mathfuncs.cpp
+++ b/modules/core/src/mathfuncs.cpp
@@ -1748,7 +1748,40 @@ double cv::solvePoly( InputArray _coeffs0, OutputArray _roots0, int maxIters )
             break;
     }
 
-    C p(1, 0), r(1, 1);
+    // Related issue: https://github.com/opencv/opencv/issues/23644,
+    // This the initialization scheme of "Initial approximations in Durand-Kerner's root finding method" by Guggenheimer.
+    // https://link.springer.com/article/10.1007/BF01935059
+    // We put the initial points equidistantly on a circle on the complex plane. This code computes the circle radius as in the paper.
+    Mat absCoeffs(n + 1, 1, CV_64F);
+    for( i = 0; i <= n; i++ )
+        absCoeffs.at<double>(i) = abs(coeffs[i]);
+
+    int nonZeroCoeffs = 0;
+    Mat u(n, 1, CV_64F, Scalar(0)), v(n, 1, CV_64F, Scalar(0));
+    for( i = 0; i <= n; i++ )
+    {
+        double coeff = absCoeffs.at<double>(i);
+        if( coeff > DBL_EPSILON )
+        {
+            if( i != n )
+                u.at<double>(i) = 2.0 * pow(coeff / absCoeffs.at<double>(n), 1.0 / (n - i));
+            if( i != 0 )
+                v.at<double>(i - 1) = 0.5 * pow(absCoeffs.at<double>(0) / coeff, 1.0 / i);
+            nonZeroCoeffs++;
+        }
+    }
+    double scale = 1;
+    if( nonZeroCoeffs > 2 )
+    {
+        Point maxU, minV;
+        minMaxLoc(u, nullptr, nullptr, nullptr, &maxU);
+        minMaxLoc(v, nullptr, nullptr, &minV);
+        u.at<double>(maxU) = 0;
+        v.at<double>(minV) = 0;
+        scale = (sum(u).val[0] + sum(v).val[0]) / (2 * nonZeroCoeffs - 2);
+    }
+
+    C p(scale, 0), r(cos(CV_2PI / n), sin(CV_2PI / n));
 
     for( i = 0; i < n; i++ )
     {
diff --git a/modules/core/test/test_math.cpp b/modules/core/test/test_math.cpp
index 1c9aa81d32..70c244d351 100644
--- a/modules/core/test/test_math.cpp
+++ b/modules/core/test/test_math.cpp
@@ -2679,6 +2679,93 @@ TEST(Core_SolvePoly, regression_5599)
     }
 }
 
+TEST(Core_SolvePoly, regression_23644)
+{
+    // x^2 - 2x - 3 = 0,  roots: 3, -1
+    cv::Mat coefs = (cv::Mat_<float>(1,3) << -3, -2, 1 );
+    cv::Mat r;
+    double prec;
+    prec = cv::solvePoly(coefs, r);
+    EXPECT_LE(prec, 1e-6);
+    EXPECT_EQ(2u, r.total());
+    ASSERT_EQ(CV_32FC2, r.type());
+    checkRoot<float>(r, 3, 0);
+    checkRoot<float>(r, -1, 0);
+}
+
+TEST(Core_SolvePoly, degree_2_polynomials)
+{
+    cv::Mat_<float> coefs(1,3);
+    cv::Mat r;
+    double prec;
+    for (int c0 = -20; c0 <= 20; c0++)
+    {
+        coefs.at<float>(0) = c0;
+        for (int c1 = -20; c1 <= 20; c1++)
+        {
+            coefs.at<float>(1) = c1;
+            for (int c2 = -20; c2 <= 20; c2++)
+            {
+                coefs.at<float>(2) = c2;
+                prec = cv::solvePoly(coefs, r);
+                EXPECT_LE(prec, 1e-6);
+            }
+        }
+    }
+}
+
+TEST(Core_SolvePoly, degree_4_polynomials)
+{
+    applyTestTag(CV_TEST_TAG_VERYLONG);
+
+    cv::Mat_<float> coefs(1,5);
+    cv::Mat r;
+    double prec;
+    for (int c0 = -10; c0 <= 10; c0++)
+    {
+        coefs.at<float>(0) = c0;
+        for (int c1 = -10; c1 <= 10; c1++)
+        {
+            coefs.at<float>(1) = c1;
+            for (int c2 = -10; c2 <= 10; c2++)
+            {
+                coefs.at<float>(2) = c2;
+                for (int c3 = -10; c3 <= 10; c3++)
+                {
+                    coefs.at<float>(3) = c3;
+                    for (int c4 = -10; c4 <= 10; c4++)
+                    {
+                        coefs.at<float>(4) = c4;
+                        prec = cv::solvePoly(coefs, r);
+                        EXPECT_LE(prec, 1e-3);
+                    }
+                }
+            }
+        }
+    }
+}
+
+TEST(Core_SolvePoly, different_magnitudes_polynomials)
+{
+    cv::Mat_<float> coefs(1,3);
+    cv::Mat r;
+    double prec;
+    for (int i = -10; i < 10; i++)
+    {
+        coefs.at<float>(0) = pow(2, i);
+        for (int j = -10; j < 10; j++)
+        {
+            coefs.at<float>(1) = pow(2, j);
+            for (int k = -10; k < 10; k++)
+            {
+                coefs.at<float>(2) = pow(2, k);
+                prec = cv::solvePoly(coefs, r);
+                EXPECT_LE(prec, 1e-6);
+            }
+        }
+    }
+}
+
 class Core_PhaseTest : public cvtest::BaseTest
 {
     int t;