remove acb_poly_refine_roots_durand_kerner and base acb_poly_find_roo…

…ts on generics with nfloat when possible

doc/source/acb_poly.rst

@@ -1103,13 +1103,6 @@ Root-finding
     it is possible that not all of the polynomial's roots are contained
     among them.
 
-.. function:: void _acb_poly_refine_roots_durand_kerner(acb_ptr roots, acb_srcptr poly, slong len, slong prec)
-
-    Refines the given roots simultaneously using a single iteration
-    of the Durand-Kerner method. The radius of each root is set to an
-    approximation of the correction, giving a rough estimate of its error (not
-    a rigorous bound).
-
 .. function:: slong _acb_poly_find_roots(acb_ptr roots, acb_srcptr poly, acb_srcptr initial, slong len, slong maxiter, slong prec)
 
 .. function:: slong acb_poly_find_roots(acb_ptr roots, const acb_poly_t poly, acb_srcptr initial, slong maxiter, slong prec)

src/acb_poly.h

@@ -465,9 +465,6 @@ void _acb_poly_root_inclusion(acb_t r, const acb_t m,
 slong _acb_poly_validate_roots(acb_ptr roots,
         acb_srcptr poly, slong len, slong prec);
 
-void _acb_poly_refine_roots_durand_kerner(acb_ptr roots,
-        acb_srcptr poly, slong len, slong prec);
-
 slong _acb_poly_find_roots(acb_ptr roots,
     acb_srcptr poly,
     acb_srcptr initial, slong len, slong maxiter, slong prec);

src/acb_poly/find_roots.c

@@ -14,6 +14,10 @@
 #include "arf.h"
 #include "fmpq.h"
 #include "acb_poly.h"
+#include "gr.h"
+#include "gr_vec.h"
+#include "gr_poly.h"
+#include "nfloat.h"
 
 slong
 _acb_get_mid_mag(const acb_t z)
@@ -144,13 +148,81 @@ _acb_poly_roots_initial_values(acb_ptr roots, acb_srcptr poly, slong deg, slong
     arf_clear(ar);
 }
 
+#define USE_ABERTH 0
+
+/* todo: this method often gets called with degree 2 or 3 polynomials;
+   consider doing something direct there */
+int
+_acb_poly_find_roots_iter(gr_ptr w, gr_ptr z, gr_srcptr f, gr_srcptr f_prime, slong deg, slong maxiter, gr_ctx_t fp_ctx, gr_ctx_t acb_ctx, slong prec)
+{
+    slong iter, i;
+    slong rootmag, max_rootmag, correction, max_correction;
+    slong sz = fp_ctx->sizeof_elem;
+    int status = GR_SUCCESS;
+
+    acb_t t;
+    acb_init(t);
+
+    for (iter = 0; iter < maxiter; iter++)
+    {
+        /* todo: support magnitude extraction in the fp_ctx */
+        {
+            max_rootmag = -WORD_MAX;
+            for (i = 0; i < deg; i++)
+            {
+                GR_MUST_SUCCEED(gr_set_other(t, GR_ENTRY(z, i, sz), fp_ctx, acb_ctx));
+                rootmag = _acb_get_mid_mag(t);
+                max_rootmag = FLINT_MAX(rootmag, max_rootmag);
+            }
+        }
+
+#if USE_ABERTH
+            status |= _gr_poly_refine_roots_aberth(w, f, fprime, deg, z, 1, fp_ctx);
+#else
+            status |= _gr_poly_refine_roots_wdk(w, f, deg, z, 1, fp_ctx);
+#endif
+
+        /* read error estimates */
+        max_correction = -ARF_PREC_EXACT;
+        for (i = 0; i < deg; i++)
+        {
+            GR_MUST_SUCCEED(gr_set_other(t, GR_ENTRY(w, i, sz), fp_ctx, acb_ctx));
+            correction = _acb_get_mid_mag(t);
+            max_correction = FLINT_MAX(correction, max_correction);
+        }
+
+        status |= _gr_vec_sub(z, z, w, deg, fp_ctx);
+
+        /* estimate the correction relative to the whole set of roots */
+        max_correction -= max_rootmag;
+        /* flint_printf("ITER %wd MAX CORRECTION: %wd\n", iter, max_correction); */
+
+        if (max_correction < -prec / 2)
+            maxiter = FLINT_MIN(maxiter, iter + 2);
+        else if (max_correction < -prec / 3)
+            maxiter = FLINT_MIN(maxiter, iter + 3);
+        else if (max_correction < -prec / 4)
+            maxiter = FLINT_MIN(maxiter, iter + 4);
+    }
+
+    acb_clear(t);
+
+    return status;
+}
+
+
 slong
 _acb_poly_find_roots(acb_ptr roots,
     acb_srcptr poly,
     acb_srcptr initial, slong len, slong maxiter, slong prec)
 {
-    slong iter, i, deg;
-    slong rootmag, max_rootmag, correction, max_correction;
+    slong i, deg;
+    gr_ptr w, z, f, fprime;
+    gr_ctx_t acb_ctx, fp_ctx;
+    slong sz;
+    acb_t t;
+    int status = GR_SUCCESS;
+    int attempt;
 
     deg = len - 1;
 
@@ -184,41 +256,69 @@ _acb_poly_find_roots(acb_ptr roots,
     if (maxiter == 0)
         maxiter = 2 * deg + n_sqrt(prec);
 
-    for (iter = 0; iter < maxiter; iter++)
+    gr_ctx_init_complex_acb(acb_ctx, prec + 64);
+    acb_init(t);
+
+    for (attempt = 0; attempt <= 1; attempt++)
     {
-        max_rootmag = -ARF_PREC_EXACT;
-        for (i = 0; i < deg; i++)
+        /* First try with nfloat if possible, otherwise fall back to acf */
+        if (attempt == 0)
         {
-            rootmag = _acb_get_mid_mag(roots + i);
-            max_rootmag = FLINT_MAX(rootmag, max_rootmag);
+            if (nfloat_complex_ctx_init(fp_ctx, prec, 0) != GR_SUCCESS)
+                continue;
+        }
+        else
+        {
+#if 0
+            flint_printf("second try: %wd\n", prec);
+#endif
+            gr_ctx_init_complex_float_acf(fp_ctx, prec);
         }
 
-        _acb_poly_refine_roots_durand_kerner(roots, poly, len, prec);
+        status = GR_SUCCESS;
 
-        max_correction = -ARF_PREC_EXACT;
+        sz = fp_ctx->sizeof_elem;
+        z = gr_heap_init_vec(4 * deg + 1, fp_ctx);
+        w = GR_ENTRY(z, deg, sz);
+        fprime = GR_ENTRY(z, 2 * deg, sz);
+        f = GR_ENTRY(z, 3 * deg, sz);
+
+        for (i = 0; i <= deg; i++)
+            status |= gr_set_other(GR_ENTRY(f, i, sz), poly + i, acb_ctx, fp_ctx);
         for (i = 0; i < deg; i++)
+            status |= gr_set_other(GR_ENTRY(z, i, sz), roots + i, acb_ctx, fp_ctx);
+
+#if USE_ABERTH
+        status |= _gr_poly_derivative(fprime, f, deg + 1, fp_ctx);
+#endif
+
+        if (status == GR_SUCCESS)
+            status = _acb_poly_find_roots_iter(w, z, f, fprime, deg, maxiter, fp_ctx, acb_ctx, prec);
+
+        if (status == GR_SUCCESS)
         {
-            correction = _acb_get_rad_mag(roots + i);
-            max_correction = FLINT_MAX(correction, max_correction);
+            /* convert back to acb */
+            for (i = 0; i < deg; i++)
+            {
+                GR_MUST_SUCCEED(gr_set_other(roots + i, GR_ENTRY(z, i, sz), fp_ctx, acb_ctx));
+                mag_zero(arb_radref(acb_realref(roots + i)));
+                mag_zero(arb_radref(acb_imagref(roots + i)));
+            }
         }
 
-        /* estimate the correction relative to the whole set of roots */
-        max_correction -= max_rootmag;
-
-        /* flint_printf("ITER %wd MAX CORRECTION: %wd\n", iter, max_correction); */
+        gr_heap_clear_vec(z, 4 * deg + 1, fp_ctx);
+        gr_ctx_clear(fp_ctx);
 
-        if (max_correction < -prec / 2)
-            maxiter = FLINT_MIN(maxiter, iter + 2);
-        else if (max_correction < -prec / 3)
-            maxiter = FLINT_MIN(maxiter, iter + 3);
-        else if (max_correction < -prec / 4)
-            maxiter = FLINT_MIN(maxiter, iter + 4);
+        if (status == GR_SUCCESS)
+            break;
     }
 
+    acb_clear(t);
+    gr_ctx_clear(acb_ctx);
+
     return _acb_poly_validate_roots(roots, poly, len, prec);
 }
 
-
 slong
 acb_poly_find_roots(acb_ptr roots,
     const acb_poly_t poly, acb_srcptr initial,
@@ -231,6 +331,5 @@ acb_poly_find_roots(acb_ptr roots,
         flint_throw(FLINT_ERROR, "find_roots: expected a nonzero polynomial");
     }
 
-    return _acb_poly_find_roots(roots, poly->coeffs, initial,
-                len, maxiter, prec);
+    return _acb_poly_find_roots(roots, poly->coeffs, initial, len, maxiter, prec);
 }

src/python/flint_ctypes.py

@@ -5059,9 +5059,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.89e-17], [1.414213562373095 +/- 4.89e-17], -1.500000000000000], [1, 1, 2])
+            ([[-1.414213562373095 +/- 4.91e-17], [1.414213562373095 +/- 4.91e-17], -1.500000000000000], [1, 1, 2])
             >>> f.roots(domain=CC)      # complex ball roots
-            ([[-1.414213562373095 +/- 4.89e-17], [1.414213562373095 +/- 4.89e-17], ([+/- 7.23e-20] + [1.000000000000000 +/- 8e-20]*I), ([+/- 7.23e-20] + [-1.000000000000000 +/- 8e-20]*I), -1.500000000000000], [1, 1, 1, 1, 2])
+            ([[-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])
             >>> f.roots(RF)     # real floating-point roots
             ([-1.414213562373095, 1.414213562373095, -1.500000000000000], [1, 1, 2])
             >>> f.roots(CF)     # complex floating-point roots
@@ -5937,7 +5937,7 @@ def eigenvalues(self, domain=None):
             >>> Mat(ZZ)([[1,2],[3,4]]).eigenvalues(domain=QQbar)
             ([Root a = 5.37228 of a^2-5*a-2, Root a = -0.372281 of a^2-5*a-2], [1, 1])
             >>> Mat(ZZ)([[1,2],[3,4]]).eigenvalues(domain=RR)
-            ([[-0.3722813232690143 +/- 3.01e-17], [5.372281323269014 +/- 3.31e-16]], [1, 1])
+            ([[-0.3722813232690143 +/- 3.01e-17], [5.372281323269014 +/- 3.32e-16]], [1, 1])
             >>> Mat(QQbar)([[1, 0, QQbar.i()], [0, 0, 1], [1, 1, 1]]).eigenvalues()
             ([Root a = 1.94721 + 0.604643*I of a^6-4*a^5+4*a^4+2*a^3-3*a^2+1, Root a = 0.654260 - 0.430857*I of a^6-4*a^5+4*a^4+2*a^3-3*a^2+1, Root a = -0.601467 - 0.173786*I of a^6-4*a^5+4*a^4+2*a^3-3*a^2+1], [1, 1, 1])
             >>> Mat(ZZi)([[1, 0, ZZi.i()], [0, 0, 1], [1, 1, 1]]).eigenvalues(domain=QQbar)