acb_poly_validate_real_roots

acb_poly.h

@@ -461,6 +461,10 @@ 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);
 
+int _acb_poly_validate_real_roots(acb_srcptr roots, acb_srcptr poly, long len, long prec);
+
+int acb_poly_validate_real_roots(acb_srcptr roots, const acb_poly_t poly, long prec);
+
 /* Special functions */
 
 void _acb_poly_pow_ui_trunc_binexp(acb_ptr res,

acb_poly/validate_real_roots.c

@@ -0,0 +1,138 @@
+/*=============================================================================
+
+    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"
+
+#ifndef __compar_fn_t
+typedef int (*__compar_fn_t) (__const void *, __const void *);
+#endif
+
+int arb_cmp_mid(const arb_t a, const arb_t b)
+{
+    return arf_cmp(arb_midref(a), arb_midref(b));
+}
+
+void _arb_vec_sort_mid(arb_ptr vec, long len)
+{
+    qsort(vec, len, sizeof(arb_struct), (__compar_fn_t) arb_cmp_mid);
+}
+
+int
+_acb_poly_validate_real_roots(acb_srcptr roots, acb_srcptr poly, long len, long prec)
+{
+    long i, deg, num_real;
+    arb_ptr real;
+    int result;
+
+    deg = len - 1;
+    num_real = 0;
+    result = 1;
+
+    if (deg <= 1)
+        return 1;
+
+    real = _arb_vec_init(deg);
+
+    /* pick out the candidate real roots */
+    for (i = 0; i < deg; i++)
+    {
+        if (arb_contains_zero(acb_imagref(roots + i)))
+        {
+            arb_set(real + num_real, acb_realref(roots + i));
+            num_real++;
+        }
+    }
+
+    /* number of real roots must be even if the polynomial is even,
+       and odd if the polynomial is odd (unless there are repeated roots...
+       in which case the input is invalid) */
+    if ((num_real % 2) != (deg % 2))
+    {
+        result = 0;
+    }
+    else if (num_real > 0)
+    {
+        int sign_neg_inf, sign_pos_inf, prev_sign;
+
+        acb_t t;
+        acb_init(t);
+
+        /* by assumption that the roots are real and isolated, the lead
+           coefficient really must be known to be either positive or negative */
+        sign_pos_inf = arb_is_positive(acb_realref(poly + deg)) ? 1 : -1;
+        sign_neg_inf = (deg % 2) ? -sign_pos_inf : sign_pos_inf;
+
+        /* now we check that there's a sign change between each root */
+        _arb_vec_sort_mid(real, num_real);
+
+        prev_sign = sign_neg_inf;
+
+        for (i = 0; i < num_real - 1; i++)
+        {
+            /* set t to the midpoint between the midpoints */
+            arb_zero(acb_imagref(t));
+            arf_add(arb_midref(acb_realref(t)),
+                arb_midref(real + i), arb_midref(real + i + 1), prec, ARF_RND_DOWN);
+            arf_mul_2exp_si(arb_midref(acb_realref(t)), arb_midref(acb_realref(t)), -1);
+            mag_zero(arb_radref(acb_realref(t)));
+
+            /* check that this point really is between both intervals (one interval
+               could be much wider than the other */
+            if (arb_lt(real + i, acb_realref(t)) && arb_lt(acb_realref(t), real + i + 1))
+            {
+                /* check sign change */
+                _acb_poly_evaluate(t, poly, len, t, prec);
+
+                if (prev_sign == 1)
+                    result = arb_is_negative(acb_realref(t));
+                else
+                    result = arb_is_positive(acb_realref(t));
+
+                if (!result)
+                    break;
+
+                prev_sign = -prev_sign;
+            }
+            else
+            {
+                result = 0;
+                break;
+            }
+        }
+
+        acb_clear(t);
+    }
+
+    _arb_vec_clear(real, deg);
+
+    return result;
+}
+
+int
+acb_poly_validate_real_roots(acb_srcptr roots, const acb_poly_t poly, long prec)
+{
+    return _acb_poly_validate_real_roots(roots, poly->coeffs, poly->length, prec);
+}
+