add arb_can_round_mpfr

arb.h

@@ -529,6 +529,9 @@ void arb_get_interval_mpfr(mpfr_t a, mpfr_t b, const arb_t x);
 void arb_union(arb_t z, const arb_t x, const arb_t y, slong prec);
 void arb_get_rand_fmpq(fmpq_t q, flint_rand_t state, const arb_t x, slong bits);
 
+int arb_can_round_arf(const arb_t x, slong prec, arf_rnd_t rnd);
+int arb_can_round_mpfr(const arb_t x, slong prec, mpfr_rnd_t rnd);
+
 void arb_add(arb_t z, const arb_t x, const arb_t y, slong prec);
 void arb_add_arf(arb_t z, const arb_t x, const arf_t y, slong prec);
 void arb_add_ui(arb_t z, const arb_t x, ulong y, slong prec);

arb/can_round_mpfr.c

@@ -0,0 +1,72 @@
+/*=============================================================================
+
+    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) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb.h"
+
+int mpfr_round_p(mp_srcptr, mp_size_t, mpfr_exp_t, mpfr_prec_t);
+
+int
+arb_can_round_arf(const arb_t x, slong prec, arf_rnd_t rnd)
+{
+    return arb_can_round_mpfr(x, prec, arf_rnd_to_mpfr(rnd));
+}
+
+int
+arb_can_round_mpfr(const arb_t x, slong prec, mpfr_rnd_t rnd)
+{
+    if (!arb_is_finite(x))
+    {
+        return 0;
+    }
+    else if (mag_is_zero(arb_radref(x)))
+    {
+        return 1;
+    }
+    else if (arf_is_zero(arb_midref(x)))
+    {
+        return 0;
+    }
+    else
+    {
+        slong e, bits;
+        mp_size_t n;
+        mp_srcptr d;
+
+        e = _fmpz_sub_small(ARF_EXPREF(arb_midref(x)), MAG_EXPREF(arb_radref(x)));
+
+        if (e < prec)
+            return 0;
+
+        /* The relative exponent could be tiny (in which case _fmpz_sub_small
+           has clamped it). Looking just past the end will be enough. */
+        bits = arb_bits(x);
+        e = FLINT_MIN(e, FLINT_MAX(bits, prec) + 10);
+
+        ARF_GET_MPN_READONLY(d, n, arb_midref(x));
+
+        return mpfr_round_p(d, n, e, prec + (rnd == MPFR_RNDN));
+    }
+}
+

arb/test/t-can_round_mpfr.c

@@ -0,0 +1,94 @@
+/*=============================================================================
+
+    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) 2016 Fredrik Johansson
+
+******************************************************************************/
+
+#include "arb.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("can_round_mpfr....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000000; iter++)
+    {
+        mpfr_t x, y1, y2;
+        int r1, r2;
+        arb_t t;
+        slong prec;
+        mpfr_rnd_t rnd;
+
+        prec = 2 + n_randint(state, 300);
+
+        mpfr_init2(x, 2 + n_randint(state, 300));
+        mpfr_init2(y1, prec);
+        mpfr_init2(y2, prec);
+
+        arb_init(t);
+
+        switch (n_randint(state, 5))
+        {
+            case 0: rnd = MPFR_RNDN; break;
+            case 1: rnd = MPFR_RNDZ; break;
+            case 2: rnd = MPFR_RNDU; break;
+            case 3: rnd = MPFR_RNDD; break;
+            default: rnd = MPFR_RNDA;
+        }
+
+        arf_randtest(arb_midref(t), state, mpfr_get_prec(x), 1 + n_randint(state, 10));
+        arf_abs(arb_midref(t), arb_midref(t));
+        arf_get_mpfr(x, arb_midref(t), MPFR_RNDN);
+        arb_root_ui(t, t, 4, 2 + n_randint(state, 300));
+
+        if (arb_can_round_mpfr(t, prec, rnd))
+        {
+            r1 = mpfr_root(y1, x, 4, rnd);
+            r2 = arf_get_mpfr(y2, arb_midref(t), rnd);
+
+            if (r1 != r2 || !mpfr_equal_p(y1, y2))
+            {
+                flint_printf("FAIL! %ld\n");
+                flint_printf("r1 = %d, r2 = %d, prec = %wd\n", r1, r2, prec);
+                flint_printf("x = "); mpfr_dump(x); flint_printf("\n");
+                flint_printf("y1 = "); mpfr_dump(y1); flint_printf("\n");
+                flint_printf("y2 = "); mpfr_dump(y2); flint_printf("\n");
+                abort();
+            }
+        }
+
+        arb_clear(t);
+        mpfr_clear(x);
+        mpfr_clear(y1);
+        mpfr_clear(y2);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

doc/source/arb.rst

@@ -425,6 +425,32 @@ Radius and interval operations
     digits. If *x* is infinite or exactly zero, the outputs are all set
     to zero.
 
+.. function:: int arb_can_round_arf(const arb_t x, slong prec, arf_rnd_t rnd)
+
+.. function:: int arb_can_round_mpfr(const arb_t x, slong prec, mpfr_rnd_t rnd)
+
+    Returns nonzero if rounding the midpoint of *x* to *prec* bits in
+    the direction *rnd* is guaranteed to give the unique correctly
+    rounded floating-point approximation for the real number represented by *x*.
+
+    In other words, if this function returns nonzero, applying
+    :func:`arf_set_round`, or :func:`arf_get_mpfr`, or :func:`arf_get_d`
+    to the midpoint of *x* is guaranteed to return a correctly rounded *arf_t*,
+    *mpfr_t* (provided that *prec* is the precision of the output variable),
+    or *double* (provided that *prec* is 53).
+    Moreover, :func:`arf_get_mpfr` is guaranteed to return the correct ternary
+    value according to MPFR semantics.
+
+    Note that the *mpfr* version of this function takes an MPFR rounding mode
+    symbol as input, while the *arf* version takes an *arf* rounding mode
+    symbol. Otherwise, the functions are identical.
+
+    This function may perform a fast, inexact test; that is, it may return
+    zero in some cases even when correct rounding actually is possible.
+
+    To be conservative, zero is returned when *x* is non-finite, even if it
+    is an "exact" infinity.
+
 Comparisons
 -------------------------------------------------------------------------------