add fmprb_get_unique_fmpz

doc/source/fmpr.rst

@@ -198,6 +198,12 @@ Comparisons
 
     Compares the absolute values of *x* and *y*.
 
+.. function:: int fmpr_cmp_2exp_si(const fmpr_t x, long e)
+
+.. function:: int fmpr_cmpabs_2exp_si(const fmpr_t x, long e)
+
+    Compares *x* (respectively its absolute value) with `2^e`.
+
 .. function:: int fmpr_sgn(const fmpr_t x)
 
     Returns `-1`, `0` or `+1` according to the sign of *x*. The sign

doc/source/fmprb.rst

@@ -271,6 +271,12 @@ Precision and comparisons
     or if the difference in magnitude between the midpoint and radius
     is so large that representing the endpoints exactly would cause overflows.
 
+.. function:: int fmprb_get_unique_fmpz(fmpz_t z, const fmprb_t x)
+
+    If *x* contains a unique integer, sets *z* to that value and returns
+    nonzero. Otherwise (if *x* represents no integers or more than one integer),
+    returns zero.
+
 .. function:: long fmprb_rel_error_bits(const fmprb_t x)
 
     Returns the effective relative error of *x* measured in bits, defined as

fmpr.h

@@ -207,7 +207,6 @@ int fmpr_cmp(const fmpr_t x, const fmpr_t y);
 
 int fmpr_cmpabs(const fmpr_t x, const fmpr_t y);
 
-
 void fmpr_randtest(fmpr_t x, flint_rand_t state, long bits, long exp_bits);
 
 void fmpr_randtest_not_zero(fmpr_t x, flint_rand_t state, long bits, long exp_bits);
@@ -580,6 +579,28 @@ fmpr_set_ui_2exp_si(fmpr_t x, ulong man, long exp)
     fmpr_mul_2exp_si(x, x, exp);
 }
 
+static __inline__ int fmpr_cmp_2exp_si(const fmpr_t x, long e)
+{
+    fmpr_t t;
+    int res;
+    fmpr_init(t);
+    fmpr_set_ui_2exp_si(t, 1, e);
+    res = fmpr_cmp(x, t);
+    fmpr_clear(t);
+    return res;
+}
+
+static __inline__ int fmpr_cmpabs_2exp_si(const fmpr_t x, long e)
+{
+    fmpr_t t;
+    int res;
+    fmpr_init(t);
+    fmpr_set_ui_2exp_si(t, 1, e);
+    res = fmpr_cmpabs(x, t);
+    fmpr_clear(t);
+    return res;
+}
+
 #define CALL_MPFR_FUNC(r, func, y, x, prec, rnd) \
     do { \
         mpfr_t __t, __u; \

fmprb.h

@@ -452,6 +452,8 @@ fmprb_get_abs_lbound_fmpr(fmpr_t u, const fmprb_t x, long prec)
 
 void fmprb_get_interval_fmpz_2exp(fmpz_t a, fmpz_t b, fmpz_t exp, const fmprb_t x);
 
+int fmprb_get_unique_fmpz(fmpz_t z, const fmprb_t x);
+
 static __inline__ long
 fmprb_rel_error_bits(const fmprb_t x)
 {

fmprb/get_unique_fmpz.c

@@ -0,0 +1,110 @@
+/*=============================================================================
+
+    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 "fmprb.h"
+
+int
+fmprb_get_unique_fmpz(fmpz_t z, const fmprb_t x)
+{
+    if (fmpr_is_inf(fmprb_radref(x)) || fmpr_is_nan(fmprb_radref(x)))
+    {
+        return 0;
+    }
+    else if (fmpr_is_zero(fmprb_radref(x)))
+    {
+        /* x = b*2^e, e >= 0 */
+        if (fmpz_sgn(fmpr_expref(fmprb_midref(x))) >= 0)
+        {
+            if (!fmpz_fits_si(fmpr_expref(fmprb_midref(x))))
+            {
+                printf("fmprb_get_unique_fmpz: too large shift\n");
+                abort();
+            }
+
+            fmpz_mul_2exp(z, fmpr_manref(fmprb_midref(x)),
+                fmpz_get_si(fmpr_expref(fmprb_midref(x))));
+
+            return 1;
+        }
+        else
+        {
+            return 0;
+        }
+    }
+    /* if the radius is >= 1, there are at least two integers */
+    else if (fmpr_cmp_2exp_si(fmprb_radref(x), 0) >= 0)
+    {
+        return 0;
+    }
+    /* there are 0 or 1 integers if the radius is < 1 */
+    else
+    {
+        /* if the midpoint is exactly an integer, it is what we want */
+        if (fmpz_sgn(fmpr_expref(fmprb_midref(x))) >= 0)
+        {
+            if (!fmpz_fits_si(fmpr_expref(fmprb_midref(x))))
+            {
+                printf("fmprb_get_unique_fmpz: too large shift\n");
+                abort();
+            }
+
+            fmpz_mul_2exp(z, fmpr_manref(fmprb_midref(x)),
+                fmpz_get_si(fmpr_expref(fmprb_midref(x))));
+
+            return 1;
+        }
+        /* if the radius is tiny, it can't be an integer */
+        else if (!fmpz_fits_si(fmpr_expref(fmprb_radref(x))))
+        {
+            return 0;
+        }
+        else
+        {
+            fmpz_t a, b, exp;
+            int res;
+
+            fmpz_init(a);
+            fmpz_init(b);
+            fmpz_init(exp);
+
+            fmprb_get_interval_fmpz_2exp(a, b, exp, x);
+
+            fmpz_cdiv_q_2exp(a, a, -fmpz_get_si(exp));
+            fmpz_fdiv_q_2exp(b, b, -fmpz_get_si(exp));
+
+            res = fmpz_equal(a, b);
+
+            if (res)
+                fmpz_set(z, a);
+
+            fmpz_clear(a);
+            fmpz_clear(b);
+            fmpz_clear(exp);
+
+            return res;
+        }
+    }
+}
+

fmprb/test/t-get_unique_fmpz.c

@@ -0,0 +1,124 @@
+/*=============================================================================
+
+    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 "fmprb.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("get_unique_fmpz....");
+    fflush(stdout);
+    flint_randinit(state);
+
+    for (iter = 0; iter < 100000; iter++)
+    {
+        fmprb_t x, z;
+        fmpz_t y, a, b, exp;
+        int unique, unique2;
+
+        fmprb_init(x);
+        fmprb_init(z);
+        fmpz_init(y);
+        fmpz_init(a);
+        fmpz_init(b);
+        fmpz_init(exp);
+
+        fmprb_randtest(x, state, 2000, 10);
+
+        unique = fmprb_get_unique_fmpz(y, x);
+
+        fmprb_get_interval_fmpz_2exp(a, b, exp, x);
+        if (fmpz_sgn(exp) >= 0)
+        {
+            fmpz_mul_2exp(a, a, fmpz_get_si(exp));
+            fmpz_mul_2exp(b, b, fmpz_get_si(exp));
+        }
+        else
+        {
+            fmpz_cdiv_q_2exp(a, a, -fmpz_get_si(exp));
+            fmpz_fdiv_q_2exp(b, b, -fmpz_get_si(exp));
+        }
+        unique2 = fmpz_equal(a, b);
+
+        if ((unique != unique2) || (unique && !fmpz_equal(y, a)))
+        {
+            printf("FAIL:\n\n");
+            printf("x = "); fmprb_print(x); printf("\n\n");
+            printf("unique = %d, unique2 = %d\n\n", unique, unique2);
+            printf("y = "); fmpz_print(y); printf("\n\n");
+            printf("a = "); fmpz_print(a); printf("\n\n");
+            printf("b = "); fmpz_print(b); printf("\n\n");
+            printf(" exp = "); fmpz_print(exp); printf("\n\n");
+            abort();
+        }
+
+        if (unique)
+        {
+            fmprb_set_fmpz(z, y);
+            fmprb_set_round(z, z, 2 + n_randint(state, 1000));
+
+            if (!fmprb_overlaps(x, z))
+            {
+                printf("FAIL (overlap):\n\n");
+                printf("x = "); fmprb_print(x); printf("\n\n");
+                printf("y = "); fmpz_print(y); printf("\n\n");
+                printf("z = "); fmprb_print(z); printf("\n\n");
+                abort();
+            }
+
+            fmpz_add_ui(b, y, 1);
+            if (fmprb_contains_fmpz(x, b))
+            {
+                printf("FAIL (contains a + 1):\n\n");
+                printf("x = "); fmprb_print(x); printf("\n\n");
+                printf("y = "); fmpz_print(y); printf("\n\n");
+                abort();
+            }
+
+            fmpz_sub_ui(b, y, 1);
+            if (fmprb_contains_fmpz(x, b))
+            {
+                printf("FAIL (contains a - 1):\n\n");
+                printf("x = "); fmprb_print(x); printf("\n\n");
+                printf("y = "); fmpz_print(y); printf("\n\n");
+                abort();
+            }
+        }
+
+        fmprb_clear(x);
+        fmprb_clear(z);
+        fmpz_clear(y);
+        fmpz_clear(a);
+        fmpz_clear(b);
+        fmpz_clear(exp);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}