From 59203d90ef6341f32a83dca9e1b43508642bc5f9 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Wed, 22 Jul 2015 13:24:08 +0200 Subject: [PATCH] implement the Fujiwara root bound --- acb_poly.h | 4 ++ acb_poly/root_bound_fujiwara.c | 75 +++++++++++++++++++++ acb_poly/test/t-root_bound_fujiwara.c | 97 +++++++++++++++++++++++++++ arb_poly.h | 4 ++ arb_poly/root_bound_fujiwara.c | 75 +++++++++++++++++++++ arb_poly/test/t-root_bound_fujiwara.c | 97 +++++++++++++++++++++++++++ doc/source/acb_poly.rst | 18 +++++ doc/source/arb_poly.rst | 18 +++++ 8 files changed, 388 insertions(+) create mode 100644 acb_poly/root_bound_fujiwara.c create mode 100644 acb_poly/test/t-root_bound_fujiwara.c create mode 100644 arb_poly/root_bound_fujiwara.c create mode 100644 arb_poly/test/t-root_bound_fujiwara.c diff --git a/acb_poly.h b/acb_poly.h index de445a2a..28f482e5 100644 --- a/acb_poly.h +++ b/acb_poly.h @@ -457,6 +457,10 @@ long acb_poly_find_roots(acb_ptr roots, const acb_poly_t poly, acb_srcptr initial, long maxiter, long prec); +void _acb_poly_root_bound_fujiwara(mag_t bound, acb_srcptr poly, long len); + +void acb_poly_root_bound_fujiwara(mag_t bound, acb_poly_t poly); + /* Special functions */ void _acb_poly_pow_ui_trunc_binexp(acb_ptr res, diff --git a/acb_poly/root_bound_fujiwara.c b/acb_poly/root_bound_fujiwara.c new file mode 100644 index 00000000..608faf39 --- /dev/null +++ b/acb_poly/root_bound_fujiwara.c @@ -0,0 +1,75 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2015 Fredrik Johansson + +******************************************************************************/ + +#include "acb_poly.h" + +void +_acb_poly_root_bound_fujiwara(mag_t bound, acb_srcptr poly, long len) +{ + mag_t t, u, v; + long i; + + if (len <= 1) + { + mag_inf(bound); + return; + } + + mag_init(t); + mag_init(u); + mag_init(v); + + /* u = 1/leading */ + acb_get_mag_lower(t, poly + len - 1); + mag_one(u); + mag_div(u, u, t); + + mag_zero(v); + + for (i = 0; i < len - 1; i++) + { + acb_get_mag(t, poly + len - 2 - i); + mag_mul(t, t, u); + + if (i == len - 2) + mag_mul_2exp_si(t, t, -1); + + mag_root(t, t, i + 1); + mag_max(v, v, t); + } + + mag_mul_2exp_si(bound, v, 1); + + mag_clear(t); + mag_clear(u); + mag_clear(v); +} + +void +acb_poly_root_bound_fujiwara(mag_t bound, acb_poly_t poly) +{ + _acb_poly_root_bound_fujiwara(bound, poly->coeffs, poly->length); +} + diff --git a/acb_poly/test/t-root_bound_fujiwara.c b/acb_poly/test/t-root_bound_fujiwara.c new file mode 100644 index 00000000..ab1c3add --- /dev/null +++ b/acb_poly/test/t-root_bound_fujiwara.c @@ -0,0 +1,97 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2015 Fredrik Johansson + +******************************************************************************/ + +#include "acb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("root_bound_fujiwara...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + acb_poly_t a; + acb_ptr roots; + acb_t t; + mag_t mag1, mag2; + long i, deg, prec; + + prec = 10 + n_randint(state, 400); + deg = n_randint(state, 10); + + acb_init(t); + acb_poly_init(a); + mag_init(mag1); + mag_init(mag2); + roots = _acb_vec_init(deg); + + for (i = 0; i < deg; i++) + acb_randtest(roots + i, state, prec, 1 + n_randint(state, 20)); + + acb_poly_product_roots(a, roots, deg, prec); + acb_randtest(t, state, prec, 1 + n_randint(state, 20)); + _acb_vec_scalar_mul(a->coeffs, a->coeffs, a->length, t, prec); + + acb_poly_root_bound_fujiwara(mag1, a); + + for (i = 0; i < deg; i++) + { + acb_get_mag(mag2, roots + i); + + /* acb_get_mag gives an upper bound which due to rounding + could be larger than mag1, so we pick a slightly + smaller number */ + mag_mul_ui(mag2, mag2, 10000); + mag_div_ui(mag2, mag2, 10001); + + if (mag_cmp(mag2, mag1) > 0) + { + printf("FAIL\n"); + printf("a = "); acb_poly_printd(a, 15); printf("\n\n"); + printf("root = "); acb_printd(roots + i, 15); printf("\n\n"); + printf("mag1 = "); mag_printd(mag1, 10); printf("\n\n"); + printf("mag2 = "); mag_printd(mag2, 10); printf("\n\n"); + abort(); + } + } + + _acb_vec_clear(roots, deg); + acb_clear(t); + acb_poly_clear(a); + mag_clear(mag1); + mag_clear(mag2); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arb_poly.h b/arb_poly.h index 3b87947c..7a2111d6 100644 --- a/arb_poly.h +++ b/arb_poly.h @@ -630,6 +630,10 @@ void _arb_poly_newton_refine_root(arb_t r, arb_srcptr poly, long eval_extra_prec, long prec); +void _arb_poly_root_bound_fujiwara(mag_t bound, arb_srcptr poly, long len); + +void arb_poly_root_bound_fujiwara(mag_t bound, arb_poly_t poly); + /* Macros */ diff --git a/arb_poly/root_bound_fujiwara.c b/arb_poly/root_bound_fujiwara.c new file mode 100644 index 00000000..871fc24c --- /dev/null +++ b/arb_poly/root_bound_fujiwara.c @@ -0,0 +1,75 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2015 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +void +_arb_poly_root_bound_fujiwara(mag_t bound, arb_srcptr poly, long len) +{ + mag_t t, u, v; + long i; + + if (len <= 1) + { + mag_inf(bound); + return; + } + + mag_init(t); + mag_init(u); + mag_init(v); + + /* u = 1/leading */ + arb_get_mag_lower(t, poly + len - 1); + mag_one(u); + mag_div(u, u, t); + + mag_zero(v); + + for (i = 0; i < len - 1; i++) + { + arb_get_mag(t, poly + len - 2 - i); + mag_mul(t, t, u); + + if (i == len - 2) + mag_mul_2exp_si(t, t, -1); + + mag_root(t, t, i + 1); + mag_max(v, v, t); + } + + mag_mul_2exp_si(bound, v, 1); + + mag_clear(t); + mag_clear(u); + mag_clear(v); +} + +void +arb_poly_root_bound_fujiwara(mag_t bound, arb_poly_t poly) +{ + _arb_poly_root_bound_fujiwara(bound, poly->coeffs, poly->length); +} + diff --git a/arb_poly/test/t-root_bound_fujiwara.c b/arb_poly/test/t-root_bound_fujiwara.c new file mode 100644 index 00000000..37cdcbad --- /dev/null +++ b/arb_poly/test/t-root_bound_fujiwara.c @@ -0,0 +1,97 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2015 Fredrik Johansson + +******************************************************************************/ + +#include "arb_poly.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("root_bound_fujiwara...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + arb_poly_t a; + arb_ptr roots; + arb_t t; + mag_t mag1, mag2; + long i, deg, prec; + + prec = 10 + n_randint(state, 400); + deg = n_randint(state, 10); + + arb_init(t); + arb_poly_init(a); + mag_init(mag1); + mag_init(mag2); + roots = _arb_vec_init(deg); + + for (i = 0; i < deg; i++) + arb_randtest(roots + i, state, prec, 1 + n_randint(state, 20)); + + arb_poly_product_roots(a, roots, deg, prec); + arb_randtest(t, state, prec, 1 + n_randint(state, 20)); + _arb_vec_scalar_mul(a->coeffs, a->coeffs, a->length, t, prec); + + arb_poly_root_bound_fujiwara(mag1, a); + + for (i = 0; i < deg; i++) + { + arb_get_mag(mag2, roots + i); + + /* arb_get_mag gives an upper bound which due to rounding + could be larger than mag1, so we pick a slightly + smaller number */ + mag_mul_ui(mag2, mag2, 10000); + mag_div_ui(mag2, mag2, 10001); + + if (mag_cmp(mag2, mag1) > 0) + { + printf("FAIL\n"); + printf("a = "); arb_poly_printd(a, 15); printf("\n\n"); + printf("root = "); arb_printd(roots + i, 15); printf("\n\n"); + printf("mag1 = "); mag_printd(mag1, 10); printf("\n\n"); + printf("mag2 = "); mag_printd(mag2, 10); printf("\n\n"); + abort(); + } + } + + _arb_vec_clear(roots, deg); + arb_clear(t); + arb_poly_clear(a); + mag_clear(mag1); + mag_clear(mag2); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/doc/source/acb_poly.rst b/doc/source/acb_poly.rst index 82f765c4..81fefeeb 100644 --- a/doc/source/acb_poly.rst +++ b/doc/source/acb_poly.rst @@ -954,6 +954,24 @@ Other special functions Root-finding ------------------------------------------------------------------------------- +.. function:: void _acb_poly_root_bound_fujiwara(mag_t bound, acb_srcptr poly, long len) + +.. function:: void acb_poly_root_bound_fujiwara(mag_t bound, acb_poly_t poly) + + Sets *bound* to an upper bound for the magnitude of all the complex + roots of *poly*. Uses Fujiwara's bound + + .. math :: + + 2 \max \left\{\left|\frac{a_{n-1}}{a_n}\right|, + \left|\frac{a_{n-2}}{a_n}\right|^{1/2}, + \cdots, + \left|\frac{a_1}{a_n}\right|^{1/(n-1)}, + \left|\frac{a_0}{2a_n}\right|^{1/n} + \right\} + + where `a_0, \ldots, a_n` are the coefficients of *poly*. + .. function:: void _acb_poly_root_inclusion(acb_t r, const acb_t m, acb_srcptr poly, acb_srcptr polyder, long len, long prec) Given any complex number `m`, and a nonconstant polynomial `f` and its diff --git a/doc/source/arb_poly.rst b/doc/source/arb_poly.rst index d3f90c19..61d85ccf 100644 --- a/doc/source/arb_poly.rst +++ b/doc/source/arb_poly.rst @@ -947,6 +947,24 @@ Zeta function Root-finding ------------------------------------------------------------------------------- +.. function:: void _arb_poly_root_bound_fujiwara(mag_t bound, arb_srcptr poly, long len) + +.. function:: void arb_poly_root_bound_fujiwara(mag_t bound, arb_poly_t poly) + + Sets *bound* to an upper bound for the magnitude of all the complex + roots of *poly*. Uses Fujiwara's bound + + .. math :: + + 2 \max \left\{\left|\frac{a_{n-1}}{a_n}\right|, + \left|\frac{a_{n-2}}{a_n}\right|^{1/2}, + \cdots, + \left|\frac{a_1}{a_n}\right|^{1/(n-1)}, + \left|\frac{a_0}{2a_n}\right|^{1/n} + \right\} + + where `a_0, \ldots, a_n` are the coefficients of *poly*. + .. function:: void _arb_poly_newton_convergence_factor(arf_t convergence_factor, arb_srcptr poly, long len, const arb_t convergence_interval, long prec) Given an interval `I` specified by *convergence_interval*, evaluates a bound