From 80bca8a04a3f9fb6618bd2d069ee30ec4b88a5db Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Wed, 5 Jun 2024 21:28:01 +0200 Subject: [PATCH] Use Bini's algorithm to select initial values in acb_poly_find_roots --- doc/source/acb_poly.rst | 8 ++- doc/source/references.rst | 4 +- src/acb_poly/find_roots.c | 111 +++++++++++++++++++++++++++++++------ src/python/flint_ctypes.py | 2 +- 4 files changed, 104 insertions(+), 21 deletions(-) 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