diff options
author | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2002-01-30 04:57:10 +0000 |
---|---|---|
committer | vlefevre <vlefevre@280ebfd0-de03-0410-8827-d642c229c3f4> | 2002-01-30 04:57:10 +0000 |
commit | dd9723f715d6c1a18616404cfcf302b082b42672 (patch) | |
tree | ba593e7c0b3ef78f35c2c7d370c616cf4e41e809 /rint.c | |
parent | 7bfdb14a6f5e2390e5028e9ba75175f8dbc00a5b (diff) | |
download | mpfr-dd9723f715d6c1a18616404cfcf302b082b42672.tar.gz |
mpfr_rint, mpfr_trunc, mpfr_floor, mpfr_ceil, mpfr_round rewritten.
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@1672 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'rint.c')
-rw-r--r-- | rint.c | 261 |
1 files changed, 261 insertions, 0 deletions
@@ -0,0 +1,261 @@ +/* mpfr_rint -- Round to an integer. + +Copyright (C) 1999-2002 Free Software Foundation, Inc. + +This file is part of the MPFR Library. + +The 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 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 MPFR Library; see the file COPYING.LIB. If not, write to +the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, +MA 02111-1307, USA. */ + +#include "gmp.h" +#include "gmp-impl.h" +#include "mpfr.h" +#include "mpfr-impl.h" + +int +mpfr_rint (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode) +{ + int sign; + mp_exp_t exp; + + if (MPFR_IS_NAN(u)) + { + MPFR_SET_NAN(r); + MPFR_RET_NAN; + } + + MPFR_CLEAR_NAN(r); + MPFR_SET_SAME_SIGN(r, u); + + if (MPFR_IS_INF(u)) + { + MPFR_SET_INF(r); + MPFR_RET(0); /* infinity is exact */ + } + + MPFR_CLEAR_INF(r); + + if (MPFR_IS_ZERO(u)) + { + MPFR_SET_ZERO(r); + MPFR_RET(0); /* zero is exact */ + } + + sign = MPFR_SIGN(u); + exp = MPFR_EXP(u); + + /* Single out the case where |u| < 1. */ + if (exp <= 0) /* 0 < |u| < 1 */ + { + if ((rnd_mode == GMP_RNDD && sign < 0) || + (rnd_mode == GMP_RNDU && sign > 0) || + (rnd_mode == GMP_RNDN && exp == 0)) + { + mp_limb_t *rp; + mp_size_t rm; + + rp = MPFR_MANT(r); + rm = (MPFR_PREC(r) - 1) / BITS_PER_MP_LIMB; + rp[rm] = MP_LIMB_T_HIGHBIT; + MPN_ZERO(rp, rm); + MPFR_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_ESIZE(u); + rn = MPFR_ESIZE(r); + sh = rn * BITS_PER_MP_LIMB - MPFR_PREC(r); + + MPFR_EXP(r) = exp; + + if ((exp - 1) / BITS_PER_MP_LIMB >= un) + { + ui = un; + idiff = 0; + uflags = 0; /* u is an integer */ + } + else + { + mp_size_t uj; + + ui = (exp - 1) / BITS_PER_MP_LIMB + 1; /* #limbs of the int part */ + 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) + { + MPFR_ASSERTN(rp != up && un > rn); + MPN_COPY(rp, up + (un - rn), rn); + /* In the rounding to the nearest mode, if the rounding bit + is 0, change the rounding mode to GMP_RNDZ. */ + if (rnd_mode == GMP_RNDN && + ((sh != 0 && (rp[0] & (MP_LIMB_T_ONE << (sh - 1))) == 0) || + (sh == 0 && (up[un - rn - 1] & MP_LIMB_T_HIGHBIT) == 0))) + rnd_mode = GMP_RNDZ; + if (uflags == 0) + { /* u is an integer; determine if it is representable */ + 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; + } + } + } + if (sh != 0) + rp[0] &= MP_LIMB_T_MAX << sh; + } + else + { + 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 */ + + MPN_ZERO(rp, rj); + + if (rp != up) + MPN_COPY(rp + rj, up + uj, ui); + + rp += rj; + rn = ui; + + ush = idiff == 0 ? 0 : BITS_PER_MP_LIMB - idiff; + if (rj == 0 && ush < sh) + { + /* In the rounding to the nearest mode, if the rounding bit + is 0, change the rounding mode to GMP_RNDZ. */ + if (rnd_mode == GMP_RNDN && + (rp[0] & (MP_LIMB_T_ONE << (sh - 1))) == 0) + rnd_mode = GMP_RNDZ; /* rounding bit is 0 */ + if (uflags == 0) + { /* u is an integer; determine if it is representable */ + mp_limb_t mask; + mask = (MP_LIMB_T_ONE << sh) - (MP_LIMB_T_ONE << ush); + if ((rp[0] & mask) != 0) + uflags = 1; /* u is not representable in r */ + } + } + else + { + sh = ush; + if (rnd_mode == GMP_RNDN && + ((ush != 0 && + (up[uj] & (MP_LIMB_T_ONE << (ush - 1))) == 0) || + (ush == 0 && + (uj == 0 || (up[uj - 1] & MP_LIMB_T_HIGHBIT) == 0)))) + rnd_mode = GMP_RNDZ; /* rounding bit is 0 */ + } + if (sh != 0) + rp[0] &= MP_LIMB_T_MAX << sh; + } + + if (uflags == 0) + MPFR_RET(0); + + /* Note: if rnd_mode == GMP_RNDN, then round away from 0 (if + the rounding bit was 0 and rnd_mode == GMP_RNDN, rnd_mode + has been changed to GMP_RNDZ). */ + + if ((rnd_mode == GMP_RNDN || + (rnd_mode == GMP_RNDD && sign < 0) || + (rnd_mode == GMP_RNDU && sign > 0)) + && mpn_add_1(rp, rp, rn, MP_LIMB_T_ONE << sh)) + { + if (exp == __mpfr_emax) + return mpfr_set_overflow(r, rnd_mode, MPFR_SIGN(r)) >= 0 ? + uflags : -uflags; + else + { + MPFR_EXP(r)++; + rp[rn-1] = MP_LIMB_T_HIGHBIT; + } + } + + MPFR_RET(rnd_mode == GMP_RNDU || + (rnd_mode == GMP_RNDZ && sign < 0) || + (rnd_mode == GMP_RNDN && 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_RNDN); +} + +#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); +} |