hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-04 17:01:40 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

implement the Fujiwara root bound

Browse source
This commit is contained in:
Fredrik Johansson 2015-07-22 13:24:08 +02:00
parent 6c81cba454
commit 59203d90ef
8 changed files with 388 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
acb_poly.h
View file

@ -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,

75
acb_poly/root_bound_fujiwara.c Normal file
View file

@ -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);
}

97
acb_poly/test/t-root_bound_fujiwara.c Normal file
View file

@ -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;
}

4
arb_poly.h
View file

@ -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 */

75
arb_poly/root_bound_fujiwara.c Normal file
View file

@ -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);
}

97
arb_poly/test/t-root_bound_fujiwara.c Normal file
View file

@ -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;
}

18
doc/source/acb_poly.rst
View file

@ -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

18
doc/source/arb_poly.rst
View file

@ -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