diff --git a/arf/set_round_ui.c b/arf/set_round_ui.c index ca435e6e..68d7b2f0 100644 --- a/arf/set_round_ui.c +++ b/arf/set_round_ui.c @@ -11,52 +11,65 @@ #include "arf.h" -/* Top-aligns the nonzero limb v and rounds it to prec bit. */ -#define ARF_NORMALISE_ROUND_LIMB(inexact, exp, v, \ - sgnbit, prec, rnd) \ - do { \ - count_leading_zeros(exp, v); \ - v <<= exp; \ - exp = FLINT_BITS - exp; \ - if (prec >= exp) \ - { \ - inexact = 0; \ - } \ - else \ - { \ - mp_limb_t __t = v; \ - v = MASK_LIMB(v, FLINT_BITS - prec); \ - inexact = (__t != v); \ - if (inexact && arf_rounds_up(rnd, sgnbit)) \ - { \ - v += (LIMB_ONE << (FLINT_BITS - prec)); \ - if (v == 0) \ - { \ - v = LIMB_TOP; \ - exp++; \ - } \ - } \ - } \ - } \ - while (0) +/* Top-aligns the single-limb integer value v and rounds it to prec bits. + Writes inexact, v, exp. Warning: macro without parentheses. */ +#define ARF_NORMALISE_ROUND_LIMB(inexact, exp, v, sgnbit, prec, rnd) \ + do { \ + count_leading_zeros(exp, v); \ + v <<= exp; \ + exp = FLINT_BITS - exp; \ + if (prec >= exp) \ + { \ + inexact = 0; \ + } \ + else \ + { \ + mp_limb_t hi_mask, lo_mask, rndn_mask, __t, __u; \ + hi_mask = LIMB_ONES << (FLINT_BITS - prec); \ + __t = v & hi_mask; \ + inexact = (__t != v); \ + if (inexact && rnd != ARF_RND_DOWN) \ + { \ + if (rnd == ARF_RND_NEAR) \ + { \ + lo_mask = LIMB_ONES >> prec; \ + rndn_mask = LIMB_ONE << (FLINT_BITS - prec - 1); \ + __u = v & lo_mask; \ + if (__u > rndn_mask || (__u == rndn_mask && \ + (__t << (prec - 1)))) \ + __t += (LIMB_ONE << (FLINT_BITS - prec)); \ + } \ + else if (arf_rounds_up(rnd, sgnbit)) \ + { \ + __t += (LIMB_ONE << (FLINT_BITS - prec)); \ + } \ + if (__t == 0) \ + { \ + __t = LIMB_TOP; \ + exp++; \ + } \ + } \ + v = __t; \ + } \ + } while (0) int _arf_set_round_ui(arf_t x, ulong v, int sgnbit, slong prec, arf_rnd_t rnd) { + _fmpz_demote(ARF_EXPREF(x)); + ARF_DEMOTE(x); + if (v == 0) { - arf_zero(x); + ARF_EXP(x) = ARF_EXP_ZERO; + ARF_XSIZE(x) = 0; return 0; } else { int exp, inexact; - ARF_NORMALISE_ROUND_LIMB(inexact, exp, v, sgnbit, prec, rnd); - - _fmpz_demote(ARF_EXPREF(x)); ARF_EXP(x) = exp; - ARF_DEMOTE(x); ARF_XSIZE(x) = ARF_MAKE_XSIZE(1, sgnbit); ARF_NOPTR_D(x)[0] = v; return inexact; diff --git a/arf/set_round_uiui.c b/arf/set_round_uiui.c index 0bd60f50..958e768a 100644 --- a/arf/set_round_uiui.c +++ b/arf/set_round_uiui.c @@ -11,34 +11,47 @@ #include "arf.h" -/* Top-aligns the nonzero limb v and rounds it to prec bit. */ -#define ARF_NORMALISE_ROUND_LIMB(inexact, exp, v, \ - sgnbit, prec, rnd) \ - do { \ - count_leading_zeros(exp, v); \ - v <<= exp; \ - exp = FLINT_BITS - exp; \ - if (prec >= exp) \ - { \ - inexact = 0; \ - } \ - else \ - { \ - mp_limb_t __t = v; \ - v = MASK_LIMB(v, FLINT_BITS - prec); \ - inexact = (__t != v); \ - if (inexact && arf_rounds_up(rnd, sgnbit)) \ - { \ - v += (LIMB_ONE << (FLINT_BITS - prec)); \ - if (v == 0) \ - { \ - v = LIMB_TOP; \ - exp++; \ - } \ - } \ - } \ - } \ - while (0) +/* Top-aligns the single-limb integer value v and rounds it to prec bits. + Writes inexact, v, exp. Warning: macro without parentheses. */ +#define ARF_NORMALISE_ROUND_LIMB(inexact, exp, v, sgnbit, prec, rnd) \ + do { \ + count_leading_zeros(exp, v); \ + v <<= exp; \ + exp = FLINT_BITS - exp; \ + if (prec >= exp) \ + { \ + inexact = 0; \ + } \ + else \ + { \ + mp_limb_t hi_mask, lo_mask, rndn_mask, __t, __u; \ + hi_mask = LIMB_ONES << (FLINT_BITS - prec); \ + __t = v & hi_mask; \ + inexact = (__t != v); \ + if (inexact && rnd != ARF_RND_DOWN) \ + { \ + if (rnd == ARF_RND_NEAR) \ + { \ + lo_mask = LIMB_ONES >> prec; \ + rndn_mask = LIMB_ONE << (FLINT_BITS - prec - 1); \ + __u = v & lo_mask; \ + if (__u > rndn_mask || (__u == rndn_mask && \ + (__t << (prec - 1)))) \ + __t += (LIMB_ONE << (FLINT_BITS - prec)); \ + } \ + else if (arf_rounds_up(rnd, sgnbit)) \ + { \ + __t += (LIMB_ONE << (FLINT_BITS - prec)); \ + } \ + if (__t == 0) \ + { \ + __t = LIMB_TOP; \ + exp++; \ + } \ + } \ + v = __t; \ + } \ + } while (0) int _arf_set_round_uiui(arf_t z, slong * fix, mp_limb_t hi, mp_limb_t lo, int sgnbit, slong prec, arf_rnd_t rnd) @@ -87,7 +100,36 @@ _arf_set_round_uiui(arf_t z, slong * fix, mp_limb_t hi, mp_limb_t lo, int sgnbit else { inexact = 1; - up = arf_rounds_up(rnd, sgnbit); + + if (rnd == ARF_RND_DOWN) + { + up = 0; + } + else if (rnd == ARF_RND_NEAR) + { + if (bc == prec + 1) + { + /* exactly one excess bit; check parity */ + if (trailing == FLINT_BITS - 1) + up = (hi != 0); + else + up = (lo >> (trailing + 1)) & 1; + } + else + { + /* two or more excess bits; test the first excess bit */ + mp_bitcnt_t pos = 2 * FLINT_BITS - leading - prec - 1; + + if (pos < FLINT_BITS) + up = (lo >> pos) & 1; + else + up = (hi >> (pos - FLINT_BITS)) & 1; + } + } + else + { + up = arf_rounds_up(rnd, sgnbit); + } if (prec <= FLINT_BITS) { diff --git a/arf/test/t-set_round_ui.c b/arf/test/t-set_round_ui.c new file mode 100644 index 00000000..f3d66e89 --- /dev/null +++ b/arf/test/t-set_round_ui.c @@ -0,0 +1,84 @@ +/* + Copyright (C) 2014 Fredrik Johansson + + This file is part of Arb. + + Arb is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 2.1 of the License, or + (at your option) any later version. See . +*/ + +#include "arf.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("set_round_ui...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 1000000 * arb_test_multiplier(); iter++) + { + arf_t x, y; + slong prec, fix1; + int ret1, ret2; + mp_limb_t t; + arf_rnd_t rnd; + + prec = 2 + n_randint(state, 1000); + + arf_init(x); + arf_init(y); + + arf_randtest_special(x, state, 1 + n_randint(state, 200), 1 + n_randint(state, 100)); + arf_randtest_special(y, state, 1 + n_randint(state, 200), 1 + n_randint(state, 100)); + + t = n_randtest(state); + + switch (n_randint(state, 5)) + { + case 0: rnd = ARF_RND_DOWN; break; + case 1: rnd = ARF_RND_UP; break; + case 2: rnd = ARF_RND_FLOOR; break; + case 3: rnd = ARF_RND_CEIL; break; + default: rnd = ARF_RND_NEAR; break; + } + + if (t == 0) + { + arf_zero(x); + ret1 = 0; + } + else + { + ret1 = _arf_set_round_mpn(x, &fix1, &t, 1, 0, prec, rnd); + fmpz_set_si(ARF_EXPREF(x), FLINT_BITS + fix1); + } + + ret2 = arf_set_round_ui(y, t, prec, rnd); + + if (!arf_equal(x, y) || (ret1 != ret2)) + { + flint_printf("FAIL\n\n"); + flint_printf("prec = %wd", prec); flint_printf("\n\n"); + flint_printf("t = %wu\n\n", t); + flint_printf("x = "); arf_print(x); flint_printf("\n\n"); + flint_printf("y = "); arf_print(y); flint_printf("\n\n"); + flint_printf("ret1 = %d, ret2 = %d\n\n", ret1, ret2); + abort(); + } + + arf_clear(x); + arf_clear(y); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/arf/test/t-set_round_uiui.c b/arf/test/t-set_round_uiui.c index c2a2e4cd..1fbee12c 100644 --- a/arf/test/t-set_round_uiui.c +++ b/arf/test/t-set_round_uiui.c @@ -44,12 +44,13 @@ int main() sgnbit = n_randint(state, 2); - switch (n_randint(state, 4)) + switch (n_randint(state, 5)) { case 0: rnd = ARF_RND_DOWN; break; case 1: rnd = ARF_RND_UP; break; case 2: rnd = ARF_RND_FLOOR; break; - default: rnd = ARF_RND_CEIL; break; + case 3: rnd = ARF_RND_CEIL; break; + default: rnd = ARF_RND_NEAR; break; } if (t[1] != 0)