From d9a6e43aba2a6baaf38069e18a374ba403722f14 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Sat, 26 Apr 2014 23:16:23 +0200 Subject: [PATCH] add more functions --- arf.h | 172 +++++++++++++++++++++++++++++------- arf/cmp_2exp_si.c | 71 +++++++++++++++ arf/cmpabs_2exp_si.c | 70 +++++++++++++++ arf/neg_round.c | 60 +++++++++++++ arf/set_round.c | 60 +++++++++++++ arf/test/t-cmp_2exp_si.c | 78 ++++++++++++++++ arf/test/t-cmpabs_2exp_si.c | 78 ++++++++++++++++ 7 files changed, 558 insertions(+), 31 deletions(-) create mode 100644 arf/cmp_2exp_si.c create mode 100644 arf/cmpabs_2exp_si.c create mode 100644 arf/neg_round.c create mode 100644 arf/set_round.c create mode 100644 arf/test/t-cmp_2exp_si.c create mode 100644 arf/test/t-cmpabs_2exp_si.c diff --git a/arf.h b/arf.h index 4fa2f08c..6e80498a 100644 --- a/arf.h +++ b/arf.h @@ -488,38 +488,9 @@ arf_set_round_fmpz(arf_t y, const fmpz_t x, long prec, arf_rnd_t rnd) return arf_set_round_mpz(y, COEFF_TO_PTR(*x), prec, rnd); } -static __inline__ int -arf_set_round(arf_t y, const arf_t x, long prec, arf_rnd_t rnd) -{ - if (arf_is_special(x)) - { - arf_set(y, x); - return 0; - } - else - { - /* XXX: fixme */ - if (y == x) - { - int inexact; - arf_t t; - arf_init(t); - arf_set(t, x); - inexact = arf_set_round(y, t, prec, rnd); - arf_clear(t); - return inexact; - } - else - { - mp_srcptr xptr; - mp_size_t xn; +int arf_set_round(arf_t y, const arf_t x, long prec, arf_rnd_t rnd); - ARF_GET_MPN_READONLY(xptr, xn, x); - return arf_set_round_mpn(y, xptr, xn, ARF_SGNBIT(x), - ARF_EXPREF(x), prec, rnd); - } - } -} +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); @@ -529,6 +500,145 @@ void arf_set_mpfr(arf_t x, const mpfr_t y); int arf_equal(const arf_t x, const arf_t y); +static __inline__ void +arf_min(arf_t z, const arf_t a, const arf_t b) +{ + if (arf_cmp(a, b) <= 0) + arf_set(z, a); + else + arf_set(z, b); +} + +static __inline__ void +arf_max(arf_t z, const arf_t a, const arf_t b) +{ + if (arf_cmp(a, b) > 0) + arf_set(z, a); + else + arf_set(z, b); +} + +static __inline__ void +arf_abs(arf_t y, const arf_t x) +{ + if (arf_sgn(x) < 0) + arf_neg(y, x); + else + arf_set(y, x); +} + +static __inline__ long +arf_bits(const arf_t x) +{ + if (arf_is_special(x)) + return 0; + else + { + mp_srcptr xp; + mp_size_t xn; + long c; + + ARF_GET_MPN_READONLY(xp, xn, x); + count_trailing_zeros(c, xp[0]); + return xn * FLINT_BITS - c; + } +} + +static __inline__ void +arf_bot(fmpz_t e, const arf_t x) +{ + if (arf_is_special(x)) + fmpz_zero(e); + else + fmpz_sub_si(e, ARF_EXPREF(x), arf_bits(x)); +} + +static __inline__ int +arf_is_int(const arf_t x) +{ + if (arf_is_special(x)) + return arf_is_zero(x); + else + { + fmpz_t t; + int r; + fmpz_init(t); + arf_bot(t, x); + r = fmpz_sgn(t) >= 0; + fmpz_clear(t); + return r; + } +} + +static __inline__ int +arf_is_int_2exp_si(const arf_t x, long e) +{ + if (arf_is_special(x)) + return arf_is_zero(x); + else + { + fmpz_t t; + int r; + fmpz_init(t); + arf_bot(t, x); + r = fmpz_cmp_si(t, e) >= 0; + fmpz_clear(t); + return r; + } +} + +int arf_cmp_2exp_si(const arf_t x, long e); + +int arf_cmpabs_2exp_si(const arf_t x, long e); + +static __inline__ void +arf_set_si_2exp_si(arf_t x, long man, long exp) +{ + arf_set_si(x, man); + if (man != 0) + fmpz_add_si_inline(ARF_EXPREF(x), ARF_EXPREF(x), exp); +} + +static __inline__ void +arf_set_ui_2exp_si(arf_t x, ulong man, long exp) +{ + arf_set_ui(x, man); + if (man != 0) + fmpz_add_si_inline(ARF_EXPREF(x), ARF_EXPREF(x), exp); +} + +static __inline__ void +arf_mul_2exp_si(arf_t y, const arf_t x, long e) +{ + arf_set(y, x); + if (!arf_is_special(y)) + fmpz_add_si_inline(ARF_EXPREF(y), ARF_EXPREF(y), e); +} + +static __inline__ void +arf_mul_2exp_fmpz(arf_t y, const arf_t x, const fmpz_t e) +{ + arf_set(y, x); + if (!arf_is_special(y)) + fmpz_add_inline(ARF_EXPREF(y), ARF_EXPREF(y), e); +} + +static __inline__ int +arf_set_round_fmpz_2exp(arf_t y, const fmpz_t x, const fmpz_t exp, long prec, arf_rnd_t rnd) +{ + if (fmpz_is_zero(x)) + { + arf_zero(y); + return 0; + } + else + { + int r = arf_set_round_fmpz(y, x, prec, rnd); + fmpz_add_inline(ARF_EXPREF(y), ARF_EXPREF(y), exp); + return r; + } +} + void arf_debug(const arf_t x); #define arf_print arf_debug diff --git a/arf/cmp_2exp_si.c b/arf/cmp_2exp_si.c new file mode 100644 index 00000000..c648191d --- /dev/null +++ b/arf/cmp_2exp_si.c @@ -0,0 +1,71 @@ +/*============================================================================= + + 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" + +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; + if (arf_is_pos_inf(x)) return 1; + if (arf_is_neg_inf(x)) return -1; + return 0; + } + + 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)) + return 0; + else + return (ARF_EXP(x) <= e) ? -1 : 1; + } + + if (pow2) + { + fmpz_t t; + fmpz_init(t); + fmpz_one(t); + fmpz_add_si(t, t, e); + if (fmpz_equal(ARF_EXPREF(x), t)) + { + fmpz_clear(t); + return 0; + } + fmpz_clear(t); + } + + return (fmpz_cmp_si(ARF_EXPREF(x), e) <= 0) ? -1 : 1; +} + diff --git a/arf/cmpabs_2exp_si.c b/arf/cmpabs_2exp_si.c new file mode 100644 index 00000000..26198424 --- /dev/null +++ b/arf/cmpabs_2exp_si.c @@ -0,0 +1,70 @@ +/*============================================================================= + + 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" + +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; + if (arf_is_inf(x)) return 1; + 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)) + return 0; + else + return (ARF_EXP(x) <= e) ? -1 : 1; + } + + if (pow2) + { + fmpz_t t; + fmpz_init(t); + + fmpz_one(t); + fmpz_add_si(t, t, e); + + if (fmpz_equal(ARF_EXPREF(x), t)) + { + fmpz_clear(t); + return 0; + } + + fmpz_clear(t); + } + + return (fmpz_cmp_si(ARF_EXPREF(x), e) <= 0) ? -1 : 1; +} + diff --git a/arf/neg_round.c b/arf/neg_round.c new file mode 100644 index 00000000..55884a71 --- /dev/null +++ b/arf/neg_round.c @@ -0,0 +1,60 @@ +/*============================================================================= + + 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" + +int +arf_neg_round(arf_t y, const arf_t x, long prec, arf_rnd_t rnd) +{ + if (arf_is_special(x)) + { + arf_neg(y, x); + return 0; + } + else + { + /* XXX: fixme */ + if (y == x) + { + int inexact; + arf_t t; + arf_init(t); + arf_set(t, x); + inexact = arf_neg_round(y, t, prec, rnd); + arf_clear(t); + return inexact; + } + else + { + mp_srcptr xptr; + mp_size_t xn; + + ARF_GET_MPN_READONLY(xptr, xn, x); + return arf_set_round_mpn(y, xptr, xn, ARF_SGNBIT(x) ^ 1, + ARF_EXPREF(x), prec, rnd); + } + } +} + diff --git a/arf/set_round.c b/arf/set_round.c new file mode 100644 index 00000000..aa32b32c --- /dev/null +++ b/arf/set_round.c @@ -0,0 +1,60 @@ +/*============================================================================= + + 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" + +int +arf_set_round(arf_t y, const arf_t x, long prec, arf_rnd_t rnd) +{ + if (arf_is_special(x)) + { + arf_set(y, x); + return 0; + } + else + { + /* XXX: fixme */ + if (y == x) + { + int inexact; + arf_t t; + arf_init(t); + arf_set(t, x); + inexact = arf_set_round(y, t, prec, rnd); + arf_clear(t); + return inexact; + } + else + { + mp_srcptr xptr; + mp_size_t xn; + + ARF_GET_MPN_READONLY(xptr, xn, x); + return arf_set_round_mpn(y, xptr, xn, ARF_SGNBIT(x), + ARF_EXPREF(x), prec, rnd); + } + } +} + diff --git a/arf/test/t-cmp_2exp_si.c b/arf/test/t-cmp_2exp_si.c new file mode 100644 index 00000000..35fb4954 --- /dev/null +++ b/arf/test/t-cmp_2exp_si.c @@ -0,0 +1,78 @@ +/*============================================================================= + + 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("cmp_2exp_si...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000; iter++) + { + long bits, e; + arf_t x, y; + int cmp1, cmp2; + + bits = 2 + n_randint(state, 1000); + e = n_randtest(state); + + arf_init(x); + arf_init(y); + + if (iter % 10 == 0) + arf_set_ui_2exp_si(x, 1, e); + else + arf_randtest_special(x, state, bits, 100); + + arf_set_ui_2exp_si(y, 1, e); + + cmp1 = arf_cmp(x, y); + cmp2 = arf_cmp_2exp_si(x, e); + + if (cmp1 != cmp2) + { + printf("FAIL\n\n"); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = "); arf_print(y); printf("\n\n"); + printf("cmp1 = %d, cmp2 = %d\n\n", cmp1, cmp2); + abort(); + } + + arf_clear(x); + arf_clear(y); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arf/test/t-cmpabs_2exp_si.c b/arf/test/t-cmpabs_2exp_si.c new file mode 100644 index 00000000..e8a8090d --- /dev/null +++ b/arf/test/t-cmpabs_2exp_si.c @@ -0,0 +1,78 @@ +/*============================================================================= + + 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("cmpabs_2exp_si...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 100000; iter++) + { + long bits, e; + arf_t x, y; + int cmp1, cmp2; + + bits = 2 + n_randint(state, 1000); + e = n_randtest(state); + + arf_init(x); + arf_init(y); + + if (iter % 10 == 0) + arf_set_ui_2exp_si(x, 1, e); + else + arf_randtest_special(x, state, bits, 100); + + arf_set_ui_2exp_si(y, 1, e); + + cmp1 = arf_cmpabs(x, y); + cmp2 = arf_cmpabs_2exp_si(x, e); + + if (cmp1 != cmp2) + { + printf("FAIL\n\n"); + printf("x = "); arf_print(x); printf("\n\n"); + printf("y = "); arf_print(y); printf("\n\n"); + printf("cmp1 = %d, cmp2 = %d\n\n", cmp1, cmp2); + abort(); + } + + arf_clear(x); + arf_clear(y); + } + + flint_randclear(state); + flint_cleanup(); + printf("PASS\n"); + return EXIT_SUCCESS; +} +