--- /dev/null
+/* mpfr_rint -- Round to an integer.
+
+Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
+Contributed by the Arenaire and Cacao projects, INRIA.
+
+This file is part of the GNU MPFR Library.
+
+The GNU MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
+option) any later version.
+
+The GNU MPFR Library 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 Lesser General Public
+License for more details.
+
+You should have received a copy of the GNU Lesser General Public License
+along with the GNU MPFR Library; see the file COPYING.LIB. If not, write to
+the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
+MA 02110-1301, USA. */
+
+#include "mpfr-impl.h"
+
+/* Merge the following mpfr_rint code with mpfr_round_raw_generic? */
+
+int
+mpfr_rint (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
+{
+ int sign;
+ int rnd_away;
+ mp_exp_t exp;
+
+ if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ))
+ {
+ if (MPFR_IS_NAN(u))
+ {
+ MPFR_SET_NAN(r);
+ MPFR_RET_NAN;
+ }
+ MPFR_SET_SAME_SIGN(r, u);
+ if (MPFR_IS_INF(u))
+ {
+ MPFR_SET_INF(r);
+ MPFR_RET(0); /* infinity is exact */
+ }
+ else /* now u is zero */
+ {
+ MPFR_ASSERTD(MPFR_IS_ZERO(u));
+ MPFR_SET_ZERO(r);
+ MPFR_RET(0); /* zero is exact */
+ }
+ }
+ MPFR_SET_SAME_SIGN (r, u); /* Does nothing if r==u */
+
+ sign = MPFR_INT_SIGN (u);
+ exp = MPFR_GET_EXP (u);
+
+ rnd_away =
+ rnd_mode == GMP_RNDD ? sign < 0 :
+ rnd_mode == GMP_RNDU ? sign > 0 :
+ rnd_mode == GMP_RNDZ ? 0 : -1;
+
+ /* rnd_away:
+ 1 if round away from zero,
+ 0 if round to zero,
+ -1 if not decided yet.
+ */
+
+ if (MPFR_UNLIKELY (exp <= 0)) /* 0 < |u| < 1 ==> round |u| to 0 or 1 */
+ {
+ /* Note: in the GMP_RNDN mode, 0.5 must be rounded to 0. */
+ if (rnd_away != 0 &&
+ (rnd_away > 0 ||
+ (exp == 0 && (rnd_mode == GMP_RNDNA ||
+ !mpfr_powerof2_raw (u)))))
+ {
+ mp_limb_t *rp;
+ mp_size_t rm;
+
+ rp = MPFR_MANT(r);
+ rm = (MPFR_PREC(r) - 1) / BITS_PER_MP_LIMB;
+ rp[rm] = MPFR_LIMB_HIGHBIT;
+ MPN_ZERO(rp, rm);
+ MPFR_SET_EXP (r, 1); /* |r| = 1 */
+ MPFR_RET(sign > 0 ? 2 : -2);
+ }
+ else
+ {
+ MPFR_SET_ZERO(r); /* r = 0 */
+ MPFR_RET(sign > 0 ? -2 : 2);
+ }
+ }
+ else /* exp > 0, |u| >= 1 */
+ {
+ mp_limb_t *up, *rp;
+ mp_size_t un, rn, ui;
+ int sh, idiff;
+ int uflags;
+
+ /*
+ * uflags will contain:
+ * _ 0 if u is an integer representable in r,
+ * _ 1 if u is an integer not representable in r,
+ * _ 2 if u is not an integer.
+ */
+
+ up = MPFR_MANT(u);
+ rp = MPFR_MANT(r);
+
+ un = MPFR_LIMB_SIZE(u);
+ rn = MPFR_LIMB_SIZE(r);
+ MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (r));
+
+ MPFR_SET_EXP (r, exp); /* Does nothing if r==u */
+
+ if ((exp - 1) / BITS_PER_MP_LIMB >= un)
+ {
+ ui = un;
+ idiff = 0;
+ uflags = 0; /* u is an integer, representable or not in r */
+ }
+ else
+ {
+ mp_size_t uj;
+
+ ui = (exp - 1) / BITS_PER_MP_LIMB + 1; /* #limbs of the int part */
+ MPFR_ASSERTD (un >= ui);
+ uj = un - ui; /* lowest limb of the integer part */
+ idiff = exp % BITS_PER_MP_LIMB; /* #int-part bits in up[uj] or 0 */
+
+ uflags = idiff == 0 || (up[uj] << idiff) == 0 ? 0 : 2;
+ if (uflags == 0)
+ while (uj > 0)
+ if (up[--uj] != 0)
+ {
+ uflags = 2;
+ break;
+ }
+ }
+
+ if (ui > rn)
+ {
+ /* More limbs in the integer part of u than in r.
+ Just round u with the precision of r. */
+ MPFR_ASSERTD (rp != up && un > rn);
+ MPN_COPY (rp, up + (un - rn), rn); /* r != u */
+ if (rnd_away < 0)
+ {
+ /* This is a rounding to nearest mode (GMP_RNDN or GMP_RNDNA).
+ Decide the rounding direction here. */
+ if (rnd_mode == GMP_RNDN &&
+ (rp[0] & (MPFR_LIMB_ONE << sh)) == 0)
+ { /* halfway cases rounded towards zero */
+ mp_limb_t a, b;
+ /* a: rounding bit and some of the following bits */
+ /* b: boundary for a (weight of the rounding bit in a) */
+ if (sh != 0)
+ {
+ a = rp[0] & ((MPFR_LIMB_ONE << sh) - 1);
+ b = MPFR_LIMB_ONE << (sh - 1);
+ }
+ else
+ {
+ a = up[un - rn - 1];
+ b = MPFR_LIMB_HIGHBIT;
+ }
+ rnd_away = a > b;
+ if (a == b)
+ {
+ mp_size_t i;
+ for (i = un - rn - 1 - (sh == 0); i >= 0; i--)
+ if (up[i] != 0)
+ {
+ rnd_away = 1;
+ break;
+ }
+ }
+ }
+ else /* halfway cases rounded away from zero */
+ rnd_away = /* rounding bit */
+ ((sh != 0 && (rp[0] & (MPFR_LIMB_ONE << (sh - 1))) != 0) ||
+ (sh == 0 && (up[un - rn - 1] & MPFR_LIMB_HIGHBIT) != 0));
+ }
+ if (uflags == 0)
+ { /* u is an integer; determine if it is representable in r */
+ if (sh != 0 && rp[0] << (BITS_PER_MP_LIMB - sh) != 0)
+ uflags = 1; /* u is not representable in r */
+ else
+ {
+ mp_size_t i;
+ for (i = un - rn - 1; i >= 0; i--)
+ if (up[i] != 0)
+ {
+ uflags = 1; /* u is not representable in r */
+ break;
+ }
+ }
+ }
+ }
+ else /* ui <= rn */
+ {
+ mp_size_t uj, rj;
+ int ush;
+
+ uj = un - ui; /* lowest limb of the integer part in u */
+ rj = rn - ui; /* lowest limb of the integer part in r */
+
+ if (MPFR_LIKELY (rp != up))
+ MPN_COPY(rp + rj, up + uj, ui);
+
+ /* Ignore the lowest rj limbs, all equal to zero. */
+ rp += rj;
+ rn = ui;
+
+ /* number of fractional bits in whole rp[0] */
+ ush = idiff == 0 ? 0 : BITS_PER_MP_LIMB - idiff;
+
+ if (rj == 0 && ush < sh)
+ {
+ /* If u is an integer (uflags == 0), we need to determine
+ if it is representable in r, i.e. if its sh - ush bits
+ in the non-significant part of r are all 0. */
+ if (uflags == 0 && (rp[0] & ((MPFR_LIMB_ONE << sh) -
+ (MPFR_LIMB_ONE << ush))) != 0)
+ uflags = 1; /* u is an integer not representable in r */
+ }
+ else /* The integer part of u fits in r, we'll round to it. */
+ sh = ush;
+
+ if (rnd_away < 0)
+ {
+ /* This is a rounding to nearest mode.
+ Decide the rounding direction here. */
+ if (uj == 0 && sh == 0)
+ rnd_away = 0; /* rounding bit = 0 (not represented in u) */
+ else if (rnd_mode == GMP_RNDN &&
+ (rp[0] & (MPFR_LIMB_ONE << sh)) == 0)
+ { /* halfway cases rounded towards zero */
+ mp_limb_t a, b;
+ /* a: rounding bit and some of the following bits */
+ /* b: boundary for a (weight of the rounding bit in a) */
+ if (sh != 0)
+ {
+ a = rp[0] & ((MPFR_LIMB_ONE << sh) - 1);
+ b = MPFR_LIMB_ONE << (sh - 1);
+ }
+ else
+ {
+ MPFR_ASSERTD (uj >= 1); /* see above */
+ a = up[uj - 1];
+ b = MPFR_LIMB_HIGHBIT;
+ }
+ rnd_away = a > b;
+ if (a == b)
+ {
+ mp_size_t i;
+ for (i = uj - 1 - (sh == 0); i >= 0; i--)
+ if (up[i] != 0)
+ {
+ rnd_away = 1;
+ break;
+ }
+ }
+ }
+ else /* halfway cases rounded away from zero */
+ rnd_away = /* rounding bit */
+ ((sh != 0 && (rp[0] & (MPFR_LIMB_ONE << (sh - 1))) != 0) ||
+ (sh == 0 && (MPFR_ASSERTD (uj >= 1),
+ up[uj - 1] & MPFR_LIMB_HIGHBIT) != 0));
+ }
+ /* Now we can make the low rj limbs to 0 */
+ MPN_ZERO (rp-rj, rj);
+ }
+
+ if (sh != 0)
+ rp[0] &= MP_LIMB_T_MAX << sh;
+
+ /* If u is a representable integer, there is no rounding. */
+ if (uflags == 0)
+ MPFR_RET(0);
+
+ MPFR_ASSERTD (rnd_away >= 0); /* rounding direction is defined */
+ if (rnd_away && mpn_add_1(rp, rp, rn, MPFR_LIMB_ONE << sh))
+ {
+ if (exp == __gmpfr_emax)
+ return mpfr_overflow(r, rnd_mode, MPFR_SIGN(r)) >= 0 ?
+ uflags : -uflags;
+ else
+ {
+ MPFR_SET_EXP(r, exp + 1);
+ rp[rn-1] = MPFR_LIMB_HIGHBIT;
+ }
+ }
+
+ MPFR_RET (rnd_away ^ (sign < 0) ? uflags : -uflags);
+ } /* exp > 0, |u| >= 1 */
+}
+
+#undef mpfr_round
+
+int
+mpfr_round (mpfr_ptr r, mpfr_srcptr u)
+{
+ return mpfr_rint(r, u, GMP_RNDNA);
+}
+
+#undef mpfr_trunc
+
+int
+mpfr_trunc (mpfr_ptr r, mpfr_srcptr u)
+{
+ return mpfr_rint(r, u, GMP_RNDZ);
+}
+
+#undef mpfr_ceil
+
+int
+mpfr_ceil (mpfr_ptr r, mpfr_srcptr u)
+{
+ return mpfr_rint(r, u, GMP_RNDU);
+}
+
+#undef mpfr_floor
+
+int
+mpfr_floor (mpfr_ptr r, mpfr_srcptr u)
+{
+ return mpfr_rint(r, u, GMP_RNDD);
+}
+
+#undef mpfr_rint_round
+
+int
+mpfr_rint_round (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
+{
+ if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
+ return mpfr_set (r, u, rnd_mode);
+ else
+ {
+ mpfr_t tmp;
+ int inex;
+ MPFR_SAVE_EXPO_DECL (expo);
+ MPFR_BLOCK_DECL (flags);
+
+ MPFR_SAVE_EXPO_MARK (expo);
+ mpfr_init2 (tmp, MPFR_PREC (u));
+ /* round(u) is representable in tmp unless an overflow occurs */
+ MPFR_BLOCK (flags, mpfr_round (tmp, u));
+ inex = (MPFR_OVERFLOW (flags)
+ ? mpfr_overflow (r, rnd_mode, MPFR_SIGN (u))
+ : mpfr_set (r, tmp, rnd_mode));
+ mpfr_clear (tmp);
+ MPFR_SAVE_EXPO_FREE (expo);
+ return mpfr_check_range (r, inex, rnd_mode);
+ }
+}
+
+#undef mpfr_rint_trunc
+
+int
+mpfr_rint_trunc (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
+{
+ if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
+ return mpfr_set (r, u, rnd_mode);
+ else
+ {
+ mpfr_t tmp;
+ int inex;
+ MPFR_SAVE_EXPO_DECL (expo);
+
+ MPFR_SAVE_EXPO_MARK (expo);
+ mpfr_init2 (tmp, MPFR_PREC (u));
+ /* trunc(u) is always representable in tmp */
+ mpfr_trunc (tmp, u);
+ inex = mpfr_set (r, tmp, rnd_mode);
+ mpfr_clear (tmp);
+ MPFR_SAVE_EXPO_FREE (expo);
+ return mpfr_check_range (r, inex, rnd_mode);
+ }
+}
+
+#undef mpfr_rint_ceil
+
+int
+mpfr_rint_ceil (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
+{
+ if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
+ return mpfr_set (r, u, rnd_mode);
+ else
+ {
+ mpfr_t tmp;
+ int inex;
+ MPFR_SAVE_EXPO_DECL (expo);
+ MPFR_BLOCK_DECL (flags);
+
+ MPFR_SAVE_EXPO_MARK (expo);
+ mpfr_init2 (tmp, MPFR_PREC (u));
+ /* ceil(u) is representable in tmp unless an overflow occurs */
+ MPFR_BLOCK (flags, mpfr_ceil (tmp, u));
+ inex = (MPFR_OVERFLOW (flags)
+ ? mpfr_overflow (r, rnd_mode, MPFR_SIGN_POS)
+ : mpfr_set (r, tmp, rnd_mode));
+ mpfr_clear (tmp);
+ MPFR_SAVE_EXPO_FREE (expo);
+ return mpfr_check_range (r, inex, rnd_mode);
+ }
+}
+
+#undef mpfr_rint_floor
+
+int
+mpfr_rint_floor (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
+{
+ if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
+ return mpfr_set (r, u, rnd_mode);
+ else
+ {
+ mpfr_t tmp;
+ int inex;
+ MPFR_SAVE_EXPO_DECL (expo);
+ MPFR_BLOCK_DECL (flags);
+
+ MPFR_SAVE_EXPO_MARK (expo);
+ mpfr_init2 (tmp, MPFR_PREC (u));
+ /* floor(u) is representable in tmp unless an overflow occurs */
+ MPFR_BLOCK (flags, mpfr_floor (tmp, u));
+ inex = (MPFR_OVERFLOW (flags)
+ ? mpfr_overflow (r, rnd_mode, MPFR_SIGN_NEG)
+ : mpfr_set (r, tmp, rnd_mode));
+ mpfr_clear (tmp);
+ MPFR_SAVE_EXPO_FREE (expo);
+ return mpfr_check_range (r, inex, rnd_mode);
+ }
+}
projects / msp430-gcc.git / blobdiff