Skip to content

Commit

Permalink
Merge pull request #2013 from fredrik-johansson/initial
Browse files Browse the repository at this point in the history 
Use Bini's algorithm to select initial values in acb_poly_find_roots
  • Loading branch information
@fredrik-johansson
fredrik-johansson authored Jun 6, 2024
2 parents 187b9bf + 80bca8a commit a559e2f
Show file tree
Hide file tree
Showing 4 changed files with 104 additions and 21 deletions.
8 changes: 7 additions & 1 deletion doc/source/acb_poly.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand Down
4 changes: 3 additions & 1 deletion doc/source/references.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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]_
111 changes: 93 additions & 18 deletions src/acb_poly/find_roots.c
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
/*
Copyright (C) 2012 Fredrik Johansson
Copyright (C) 2012, 2024 Fredrik Johansson
This file is part of FLINT.
Expand All @@ -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"

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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);

Expand Down
2 changes: 1 addition & 1 deletion src/python/flint_ctypes.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit a559e2f

Please sign in to comment.