1/* mpfr_rint -- Round to an integer. 2 3Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#include "mpfr-impl.h" 24 25/* Merge the following mpfr_rint code with mpfr_round_raw_generic? */ 26 27int 28mpfr_rint (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode) 29{ 30 int sign; 31 int rnd_away; 32 mpfr_exp_t exp; 33 34 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) )) 35 { 36 if (MPFR_IS_NAN(u)) 37 { 38 MPFR_SET_NAN(r); 39 MPFR_RET_NAN; 40 } 41 MPFR_SET_SAME_SIGN(r, u); 42 if (MPFR_IS_INF(u)) 43 { 44 MPFR_SET_INF(r); 45 MPFR_RET(0); /* infinity is exact */ 46 } 47 else /* now u is zero */ 48 { 49 MPFR_ASSERTD(MPFR_IS_ZERO(u)); 50 MPFR_SET_ZERO(r); 51 MPFR_RET(0); /* zero is exact */ 52 } 53 } 54 MPFR_SET_SAME_SIGN (r, u); /* Does nothing if r==u */ 55 56 sign = MPFR_INT_SIGN (u); 57 exp = MPFR_GET_EXP (u); 58 59 rnd_away = 60 rnd_mode == MPFR_RNDD ? sign < 0 : 61 rnd_mode == MPFR_RNDU ? sign > 0 : 62 rnd_mode == MPFR_RNDZ ? 0 : 63 rnd_mode == MPFR_RNDA ? 1 : 64 -1; /* round to nearest-even (RNDN) or nearest-away (RNDNA) */ 65 66 /* rnd_away: 67 1 if round away from zero, 68 0 if round to zero, 69 -1 if not decided yet. 70 */ 71 72 if (MPFR_UNLIKELY (exp <= 0)) /* 0 < |u| < 1 ==> round |u| to 0 or 1 */ 73 { 74 /* Note: in the MPFR_RNDN mode, 0.5 must be rounded to 0. */ 75 if (rnd_away != 0 && 76 (rnd_away > 0 || 77 (exp == 0 && (rnd_mode == MPFR_RNDNA || 78 !mpfr_powerof2_raw (u))))) 79 { 80 mp_limb_t *rp; 81 mp_size_t rm; 82 83 rp = MPFR_MANT(r); 84 rm = (MPFR_PREC(r) - 1) / GMP_NUMB_BITS; 85 rp[rm] = MPFR_LIMB_HIGHBIT; 86 MPN_ZERO(rp, rm); 87 MPFR_SET_EXP (r, 1); /* |r| = 1 */ 88 MPFR_RET(sign > 0 ? 2 : -2); 89 } 90 else 91 { 92 MPFR_SET_ZERO(r); /* r = 0 */ 93 MPFR_RET(sign > 0 ? -2 : 2); 94 } 95 } 96 else /* exp > 0, |u| >= 1 */ 97 { 98 mp_limb_t *up, *rp; 99 mp_size_t un, rn, ui; 100 int sh, idiff; 101 int uflags; 102 103 /* 104 * uflags will contain: 105 * _ 0 if u is an integer representable in r, 106 * _ 1 if u is an integer not representable in r, 107 * _ 2 if u is not an integer. 108 */ 109 110 up = MPFR_MANT(u); 111 rp = MPFR_MANT(r); 112 113 un = MPFR_LIMB_SIZE(u); 114 rn = MPFR_LIMB_SIZE(r); 115 MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (r)); 116 117 MPFR_SET_EXP (r, exp); /* Does nothing if r==u */ 118 119 if ((exp - 1) / GMP_NUMB_BITS >= un) 120 { 121 ui = un; 122 idiff = 0; 123 uflags = 0; /* u is an integer, representable or not in r */ 124 } 125 else 126 { 127 mp_size_t uj; 128 129 ui = (exp - 1) / GMP_NUMB_BITS + 1; /* #limbs of the int part */ 130 MPFR_ASSERTD (un >= ui); 131 uj = un - ui; /* lowest limb of the integer part */ 132 idiff = exp % GMP_NUMB_BITS; /* #int-part bits in up[uj] or 0 */ 133 134 uflags = idiff == 0 || (up[uj] << idiff) == 0 ? 0 : 2; 135 if (uflags == 0) 136 while (uj > 0) 137 if (up[--uj] != 0) 138 { 139 uflags = 2; 140 break; 141 } 142 } 143 144 if (ui > rn) 145 { 146 /* More limbs in the integer part of u than in r. 147 Just round u with the precision of r. */ 148 MPFR_ASSERTD (rp != up && un > rn); 149 MPN_COPY (rp, up + (un - rn), rn); /* r != u */ 150 if (rnd_away < 0) 151 { 152 /* This is a rounding to nearest mode (MPFR_RNDN or MPFR_RNDNA). 153 Decide the rounding direction here. */ 154 if (rnd_mode == MPFR_RNDN && 155 (rp[0] & (MPFR_LIMB_ONE << sh)) == 0) 156 { /* halfway cases rounded toward zero */ 157 mp_limb_t a, b; 158 /* a: rounding bit and some of the following bits */ 159 /* b: boundary for a (weight of the rounding bit in a) */ 160 if (sh != 0) 161 { 162 a = rp[0] & ((MPFR_LIMB_ONE << sh) - 1); 163 b = MPFR_LIMB_ONE << (sh - 1); 164 } 165 else 166 { 167 a = up[un - rn - 1]; 168 b = MPFR_LIMB_HIGHBIT; 169 } 170 rnd_away = a > b; 171 if (a == b) 172 { 173 mp_size_t i; 174 for (i = un - rn - 1 - (sh == 0); i >= 0; i--) 175 if (up[i] != 0) 176 { 177 rnd_away = 1; 178 break; 179 } 180 } 181 } 182 else /* halfway cases rounded away from zero */ 183 rnd_away = /* rounding bit */ 184 ((sh != 0 && (rp[0] & (MPFR_LIMB_ONE << (sh - 1))) != 0) || 185 (sh == 0 && (up[un - rn - 1] & MPFR_LIMB_HIGHBIT) != 0)); 186 } 187 if (uflags == 0) 188 { /* u is an integer; determine if it is representable in r */ 189 if (sh != 0 && rp[0] << (GMP_NUMB_BITS - sh) != 0) 190 uflags = 1; /* u is not representable in r */ 191 else 192 { 193 mp_size_t i; 194 for (i = un - rn - 1; i >= 0; i--) 195 if (up[i] != 0) 196 { 197 uflags = 1; /* u is not representable in r */ 198 break; 199 } 200 } 201 } 202 } 203 else /* ui <= rn */ 204 { 205 mp_size_t uj, rj; 206 int ush; 207 208 uj = un - ui; /* lowest limb of the integer part in u */ 209 rj = rn - ui; /* lowest limb of the integer part in r */ 210 211 if (MPFR_LIKELY (rp != up)) 212 MPN_COPY(rp + rj, up + uj, ui); 213 214 /* Ignore the lowest rj limbs, all equal to zero. */ 215 rp += rj; 216 rn = ui; 217 218 /* number of fractional bits in whole rp[0] */ 219 ush = idiff == 0 ? 0 : GMP_NUMB_BITS - idiff; 220 221 if (rj == 0 && ush < sh) 222 { 223 /* If u is an integer (uflags == 0), we need to determine 224 if it is representable in r, i.e. if its sh - ush bits 225 in the non-significant part of r are all 0. */ 226 if (uflags == 0 && (rp[0] & ((MPFR_LIMB_ONE << sh) - 227 (MPFR_LIMB_ONE << ush))) != 0) 228 uflags = 1; /* u is an integer not representable in r */ 229 } 230 else /* The integer part of u fits in r, we'll round to it. */ 231 sh = ush; 232 233 if (rnd_away < 0) 234 { 235 /* This is a rounding to nearest mode. 236 Decide the rounding direction here. */ 237 if (uj == 0 && sh == 0) 238 rnd_away = 0; /* rounding bit = 0 (not represented in u) */ 239 else if (rnd_mode == MPFR_RNDN && 240 (rp[0] & (MPFR_LIMB_ONE << sh)) == 0) 241 { /* halfway cases rounded toward zero */ 242 mp_limb_t a, b; 243 /* a: rounding bit and some of the following bits */ 244 /* b: boundary for a (weight of the rounding bit in a) */ 245 if (sh != 0) 246 { 247 a = rp[0] & ((MPFR_LIMB_ONE << sh) - 1); 248 b = MPFR_LIMB_ONE << (sh - 1); 249 } 250 else 251 { 252 MPFR_ASSERTD (uj >= 1); /* see above */ 253 a = up[uj - 1]; 254 b = MPFR_LIMB_HIGHBIT; 255 } 256 rnd_away = a > b; 257 if (a == b) 258 { 259 mp_size_t i; 260 for (i = uj - 1 - (sh == 0); i >= 0; i--) 261 if (up[i] != 0) 262 { 263 rnd_away = 1; 264 break; 265 } 266 } 267 } 268 else /* halfway cases rounded away from zero */ 269 rnd_away = /* rounding bit */ 270 ((sh != 0 && (rp[0] & (MPFR_LIMB_ONE << (sh - 1))) != 0) || 271 (sh == 0 && (MPFR_ASSERTD (uj >= 1), 272 up[uj - 1] & MPFR_LIMB_HIGHBIT) != 0)); 273 } 274 /* Now we can make the low rj limbs to 0 */ 275 MPN_ZERO (rp-rj, rj); 276 } 277 278 if (sh != 0) 279 rp[0] &= MP_LIMB_T_MAX << sh; 280 281 /* If u is a representable integer, there is no rounding. */ 282 if (uflags == 0) 283 MPFR_RET(0); 284 285 MPFR_ASSERTD (rnd_away >= 0); /* rounding direction is defined */ 286 if (rnd_away && mpn_add_1(rp, rp, rn, MPFR_LIMB_ONE << sh)) 287 { 288 if (exp == __gmpfr_emax) 289 return mpfr_overflow(r, rnd_mode, MPFR_SIGN(r)) >= 0 ? 290 uflags : -uflags; 291 else 292 { 293 MPFR_SET_EXP(r, exp + 1); 294 rp[rn-1] = MPFR_LIMB_HIGHBIT; 295 } 296 } 297 298 MPFR_RET (rnd_away ^ (sign < 0) ? uflags : -uflags); 299 } /* exp > 0, |u| >= 1 */ 300} 301 302#undef mpfr_round 303 304int 305mpfr_round (mpfr_ptr r, mpfr_srcptr u) 306{ 307 return mpfr_rint (r, u, MPFR_RNDNA); 308} 309 310#undef mpfr_trunc 311 312int 313mpfr_trunc (mpfr_ptr r, mpfr_srcptr u) 314{ 315 return mpfr_rint (r, u, MPFR_RNDZ); 316} 317 318#undef mpfr_ceil 319 320int 321mpfr_ceil (mpfr_ptr r, mpfr_srcptr u) 322{ 323 return mpfr_rint (r, u, MPFR_RNDU); 324} 325 326#undef mpfr_floor 327 328int 329mpfr_floor (mpfr_ptr r, mpfr_srcptr u) 330{ 331 return mpfr_rint (r, u, MPFR_RNDD); 332} 333 334#undef mpfr_rint_round 335 336int 337mpfr_rint_round (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode) 338{ 339 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u)) 340 return mpfr_set (r, u, rnd_mode); 341 else 342 { 343 mpfr_t tmp; 344 int inex; 345 MPFR_SAVE_EXPO_DECL (expo); 346 MPFR_BLOCK_DECL (flags); 347 348 MPFR_SAVE_EXPO_MARK (expo); 349 mpfr_init2 (tmp, MPFR_PREC (u)); 350 /* round(u) is representable in tmp unless an overflow occurs */ 351 MPFR_BLOCK (flags, mpfr_round (tmp, u)); 352 inex = (MPFR_OVERFLOW (flags) 353 ? mpfr_overflow (r, rnd_mode, MPFR_SIGN (u)) 354 : mpfr_set (r, tmp, rnd_mode)); 355 mpfr_clear (tmp); 356 MPFR_SAVE_EXPO_FREE (expo); 357 return mpfr_check_range (r, inex, rnd_mode); 358 } 359} 360 361#undef mpfr_rint_trunc 362 363int 364mpfr_rint_trunc (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode) 365{ 366 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u)) 367 return mpfr_set (r, u, rnd_mode); 368 else 369 { 370 mpfr_t tmp; 371 int inex; 372 MPFR_SAVE_EXPO_DECL (expo); 373 374 MPFR_SAVE_EXPO_MARK (expo); 375 mpfr_init2 (tmp, MPFR_PREC (u)); 376 /* trunc(u) is always representable in tmp */ 377 mpfr_trunc (tmp, u); 378 inex = mpfr_set (r, tmp, rnd_mode); 379 mpfr_clear (tmp); 380 MPFR_SAVE_EXPO_FREE (expo); 381 return mpfr_check_range (r, inex, rnd_mode); 382 } 383} 384 385#undef mpfr_rint_ceil 386 387int 388mpfr_rint_ceil (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode) 389{ 390 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u)) 391 return mpfr_set (r, u, rnd_mode); 392 else 393 { 394 mpfr_t tmp; 395 int inex; 396 MPFR_SAVE_EXPO_DECL (expo); 397 MPFR_BLOCK_DECL (flags); 398 399 MPFR_SAVE_EXPO_MARK (expo); 400 mpfr_init2 (tmp, MPFR_PREC (u)); 401 /* ceil(u) is representable in tmp unless an overflow occurs */ 402 MPFR_BLOCK (flags, mpfr_ceil (tmp, u)); 403 inex = (MPFR_OVERFLOW (flags) 404 ? mpfr_overflow (r, rnd_mode, MPFR_SIGN_POS) 405 : mpfr_set (r, tmp, rnd_mode)); 406 mpfr_clear (tmp); 407 MPFR_SAVE_EXPO_FREE (expo); 408 return mpfr_check_range (r, inex, rnd_mode); 409 } 410} 411 412#undef mpfr_rint_floor 413 414int 415mpfr_rint_floor (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode) 416{ 417 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u)) 418 return mpfr_set (r, u, rnd_mode); 419 else 420 { 421 mpfr_t tmp; 422 int inex; 423 MPFR_SAVE_EXPO_DECL (expo); 424 MPFR_BLOCK_DECL (flags); 425 426 MPFR_SAVE_EXPO_MARK (expo); 427 mpfr_init2 (tmp, MPFR_PREC (u)); 428 /* floor(u) is representable in tmp unless an overflow occurs */ 429 MPFR_BLOCK (flags, mpfr_floor (tmp, u)); 430 inex = (MPFR_OVERFLOW (flags) 431 ? mpfr_overflow (r, rnd_mode, MPFR_SIGN_NEG) 432 : mpfr_set (r, tmp, rnd_mode)); 433 mpfr_clear (tmp); 434 MPFR_SAVE_EXPO_FREE (expo); 435 return mpfr_check_range (r, inex, rnd_mode); 436 } 437} 438