From 30e004ecb1e422fdf7ea70e18a4e6d572d865cb7 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Mon, 15 Oct 2012 17:31:04 +0200 Subject: [PATCH] add helper functions for root isolation --- fmpcb.h | 2 +- fmpcb_poly.h | 6 + fmpcb_poly/refine_roots_durand_kerner.c | 191 ++++++++++++++++++++++++ fmpcb_poly/validate_roots.c | 109 ++++++++++++++ fmprb.h | 2 +- 5 files changed, 308 insertions(+), 2 deletions(-) create mode 100644 fmpcb_poly/refine_roots_durand_kerner.c create mode 100644 fmpcb_poly/validate_roots.c diff --git a/fmpcb.h b/fmpcb.h index c8bea00f..b645d903 100644 --- a/fmpcb.h +++ b/fmpcb.h @@ -132,7 +132,7 @@ fmpcb_overlaps(const fmpcb_t x, const fmpcb_t y) fmprb_overlaps(fmpcb_imagref(x), fmpcb_imagref(y)); } -void +static __inline__ void fmpcb_get_abs_ubound_fmpr(fmpr_t u, const fmpcb_t z, long prec) { fmpr_t v; diff --git a/fmpcb_poly.h b/fmpcb_poly.h index 8e4499c9..bda678f4 100644 --- a/fmpcb_poly.h +++ b/fmpcb_poly.h @@ -98,4 +98,10 @@ void _fmpcb_poly_root_inclusion(fmpcb_t r, const fmpcb_t m, const fmpcb_struct * poly, const fmpcb_struct * polyder, long len, long prec); +long _fmpcb_poly_validate_roots(fmpcb_struct * roots, + const fmpcb_struct * poly, long len, long prec); + +void _fmpcb_poly_refine_roots_durand_kerner(fmpcb_struct * roots, + const fmpcb_struct * poly, long len, long prec); + #endif diff --git a/fmpcb_poly/refine_roots_durand_kerner.c b/fmpcb_poly/refine_roots_durand_kerner.c new file mode 100644 index 00000000..acd39010 --- /dev/null +++ b/fmpcb_poly/refine_roots_durand_kerner.c @@ -0,0 +1,191 @@ +/*============================================================================= + + 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) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + +/* we don't need any error bounding, so we define a few helper + functions that ignore the radii */ + +static __inline__ void +fmpcb_sub_mid(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec) +{ + fmpr_sub(fmprb_midref(fmpcb_realref(z)), + fmprb_midref(fmpcb_realref(x)), + fmprb_midref(fmpcb_realref(y)), prec, FMPR_RND_DOWN); + fmpr_sub(fmprb_midref(fmpcb_imagref(z)), + fmprb_midref(fmpcb_imagref(x)), + fmprb_midref(fmpcb_imagref(y)), prec, FMPR_RND_DOWN); +} + +static __inline__ void +fmpcb_add_mid(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec) +{ + fmpr_add(fmprb_midref(fmpcb_realref(z)), + fmprb_midref(fmpcb_realref(x)), + fmprb_midref(fmpcb_realref(y)), prec, FMPR_RND_DOWN); + fmpr_add(fmprb_midref(fmpcb_imagref(z)), + fmprb_midref(fmpcb_imagref(x)), + fmprb_midref(fmpcb_imagref(y)), prec, FMPR_RND_DOWN); +} + +static __inline__ void +fmpcb_mul_mid(fmpcb_t z, const fmpcb_t x, const fmpcb_t y, long prec) +{ +#define a fmprb_midref(fmpcb_realref(x)) +#define b fmprb_midref(fmpcb_imagref(x)) +#define c fmprb_midref(fmpcb_realref(y)) +#define d fmprb_midref(fmpcb_imagref(y)) +#define e fmprb_midref(fmpcb_realref(z)) +#define f fmprb_midref(fmpcb_imagref(z)) + + fmpr_t t, u, v; + + fmpr_init(t); + fmpr_init(u); + fmpr_init(v); + + fmpr_add(t, a, b, prec, FMPR_RND_DOWN); + fmpr_add(u, c, d, prec, FMPR_RND_DOWN); + fmpr_mul(v, t, u, prec, FMPR_RND_DOWN); + + fmpr_mul(t, a, c, prec, FMPR_RND_DOWN); + fmpr_mul(u, b, d, prec, FMPR_RND_DOWN); + + fmpr_sub(e, t, u, prec, FMPR_RND_DOWN); + fmpr_sub(f, v, t, prec, FMPR_RND_DOWN); + fmpr_sub(f, f, u, prec, FMPR_RND_DOWN); + + fmpr_clear(t); + fmpr_clear(u); + fmpr_clear(v); + +#undef a +#undef b +#undef c +#undef d +#undef e +#undef f +} + +static __inline__ void +fmpcb_inv_mid(fmpcb_t z, const fmpcb_t x, long prec) +{ + fmpr_t t; + fmpr_init(t); + +#define a fmprb_midref(fmpcb_realref(x)) +#define b fmprb_midref(fmpcb_imagref(x)) +#define e fmprb_midref(fmpcb_realref(z)) +#define f fmprb_midref(fmpcb_imagref(z)) + + fmpr_mul(t, a, a, prec, FMPR_RND_DOWN); + fmpr_addmul(t, b, b, prec, FMPR_RND_DOWN); + + fmpr_div(e, a, t, prec, FMPR_RND_DOWN); + fmpr_div(f, b, t, prec, FMPR_RND_DOWN); + + fmpr_neg(f, f); + +#undef a +#undef b +#undef e +#undef f + + fmpr_clear(t); +} + +void +_fmpcb_poly_evaluate_mid(fmpcb_t res, const fmpcb_struct * f, long len, + const fmpcb_t a, long prec) +{ + long i = len - 1; + fmpcb_t t; + + fmpcb_init(t); + fmpcb_set(res, f + i); + + for (i = len - 2; i >= 0; i--) + { + fmpcb_mul_mid(t, res, a, prec); + fmpcb_add_mid(res, f + i, t, prec); + } + + fmpcb_clear(t); +} + + +/* +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). +*/ + +void +_fmpcb_poly_refine_roots_durand_kerner(fmpcb_struct * roots, + const fmpcb_struct * poly, long len, long prec) +{ + long i, j; + + fmpcb_t x, y, t; + + fmpcb_init(x); + fmpcb_init(y); + fmpcb_init(t); + + for (i = 0; i < len - 1; i++) + { + _fmpcb_poly_evaluate_mid(x, poly, len, roots + i, prec); + + fmpcb_set(y, poly + len - 1); + + for (j = 0; j < len - 1; j++) + { + if (i != j) + { + fmpcb_sub_mid(t, roots + i, roots + j, prec); + fmpcb_mul_mid(y, y, t, prec); + } + } + + fmpr_zero(fmprb_radref(fmpcb_realref(y))); + fmpr_zero(fmprb_radref(fmpcb_imagref(y))); + + fmpcb_inv_mid(t, y, prec); + fmpcb_mul_mid(t, t, x, prec); + + fmpcb_sub_mid(roots + i, roots + i, t, prec); + + fmpr_set_round(fmprb_radref(fmpcb_realref(roots + i)), + fmprb_midref(fmpcb_realref(t)), FMPRB_RAD_PREC, FMPR_RND_UP); + fmpr_set_round(fmprb_radref(fmpcb_imagref(roots + i)), + fmprb_midref(fmpcb_imagref(t)), FMPRB_RAD_PREC, FMPR_RND_UP); + } + + fmpcb_clear(x); + fmpcb_clear(y); + fmpcb_clear(t); +} + diff --git a/fmpcb_poly/validate_roots.c b/fmpcb_poly/validate_roots.c new file mode 100644 index 00000000..4a5d41b2 --- /dev/null +++ b/fmpcb_poly/validate_roots.c @@ -0,0 +1,109 @@ +/*============================================================================= + + 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) 2012 Fredrik Johansson + +******************************************************************************/ + +#include "fmpcb_poly.h" + +/* +Given a list of approximate roots of the input polynomial, this +function sets a rigorous bounding interval for each root, and determines +which roots are isolated from all the other roots. +It then rearranges the list of roots so that the isolated roots +are at the front of the list, and returns the count of isolated roots. + +If the return value equals the degree of the polynomial, then all +roots have been separated. If the return value is smaller, all the +remaining intervals are guaranteed to contain roots, but +it is possible that not all of the polynomial's roots are contained +among them. +*/ + +long +_fmpcb_poly_validate_roots(fmpcb_struct * roots, + const fmpcb_struct * poly, long len, long prec) +{ + long i, j, deg; + long isolated, nonisolated, total_isolated; + fmpcb_struct * deriv; + fmpcb_struct * tmp; + int *overlap; + + deg = len - 1; + + deriv = _fmpcb_vec_init(deg); + overlap = flint_calloc(deg, sizeof(int)); + tmp = flint_malloc(sizeof(fmpcb_struct) * deg); + + _fmpcb_poly_derivative(deriv, poly, len, prec); + + /* compute an inclusion interval for each point */ + for (i = 0; i < deg; i++) + { + _fmpcb_poly_root_inclusion(roots + i, roots + i, + poly, deriv, len, prec); + } + + /* find which points do not overlap with any other points */ + for (i = 0; i < deg; i++) + { + for (j = i + 1; j < deg; j++) + { + if (fmpcb_overlaps(roots + i, roots + j)) + { + overlap[i] = overlap[j] = 1; + } + } + } + + /* count and move all isolated roots to the front of the array */ + total_isolated = 0; + for (i = 0; i < deg; i++) + total_isolated += (overlap[i] == 0); + + for (i = 0; i < deg; i++) + tmp[i] = roots[i]; + + isolated = 0; + nonisolated = 0; + for (i = 0; i < deg; i++) + { + if (overlap[i] == 0) + { + roots[isolated] = tmp[i]; + isolated++; + } + else + { + roots[total_isolated + nonisolated] = tmp[i]; + nonisolated++; + } + } + + _fmpcb_vec_clear(deriv, deg); + flint_free(tmp); + flint_free(overlap); + + return isolated; +} + diff --git a/fmprb.h b/fmprb.h index 295af983..06a8bfbc 100644 --- a/fmprb.h +++ b/fmprb.h @@ -350,7 +350,7 @@ int fmprb_contains_zero(const fmprb_t x); int fmprb_overlaps(const fmprb_t x, const fmprb_t y); -void +static __inline__ void fmprb_get_abs_ubound_fmpr(fmpr_t u, const fmprb_t x, long prec) { if (fmpr_sgn(fmprb_midref(x)) < 0)