contrib/mpfr/rint.c

   1 /* mpfr_rint -- Round to an integer.
   2
   3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Cacao 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 2.1 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.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  21 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   mp_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 == GMP_RNDD ? sign < 0 :
  61     rnd_mode == GMP_RNDU ? sign > 0 :
  62     rnd_mode == GMP_RNDZ ? 0 : -1;
  63
  64   /* rnd_away:
  65      1 if round away from zero,
  66      0 if round to zero,
  67      -1 if not decided yet.
  68    */
  69
  70   if (MPFR_UNLIKELY (exp <= 0))  /* 0 < |u| < 1 ==> round |u| to 0 or 1 */
  71     {
  72       /* Note: in the GMP_RNDN mode, 0.5 must be rounded to 0. */
  73       if (rnd_away != 0 &&
  74           (rnd_away > 0 ||
  75            (exp == 0 && (rnd_mode == GMP_RNDNA ||
  76                          !mpfr_powerof2_raw (u)))))
  77         {
  78           mp_limb_t *rp;
  79           mp_size_t rm;
  80
  81           rp = MPFR_MANT(r);
  82           rm = (MPFR_PREC(r) - 1) / BITS_PER_MP_LIMB;
  83           rp[rm] = MPFR_LIMB_HIGHBIT;
  84           MPN_ZERO(rp, rm);
  85           MPFR_SET_EXP (r, 1);  /* |r| = 1 */
  86           MPFR_RET(sign > 0 ? 2 : -2);
  87         }
  88       else
  89         {
  90           MPFR_SET_ZERO(r);  /* r = 0 */
  91           MPFR_RET(sign > 0 ? -2 : 2);
  92         }
  93     }
  94   else  /* exp > 0, |u| >= 1 */
  95     {
  96       mp_limb_t *up, *rp;
  97       mp_size_t un, rn, ui;
  98       int sh, idiff;
  99       int uflags;
 100
 101       /*
 102        * uflags will contain:
 103        *   _ 0 if u is an integer representable in r,
 104        *   _ 1 if u is an integer not representable in r,
 105        *   _ 2 if u is not an integer.
 106        */
 107
 108       up = MPFR_MANT(u);
 109       rp = MPFR_MANT(r);
 110
 111       un = MPFR_LIMB_SIZE(u);
 112       rn = MPFR_LIMB_SIZE(r);
 113       MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (r));
 114
 115       MPFR_SET_EXP (r, exp); /* Does nothing if r==u */
 116
 117       if ((exp - 1) / BITS_PER_MP_LIMB >= un)
 118         {
 119           ui = un;
 120           idiff = 0;
 121           uflags = 0;  /* u is an integer, representable or not in r */
 122         }
 123       else
 124         {
 125           mp_size_t uj;
 126
 127           ui = (exp - 1) / BITS_PER_MP_LIMB + 1;  /* #limbs of the int part */
 128           MPFR_ASSERTD (un >= ui);
 129           uj = un - ui;  /* lowest limb of the integer part */
 130           idiff = exp % BITS_PER_MP_LIMB;  /* #int-part bits in up[uj] or 0 */
 131
 132           uflags = idiff == 0 || (up[uj] << idiff) == 0 ? 0 : 2;
 133           if (uflags == 0)
 134             while (uj > 0)
 135               if (up[--uj] != 0)
 136                 {
 137                   uflags = 2;
 138                   break;
 139                 }
 140         }
 141
 142       if (ui > rn)
 143         {
 144           /* More limbs in the integer part of u than in r.
 145              Just round u with the precision of r. */
 146           MPFR_ASSERTD (rp != up && un > rn);
 147           MPN_COPY (rp, up + (un - rn), rn); /* r != u */
 148           if (rnd_away < 0)
 149             {
 150               /* This is a rounding to nearest mode (GMP_RNDN or GMP_RNDNA).
 151                  Decide the rounding direction here. */
 152               if (rnd_mode == GMP_RNDN &&
 153                   (rp[0] & (MPFR_LIMB_ONE << sh)) == 0)
 154                 { /* halfway cases rounded toward zero */
 155                   mp_limb_t a, b;
 156                   /* a: rounding bit and some of the following bits */
 157                   /* b: boundary for a (weight of the rounding bit in a) */
 158                   if (sh != 0)
 159                     {
 160                       a = rp[0] & ((MPFR_LIMB_ONE << sh) - 1);
 161                       b = MPFR_LIMB_ONE << (sh - 1);
 162                     }
 163                   else
 164                     {
 165                       a = up[un - rn - 1];
 166                       b = MPFR_LIMB_HIGHBIT;
 167                     }
 168                   rnd_away = a > b;
 169                   if (a == b)
 170                     {
 171                       mp_size_t i;
 172                       for (i = un - rn - 1 - (sh == 0); i >= 0; i--)
 173                         if (up[i] != 0)
 174                           {
 175                             rnd_away = 1;
 176                             break;
 177                           }
 178                     }
 179                 }
 180               else  /* halfway cases rounded away from zero */
 181                 rnd_away =  /* rounding bit */
 182                   ((sh != 0 && (rp[0] & (MPFR_LIMB_ONE << (sh - 1))) != 0) ||
 183                    (sh == 0 && (up[un - rn - 1] & MPFR_LIMB_HIGHBIT) != 0));
 184             }
 185           if (uflags == 0)
 186             { /* u is an integer; determine if it is representable in r */
 187               if (sh != 0 && rp[0] << (BITS_PER_MP_LIMB - sh) != 0)
 188                 uflags = 1;  /* u is not representable in r */
 189               else
 190                 {
 191                   mp_size_t i;
 192                   for (i = un - rn - 1; i >= 0; i--)
 193                     if (up[i] != 0)
 194                       {
 195                         uflags = 1;  /* u is not representable in r */
 196                         break;
 197                       }
 198                 }
 199             }
 200         }
 201       else  /* ui <= rn */
 202         {
 203           mp_size_t uj, rj;
 204           int ush;
 205
 206           uj = un - ui;  /* lowest limb of the integer part in u */
 207           rj = rn - ui;  /* lowest limb of the integer part in r */
 208
 209           if (MPFR_LIKELY (rp != up))
 210             MPN_COPY(rp + rj, up + uj, ui);
 211
 212           /* Ignore the lowest rj limbs, all equal to zero. */
 213           rp += rj;
 214           rn = ui;
 215
 216           /* number of fractional bits in whole rp[0] */
 217           ush = idiff == 0 ? 0 : BITS_PER_MP_LIMB - idiff;
 218
 219           if (rj == 0 && ush < sh)
 220             {
 221               /* If u is an integer (uflags == 0), we need to determine
 222                  if it is representable in r, i.e. if its sh - ush bits
 223                  in the non-significant part of r are all 0. */
 224               if (uflags == 0 && (rp[0] & ((MPFR_LIMB_ONE << sh) -
 225                                            (MPFR_LIMB_ONE << ush))) != 0)
 226                 uflags = 1;  /* u is an integer not representable in r */
 227             }
 228           else  /* The integer part of u fits in r, we'll round to it. */
 229             sh = ush;
 230
 231           if (rnd_away < 0)
 232             {
 233               /* This is a rounding to nearest mode.
 234                  Decide the rounding direction here. */
 235               if (uj == 0 && sh == 0)
 236                 rnd_away = 0; /* rounding bit = 0 (not represented in u) */
 237               else if (rnd_mode == GMP_RNDN &&
 238                        (rp[0] & (MPFR_LIMB_ONE << sh)) == 0)
 239                 { /* halfway cases rounded toward zero */
 240                   mp_limb_t a, b;
 241                   /* a: rounding bit and some of the following bits */
 242                   /* b: boundary for a (weight of the rounding bit in a) */
 243                   if (sh != 0)
 244                     {
 245                       a = rp[0] & ((MPFR_LIMB_ONE << sh) - 1);
 246                       b = MPFR_LIMB_ONE << (sh - 1);
 247                     }
 248                   else
 249                     {
 250                       MPFR_ASSERTD (uj >= 1);  /* see above */
 251                       a = up[uj - 1];
 252                       b = MPFR_LIMB_HIGHBIT;
 253                     }
 254                   rnd_away = a > b;
 255                   if (a == b)
 256                     {
 257                       mp_size_t i;
 258                       for (i = uj - 1 - (sh == 0); i >= 0; i--)
 259                         if (up[i] != 0)
 260                           {
 261                             rnd_away = 1;
 262                             break;
 263                           }
 264                     }
 265                 }
 266               else  /* halfway cases rounded away from zero */
 267                 rnd_away =  /* rounding bit */
 268                   ((sh != 0 && (rp[0] & (MPFR_LIMB_ONE << (sh - 1))) != 0) ||
 269                    (sh == 0 && (MPFR_ASSERTD (uj >= 1),
 270                                 up[uj - 1] & MPFR_LIMB_HIGHBIT) != 0));
 271             }
 272           /* Now we can make the low rj limbs to 0 */
 273           MPN_ZERO (rp-rj, rj);
 274         }
 275
 276       if (sh != 0)
 277         rp[0] &= MP_LIMB_T_MAX << sh;
 278
 279       /* If u is a representable integer, there is no rounding. */
 280       if (uflags == 0)
 281         MPFR_RET(0);
 282
 283       MPFR_ASSERTD (rnd_away >= 0);  /* rounding direction is defined */
 284       if (rnd_away && mpn_add_1(rp, rp, rn, MPFR_LIMB_ONE << sh))
 285         {
 286           if (exp == __gmpfr_emax)
 287             return mpfr_overflow(r, rnd_mode, MPFR_SIGN(r)) >= 0 ?
 288               uflags : -uflags;
 289           else
 290             {
 291               MPFR_SET_EXP(r, exp + 1);
 292               rp[rn-1] = MPFR_LIMB_HIGHBIT;
 293             }
 294         }
 295
 296       MPFR_RET (rnd_away ^ (sign < 0) ? uflags : -uflags);
 297     }  /* exp > 0, |u| >= 1 */
 298 }
 299
 300 #undef mpfr_round
 301
 302 int
 303 mpfr_round (mpfr_ptr r, mpfr_srcptr u)
 304 {
 305   return mpfr_rint(r, u, GMP_RNDNA);
 306 }
 307
 308 #undef mpfr_trunc
 309
 310 int
 311 mpfr_trunc (mpfr_ptr r, mpfr_srcptr u)
 312 {
 313   return mpfr_rint(r, u, GMP_RNDZ);
 314 }
 315
 316 #undef mpfr_ceil
 317
 318 int
 319 mpfr_ceil (mpfr_ptr r, mpfr_srcptr u)
 320 {
 321   return mpfr_rint(r, u, GMP_RNDU);
 322 }
 323
 324 #undef mpfr_floor
 325
 326 int
 327 mpfr_floor (mpfr_ptr r, mpfr_srcptr u)
 328 {
 329   return mpfr_rint(r, u, GMP_RNDD);
 330 }
 331
 332 #undef mpfr_rint_round
 333
 334 int
 335 mpfr_rint_round (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
 336 {
 337   if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
 338     return mpfr_set (r, u, rnd_mode);
 339   else
 340     {
 341       mpfr_t tmp;
 342       int inex;
 343       MPFR_SAVE_EXPO_DECL (expo);
 344       MPFR_BLOCK_DECL (flags);
 345
 346       MPFR_SAVE_EXPO_MARK (expo);
 347       mpfr_init2 (tmp, MPFR_PREC (u));
 348       /* round(u) is representable in tmp unless an overflow occurs */
 349       MPFR_BLOCK (flags, mpfr_round (tmp, u));
 350       inex = (MPFR_OVERFLOW (flags)
 351               ? mpfr_overflow (r, rnd_mode, MPFR_SIGN (u))
 352               : mpfr_set (r, tmp, rnd_mode));
 353       mpfr_clear (tmp);
 354       MPFR_SAVE_EXPO_FREE (expo);
 355       return mpfr_check_range (r, inex, rnd_mode);
 356     }
 357 }
 358
 359 #undef mpfr_rint_trunc
 360
 361 int
 362 mpfr_rint_trunc (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
 363 {
 364   if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
 365     return mpfr_set (r, u, rnd_mode);
 366   else
 367     {
 368       mpfr_t tmp;
 369       int inex;
 370       MPFR_SAVE_EXPO_DECL (expo);
 371
 372       MPFR_SAVE_EXPO_MARK (expo);
 373       mpfr_init2 (tmp, MPFR_PREC (u));
 374       /* trunc(u) is always representable in tmp */
 375       mpfr_trunc (tmp, u);
 376       inex = mpfr_set (r, tmp, rnd_mode);
 377       mpfr_clear (tmp);
 378       MPFR_SAVE_EXPO_FREE (expo);
 379       return mpfr_check_range (r, inex, rnd_mode);
 380     }
 381 }
 382
 383 #undef mpfr_rint_ceil
 384
 385 int
 386 mpfr_rint_ceil (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
 387 {
 388   if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
 389     return mpfr_set (r, u, rnd_mode);
 390   else
 391     {
 392       mpfr_t tmp;
 393       int inex;
 394       MPFR_SAVE_EXPO_DECL (expo);
 395       MPFR_BLOCK_DECL (flags);
 396
 397       MPFR_SAVE_EXPO_MARK (expo);
 398       mpfr_init2 (tmp, MPFR_PREC (u));
 399       /* ceil(u) is representable in tmp unless an overflow occurs */
 400       MPFR_BLOCK (flags, mpfr_ceil (tmp, u));
 401       inex = (MPFR_OVERFLOW (flags)
 402               ? mpfr_overflow (r, rnd_mode, MPFR_SIGN_POS)
 403               : mpfr_set (r, tmp, rnd_mode));
 404       mpfr_clear (tmp);
 405       MPFR_SAVE_EXPO_FREE (expo);
 406       return mpfr_check_range (r, inex, rnd_mode);
 407     }
 408 }
 409
 410 #undef mpfr_rint_floor
 411
 412 int
 413 mpfr_rint_floor (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
 414 {
 415   if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
 416     return mpfr_set (r, u, rnd_mode);
 417   else
 418     {
 419       mpfr_t tmp;
 420       int inex;
 421       MPFR_SAVE_EXPO_DECL (expo);
 422       MPFR_BLOCK_DECL (flags);
 423
 424       MPFR_SAVE_EXPO_MARK (expo);
 425       mpfr_init2 (tmp, MPFR_PREC (u));
 426       /* floor(u) is representable in tmp unless an overflow occurs */
 427       MPFR_BLOCK (flags, mpfr_floor (tmp, u));
 428       inex = (MPFR_OVERFLOW (flags)
 429               ? mpfr_overflow (r, rnd_mode, MPFR_SIGN_NEG)
 430               : mpfr_set (r, tmp, rnd_mode));
 431       mpfr_clear (tmp);
 432       MPFR_SAVE_EXPO_FREE (expo);
 433       return mpfr_check_range (r, inex, rnd_mode);
 434     }
 435 }