1/* mpfr_rint -- Round to an integer. 2 3Copyright 1999-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba 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 20https://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 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 36int 37mpfr_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 || MPFR_LIMB_LSHIFT(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 && MPFR_LIMB_LSHIFT(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 309int 310mpfr_roundeven (mpfr_ptr r, mpfr_srcptr u) 311{ 312 return mpfr_rint (r, u, MPFR_RNDN); 313} 314 315#undef mpfr_round 316 317int 318mpfr_round (mpfr_ptr r, mpfr_srcptr u) 319{ 320 return mpfr_rint (r, u, MPFR_RNDNA); 321} 322 323#undef mpfr_trunc 324 325int 326mpfr_trunc (mpfr_ptr r, mpfr_srcptr u) 327{ 328 return mpfr_rint (r, u, MPFR_RNDZ); 329} 330 331#undef mpfr_ceil 332 333int 334mpfr_ceil (mpfr_ptr r, mpfr_srcptr u) 335{ 336 return mpfr_rint (r, u, MPFR_RNDU); 337} 338 339#undef mpfr_floor 340 341int 342mpfr_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 354int 355mpfr_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 380int 381mpfr_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 406int 407mpfr_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 429int 430mpfr_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 455int 456mpfr_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