support ARF_RND_NEAR in _arf_set_round_ui/uiui

This commit is contained in:
Fredrik Johansson 2016-04-27 23:20:35 +02:00
parent 8e6bee9824
commit b6a273370f
4 changed files with 204 additions and 64 deletions

View file

@ -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;

View file

@ -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)
{

84
arf/test/t-set_round_ui.c Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
*/
#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;
}

View file

@ -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)