path: root/rint.c

/* mpfr_rint -- Round to an integer.

Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation, Inc.

This file is part of the MPFR Library.

The MPFR Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.

The MPFR Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the MPFR Library; see the file COPYING.LIB.  If not, write to
the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
MA 02111-1307, USA. */

#include "gmp.h"
#include "gmp-impl.h"
#include "mpfr.h"
#include "mpfr-impl.h"

/* This may be a standard rounding mode in the future. */
#define NEAREST_AWAY ((mp_rnd_t) -1)

/* Merge the following mpfr_rint code with mpfr_round_raw_generic? */

int
mpfr_rint (mpfr_ptr r, mpfr_srcptr u, mp_rnd_t rnd_mode)
{
  int sign;
  int rnd_away;
  mp_exp_t exp;

  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ))
    {
      if (MPFR_IS_NAN(u))
	{
	  MPFR_SET_NAN(r);
	  MPFR_RET_NAN;
	}
      MPFR_SET_SAME_SIGN(r, u);
      if (MPFR_IS_INF(u))
	{
	  MPFR_SET_INF(r);
	  MPFR_RET(0);  /* infinity is exact */
	}
      if (MPFR_IS_ZERO(u))
	{
	  MPFR_SET_ZERO(r);
	  MPFR_RET(0);  /* zero is exact */
	}
      MPFR_ASSERTN(0);
    }
  MPFR_SET_SAME_SIGN(r, u);
 
  sign = MPFR_INT_SIGN(u);
  exp = MPFR_GET_EXP (u);

  rnd_away =
    rnd_mode == GMP_RNDD ? sign < 0 :
    rnd_mode == GMP_RNDU ? sign > 0 :
    rnd_mode == GMP_RNDZ ? 0 : -1;

  /* rnd_away:
     1 if round away from zero,
     0 if round to zero,
     -1 if not decided yet.
   */

  if (exp <= 0)  /* 0 < |u| < 1 ==> round |u| to 0 or 1 */
    {
      /* Note: in the GMP_RNDN mode, 0.5 must be rounded to 0. */
      if (rnd_away != 0 &&
          (rnd_away > 0 ||
           (exp == 0 && (rnd_mode == NEAREST_AWAY ||
                         !mpfr_powerof2_raw (u)))))
        {
          mp_limb_t *rp;
          mp_size_t rm;

          rp = MPFR_MANT(r);
          rm = (MPFR_PREC(r) - 1) / BITS_PER_MP_LIMB;
          rp[rm] = MPFR_LIMB_HIGHBIT;
          MPN_ZERO(rp, rm);
          MPFR_SET_EXP (r, 1);  /* |r| = 1 */
          MPFR_RET(sign > 0 ? 2 : -2);
        }
      else
        {
          MPFR_SET_ZERO(r);  /* r = 0 */
          MPFR_RET(sign > 0 ? -2 : 2);
        }
    }
  else  /* exp > 0, |u| >= 1 */
    {
      mp_limb_t *up, *rp;
      mp_size_t un, rn, ui;
      int sh, idiff;
      int uflags;

      /*
       * uflags will contain:
       *   _ 0 if u is an integer representable in r,
       *   _ 1 if u is an integer not representable in r,
       *   _ 2 if u is not an integer.
       */

      up = MPFR_MANT(u);
      rp = MPFR_MANT(r);

      un = MPFR_LIMB_SIZE(u);
      rn = MPFR_LIMB_SIZE(r);
      sh = (mp_prec_t) rn * BITS_PER_MP_LIMB - MPFR_PREC(r);

      MPFR_SET_EXP (r, exp);

      if ((exp - 1) / BITS_PER_MP_LIMB >= un)
        {
          ui = un;
          idiff = 0;
          uflags = 0;  /* u is an integer, representable or not in r */
        }
      else
        {
          mp_size_t uj;

          ui = (exp - 1) / BITS_PER_MP_LIMB + 1;  /* #limbs of the int part */
          MPFR_ASSERTN (un >= ui);
          uj = un - ui;  /* lowest limb of the integer part */
          idiff = exp % BITS_PER_MP_LIMB;  /* #int-part bits in up[uj] or 0 */

          uflags = idiff == 0 || (up[uj] << idiff) == 0 ? 0 : 2;
          if (uflags == 0)
            while (uj > 0)
              if (up[--uj] != 0)
                {
                  uflags = 2;
                  break;
                }
        }

      if (ui > rn)
        {
          /* More limbs in the integer part of u than in r.
             Just round u with the precision of r. */
          MPFR_ASSERTN(rp != up && un > rn);
          MPN_COPY(rp, up + (un - rn), rn);
          if (rnd_away < 0)
            {
              /* This is a rounding to nearest mode.
                 Decide the rounding direction here. */
              if (rnd_mode == GMP_RNDN &&
                  (rp[0] & (MP_LIMB_T_ONE << sh)) == 0)
                { /* halfway cases rounded towards zero */
                  mp_limb_t a, b;
                  /* a: rounding bit and some of the following bits */
                  /* b: boundary for a (weight of the rounding bit in a) */
                  if (sh != 0)
                    {
                      a = rp[0] & ((MP_LIMB_T_ONE << sh) - 1);
                      b = MP_LIMB_T_ONE << (sh - 1);
                    }
                  else
                    {
                      a = up[un - rn - 1];
                      b = MPFR_LIMB_HIGHBIT;
                    }
                  rnd_away = a > b;
                  if (a == b)
                    {
                      mp_size_t i;
                      for (i = un - rn - 1 - (sh == 0); i >= 0; i--)
                        if (up[i] != 0)
                          {
                            rnd_away = 1;
                            break;
                          }
                    }
                }
              else  /* halfway cases rounded away from zero */
                rnd_away =  /* rounding bit */
                  ((sh != 0 && (rp[0] & (MP_LIMB_T_ONE << (sh - 1))) != 0) ||
                   (sh == 0 && (up[un - rn - 1] & MPFR_LIMB_HIGHBIT) != 0));
            }
          if (uflags == 0)
            { /* u is an integer; determine if it is representable in r */
              if (sh != 0 && rp[0] << (BITS_PER_MP_LIMB - sh) != 0)
                uflags = 1;  /* u is not representable in r */
              else
                {
                  mp_size_t i;
                  for (i = un - rn - 1; i >= 0; i--)
                    if (up[i] != 0)
                      {
                        uflags = 1;  /* u is not representable in r */
                        break;
                      }
                }
            }
        }
      else  /* ui <= rn */
        {
          mp_size_t uj, rj;
          int ush;

          uj = un - ui;  /* lowest limb of the integer part in u */
          rj = rn - ui;  /* lowest limb of the integer part in r */

          MPN_ZERO(rp, rj);

          if (rp != up)
            MPN_COPY(rp + rj, up + uj, ui);

          /* Ignore the lowest rj limbs, all equal to zero. */
          rp += rj;
          rn = ui;

          /* number of fractional bits in whole rp[0] */
          ush = idiff == 0 ? 0 : BITS_PER_MP_LIMB - idiff;

          if (rj == 0 && ush < sh)
            {
              /* If u is an integer (uflags == 0), we need to determine
                 if it is representable in r, i.e. if its sh - ush bits
                 in the non-significant part of r are all 0. */
              if (uflags == 0 && (rp[0] & ((MP_LIMB_T_ONE << sh) -
                                           (MP_LIMB_T_ONE << ush))) != 0)
                uflags = 1;  /* u is an integer not representable in r */
            }
          else  /* The integer part of u fits in r, we'll round to it. */
            sh = ush;

          if (rnd_away < 0)
            {
              /* This is a rounding to nearest mode.
                 Decide the rounding direction here. */
              if (uj == 0 && sh == 0)
                rnd_away = 0; /* rounding bit = 0 (not represented in u) */
              else if (rnd_mode == GMP_RNDN &&
                       (rp[0] & (MP_LIMB_T_ONE << sh)) == 0)
                { /* halfway cases rounded towards zero */
                  mp_limb_t a, b;
                  /* a: rounding bit and some of the following bits */
                  /* b: boundary for a (weight of the rounding bit in a) */
                  if (sh != 0)
                    {
                      a = rp[0] & ((MP_LIMB_T_ONE << sh) - 1);
                      b = MP_LIMB_T_ONE << (sh - 1);
                    }
                  else
                    {
                      MPFR_ASSERTN (uj >= 1);  /* see above */
                      a = up[uj - 1];
                      b = MPFR_LIMB_HIGHBIT;
                    }
                  rnd_away = a > b;
                  if (a == b)
                    {
                      mp_size_t i;
                      for (i = uj - 1 - (sh == 0); i >= 0; i--)
                        if (up[i] != 0)
                          {
                            rnd_away = 1;
                            break;
                          }
                    }
                }
              else  /* halfway cases rounded away from zero */
                rnd_away =  /* rounding bit */
                  ((sh != 0 && (rp[0] & (MP_LIMB_T_ONE << (sh - 1))) != 0) ||
                   (sh == 0 && (MPFR_ASSERTN (uj >= 1),
                                up[uj - 1] & MPFR_LIMB_HIGHBIT) != 0));
            }
        }

      if (sh != 0)
        rp[0] &= MP_LIMB_T_MAX << sh;

      /* If u is a representable integer, there is no rounding. */
      if (uflags == 0)
        MPFR_RET(0);

      MPFR_ASSERTN (rnd_away >= 0);  /* rounding direction is defined */
      if (rnd_away && mpn_add_1(rp, rp, rn, MP_LIMB_T_ONE << sh))
        {
          if (exp == __gmpfr_emax)
            return mpfr_set_overflow(r, rnd_mode, MPFR_SIGN(r)) >= 0 ?
              uflags : -uflags;
          else
            {
              MPFR_SET_EXP(r, exp + 1);
              rp[rn-1] = MPFR_LIMB_HIGHBIT;
            }
        }

      MPFR_RET (rnd_away ^ (sign < 0) ? uflags : -uflags);
    }  /* exp > 0, |u| >= 1 */
}

#undef mpfr_round

int
mpfr_round (mpfr_ptr r, mpfr_srcptr u)
{
  return mpfr_rint(r, u, NEAREST_AWAY);
}

#undef mpfr_trunc

int
mpfr_trunc (mpfr_ptr r, mpfr_srcptr u)
{
  return mpfr_rint(r, u, GMP_RNDZ);
}

#undef mpfr_ceil

int
mpfr_ceil (mpfr_ptr r, mpfr_srcptr u)
{
  return mpfr_rint(r, u, GMP_RNDU);
}

#undef mpfr_floor

int
mpfr_floor (mpfr_ptr r, mpfr_srcptr u)
{
  return mpfr_rint(r, u, GMP_RNDD);
}