1 /* mpfr_rint -- Round to an integer. 2 3 Copyright 1999-2018 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramba 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 /* For all the round-to-integer functions, we don't need to extend the 28 * exponent range. And it is better not to do so, so that we can test 29 * the flag setting for intermediate overflow in the test suite without 30 * involving huge non-integer numbers (thus in huge precision). This 31 * should also be faster. 32 * 33 * We also need to be careful with the flags. 34 */ 35 36 int 37 mpfr_rint (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode) 38 { 39 int sign; 40 int rnd_away; 41 mpfr_exp_t exp; 42 43 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) )) 44 { 45 if (MPFR_IS_NAN(u)) 46 { 47 MPFR_SET_NAN(r); 48 MPFR_RET_NAN; 49 } 50 MPFR_SET_SAME_SIGN(r, u); 51 if (MPFR_IS_INF(u)) 52 { 53 MPFR_SET_INF(r); 54 MPFR_RET(0); /* infinity is exact */ 55 } 56 else /* now u is zero */ 57 { 58 MPFR_ASSERTD(MPFR_IS_ZERO(u)); 59 MPFR_SET_ZERO(r); 60 MPFR_RET(0); /* zero is exact */ 61 } 62 } 63 MPFR_SET_SAME_SIGN (r, u); /* Does nothing if r==u */ 64 65 sign = MPFR_INT_SIGN (u); 66 exp = MPFR_GET_EXP (u); 67 68 rnd_away = 69 rnd_mode == MPFR_RNDD ? sign < 0 : 70 rnd_mode == MPFR_RNDU ? sign > 0 : 71 rnd_mode == MPFR_RNDZ ? 0 : 72 rnd_mode == MPFR_RNDA ? 1 : 73 -1; /* round to nearest-even (RNDN) or nearest-away (RNDNA) */ 74 75 /* rnd_away: 76 1 if round away from zero, 77 0 if round to zero, 78 -1 if not decided yet. 79 */ 80 81 if (MPFR_UNLIKELY (exp <= 0)) /* 0 < |u| < 1 ==> round |u| to 0 or 1 */ 82 { 83 /* Note: in the MPFR_RNDN mode, 0.5 must be rounded to 0. */ 84 if (rnd_away != 0 && 85 (rnd_away > 0 || 86 (exp == 0 && (rnd_mode == MPFR_RNDNA || 87 !mpfr_powerof2_raw (u))))) 88 { 89 /* The flags will correctly be set and overflow will correctly 90 be handled by mpfr_set_si. */ 91 mpfr_set_si (r, sign, rnd_mode); 92 MPFR_RET(sign > 0 ? 2 : -2); 93 } 94 else 95 { 96 MPFR_SET_ZERO(r); /* r = 0 */ 97 MPFR_RET(sign > 0 ? -2 : 2); 98 } 99 } 100 else /* exp > 0, |u| >= 1 */ 101 { 102 mp_limb_t *up, *rp; 103 mp_size_t un, rn, ui; 104 int sh, idiff; 105 int uflags; 106 107 /* 108 * uflags will contain: 109 * _ 0 if u is an integer representable in r, 110 * _ 1 if u is an integer not representable in r, 111 * _ 2 if u is not an integer. 112 */ 113 114 up = MPFR_MANT(u); 115 rp = MPFR_MANT(r); 116 117 un = MPFR_LIMB_SIZE(u); 118 rn = MPFR_LIMB_SIZE(r); 119 MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (r)); 120 121 /* exp is in the current exponent range: obtained from the input. */ 122 MPFR_SET_EXP (r, exp); /* Does nothing if r==u */ 123 124 if ((exp - 1) / GMP_NUMB_BITS >= un) 125 { 126 ui = un; 127 idiff = 0; 128 uflags = 0; /* u is an integer, representable or not in r */ 129 } 130 else 131 { 132 mp_size_t uj; 133 134 ui = (exp - 1) / GMP_NUMB_BITS + 1; /* #limbs of the int part */ 135 MPFR_ASSERTD (un >= ui); 136 uj = un - ui; /* lowest limb of the integer part */ 137 idiff = exp % GMP_NUMB_BITS; /* #int-part bits in up[uj] or 0 */ 138 139 uflags = idiff == 0 || (up[uj] << idiff) == 0 ? 0 : 2; 140 if (uflags == 0) 141 while (uj > 0) 142 if (up[--uj] != 0) 143 { 144 uflags = 2; 145 break; 146 } 147 } 148 149 if (ui > rn) 150 { 151 /* More limbs in the integer part of u than in r. 152 Just round u with the precision of r. */ 153 MPFR_ASSERTD (rp != up && un > rn); 154 MPN_COPY (rp, up + (un - rn), rn); /* r != u */ 155 if (rnd_away < 0) 156 { 157 /* This is a rounding to nearest mode (MPFR_RNDN or MPFR_RNDNA). 158 Decide the rounding direction here. */ 159 if (rnd_mode == MPFR_RNDN && 160 (rp[0] & (MPFR_LIMB_ONE << sh)) == 0) 161 { /* halfway cases rounded toward zero */ 162 mp_limb_t a, b; 163 /* a: rounding bit and some of the following bits */ 164 /* b: boundary for a (weight of the rounding bit in a) */ 165 if (sh != 0) 166 { 167 a = rp[0] & ((MPFR_LIMB_ONE << sh) - 1); 168 b = MPFR_LIMB_ONE << (sh - 1); 169 } 170 else 171 { 172 a = up[un - rn - 1]; 173 b = MPFR_LIMB_HIGHBIT; 174 } 175 rnd_away = a > b; 176 if (a == b) 177 { 178 mp_size_t i; 179 for (i = un - rn - 1 - (sh == 0); i >= 0; i--) 180 if (up[i] != 0) 181 { 182 rnd_away = 1; 183 break; 184 } 185 } 186 } 187 else /* halfway cases rounded away from zero */ 188 rnd_away = /* rounding bit */ 189 ((sh != 0 && (rp[0] & (MPFR_LIMB_ONE << (sh - 1))) != 0) || 190 (sh == 0 && (up[un - rn - 1] & MPFR_LIMB_HIGHBIT) != 0)); 191 } 192 if (uflags == 0) 193 { /* u is an integer; determine if it is representable in r */ 194 if (sh != 0 && rp[0] << (GMP_NUMB_BITS - sh) != 0) 195 uflags = 1; /* u is not representable in r */ 196 else 197 { 198 mp_size_t i; 199 for (i = un - rn - 1; i >= 0; i--) 200 if (up[i] != 0) 201 { 202 uflags = 1; /* u is not representable in r */ 203 break; 204 } 205 } 206 } 207 } 208 else /* ui <= rn */ 209 { 210 mp_size_t uj, rj; 211 int ush; 212 213 uj = un - ui; /* lowest limb of the integer part in u */ 214 rj = rn - ui; /* lowest limb of the integer part in r */ 215 216 if (rp != up) 217 MPN_COPY(rp + rj, up + uj, ui); 218 219 /* Ignore the lowest rj limbs, all equal to zero. */ 220 rp += rj; 221 rn = ui; 222 223 /* number of fractional bits in whole rp[0] */ 224 ush = idiff == 0 ? 0 : GMP_NUMB_BITS - idiff; 225 226 if (rj == 0 && ush < sh) 227 { 228 /* If u is an integer (uflags == 0), we need to determine 229 if it is representable in r, i.e. if its sh - ush bits 230 in the non-significant part of r are all 0. */ 231 if (uflags == 0 && (rp[0] & ((MPFR_LIMB_ONE << sh) - 232 (MPFR_LIMB_ONE << ush))) != 0) 233 uflags = 1; /* u is an integer not representable in r */ 234 } 235 else /* The integer part of u fits in r, we'll round to it. */ 236 sh = ush; 237 238 if (rnd_away < 0) 239 { 240 /* This is a rounding to nearest mode. 241 Decide the rounding direction here. */ 242 if (uj == 0 && sh == 0) 243 rnd_away = 0; /* rounding bit = 0 (not represented in u) */ 244 else if (rnd_mode == MPFR_RNDN && 245 (rp[0] & (MPFR_LIMB_ONE << sh)) == 0) 246 { /* halfway cases rounded toward zero */ 247 mp_limb_t a, b; 248 /* a: rounding bit and some of the following bits */ 249 /* b: boundary for a (weight of the rounding bit in a) */ 250 if (sh != 0) 251 { 252 a = rp[0] & ((MPFR_LIMB_ONE << sh) - 1); 253 b = MPFR_LIMB_ONE << (sh - 1); 254 } 255 else 256 { 257 MPFR_ASSERTD (uj >= 1); /* see above */ 258 a = up[uj - 1]; 259 b = MPFR_LIMB_HIGHBIT; 260 } 261 rnd_away = a > b; 262 if (a == b) 263 { 264 mp_size_t i; 265 for (i = uj - 1 - (sh == 0); i >= 0; i--) 266 if (up[i] != 0) 267 { 268 rnd_away = 1; 269 break; 270 } 271 } 272 } 273 else /* halfway cases rounded away from zero */ 274 rnd_away = /* rounding bit */ 275 ((sh != 0 && (rp[0] & (MPFR_LIMB_ONE << (sh - 1))) != 0) || 276 (sh == 0 && (MPFR_ASSERTD (uj >= 1), 277 up[uj - 1] & MPFR_LIMB_HIGHBIT) != 0)); 278 } 279 /* Now we can make the low rj limbs to 0 */ 280 MPN_ZERO (rp-rj, rj); 281 } 282 283 if (sh != 0) 284 rp[0] &= MPFR_LIMB_MAX << sh; 285 286 /* If u is a representable integer, there is no rounding. */ 287 if (uflags == 0) 288 MPFR_RET(0); 289 290 MPFR_ASSERTD (rnd_away >= 0); /* rounding direction is defined */ 291 if (rnd_away && mpn_add_1(rp, rp, rn, MPFR_LIMB_ONE << sh)) 292 { 293 if (exp == __gmpfr_emax) 294 return mpfr_overflow (r, rnd_mode, sign) >= 0 ? 295 uflags : -uflags; 296 else /* no overflow */ 297 { 298 MPFR_SET_EXP(r, exp + 1); 299 rp[rn-1] = MPFR_LIMB_HIGHBIT; 300 } 301 } 302 303 MPFR_RET (rnd_away ^ (sign < 0) ? uflags : -uflags); 304 } /* exp > 0, |u| >= 1 */ 305 } 306 307 #undef mpfr_roundeven 308 309 int 310 mpfr_roundeven (mpfr_ptr r, mpfr_srcptr u) 311 { 312 return mpfr_rint (r, u, MPFR_RNDN); 313 } 314 315 #undef mpfr_round 316 317 int 318 mpfr_round (mpfr_ptr r, mpfr_srcptr u) 319 { 320 return mpfr_rint (r, u, MPFR_RNDNA); 321 } 322 323 #undef mpfr_trunc 324 325 int 326 mpfr_trunc (mpfr_ptr r, mpfr_srcptr u) 327 { 328 return mpfr_rint (r, u, MPFR_RNDZ); 329 } 330 331 #undef mpfr_ceil 332 333 int 334 mpfr_ceil (mpfr_ptr r, mpfr_srcptr u) 335 { 336 return mpfr_rint (r, u, MPFR_RNDU); 337 } 338 339 #undef mpfr_floor 340 341 int 342 mpfr_floor (mpfr_ptr r, mpfr_srcptr u) 343 { 344 return mpfr_rint (r, u, MPFR_RNDD); 345 } 346 347 /* We need to save the flags and restore them after calling the mpfr_round, 348 * mpfr_trunc, mpfr_ceil, mpfr_floor functions because these functions set 349 * the inexact flag when the argument is not an integer. 350 */ 351 352 #undef mpfr_rint_roundeven 353 354 int 355 mpfr_rint_roundeven (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode) 356 { 357 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u)) 358 return mpfr_set (r, u, rnd_mode); 359 else 360 { 361 mpfr_t tmp; 362 int inex; 363 mpfr_flags_t saved_flags = __gmpfr_flags; 364 MPFR_BLOCK_DECL (flags); 365 366 mpfr_init2 (tmp, MPFR_PREC (u)); 367 /* round(u) is representable in tmp unless an overflow occurs */ 368 MPFR_BLOCK (flags, mpfr_roundeven (tmp, u)); 369 __gmpfr_flags = saved_flags; 370 inex = (MPFR_OVERFLOW (flags) 371 ? mpfr_overflow (r, rnd_mode, MPFR_SIGN (u)) 372 : mpfr_set (r, tmp, rnd_mode)); 373 mpfr_clear (tmp); 374 return inex; 375 } 376 } 377 378 #undef mpfr_rint_round 379 380 int 381 mpfr_rint_round (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode) 382 { 383 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u)) 384 return mpfr_set (r, u, rnd_mode); 385 else 386 { 387 mpfr_t tmp; 388 int inex; 389 mpfr_flags_t saved_flags = __gmpfr_flags; 390 MPFR_BLOCK_DECL (flags); 391 392 mpfr_init2 (tmp, MPFR_PREC (u)); 393 /* round(u) is representable in tmp unless an overflow occurs */ 394 MPFR_BLOCK (flags, mpfr_round (tmp, u)); 395 __gmpfr_flags = saved_flags; 396 inex = (MPFR_OVERFLOW (flags) 397 ? mpfr_overflow (r, rnd_mode, MPFR_SIGN (u)) 398 : mpfr_set (r, tmp, rnd_mode)); 399 mpfr_clear (tmp); 400 return inex; 401 } 402 } 403 404 #undef mpfr_rint_trunc 405 406 int 407 mpfr_rint_trunc (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode) 408 { 409 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u)) 410 return mpfr_set (r, u, rnd_mode); 411 else 412 { 413 mpfr_t tmp; 414 int inex; 415 mpfr_flags_t saved_flags = __gmpfr_flags; 416 417 mpfr_init2 (tmp, MPFR_PREC (u)); 418 /* trunc(u) is always representable in tmp */ 419 mpfr_trunc (tmp, u); 420 __gmpfr_flags = saved_flags; 421 inex = mpfr_set (r, tmp, rnd_mode); 422 mpfr_clear (tmp); 423 return inex; 424 } 425 } 426 427 #undef mpfr_rint_ceil 428 429 int 430 mpfr_rint_ceil (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode) 431 { 432 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u)) 433 return mpfr_set (r, u, rnd_mode); 434 else 435 { 436 mpfr_t tmp; 437 int inex; 438 mpfr_flags_t saved_flags = __gmpfr_flags; 439 MPFR_BLOCK_DECL (flags); 440 441 mpfr_init2 (tmp, MPFR_PREC (u)); 442 /* ceil(u) is representable in tmp unless an overflow occurs */ 443 MPFR_BLOCK (flags, mpfr_ceil (tmp, u)); 444 __gmpfr_flags = saved_flags; 445 inex = (MPFR_OVERFLOW (flags) 446 ? mpfr_overflow (r, rnd_mode, MPFR_SIGN_POS) 447 : mpfr_set (r, tmp, rnd_mode)); 448 mpfr_clear (tmp); 449 return inex; 450 } 451 } 452 453 #undef mpfr_rint_floor 454 455 int 456 mpfr_rint_floor (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode) 457 { 458 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u)) 459 return mpfr_set (r, u, rnd_mode); 460 else 461 { 462 mpfr_t tmp; 463 int inex; 464 mpfr_flags_t saved_flags = __gmpfr_flags; 465 MPFR_BLOCK_DECL (flags); 466 467 mpfr_init2 (tmp, MPFR_PREC (u)); 468 /* floor(u) is representable in tmp unless an overflow occurs */ 469 MPFR_BLOCK (flags, mpfr_floor (tmp, u)); 470 __gmpfr_flags = saved_flags; 471 inex = (MPFR_OVERFLOW (flags) 472 ? mpfr_overflow (r, rnd_mode, MPFR_SIGN_NEG) 473 : mpfr_set (r, tmp, rnd_mode)); 474 mpfr_clear (tmp); 475 return inex; 476 } 477 } 478