From 9f063941c093031184dc2716cbc56955ffcf39b8 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Fri, 7 Jun 2024 15:33:48 +0200 Subject: [PATCH] remove acb_poly_refine_roots_durand_kerner and base acb_poly_find_roots on generics with nfloat when possible --- doc/source/acb_poly.rst | 7 - src/acb_poly.h | 3 - src/acb_poly/find_roots.c | 147 +++++++++++++++++---- src/acb_poly/refine_roots_durand_kerner.c | 153 ---------------------- src/python/flint_ctypes.py | 6 +- 5 files changed, 126 insertions(+), 190 deletions(-) delete mode 100644 src/acb_poly/refine_roots_durand_kerner.c diff --git a/doc/source/acb_poly.rst b/doc/source/acb_poly.rst index 759f7f16d3..b102bcbc3b 100644 --- a/doc/source/acb_poly.rst +++ b/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) diff --git a/src/acb_poly.h b/src/acb_poly.h index 5cafd3f05b..a7bb1d307b 100644 --- a/src/acb_poly.h +++ b/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); diff --git a/src/acb_poly/find_roots.c b/src/acb_poly/find_roots.c index 586f3293e3..1ce9e1a3c1 100644 --- a/src/acb_poly/find_roots.c +++ b/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); } diff --git a/src/acb_poly/refine_roots_durand_kerner.c b/src/acb_poly/refine_roots_durand_kerner.c deleted file mode 100644 index f9280238c5..0000000000 --- a/src/acb_poly/refine_roots_durand_kerner.c +++ /dev/null @@ -1,153 +0,0 @@ -/* - Copyright (C) 2012 Fredrik Johansson - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "acb_poly.h" - -/* we don't need any error bounding, so we define a few helper - functions that ignore the radii */ - -static inline void -acb_sub_mid(acb_t z, const acb_t x, const acb_t y, slong prec) -{ - arf_sub(arb_midref(acb_realref(z)), - arb_midref(acb_realref(x)), - arb_midref(acb_realref(y)), prec, ARF_RND_DOWN); - arf_sub(arb_midref(acb_imagref(z)), - arb_midref(acb_imagref(x)), - arb_midref(acb_imagref(y)), prec, ARF_RND_DOWN); -} - -static inline void -acb_add_mid(acb_t z, const acb_t x, const acb_t y, slong prec) -{ - arf_add(arb_midref(acb_realref(z)), - arb_midref(acb_realref(x)), - arb_midref(acb_realref(y)), prec, ARF_RND_DOWN); - arf_add(arb_midref(acb_imagref(z)), - arb_midref(acb_imagref(x)), - arb_midref(acb_imagref(y)), prec, ARF_RND_DOWN); -} - -static inline void -acb_mul_mid(acb_t z, const acb_t x, const acb_t y, slong prec) -{ -#define a arb_midref(acb_realref(x)) -#define b arb_midref(acb_imagref(x)) -#define c arb_midref(acb_realref(y)) -#define d arb_midref(acb_imagref(y)) -#define e arb_midref(acb_realref(z)) -#define f arb_midref(acb_imagref(z)) - - arf_complex_mul(e, f, a, b, c, d, prec, ARF_RND_DOWN); - -#undef a -#undef b -#undef c -#undef d -#undef e -#undef f -} - -static inline void -acb_inv_mid(acb_t z, const acb_t x, slong prec) -{ - arf_t t; - arf_init(t); - -#define a arb_midref(acb_realref(x)) -#define b arb_midref(acb_imagref(x)) -#define e arb_midref(acb_realref(z)) -#define f arb_midref(acb_imagref(z)) - - arf_mul(t, a, a, prec, ARF_RND_DOWN); - arf_addmul(t, b, b, prec, ARF_RND_DOWN); - - arf_div(e, a, t, prec, ARF_RND_DOWN); - arf_div(f, b, t, prec, ARF_RND_DOWN); - - arf_neg(f, f); - -#undef a -#undef b -#undef e -#undef f - - arf_clear(t); -} - -void -_acb_poly_evaluate_mid(acb_t res, acb_srcptr f, slong len, - const acb_t a, slong prec) -{ - slong i = len - 1; - acb_t t; - - acb_init(t); - acb_set(res, f + i); - - for (i = len - 2; i >= 0; i--) - { - acb_mul_mid(t, res, a, prec); - acb_add_mid(res, f + i, t, prec); - } - - acb_clear(t); -} - -void -_acb_poly_refine_roots_durand_kerner(acb_ptr roots, - acb_srcptr poly, slong len, slong prec) -{ - slong i, j; - - acb_t x, y, t; - - acb_init(x); - acb_init(y); - acb_init(t); - - for (i = 0; i < len - 1; i++) - { - _acb_poly_evaluate_mid(x, poly, len, roots + i, prec); - - acb_set(y, poly + len - 1); - - for (j = 0; j < len - 1; j++) - { - if (i != j) - { - acb_sub_mid(t, roots + i, roots + j, prec); - acb_mul_mid(y, y, t, prec); - } - } - - mag_zero(arb_radref(acb_realref(y))); - mag_zero(arb_radref(acb_imagref(y))); - - /* In the extremely unlikely event that y = 0, don't divide by 0. */ - if (arf_is_zero(arb_midref(acb_realref(y))) && arf_is_zero(arb_midref(acb_imagref(y)))) - { - arf_set_ui_2exp_si(arb_midref(acb_realref(y)), 1, -prec); - arf_set_ui_2exp_si(arb_midref(acb_imagref(y)), 1, -prec); - } - - acb_inv_mid(t, y, prec); - acb_mul_mid(t, t, x, prec); - acb_sub_mid(roots + i, roots + i, t, prec); - - arf_get_mag(arb_radref(acb_realref(roots + i)), arb_midref(acb_realref(t))); - arf_get_mag(arb_radref(acb_imagref(roots + i)), arb_midref(acb_imagref(t))); - } - - acb_clear(x); - acb_clear(y); - acb_clear(t); -} diff --git a/src/python/flint_ctypes.py b/src/python/flint_ctypes.py index a7f513d2eb..d2641b32e3 100644 --- a/src/python/flint_ctypes.py +++ b/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)