Merge pull request #2422 from fredrik-johansson/polyinit

Fix initial value selection in acb_poly_find_roots for overlapping radii

Author: fredrik-johansson

src/acb_poly/find_roots.c

@@ -123,10 +123,17 @@ _acb_poly_roots_initial_values(acb_ptr roots, acb_srcptr poly, slong deg, slong
         flint_printf("initial values: radius %{mag} multiplicity %wd\n", r, m);
 #endif
 
-        /* pick m points on the circle with radius r */
+        /* Pick m points on the circle with radius r */
         for (j = 0; j < m; j++)
         {
-            theta = 6.283185307179586 * ((j + 1.0) / m + (double) i / deg) + 0.577216;
+            /* Hack: the angle 2pi * ((j+1)/m + i/deg) distributes points uniformly
+               when there are many radii. However, if two radii are
+               very close, points can end up overlapping if (j+1)/m + i/deg
+               happens to give the same fraction for two different i. Ensure
+               that this doesn't happen. */
+            double irrational_factor = 0.99071828182845905;
+
+            theta = 6.283185307179586 * (irrational_factor * (j + 1.0) / m + (double) i / deg) + 0.577216;
 
             ri = roots + total;
             acb_zero(ri);

src/acb_poly/test/t-find_roots.c

@@ -85,5 +85,34 @@ TEST_FUNCTION_START(acb_poly_find_roots, state)
         acb_poly_clear(C);
     }
 
+    /* Check that #2421 is fixed */
+    {
+        acb_poly_t f;
+        acb_ptr roots;
+        acb_t c;
+
+        acb_poly_init(f);
+        acb_init(c);
+        roots = _acb_vec_init(4);
+
+        acb_poly_set_coeff_si(f, 0, 10);
+        acb_poly_set_coeff_si(f, 1, -170);
+        acb_poly_set_coeff_si(f, 2, 2890);
+        acb_poly_set_coeff_si(f, 3, -49130);
+        acb_poly_set_coeff_si(f, 4, 835210);
+        acb_set_ui(c, 83521);
+        acb_poly_scalar_div(f, f, c, 64);
+
+        if (acb_poly_find_roots(roots, f, NULL, 40, 64) != 4)
+        {
+            flint_printf("FAIL: geometric example\n");
+            flint_abort();
+        }
+
+        _acb_vec_clear(roots, 4);
+        acb_poly_clear(f);
+        acb_clear(c);
+    }
+
     TEST_FUNCTION_END(state);
 }

src/python/flint_ctypes.py

@@ -5422,9 +5422,9 @@ def roots(self, domain=None):
             >>> f.roots(domain=QQbar)     # complex algebraic roots
             ([Root a = 1.00000*I of a^2+1, Root a = -1.00000*I of a^2+1, Root a = 1.41421 of a^2-2, Root a = -1.41421 of a^2-2, -3/2], [1, 1, 1, 1, 2])
             >>> f.roots(domain=RR)      # real ball roots
-            ([[-1.414213562373095 +/- 4.91e-17], [1.414213562373095 +/- 4.91e-17], -1.500000000000000], [1, 1, 2])
+            ([[-1.414213562373095 +/- 4.89e-17], [1.414213562373095 +/- 4.89e-17], -1.500000000000000], [1, 1, 2])
             >>> f.roots(domain=CC)      # complex ball roots
-            ([[-1.414213562373095 +/- 4.91e-17], [1.414213562373095 +/- 4.91e-17], ([+/- 1.45e-19] + [1.000000000000000 +/- 2.7e-19]*I), ([+/- 1.45e-19] + [-1.000000000000000 +/- 2.7e-19]*I), -1.500000000000000], [1, 1, 1, 1, 2])
+            ([[-1.414213562373095 +/- 4.89e-17], [1.414213562373095 +/- 4.89e-17], ([+/- 1.45e-19] + [1.000000000000000 +/- 2.7e-19]*I), ([+/- 1.45e-19] + [-1.000000000000000 +/- 2.7e-19]*I), -1.500000000000000], [1, 1, 1, 1, 2])
             >>> f.roots(RF)     # real floating-point roots
             ([-1.414213562373095, 1.414213562373095, -1.500000000000000], [1, 1, 2])
             >>> f.roots(CF)     # complex floating-point roots