diff --git a/doc/source/fmpr.rst b/doc/source/fmpr.rst index 0c0b620f..692af289 100644 --- a/doc/source/fmpr.rst +++ b/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 diff --git a/doc/source/fmprb.rst b/doc/source/fmprb.rst index e35b08a3..968eee57 100644 --- a/doc/source/fmprb.rst +++ b/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 diff --git a/fmpr.h b/fmpr.h index 89d9b4d7..14a1cead 100644 --- a/fmpr.h +++ b/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; \ diff --git a/fmprb.h b/fmprb.h index ded705cc..a165a537 100644 --- a/fmprb.h +++ b/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) { diff --git a/fmprb/get_unique_fmpz.c b/fmprb/get_unique_fmpz.c new file mode 100644 index 00000000..3ba80250 --- /dev/null +++ b/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; + } + } +} + diff --git a/fmprb/test/t-get_unique_fmpz.c b/fmprb/test/t-get_unique_fmpz.c new file mode 100644 index 00000000..57df2efb --- /dev/null +++ b/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; +}