Annotation of OpenXM_contrib/gmp/mpfr/rint.c, Revision 1.1
1.1 ! ohara 1: /* mpfr_rint -- Round to an integer.
! 2:
! 3: Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
! 4:
! 5: This file is part of the MPFR Library.
! 6:
! 7: The MPFR Library is free software; you can redistribute it and/or modify
! 8: it under the terms of the GNU Lesser General Public License as published by
! 9: the Free Software Foundation; either version 2.1 of the License, or (at your
! 10: option) any later version.
! 11:
! 12: The MPFR Library is distributed in the hope that it will be useful, but
! 13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! 15: License for more details.
! 16:
! 17: You should have received a copy of the GNU Lesser General Public License
! 18: along with the MPFR Library; see the file COPYING.LIB. If not, write to
! 19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 20: MA 02111-1307, USA. */
! 21:
! 22: #include "gmp.h"
! 23: #include "gmp-impl.h"
! 24: #include "mpfr.h"
! 25: #include "mpfr-impl.h"
! 26:
! 27: int
! 28: mpfr_rint (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode)
! 29: {
! 30: int sign;
! 31: mp_exp_t exp;
! 32:
! 33: if (MPFR_IS_NAN(u))
! 34: {
! 35: MPFR_SET_NAN(r);
! 36: MPFR_RET_NAN;
! 37: }
! 38:
! 39: MPFR_CLEAR_NAN(r);
! 40: MPFR_SET_SAME_SIGN(r, u);
! 41:
! 42: if (MPFR_IS_INF(u))
! 43: {
! 44: MPFR_SET_INF(r);
! 45: MPFR_RET(0); /* infinity is exact */
! 46: }
! 47:
! 48: MPFR_CLEAR_INF(r);
! 49:
! 50: if (MPFR_IS_ZERO(u))
! 51: {
! 52: MPFR_SET_ZERO(r);
! 53: MPFR_RET(0); /* zero is exact */
! 54: }
! 55:
! 56: sign = MPFR_SIGN(u);
! 57: exp = MPFR_EXP(u);
! 58:
! 59: /* Single out the case where |u| < 1. */
! 60: if (exp <= 0) /* 0 < |u| < 1 */
! 61: {
! 62: if ((rnd_mode == GMP_RNDD && sign < 0) ||
! 63: (rnd_mode == GMP_RNDU && sign > 0) ||
! 64: (rnd_mode == GMP_RNDN && exp == 0))
! 65: {
! 66: mp_limb_t *rp;
! 67: mp_size_t rm;
! 68:
! 69: rp = MPFR_MANT(r);
! 70: rm = (MPFR_PREC(r) - 1) / BITS_PER_MP_LIMB;
! 71: rp[rm] = GMP_LIMB_HIGHBIT;
! 72: MPN_ZERO(rp, rm);
! 73: MPFR_EXP(r) = 1; /* |r| = 1 */
! 74: MPFR_RET(sign > 0 ? 2 : -2);
! 75: }
! 76: else
! 77: {
! 78: MPFR_SET_ZERO(r); /* r = 0 */
! 79: MPFR_RET(sign > 0 ? -2 : 2);
! 80: }
! 81: }
! 82: else /* exp > 0, |u| >= 1 */
! 83: {
! 84: mp_limb_t *up, *rp;
! 85: mp_size_t un, rn, ui;
! 86: int sh, idiff;
! 87: int uflags;
! 88:
! 89: /*
! 90: * uflags will contain:
! 91: * _ 0 if u is an integer representable in r,
! 92: * _ 1 if u is an integer not representable in r,
! 93: * _ 2 if u is not an integer.
! 94: */
! 95:
! 96: up = MPFR_MANT(u);
! 97: rp = MPFR_MANT(r);
! 98:
! 99: un = MPFR_ESIZE(u);
! 100: rn = MPFR_ESIZE(r);
! 101: sh = rn * BITS_PER_MP_LIMB - MPFR_PREC(r);
! 102:
! 103: MPFR_EXP(r) = exp;
! 104:
! 105: if ((exp - 1) / BITS_PER_MP_LIMB >= un)
! 106: {
! 107: ui = un;
! 108: idiff = 0;
! 109: uflags = 0; /* u is an integer */
! 110: }
! 111: else
! 112: {
! 113: mp_size_t uj;
! 114:
! 115: ui = (exp - 1) / BITS_PER_MP_LIMB + 1; /* #limbs of the int part */
! 116: uj = un - ui; /* lowest limb of the integer part */
! 117: idiff = exp % BITS_PER_MP_LIMB; /* #int-part bits in up[uj] or 0 */
! 118:
! 119: uflags = idiff == 0 || (up[uj] << idiff) == 0 ? 0 : 2;
! 120: if (uflags == 0)
! 121: while (uj > 0)
! 122: if (up[--uj] != 0)
! 123: {
! 124: uflags = 2;
! 125: break;
! 126: }
! 127: }
! 128:
! 129: if (ui > rn)
! 130: {
! 131: MPFR_ASSERTN(rp != up && un > rn);
! 132: MPN_COPY(rp, up + (un - rn), rn);
! 133: /* In the rounding to the nearest mode, if the rounding bit
! 134: is 0, change the rounding mode to GMP_RNDZ. */
! 135: if (rnd_mode == GMP_RNDN &&
! 136: ((sh != 0 && (rp[0] & (MP_LIMB_T_ONE << (sh - 1))) == 0) ||
! 137: (sh == 0 && (up[un - rn - 1] & GMP_LIMB_HIGHBIT) == 0)))
! 138: rnd_mode = GMP_RNDZ;
! 139: if (uflags == 0)
! 140: { /* u is an integer; determine if it is representable */
! 141: if (sh != 0 && rp[0] << (BITS_PER_MP_LIMB - sh) != 0)
! 142: uflags = 1; /* u is not representable in r */
! 143: else
! 144: {
! 145: mp_size_t i;
! 146: for (i = un - rn - 1; i >= 0; i--)
! 147: if (up[i] != 0)
! 148: {
! 149: uflags = 1; /* u is not representable in r */
! 150: break;
! 151: }
! 152: }
! 153: }
! 154: if (sh != 0)
! 155: rp[0] &= MP_LIMB_T_MAX << sh;
! 156: }
! 157: else
! 158: {
! 159: mp_size_t uj, rj;
! 160: int ush;
! 161:
! 162: uj = un - ui; /* lowest limb of the integer part in u */
! 163: rj = rn - ui; /* lowest limb of the integer part in r */
! 164:
! 165: MPN_ZERO(rp, rj);
! 166:
! 167: if (rp != up)
! 168: MPN_COPY(rp + rj, up + uj, ui);
! 169:
! 170: rp += rj;
! 171: rn = ui;
! 172:
! 173: ush = idiff == 0 ? 0 : BITS_PER_MP_LIMB - idiff;
! 174: if (rj == 0 && ush < sh)
! 175: {
! 176: /* In the rounding to the nearest mode, if the rounding bit
! 177: is 0, change the rounding mode to GMP_RNDZ. */
! 178: if (rnd_mode == GMP_RNDN &&
! 179: (rp[0] & (MP_LIMB_T_ONE << (sh - 1))) == 0)
! 180: rnd_mode = GMP_RNDZ; /* rounding bit is 0 */
! 181: if (uflags == 0)
! 182: { /* u is an integer; determine if it is representable */
! 183: mp_limb_t mask;
! 184: mask = (MP_LIMB_T_ONE << sh) - (MP_LIMB_T_ONE << ush);
! 185: if ((rp[0] & mask) != 0)
! 186: uflags = 1; /* u is not representable in r */
! 187: }
! 188: }
! 189: else
! 190: {
! 191: sh = ush;
! 192: if (rnd_mode == GMP_RNDN &&
! 193: ((ush != 0 &&
! 194: (up[uj] & (MP_LIMB_T_ONE << (ush - 1))) == 0) ||
! 195: (ush == 0 &&
! 196: (uj == 0 || (up[uj - 1] & GMP_LIMB_HIGHBIT) == 0))))
! 197: rnd_mode = GMP_RNDZ; /* rounding bit is 0 */
! 198: }
! 199: if (sh != 0)
! 200: rp[0] &= MP_LIMB_T_MAX << sh;
! 201: }
! 202:
! 203: if (uflags == 0)
! 204: MPFR_RET(0);
! 205:
! 206: /* Note: if rnd_mode == GMP_RNDN, then round away from 0 (if
! 207: the rounding bit was 0 and rnd_mode == GMP_RNDN, rnd_mode
! 208: has been changed to GMP_RNDZ). */
! 209:
! 210: if ((rnd_mode == GMP_RNDN ||
! 211: (rnd_mode == GMP_RNDD && sign < 0) ||
! 212: (rnd_mode == GMP_RNDU && sign > 0))
! 213: && mpn_add_1(rp, rp, rn, MP_LIMB_T_ONE << sh))
! 214: {
! 215: if (exp == __mpfr_emax)
! 216: return mpfr_set_overflow(r, rnd_mode, MPFR_SIGN(r)) >= 0 ?
! 217: uflags : -uflags;
! 218: else
! 219: {
! 220: MPFR_EXP(r)++;
! 221: rp[rn-1] = GMP_LIMB_HIGHBIT;
! 222: }
! 223: }
! 224:
! 225: MPFR_RET(rnd_mode == GMP_RNDU ||
! 226: (rnd_mode == GMP_RNDZ && sign < 0) ||
! 227: (rnd_mode == GMP_RNDN && sign > 0) ? uflags : -uflags);
! 228: } /* exp > 0, |u| >= 1 */
! 229: }
! 230:
! 231: #undef mpfr_round
! 232:
! 233: int
! 234: mpfr_round (mpfr_ptr r, mpfr_srcptr u)
! 235: {
! 236: return mpfr_rint(r, u, GMP_RNDN);
! 237: }
! 238:
! 239: #undef mpfr_trunc
! 240:
! 241: int
! 242: mpfr_trunc (mpfr_ptr r, mpfr_srcptr u)
! 243: {
! 244: return mpfr_rint(r, u, GMP_RNDZ);
! 245: }
! 246:
! 247: #undef mpfr_ceil
! 248:
! 249: int
! 250: mpfr_ceil (mpfr_ptr r, mpfr_srcptr u)
! 251: {
! 252: return mpfr_rint(r, u, GMP_RNDU);
! 253: }
! 254:
! 255: #undef mpfr_floor
! 256:
! 257: int
! 258: mpfr_floor (mpfr_ptr r, mpfr_srcptr u)
! 259: {
! 260: return mpfr_rint(r, u, GMP_RNDD);
! 261: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to rint.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr