add helper functions for root isolation

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;

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

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

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

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)