diff --git a/doc/source/acb_poly.rst b/doc/source/acb_poly.rst
index abe11d4378..759f7f16d3 100644
--- a/doc/source/acb_poly.rst
+++ b/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),
diff --git a/doc/source/references.rst b/doc/source/references.rst
index e1e5185524..83ed474d26 100644
--- a/doc/source/references.rst
+++ b/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]_
diff --git a/src/acb_poly/find_roots.c b/src/acb_poly/find_roots.c
index 0f4b2bb05f..586f3293e3 100644
--- a/src/acb_poly/find_roots.c
+++ b/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 .
*/
+#include
#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);
diff --git a/src/python/flint_ctypes.py b/src/python/flint_ctypes.py
index 4eeb4e4608..a7f513d2eb 100644
--- a/src/python/flint_ctypes.py
+++ b/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