contrib/mpfr/src/rint.c

   1 /* mpfr_rint -- Round to an integer.
   2
   3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
   4 Contributed by the Arenaire 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 }