1 /* mpfr_rint -- Round to an integer. 2 3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramel projects, INRIA. 5 6 This file is part of the GNU MPFR Library. 7 8 The GNU MPFR Library is free software; you can redistribute it and/or modify 9 it under the terms of the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or (at your 11 option) any later version. 12 13 The GNU MPFR Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 21 51 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 27 int 28 mpfr_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 304 int 305 mpfr_round (mpfr_ptr r, mpfr_srcptr u) 306 { 307 return mpfr_rint (r, u, MPFR_RNDNA); 308 } 309 310 #undef mpfr_trunc 311 312 int 313 mpfr_trunc (mpfr_ptr r, mpfr_srcptr u) 314 { 315 return mpfr_rint (r, u, MPFR_RNDZ); 316 } 317 318 #undef mpfr_ceil 319 320 int 321 mpfr_ceil (mpfr_ptr r, mpfr_srcptr u) 322 { 323 return mpfr_rint (r, u, MPFR_RNDU); 324 } 325 326 #undef mpfr_floor 327 328 int 329 mpfr_floor (mpfr_ptr r, mpfr_srcptr u) 330 { 331 return mpfr_rint (r, u, MPFR_RNDD); 332 } 333 334 #undef mpfr_rint_round 335 336 int 337 mpfr_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 363 int 364 mpfr_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 387 int 388 mpfr_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 414 int 415 mpfr_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