Use Bini's algorithm to select initial values in acb_poly_find_roots

doc/source/acb_poly.rst

@@ -1132,7 +1132,13 @@ Root-finding
     a default value. Finally, the approximate roots are validated rigorously.
 
     Initial values for the iteration can be provided as the array *initial*.
-    If *initial* is set to *NULL*, default values `(0.4+0.9i)^k` are used.
+    If *initial* is set to *NULL*, default values are chosen using the
+    Newton polygon technique of [Bin1996]_.
+    Suppose that the polynomial is `\sum_k a_k x^k` and
+    let `k_1, \ldots, k_q` denote the index coordinates of the
+    vertices in the convex hull of `(k, \log |a_k|)`.
+    For each `2 \le i \le q`, we pick `m_i = k_i - k_{i-1}` points
+    on the circle with radius `r_i = (|a_{k_{i-1}} / a_{k_i} |)^{1/m_i}`.
 
     The polynomial is assumed to be squarefree. If there are repeated
     roots, the iteration is likely to find them (with low numerical accuracy),

doc/source/references.rst

@@ -41,6 +41,8 @@ References
 
 .. [BerTas2010] \D. Berend and T. Tassa : Improved bounds on Bell numbers and on moments of sums of random variables, Probability and Mathematical Statistics vol. 30 (2010) 185--205
 
+.. [Bin1996] \D. A. Bini, "Numerical computation of polynomial zeros by means of Aberth's method". Numerical algorithms 13 (1996): 179-200. https://doi.org/10.1007/BF02207694
+
 .. [Blo2009] \R. Bloemen, "Even faster zeta(2n) calculation!", https://web.archive.org/web/20141101133659/http://xn--2-umb.com/09/11/even-faster-zeta-calculation
 
 .. [Bodrato2010] \Bodrato, Marco : A Strassen-like Matrix Multiplication Suited for Squaring and Higher Power Computation. Proceedings of the ISSAC 2010 München, Germany, 25-28 July, 2010
@@ -337,4 +339,4 @@ References
 
 .. [vdH2006] \J. van der Hoeven, "Computations with effective real numbers". Theoretical Computer Science, Volume 351, Issue 1, 14 February 2006, Pages 52-60. https://doi.org/10.1016/j.tcs.2005.09.060
 
-All referenced works: [AbbottBronsteinMulders1999]_, [Apostol1997]_, [Ari2011]_, [Ari2012]_, [Arn2010]_, [Arn2012]_, [ArnoldMonagan2011]_, [BBC1997]_, [BBC2000]_, [BBK2014]_, [BD1992]_, [BF2020]_, [BFSS2006]_, [BJ2013]_, [BM1980]_, [BZ1992]_, [BZ2011]_, [BaiWag1980]_, [BerTas2010]_, [Blo2009]_, [Bodrato2010]_, [Boe2020]_, [Bog2012]_, [Bol1887]_, [Bor1987]_, [Bor2000]_, [Bre1978]_, [Bre1979]_, [Bre2010]_, [BrentKung1978]_, [BuhlerCrandallSompolski1992]_, [CFG2017]_, [CFG2019]_, [CGHJK1996]_, [CP2005]_, [Car1995]_, [Car2004]_, [Chen2003]_, [Cho1999]_, [Coh1996]_, [Coh2000]_, [Col1971]_, [CraPom2005]_, [DHBHS2004]_, [DYF1999]_, [DelegliseNicolasZimmermann2009]_, [DomKanTro1987]_, [Dup2006]_, [Dus1999]_, [EHJ2016]_, [EM2004]_, [EK2023]_, [Fie2007]_, [FieHof2014]_, [Fil1992]_, [GCL1992]_, [GG2003]_, [GS2003]_, [GVL1996]_, [Gas2018]_, [Gos1974]_, [GowWag2008]_, [GraMol2010]_, [HM2017]_, [HS1967]_, [HZ2004]_, [HanZim2004]_, [Har2010]_, [HZ2011]_, [Har2012]_, [Har2015]_, [Har2018]_, [Hart2010]_, [Hen1956]_, [Hoe2001]_, [Hoe2009]_, [Hor1972]_, [Iliopoulos1989]_, [Igu1972]_, [Igu1979]_, [JB2018]_, [JM2018]_, [JR1999]_, [Joh2012]_, [Joh2013]_, [Joh2014a]_, [Joh2014b]_, [Joh2014c]_, [Joh2015]_, [Joh2015b]_, [Joh2016]_, [Joh2017]_, [Joh2017a]_, [Joh2017b]_, [Joh2018a]_, [Joh2018b]_, [JvdP2002]_, [Kahan1991]_, [KanBac1979]_, [Kar1998]_, [Knu1997]_, [Kob2010]_, [Kri2013]_, [LT2016]_, [Leh1970]_, [LukPatWil1996]_, [MN2019]_, [MP2006]_, [MPFR2012]_, [MasRob1996]_, [Mic2007]_, [Miy2010]_, [Mos1971]_, [Mul2000]_, [Mum1983]_, [Mum1984]_, [NIST2012]_, [NakTurWil1997]_, [Olv1997]_, [PP2010]_, [PS1973]_, [PS1991]_, [Paterson1973]_, [PernetStein2010]_, [Pet1999]_, [Pla2011]_, [Pla2017]_, [RF1994]_, [Rad1973]_, [Rademacher1937]_, [Ric1992]_, [Ric1995]_, [Ric1997]_, [Ric2007]_, [Ric2009]_, [RosSch1962]_, [Rum2010]_, [Smi2001]_, [SorWeb2016]_, [Ste2002]_, [Ste2010]_, [Stehle2010]_, [Stein2007]_, [Sut2007]_, [StoMul1998]_, [Str2014]_, [Str1997]_, [Str2012]_, [Tak2000]_, [ThullYap1990]_, [Tre2008]_, [Tru2011]_, [Tru2014]_, [Tur1953]_, [Villard2007]_, [WaktinsZeitlin1993]_, [Wei2000]_, [Whiteman1956]_, [Zip1985]_, [Zun2023]_, [Zun2023b]_, [vHP2012]_, [vdH1995]_, [vdH2006]_
+All referenced works: [AbbottBronsteinMulders1999]_, [Apostol1997]_, [Ari2011]_, [Ari2012]_, [Arn2010]_, [Arn2012]_, [ArnoldMonagan2011]_, [BBC1997]_, [BBC2000]_, [BBK2014]_, [BD1992]_, [BF2020]_, [BFSS2006]_, [BJ2013]_, [BM1980]_, [BZ1992]_, [BZ2011]_, [BaiWag1980]_, [BerTas2010]_, [Bin1996]_, [Blo2009]_, [Bodrato2010]_, [Boe2020]_, [Bog2012]_, [Bol1887]_, [Bor1987]_, [Bor2000]_, [Bre1978]_, [Bre1979]_, [Bre2010]_, [BrentKung1978]_, [BuhlerCrandallSompolski1992]_, [CFG2017]_, [CFG2019]_, [CGHJK1996]_, [CP2005]_, [Car1995]_, [Car2004]_, [Chen2003]_, [Cho1999]_, [Coh1996]_, [Coh2000]_, [Col1971]_, [CraPom2005]_, [DHBHS2004]_, [DYF1999]_, [DelegliseNicolasZimmermann2009]_, [DomKanTro1987]_, [Dup2006]_, [Dus1999]_, [EHJ2016]_, [EM2004]_, [EK2023]_, [Fie2007]_, [FieHof2014]_, [Fil1992]_, [GCL1992]_, [GG2003]_, [GS2003]_, [GVL1996]_, [Gas2018]_, [Gos1974]_, [GowWag2008]_, [GraMol2010]_, [HM2017]_, [HS1967]_, [HZ2004]_, [HanZim2004]_, [Har2010]_, [HZ2011]_, [Har2012]_, [Har2015]_, [Har2018]_, [Hart2010]_, [Hen1956]_, [Hoe2001]_, [Hoe2009]_, [Hor1972]_, [Iliopoulos1989]_, [Igu1972]_, [Igu1979]_, [JB2018]_, [JM2018]_, [JR1999]_, [Joh2012]_, [Joh2013]_, [Joh2014a]_, [Joh2014b]_, [Joh2014c]_, [Joh2015]_, [Joh2015b]_, [Joh2016]_, [Joh2017]_, [Joh2017a]_, [Joh2017b]_, [Joh2018a]_, [Joh2018b]_, [JvdP2002]_, [Kahan1991]_, [KanBac1979]_, [Kar1998]_, [Knu1997]_, [Kob2010]_, [Kri2013]_, [LT2016]_, [Leh1970]_, [LukPatWil1996]_, [MN2019]_, [MP2006]_, [MPFR2012]_, [MasRob1996]_, [Mic2007]_, [Miy2010]_, [Mos1971]_, [Mul2000]_, [Mum1983]_, [Mum1984]_, [NIST2012]_, [NakTurWil1997]_, [Olv1997]_, [PP2010]_, [PS1973]_, [PS1991]_, [Paterson1973]_, [PernetStein2010]_, [Pet1999]_, [Pla2011]_, [Pla2017]_, [RF1994]_, [Rad1973]_, [Rademacher1937]_, [Ric1992]_, [Ric1995]_, [Ric1997]_, [Ric2007]_, [Ric2009]_, [RosSch1962]_, [Rum2010]_, [Smi2001]_, [SorWeb2016]_, [Ste2002]_, [Ste2010]_, [Stehle2010]_, [Stein2007]_, [Sut2007]_, [StoMul1998]_, [Str2014]_, [Str1997]_, [Str2012]_, [Tak2000]_, [ThullYap1990]_, [Tre2008]_, [Tru2011]_, [Tru2014]_, [Tur1953]_, [Villard2007]_, [WaktinsZeitlin1993]_, [Wei2000]_, [Whiteman1956]_, [Zip1985]_, [Zun2023]_, [Zun2023b]_, [vHP2012]_, [vdH1995]_, [vdH2006]_

src/acb_poly/find_roots.c

@@ -1,5 +1,5 @@
 /*
-    Copyright (C) 2012 Fredrik Johansson
+    Copyright (C) 2012, 2024 Fredrik Johansson
 
     This file is part of FLINT.
 
@@ -9,7 +9,9 @@
     (at your option) any later version.  See <https://www.gnu.org/licenses/>.
 */
 
+#include <math.h>
 #include "ulong_extras.h"
+#include "arf.h"
 #include "fmpq.h"
 #include "acb_poly.h"
 
@@ -44,29 +46,102 @@ _acb_get_rad_mag(const acb_t z)
     return FLINT_MAX(rm, im);
 }
 
-void
-_acb_poly_roots_initial_values(acb_ptr roots, slong deg, slong prec)
+/* Compute res[0, ..., n-1] = {i} indexing the convex hull of {i, y[i]}. */
+static slong convex_hull(slong * res, const double * y, slong len)
 {
-    slong i;
+    slong i, n = 0;
 
-    fmpq_t q;
-    fmpq_init(q);
+    for (i = 0; i < len; i++)
+    {
+        while (n >= 2 && (res[n - 2] - res[n - 1]) * (y[i] - y[res[n - 1]])
+                <= (i - res[n - 1]) * (y[res[n - 2]] - y[res[n - 1]]))
+            n--;
 
-    fmpq_set_si(q, 4, 10);
-    arb_set_fmpq(acb_realref(roots + 0), q, prec);
-    mag_zero(arb_radref(acb_realref(roots + 0)));
+        res[n] = i;
+        n++;
+    }
 
-    fmpq_set_si(q, 9, 10);
-    arb_set_fmpq(acb_imagref(roots + 0), q, prec);
-    mag_zero(arb_radref(acb_imagref(roots + 0)));
-    fmpq_clear(q);
+    return n;
+}
 
-    for (i = 1; i < deg; i++)
+void
+_acb_poly_roots_initial_values(acb_ptr roots, acb_srcptr poly, slong deg, slong prec)
+{
+    double * alog;
+    mag_ptr amag;
+    mag_t r;
+    arf_t ar;
+    slong i, j;
+    slong * ki, num, m, total;
+    double theta;
+    acb_t cmid;
+    acb_ptr ri;
+
+    amag = _mag_vec_init(deg + 1);
+    alog = flint_malloc(sizeof(double) * (deg + 1));
+    ki = flint_malloc((deg + 1) * sizeof(slong));
+    mag_init(r);
+    arf_init(ar);
+    acb_init(cmid);  /* shallow only; will not be cleared */
+
+    for (i = 0; i <= deg; i++)
     {
-        acb_mul(roots + i, roots + i - 1, roots + 0, prec);
-        mag_zero(arb_radref(acb_realref(roots + i)));
-        mag_zero(arb_radref(acb_imagref(roots + i)));
+        /* shallow midpoint */
+        *arb_midref(acb_realref(cmid)) = *arb_midref(acb_realref(poly + i));
+        *arb_midref(acb_imagref(cmid)) = *arb_midref(acb_imagref(poly + i));
+        acb_get_mag(amag + i, cmid);
+
+        /* todo: mag_get_d_log2_approx is not very precise; probably
+           this does not matter though? */
+        alog[i] = mag_get_d_log2_approx(amag + i);
     }
+
+    num = convex_hull(ki, alog, deg + 1);
+    total = 0;
+
+    for (i = 1; i < num; i++)
+    {
+        /* multiplicity */
+        m = ki[i] - ki[i - 1];
+
+        /* radius */
+        mag_div(r, amag + ki[i - 1], amag + ki[i]);
+        mag_root(r, r, m);
+
+        if (mag_is_zero(r))
+            mag_set_ui_2exp_si(r, 1, -prec);
+        if (mag_is_inf(r))
+            mag_set_ui_2exp_si(r, 1, prec);
+
+        arf_set_mag(ar, r);
+
+#if 0
+        flint_printf("initial values: radius %{mag} multiplicity %wd\n", r, m);
+#endif
+
+        /* 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;
+
+            ri = roots + total;
+            acb_zero(ri);
+            arf_set_d(arb_midref(acb_realref(ri)), cos(theta));
+            arf_set_d(arb_midref(acb_imagref(ri)), sin(theta));
+            arf_mul(arb_midref(acb_realref(ri)), arb_midref(acb_realref(ri)), ar, MAG_BITS, ARF_RND_DOWN);
+            arf_mul(arb_midref(acb_imagref(ri)), arb_midref(acb_imagref(ri)), ar, MAG_BITS, ARF_RND_DOWN);
+            total++;
+        }
+    }
+
+    if (total != deg)
+        flint_abort();
+
+    _mag_vec_clear(amag, deg + 1);
+    flint_free(alog);
+    flint_free(ki);
+    mag_clear(r);
+    arf_clear(ar);
 }
 
 slong
@@ -102,7 +177,7 @@ _acb_poly_find_roots(acb_ptr roots,
     }
 
     if (initial == NULL)
-        _acb_poly_roots_initial_values(roots, deg, prec);
+        _acb_poly_roots_initial_values(roots, poly, deg, prec);
     else
         _acb_vec_set(roots, initial, deg);
 

src/python/flint_ctypes.py

@@ -5061,7 +5061,7 @@ def roots(self, domain=None):
             >>> f.roots(domain=RR)      # real ball roots
             ([[-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.89e-17], [1.414213562373095 +/- 4.89e-17], 1.000000000000000*I, -1.000000000000000*I, -1.500000000000000], [1, 1, 1, 1, 2])
+            ([[-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])
             >>> f.roots(RF)     # real floating-point roots
             ([-1.414213562373095, 1.414213562373095, -1.500000000000000], [1, 1, 2])
             >>> f.roots(CF)     # complex floating-point roots