From 24a8a6b31f374f07303ae06bebbf2efe57c82dac Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Sun, 27 Apr 2014 00:32:37 +0200 Subject: [PATCH] more functions and test functions --- arf.h | 37 ++++++++++ arf/abs_bound_lt_2exp_si.c | 56 ++++++++++++++ arf/cmp_2exp_si.c | 8 +- arf/cmpabs_2exp_si.c | 8 +- arf/get_fmpr.c | 11 +-- arf/get_fmpz_2exp.c | 72 ++++++++++++++++++ arf/test/t-abs_bound_le_2exp_fmpz.c | 82 +++++++++++++++++++++ arf/test/t-abs_bound_lt_2exp_fmpz.c | 82 +++++++++++++++++++++ arf/test/t-abs_bound_lt_2exp_si.c | 103 ++++++++++++++++++++++++++ arf/test/t-set_fmpr.c | 110 ++++++++++++++++++++++++++++ arf/test/t-set_fmpz_2exp.c | 79 ++++++++++++++++++++ 11 files changed, 626 insertions(+), 22 deletions(-) create mode 100644 arf/abs_bound_lt_2exp_si.c create mode 100644 arf/get_fmpz_2exp.c create mode 100644 arf/test/t-abs_bound_le_2exp_fmpz.c create mode 100644 arf/test/t-abs_bound_lt_2exp_fmpz.c create mode 100644 arf/test/t-abs_bound_lt_2exp_si.c create mode 100644 arf/test/t-set_fmpr.c create mode 100644 arf/test/t-set_fmpz_2exp.c diff --git a/arf.h b/arf.h index 6e80498a..6ef0f4d8 100644 --- a/arf.h +++ b/arf.h @@ -117,6 +117,9 @@ arf_rnd_to_mpfr(arf_rnd_t rnd) /* Encoding for special values. */ #define ARF_IS_SPECIAL(x) (ARF_XSIZE(x) == 0) +/* Value is +/- a power of two */ +#define ARF_IS_POW2(x) (ARF_SIZE(x) == 1) && (ARF_NOPTR_D(x)[0] == LIMB_TOP) + /* To set a special value, first call this and then set the exponent. */ #define ARF_MAKE_SPECIAL(x) \ do { \ @@ -494,6 +497,8 @@ int arf_neg_round(arf_t y, const arf_t x, long prec, arf_rnd_t rnd); void arf_get_fmpr(fmpr_t y, const arf_t x); +void arf_set_fmpr(arf_t y, const fmpr_t x); + int arf_get_mpfr(mpfr_t x, const arf_t y, mpfr_rnd_t rnd); void arf_set_mpfr(arf_t x, const mpfr_t y); @@ -639,6 +644,38 @@ arf_set_round_fmpz_2exp(arf_t y, const fmpz_t x, const fmpz_t exp, long prec, ar } } +static __inline__ void +arf_abs_bound_lt_2exp_fmpz(fmpz_t b, const arf_t x) +{ + if (arf_is_special(x)) + fmpz_zero(b); + else + fmpz_set(b, ARF_EXPREF(x)); +} + +static __inline__ void +arf_abs_bound_le_2exp_fmpz(fmpz_t b, const arf_t x) +{ + if (arf_is_special(x)) + fmpz_zero(b); + else if (ARF_IS_POW2(x)) + fmpz_sub_ui(b, ARF_EXPREF(x), 1); + else + fmpz_set(b, ARF_EXPREF(x)); +} + +long arf_abs_bound_lt_2exp_si(const arf_t x); + +void arf_get_fmpz_2exp(fmpz_t man, fmpz_t exp, const arf_t x); + +static __inline__ void +arf_set_fmpz_2exp(arf_t x, const fmpz_t man, const fmpz_t exp) +{ + arf_set_fmpz(x, man); + if (!arf_is_zero(x)) + fmpz_add_inline(ARF_EXPREF(x), ARF_EXPREF(x), exp); +} + void arf_debug(const arf_t x); #define arf_print arf_debug diff --git a/arf/abs_bound_lt_2exp_si.c b/arf/abs_bound_lt_2exp_si.c new file mode 100644 index 00000000..b3372ac2 --- /dev/null +++ b/arf/abs_bound_lt_2exp_si.c @@ -0,0 +1,56 @@ +/*============================================================================= + + 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) 2014 Fredrik Johansson + +******************************************************************************/ + +#include "arf.h" + +long +arf_abs_bound_lt_2exp_si(const arf_t x) +{ + long res; + + if (arf_is_special(x)) + { + if (arf_is_zero(x)) + return -ARF_PREC_EXACT; + else + return ARF_PREC_EXACT; + } + + if (!COEFF_IS_MPZ(ARF_EXP(x))) + return ARF_EXP(x); + + if (fmpz_fits_si(ARF_EXPREF(x))) + res = fmpz_get_si(ARF_EXPREF(x)); + else + res = fmpz_sgn(ARF_EXPREF(x)) < 0 ? -ARF_PREC_EXACT : ARF_PREC_EXACT; + + if (res < -ARF_PREC_EXACT) + res = -ARF_PREC_EXACT; + if (res > ARF_PREC_EXACT) + res = ARF_PREC_EXACT; + + return res; +} + diff --git a/arf/cmp_2exp_si.c b/arf/cmp_2exp_si.c index c648191d..a42b3123 100644 --- a/arf/cmp_2exp_si.c +++ b/arf/cmp_2exp_si.c @@ -28,8 +28,6 @@ int arf_cmp_2exp_si(const arf_t x, long e) { - int pow2; - if (arf_is_special(x)) { if (arf_is_zero(x)) return -1; @@ -41,18 +39,16 @@ arf_cmp_2exp_si(const arf_t x, long e) if (ARF_SGNBIT(x)) return -1; - pow2 = (ARF_SIZE(x) == 1) && (ARF_NOPTR_D(x)[0] == LIMB_TOP); - /* Fast path. */ if (!COEFF_IS_MPZ(ARF_EXP(x))) { - if (pow2 && (ARF_EXP(x) - 1 == e)) + if (ARF_IS_POW2(x) && (ARF_EXP(x) - 1 == e)) return 0; else return (ARF_EXP(x) <= e) ? -1 : 1; } - if (pow2) + if (ARF_IS_POW2(x)) { fmpz_t t; fmpz_init(t); diff --git a/arf/cmpabs_2exp_si.c b/arf/cmpabs_2exp_si.c index 26198424..0fc527d0 100644 --- a/arf/cmpabs_2exp_si.c +++ b/arf/cmpabs_2exp_si.c @@ -28,8 +28,6 @@ int arf_cmpabs_2exp_si(const arf_t x, long e) { - int pow2; - if (arf_is_special(x)) { if (arf_is_zero(x)) return -1; @@ -37,18 +35,16 @@ arf_cmpabs_2exp_si(const arf_t x, long e) return 0; } - pow2 = (ARF_SIZE(x) == 1) && (ARF_NOPTR_D(x)[0] == LIMB_TOP); - /* Fast path. */ if (!COEFF_IS_MPZ(ARF_EXP(x))) { - if (pow2 && (ARF_EXP(x) - 1 == e)) + if (ARF_IS_POW2(x) && (ARF_EXP(x) - 1 == e)) return 0; else return (ARF_EXP(x) <= e) ? -1 : 1; } - if (pow2) + if (ARF_IS_POW2(x)) { fmpz_t t; fmpz_init(t); diff --git a/arf/get_fmpr.c b/arf/get_fmpr.c index fb2b3f42..ee0a583a 100644 --- a/arf/get_fmpr.c +++ b/arf/get_fmpr.c @@ -41,16 +41,7 @@ arf_get_fmpr(fmpr_t y, const arf_t x) } else { - mp_srcptr xptr; - mp_size_t xn; - long shift; - - ARF_GET_MPN_READONLY(xptr, xn, x); - - _fmpr_set_round_mpn(&shift, fmpr_manref(y), - xptr, xn, ARF_SGNBIT(x), FMPR_PREC_EXACT, FMPR_RND_DOWN); - - fmpz_add_si(fmpr_expref(y), ARF_EXPREF(x), shift - xn * FLINT_BITS); + arf_get_fmpz_2exp(fmpr_manref(y), fmpr_expref(y), x); } } diff --git a/arf/get_fmpz_2exp.c b/arf/get_fmpz_2exp.c new file mode 100644 index 00000000..bc80faab --- /dev/null +++ b/arf/get_fmpz_2exp.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) 2014 Fredrik Johansson + +******************************************************************************/ + +#include "arf.h" + +void +arf_get_fmpz_2exp(fmpz_t man, fmpz_t exp, const arf_t x) +{ + if (arf_is_special(x)) + { + fmpz_zero(man); + fmpz_zero(exp); + } + else + { + mp_srcptr xptr; + mp_size_t xn; + int shift; + + ARF_GET_MPN_READONLY(xptr, xn, x); + + count_trailing_zeros(shift, xptr[0]); + + fmpz_sub_ui(exp, ARF_EXPREF(x), xn * FLINT_BITS - shift); + + if (xn == 1) + { + if (ARF_SGNBIT(x)) + fmpz_neg_ui(man, xptr[0] >> shift); + else + fmpz_set_ui(man, xptr[0] >> shift); + } + else + { + __mpz_struct * z = _fmpz_promote(man); + + if (z->_mp_alloc < xn) + mpz_realloc(z, xn); + + if (shift == 0) + flint_mpn_copyi(z->_mp_d, xptr, xn); + else + mpn_rshift(z->_mp_d, xptr, xn, shift); + + /* top limb cannot be zero */ + z->_mp_size = ARF_SGNBIT(x) ? -xn : xn; + } + } +} + diff --git a/arf/test/t-abs_bound_le_2exp_fmpz.c b/arf/test/t-abs_bound_le_2exp_fmpz.c new file mode 100644 index 00000000..9fdb2dae --- /dev/null +++ b/arf/test/t-abs_bound_le_2exp_fmpz.c @@ -0,0 +1,82 @@ +/*============================================================================= + + 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) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arf.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("abs_bound_le_2exp_fmpz...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + arf_t x, y; + fmpz_t b; + int cmp1, cmp2; + + arf_init(x); + arf_init(y); + fmpz_init(b); + + arf_randtest_not_zero(x, state, 2 + n_randint(state, 1000), 100); + arf_abs_bound_le_2exp_fmpz(b, x); + + arf_one(y); + arf_mul_2exp_fmpz(y, y, b); + + cmp1 = (arf_cmpabs(x, y) <= 0); + + arf_mul_2exp_si(y, y, -1); + + cmp2 = (arf_cmpabs(y, x) < 0); + + arf_mul_2exp_si(y, y, 1); + + if (!cmp1 || !cmp2) + { + printf("FAIL\n\n"); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = "); arf_print(y); printf("\n\n"); + printf("b = "); fmpz_print(b); printf("\n\n"); + printf("cmp1 = %d, cmp2 = %d\n\n", cmp1, cmp2); + abort(); + } + + arf_clear(x); + arf_clear(y); + fmpz_clear(b); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arf/test/t-abs_bound_lt_2exp_fmpz.c b/arf/test/t-abs_bound_lt_2exp_fmpz.c new file mode 100644 index 00000000..37e915bc --- /dev/null +++ b/arf/test/t-abs_bound_lt_2exp_fmpz.c @@ -0,0 +1,82 @@ +/*============================================================================= + + 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) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arf.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("abs_bound_lt_2exp_fmpz...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + arf_t x, y; + fmpz_t b; + int cmp1, cmp2; + + arf_init(x); + arf_init(y); + fmpz_init(b); + + arf_randtest_not_zero(x, state, 2 + n_randint(state, 1000), 100); + arf_abs_bound_lt_2exp_fmpz(b, x); + + arf_one(y); + arf_mul_2exp_fmpz(y, y, b); + + cmp1 = (arf_cmpabs(x, y) < 0); + + arf_mul_2exp_si(y, y, -1); + + cmp2 = (arf_cmpabs(y, x) <= 0); + + arf_mul_2exp_si(y, y, 1); + + if (!cmp1 || !cmp2) + { + printf("FAIL\n\n"); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = "); arf_print(y); printf("\n\n"); + printf("b = "); fmpz_print(b); printf("\n\n"); + printf("cmp1 = %d, cmp2 = %d\n\n", cmp1, cmp2); + abort(); + } + + arf_clear(x); + arf_clear(y); + fmpz_clear(b); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arf/test/t-abs_bound_lt_2exp_si.c b/arf/test/t-abs_bound_lt_2exp_si.c new file mode 100644 index 00000000..d2974bdd --- /dev/null +++ b/arf/test/t-abs_bound_lt_2exp_si.c @@ -0,0 +1,103 @@ +/*============================================================================= + + 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) 2013 Fredrik Johansson + +******************************************************************************/ + +#include "arf.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("abs_bound_lt_2exp_si...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + arf_t x; + fmpz_t b; + long c; + + arf_init(x); + fmpz_init(b); + + if (iter < 5) + { + arf_set_ui_2exp_si(x, 1, LONG_MIN + 2); + arf_mul_2exp_si(x, x, -iter); + } + else if (iter < 10) + { + arf_set_ui_2exp_si(x, 1, LONG_MAX - 2); + arf_mul_2exp_si(x, x, iter - 5); + } + else + { + arf_randtest_special(x, state, 2 + n_randint(state, 1000), 100); + } + + arf_abs_bound_lt_2exp_fmpz(b, x); + c = arf_abs_bound_lt_2exp_si(x); + + if (c == -ARF_PREC_EXACT) + { + if (!(arf_is_zero(x) || fmpz_cmp_si(b, -ARF_PREC_EXACT) <= 0)) + { + printf("FAIL (small/zero)\n\n"); + abort(); + } + } + else if (c == ARF_PREC_EXACT) + { + if (!(arf_is_inf(x) || arf_is_nan(x) || + fmpz_cmp_si(b, ARF_PREC_EXACT) >= 0)) + { + printf("FAIL (large/inf/nan)\n\n"); + abort(); + } + } + else + { + if (fmpz_cmp_si(b, c) != 0) + { + printf("FAIL (normal)\n\n"); + printf("x = "); arf_print(x); printf("\n\n"); + printf("b = "); fmpz_print(b); printf("\n\n"); + printf("c = %ld\n\n", c); + abort(); + } + } + + arf_clear(x); + fmpz_clear(b); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arf/test/t-set_fmpr.c b/arf/test/t-set_fmpr.c new file mode 100644 index 00000000..6895920a --- /dev/null +++ b/arf/test/t-set_fmpr.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 "arf.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("set_fmpr...."); + fflush(stdout); + + flint_randinit(state); + + /* test exact roundtrip R -> Q -> R */ + for (iter = 0; iter < 100000; iter++) + { + long bits; + arf_t x, z; + fmpr_t y; + + bits = 2 + n_randint(state, 1000); + + arf_init(x); + arf_init(z); + fmpr_init(y); + + arf_randtest_special(x, state, bits, 1 + n_randint(state, 100)); + arf_randtest_special(z, state, bits, 1 + n_randint(state, 100)); + + arf_get_fmpr(y, x); + arf_set_fmpr(z, y); + + if (!arf_equal(x, z)) + { + printf("FAIL\n\n"); + printf("bits: %ld\n", bits); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = "); fmpr_print(y); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + abort(); + } + + arf_clear(x); + arf_clear(z); + fmpr_clear(y); + } + + /* test exact roundtrip Q -> R -> Q */ + for (iter = 0; iter < 100000; iter++) + { + long bits; + fmpr_t x, z; + arf_t y; + + bits = 2 + n_randint(state, 1000); + + fmpr_init(x); + fmpr_init(z); + arf_init(y); + + fmpr_randtest_special(x, state, bits, 1 + n_randint(state, 100)); + fmpr_randtest_special(z, state, bits, 1 + n_randint(state, 100)); + + arf_set_fmpr(y, x); + arf_get_fmpr(z, y); + + if (!fmpr_equal(x, z)) + { + printf("FAIL (2)\n\n"); + printf("bits: %ld\n", bits); + printf("x = "); fmpr_print(x); printf("\n\n"); + printf("y = "); arf_print(y); printf("\n\n"); + printf("z = "); fmpr_print(z); printf("\n\n"); + abort(); + } + + fmpr_clear(x); + fmpr_clear(z); + arf_clear(y); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/arf/test/t-set_fmpz_2exp.c b/arf/test/t-set_fmpz_2exp.c new file mode 100644 index 00000000..063ec264 --- /dev/null +++ b/arf/test/t-set_fmpz_2exp.c @@ -0,0 +1,79 @@ +/*============================================================================= + + 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 "arf.h" + +int main() +{ + long iter; + flint_rand_t state; + + printf("set_fmpz_2exp...."); + fflush(stdout); + + flint_randinit(state); + + /* test exact roundtrip R -> Q -> R */ + for (iter = 0; iter < 100000; iter++) + { + long bits; + arf_t x, z; + fmpz_t y, e; + + bits = 2 + n_randint(state, 200); + + arf_init(x); + arf_init(z); + fmpz_init(y); + fmpz_init(e); + + arf_randtest(x, state, bits, 1 + n_randint(state, 100)); + arf_randtest(z, state, bits, 1 + n_randint(state, 100)); + + arf_get_fmpz_2exp(y, e, x); + arf_set_fmpz_2exp(z, y, e); + + if (!arf_equal(x, z)) + { + printf("FAIL\n\n"); + printf("bits: %ld\n", bits); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = "); fmpz_print(y); printf("\n\n"); + printf("e = "); fmpz_print(e); printf("\n\n"); + printf("z = "); arf_print(z); printf("\n\n"); + abort(); + } + + arf_clear(x); + arf_clear(z); + fmpz_clear(y); + fmpz_clear(e); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +}