path: root/src/round_prec.c

/* mpfr_round_raw_generic, mpfr_round_raw2, mpfr_round_raw, mpfr_prec_round,
   mpfr_can_round, mpfr_can_round_raw -- various rounding functions

Copyright 1999-2018 Free Software Foundation, Inc.
Contributed by the AriC and Caramba projects, INRIA.

This file is part of the GNU MPFR Library.

The GNU 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 3 of the License, or (at your
option) any later version.

The GNU 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 GNU MPFR Library; see the file COPYING.LESSER.  If not, see
http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */

#include "mpfr-impl.h"

#define mpfr_round_raw_generic mpfr_round_raw
#define flag 0
#define use_inexp 1
#include "round_raw_generic.c"

/* mpfr_round_raw_2 is called from mpfr_round_raw2 */
#define mpfr_round_raw_generic mpfr_round_raw_2
#define flag 1
#define use_inexp 0
#include "round_raw_generic.c"

/* Seems to be unused. Remove comment to implement it.
#define mpfr_round_raw_generic mpfr_round_raw_3
#define flag 1
#define use_inexp 1
#include "round_raw_generic.c"
*/

#define mpfr_round_raw_generic mpfr_round_raw_4
#define flag 0
#define use_inexp 0
#include "round_raw_generic.c"

/* Note: if the new prec is lower than the current one, a reallocation
   must not be done (see exp_2.c). */

int
mpfr_prec_round (mpfr_ptr x, mpfr_prec_t prec, mpfr_rnd_t rnd_mode)
{
  mp_limb_t *tmp, *xp;
  int carry, inexact;
  mpfr_prec_t nw, ow;
  MPFR_TMP_DECL(marker);

  MPFR_ASSERTN (MPFR_PREC_COND (prec));

  nw = MPFR_PREC2LIMBS (prec); /* needed allocated limbs */

  /* check if x has enough allocated space for the significand */
  /* Get the number of limbs from the precision.
     (Compatible with all allocation methods) */
  ow = MPFR_LIMB_SIZE (x);
  if (MPFR_UNLIKELY (nw > ow))
    {
      /* FIXME: Variable can't be created using custom allocation,
         MPFR_DECL_INIT or GROUP_ALLOC: How to detect? */
      ow = MPFR_GET_ALLOC_SIZE(x);
      if (nw > ow)
       {
         mpfr_size_limb_t *tmpx;

         /* Realloc significand */
         tmpx = (mpfr_size_limb_t *) mpfr_reallocate_func
           (MPFR_GET_REAL_PTR(x), MPFR_MALLOC_SIZE(ow), MPFR_MALLOC_SIZE(nw));
         MPFR_SET_MANT_PTR(x, tmpx); /* mant ptr must be set
                                        before alloc size */
         MPFR_SET_ALLOC_SIZE(x, nw); /* new number of allocated limbs */
       }
    }

  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) ))
    {
      MPFR_PREC(x) = prec; /* Special value: need to set prec */
      if (MPFR_IS_NAN(x))
        MPFR_RET_NAN;
      MPFR_ASSERTD(MPFR_IS_INF(x) || MPFR_IS_ZERO(x));
      return 0; /* infinity and zero are exact */
    }

  /* x is a non-zero real number */

  MPFR_TMP_MARK(marker);
  tmp = MPFR_TMP_LIMBS_ALLOC (nw);
  xp = MPFR_MANT(x);
  carry = mpfr_round_raw (tmp, xp, MPFR_PREC(x), MPFR_IS_NEG(x),
                          prec, rnd_mode, &inexact);
  MPFR_PREC(x) = prec;

  if (MPFR_UNLIKELY(carry))
    {
      mpfr_exp_t exp = MPFR_EXP (x);

      if (MPFR_UNLIKELY(exp == __gmpfr_emax))
        (void) mpfr_overflow(x, rnd_mode, MPFR_SIGN(x));
      else
        {
          MPFR_ASSERTD (exp < __gmpfr_emax);
          MPFR_SET_EXP (x, exp + 1);
          xp[nw - 1] = MPFR_LIMB_HIGHBIT;
          if (nw - 1 > 0)
            MPN_ZERO(xp, nw - 1);
        }
    }
  else
    MPN_COPY(xp, tmp, nw);

  MPFR_TMP_FREE(marker);
  return inexact;
}

/* assumption: GMP_NUMB_BITS is a power of 2 */

/* assuming b is an approximation to x in direction rnd1 with error at
   most 2^(MPFR_EXP(b)-err), returns 1 if one is able to round exactly
   x to precision prec with direction rnd2, and 0 otherwise.
   Side effects: none.

   rnd1 = RNDN and RNDF are similar: the sign of the error is unknown.

   rnd2 = RNDF: assume that the user will round the approximation b
   toward the direction of x, i.e. the opposite of rnd1 in directed
   rounding modes, otherwise RNDN. Some details:

                u   xinf        v xsup          w
           -----|----+----------|--+------------|-----
                     [----- x -----]
     rnd1 = RNDD     b             |
     rnd1 = RNDU                   b

   where u, v and w are consecutive machine numbers.

   * If [xinf,xsup] contains no machine numbers, then return 1.

   * If [xinf,xsup] contains 2 machine numbers, then return 0.

   * If [xinf,xsup] contains a single machine number, then return 1 iff
     the rounding of b is this machine number.
     With the above choice for the rounding of b, this will always be
     the case if rnd1 is a directed rounding mode; said otherwise, for
     rnd2 = RNDF and rnd1 being a directed rounding mode, return 1 iff
     [xinf,xsup] contains at most 1 machine number.
*/

int
mpfr_can_round (mpfr_srcptr b, mpfr_exp_t err, mpfr_rnd_t rnd1,
                mpfr_rnd_t rnd2, mpfr_prec_t prec)
{
  if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(b)))
    return 0; /* We cannot round if Zero, Nan or Inf */
  else
    return mpfr_can_round_raw (MPFR_MANT(b), MPFR_LIMB_SIZE(b),
                               MPFR_SIGN(b), err, rnd1, rnd2, prec);
}

/* TODO: mpfr_can_round_raw currently does a memory allocation and some
   mpn operations. A bit inspection like for mpfr_round_p (round_p.c) may
   be sufficient, though this would be more complex than the one done in
   mpfr_round_p, and in particular, for some rnd1/rnd2 combinations, one
   needs to take care of changes of binade when the value is close to a
   power of 2. */

int
mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err,
                    mpfr_rnd_t rnd1, mpfr_rnd_t rnd2, mpfr_prec_t prec)
{
  mpfr_prec_t prec2;
  mp_size_t k, k1, tn;
  int s, s1;
  mp_limb_t cc, cc2;
  mp_limb_t *tmp;
  mp_limb_t cy = 0, tmp_hi;
  int res;
  MPFR_TMP_DECL(marker);

  /* Since mpfr_can_round is a function in the API, use MPFR_ASSERTN.
     The specification makes sense only for prec >= 1. */
  MPFR_ASSERTN (prec >= 1);

  MPFR_ASSERTD(bp[bn - 1] & MPFR_LIMB_HIGHBIT);

  MPFR_ASSERT_SIGN(neg);
  neg = MPFR_IS_NEG_SIGN(neg);

  /* Transform RNDD and RNDU to Zero / Away */
  MPFR_ASSERTD (neg == 0 || neg == 1);
  /* transform RNDF to RNDN */
  if (rnd1 == MPFR_RNDF)
    rnd1 = MPFR_RNDN;
  if (rnd1 != MPFR_RNDN)
    rnd1 = MPFR_IS_LIKE_RNDZ(rnd1, neg) ? MPFR_RNDZ : MPFR_RNDA;
  if (rnd2 == MPFR_RNDF)
    rnd2 = (rnd1 == MPFR_RNDN) ? MPFR_RNDN :
      MPFR_IS_LIKE_RNDZ(rnd1, neg) ? MPFR_RNDA : MPFR_RNDZ;
  if (rnd2 != MPFR_RNDN)
    rnd2 = MPFR_IS_LIKE_RNDZ(rnd2, neg) ? MPFR_RNDZ : MPFR_RNDA;

  /* now rnd2 is either RNDN, RNDZ or RNDZ */

  /* For err < prec (+1 for rnd1=RNDN), we can never round correctly, since
     the error is at least 2*ulp(b) >= ulp(round(b)).
     However for err = prec (+1 for rnd1=RNDN), we can round correctly in some
     rare cases where ulp(b) = 1/2*ulp(U) [see below for the definition of U],
     which implies rnd1 = RNDZ or RNDN, and rnd2 = RNDA or RNDN. */

  if (MPFR_UNLIKELY (err < prec + (rnd1 == MPFR_RNDN) ||
                     (err == prec + (rnd1 == MPFR_RNDN) &&
                      (rnd1 == MPFR_RNDA ||
                       rnd2 == MPFR_RNDZ))))
    return 0;  /* can't round */

  /* As a consequence... */
  MPFR_ASSERTD (err >= prec);

  /* The bound c on the error |x-b| is: c = 2^(MPFR_EXP(b)-err) <= b/2.
   * So, we now know that x and b have the same sign. By symmetry,
   * assume x > 0 and b > 0. We have: L <= x <= U, where, depending
   * on rnd1:
   *   MPFR_RNDN: L = b-c, U = b+c
   *   MPFR_RNDZ: L = b,   U = b+c
   *   MPFR_RNDA: L = b-c, U = b
   *
   * We can round x iff round(L,prec,rnd2) = round(U,prec,rnd2).
   */

  if (MPFR_UNLIKELY (prec > (mpfr_prec_t) bn * GMP_NUMB_BITS))
    { /* Then prec > PREC(b): we can round:
         (i) in rounding to the nearest as long as err >= prec + 2.
             When err = prec + 1 and b is not a power
             of two (so that a change of binade cannot occur), then one
             can round to nearest thanks to the even rounding rule (in the
             target precision prec, the significand of b ends with a 0).
             When err = prec + 1 and b is a power of two, when rnd1 = RNDZ one
             can round too.
         (ii) in directed rounding mode iff rnd1 is compatible with rnd2
              and err >= prec + 1, unless b = 2^k and rnd1 = RNDA or RNDN in
              which case we need err >= prec + 2.
      */
      if ((rnd1 == rnd2 || rnd2 == MPFR_RNDN) && err >= prec + 1)
        {
          if (rnd1 != MPFR_RNDZ &&
              err == prec + 1 &&
              mpfr_powerof2_raw2 (bp, bn))
            return 0;
          else
            return 1;
        }
      return 0;
    }

  /* now prec <= bn * GMP_NUMB_BITS */

  if (MPFR_UNLIKELY (err > (mpfr_prec_t) bn * GMP_NUMB_BITS))
    {
      /* we distinguish the case where b is a power of two:
         rnd1 rnd2 can round?
         RNDZ RNDZ ok
         RNDZ RNDA no
         RNDZ RNDN ok
         RNDA RNDZ no
         RNDA RNDA ok except when err = prec + 1
         RNDA RNDN ok except when err = prec + 1
         RNDN RNDZ no
         RNDN RNDA no
         RNDN RNDN ok except when err = prec + 1 */
      if (mpfr_powerof2_raw2 (bp, bn))
        {
          if ((rnd2 == MPFR_RNDZ || rnd2 == MPFR_RNDA) && rnd1 != rnd2)
            return 0;
          else if (rnd1 == MPFR_RNDZ)
            return 1; /* RNDZ RNDZ and RNDZ RNDN */
          else
            return err > prec + 1;
        }

      /* now the general case where b is not a power of two:
         rnd1 rnd2 can round?
         RNDZ RNDZ ok
         RNDZ RNDA except when b is representable in precision 'prec'
         RNDZ RNDN except when b is the middle of two representable numbers in
                   precision 'prec' and b ends with 'xxx0[1]',
                   or b is representable in precision 'prec'
                   and err = prec + 1 and b ends with '1'.
         RNDA RNDZ except when b is representable in precision 'prec'
         RNDA RNDA ok
         RNDA RNDN except when b is the middle of two representable numbers in
                   precision 'prec' and b ends with 'xxx1[1]',
                   or b is representable in precision 'prec'
                   and err = prec + 1 and b ends with '1'.
         RNDN RNDZ except when b is representable in precision 'prec'
         RNDN RNDA except when b is representable in precision 'prec'
         RNDN RNDN except when b is the middle of two representable numbers in
                   precision 'prec', or b is representable in precision 'prec'
                   and err = prec + 1 and b ends with '1'. */
      if (rnd2 == MPFR_RNDN)
        {
          if (err == prec + 1 && (bp[0] & 1))
            return 0; /* err == prec + 1 implies prec = bn * GMP_NUMB_BITS */
          if (prec < (mpfr_prec_t) bn * GMP_NUMB_BITS)
            {
              k1 = MPFR_PREC2LIMBS (prec + 1);
              MPFR_UNSIGNED_MINUS_MODULO(s1, prec + 1);
              if (((bp[bn - k1] >> s1) & 1) &&
                  mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec + 1) == 0)
                { /* b is the middle of two representable numbers */
                  if (rnd1 == MPFR_RNDN)
                    return 0;
                  k1 = MPFR_PREC2LIMBS (prec);
                  MPFR_UNSIGNED_MINUS_MODULO(s1, prec);
                  return (rnd1 == MPFR_RNDZ) ^
                    (((bp[bn - k1] >> s1) & 1) == 0);
                }
            }
          return 1;
        }
      else if (rnd1 == rnd2) /* cases RNDZ RNDZ or RNDA RNDA: ok */
        return 1;
      else
        return mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec) != 0;
    }

  /* now err <= bn * GMP_NUMB_BITS */

  /* warning: if k = m*GMP_NUMB_BITS, consider limb m-1 and not m */
  k = (err - 1) / GMP_NUMB_BITS;
  MPFR_UNSIGNED_MINUS_MODULO(s, err);
  /* the error corresponds to bit s in limb k, the most significant limb
     being limb 0; in memory, limb k is bp[bn-1-k]. */

  k1 = (prec - 1) / GMP_NUMB_BITS;
  MPFR_UNSIGNED_MINUS_MODULO(s1, prec);
  /* the least significant bit is bit s1 in limb k1 */

  /* We don't need to consider the k1 most significant limbs.
     They will be considered later only to detect when subtracting
     the error bound yields a change of binade.
     Warning! The number with updated bn may no longer be normalized. */
  k -= k1;
  bn -= k1;
  prec2 = prec - (mpfr_prec_t) k1 * GMP_NUMB_BITS;

  /* We can decide of the correct rounding if rnd2(b-eps) and rnd2(b+eps)
     give the same result to the target precision 'prec', i.e., if when
     adding or subtracting (1 << s) in bp[bn-1-k], it does not change the
     rounding in direction 'rnd2' at ulp-position bp[bn-1] >> s1, taking also
     into account the possible change of binade. */
  MPFR_TMP_MARK(marker);
  tn = bn;
  k++; /* since we work with k+1 everywhere */
  tmp = MPFR_TMP_LIMBS_ALLOC (tn);
  if (bn > k)
    MPN_COPY (tmp, bp, bn - k);

  MPFR_ASSERTD (k > 0);

  switch (rnd1)
    {
    case MPFR_RNDZ:
      /* rnd1 = Round to Zero */
      cc = (bp[bn - 1] >> s1) & 1;
      /* mpfr_round_raw2 returns 1 if one should add 1 at ulp(b,prec),
         and 0 otherwise */
      cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec2);
      /* cc is the new value of bit s1 in bp[bn-1] after rounding 'rnd2' */

      /* now round b + 2^(MPFR_EXP(b)-err) */
      cy = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
      /* propagate carry up to most significant limb */
      for (tn = 0; tn + 1 < k1 && cy != 0; tn ++)
        cy = ~bp[bn + tn] == 0;
      if (cy == 0 && err == prec)
        {
          res = 0;
          goto end;
        }
      if (MPFR_UNLIKELY(cy))
        {
          /* when a carry occurs, we have b < 2^h <= b+c, we can round iff:
             rnd2 = RNDZ: never, since b and b+c round to different values;
             rnd2 = RNDA: when b+c is an exact power of two, and err > prec
                          (since for err = prec, b = 2^h - 1/2*ulp(2^h) is
                          exactly representable and thus rounds to itself);
             rnd2 = RNDN: whenever cc = 0, since err >= prec implies
                          c <= ulp(b) = 1/2*ulp(2^h), thus b+c rounds to 2^h,
                          and b+c >= 2^h implies that bit 'prec' of b is 1,
                          thus cc = 0 means that b is rounded to 2^h too. */
          res = (rnd2 == MPFR_RNDZ) ? 0
            : (rnd2 == MPFR_RNDA) ? (err > prec && k == bn && tmp[0] == 0)
            : cc == 0;
          goto end;
        }
      break;
    case MPFR_RNDN:
      /* rnd1 = Round to nearest */

      /* first round b+2^(MPFR_EXP(b)-err) */
      cy = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
      /* propagate carry up to most significant limb */
      for (tn = 0; tn + 1 < k1 && cy != 0; tn ++)
        cy = ~bp[bn + tn] == 0;
      cc = (tmp[bn - 1] >> s1) & 1; /* gives 0 when cc=1 */
      cc ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec2);
      /* cc is the new value of bit s1 in bp[bn-1]+eps after rounding 'rnd2' */
      if (MPFR_UNLIKELY (cy != 0))
        {
          /* when a carry occurs, we have b-c < b < 2^h <= b+c, we can round
             iff:
             rnd2 = RNDZ: never, since b-c and b+c round to different values;
             rnd2 = RNDA: when b+c is an exact power of two, and
                          err > prec + 1 (since for err <= prec + 1,
                          b-c <= 2^h - 1/2*ulp(2^h) is exactly representable
                          and thus rounds to itself);
             rnd2 = RNDN: whenever err > prec + 1, since for err = prec + 1,
                          b+c rounds to 2^h, and b-c rounds to nextbelow(2^h).
                          For err > prec + 1, c <= 1/4*ulp(b) <= 1/8*ulp(2^h),
                          thus
                          2^h - 1/4*ulp(b) <= b-c < b+c <= 2^h + 1/8*ulp(2^h),
                          therefore both b-c and b+c round to 2^h. */
          res = (rnd2 == MPFR_RNDZ) ? 0
            : (rnd2 == MPFR_RNDA) ? (err > prec + 1 && k == bn && tmp[0] == 0)
            : err > prec + 1;
          goto end;
        }
    subtract_eps:
      /* now round b-2^(MPFR_EXP(b)-err), this happens for
         rnd1 = RNDN or RNDA */
      MPFR_ASSERTD(rnd1 == MPFR_RNDN || rnd1 == MPFR_RNDA);
      cy = mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
      /* propagate the potential borrow up to the most significant limb
         (it cannot propagate further since the most significant limb is
         at least MPFR_LIMB_HIGHBIT).
         Note: we use the same limb tmp[bn-1] to subtract. */
      tmp_hi = tmp[bn - 1];
      for (tn = 0; tn < k1 && cy != 0; tn ++)
        cy = mpn_sub_1 (&tmp_hi, bp + bn + tn, 1, cy);
      /* We have an exponent decrease when tn = k1 and
         tmp[bn-1] < MPFR_LIMB_HIGHBIT:
         b-c < 2^h <= b (for RNDA) or b+c (for RNDN).
         Then we surely cannot round when rnd2 = RNDZ, since b or b+c round to
         a value >= 2^h, and b-c rounds to a value < 2^h.
         We also surely cannot round when (rnd1,rnd2) = (RNDN,RNDA), since
         b-c rounds to a value <= 2^h, and b+c > 2^h rounds to a value > 2^h.
         It thus remains:
         (rnd1,rnd2) = (RNDA,RNDA), (RNDA,RNDN) and (RNDN,RNDN).
         For (RNDA,RNDA) we can round only when b-c and b round to 2^h, which
         implies b = 2^h and err > prec (which is true in that case):
         a necessary condition is that cc = 0.
         For (RNDA,RNDN) we can round only when b-c and b round to 2^h, which
         implies b-c >= 2^h - 1/4*ulp(2^h), and b <= 2^h + 1/2*ulp(2^h);
         since ulp(2^h) = ulp(b), this implies c <= 3/4*ulp(b), thus
         err > prec.
         For (RNDN,RNDN) we can round only when b-c and b+c round to 2^h,
         which implies b-c >= 2^h - 1/4*ulp(2^h), and
         b+c <= 2^h + 1/2*ulp(2^h);
         since ulp(2^h) = ulp(b), this implies 2*c <= 3/4*ulp(b), thus
         err > prec+1.
      */
      if (tn == k1 && tmp_hi < MPFR_LIMB_HIGHBIT) /* exponent decrease */
        {
          if (rnd2 == MPFR_RNDZ || (rnd1 == MPFR_RNDN && rnd2 == MPFR_RNDA) ||
              cc != 0 /* b or b+c does not round to 2^h */)
            {
              res = 0;
              goto end;
            }
          /* in that case since the most significant bit of tmp is 0, we
             should consider one more bit; res = 0 when b-c does not round
             to 2^h. */
          res = mpfr_round_raw2 (tmp, bn, neg, rnd2, prec2 + 1) != 0;
          goto end;
        }
      if (err == prec + (rnd1 == MPFR_RNDN))
        {
          /* No exponent increase nor decrease, thus we have |U-L| = ulp(b).
             For rnd2 = RNDZ or RNDA, either [L,U] contains one representable
             number in the target precision, and then L and U round
             differently; or both L and U are representable: they round
             differently too; thus in all cases we cannot round.
             For rnd2 = RNDN, the only case where we can round is when the
             middle of [L,U] (i.e. b) is representable, and ends with a 0. */
          res = (rnd2 == MPFR_RNDN && (((bp[bn - 1] >> s1) & 1) == 0) &&
                 mpfr_round_raw2 (bp, bn, neg, MPFR_RNDZ, prec2) ==
                 mpfr_round_raw2 (bp, bn, neg, MPFR_RNDA, prec2));
          goto end;
        }
      break;
    default:
      /* rnd1 = Round away */
      MPFR_ASSERTD (rnd1 == MPFR_RNDA);
      cc = (bp[bn - 1] >> s1) & 1;
      /* the mpfr_round_raw2() call below returns whether one should add 1 or
         not for rounding */
      cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec2);
      /* cc is the new value of bit s1 in bp[bn-1]+eps after rounding 'rnd2' */

      goto subtract_eps;
    }

  cc2 = (tmp[bn - 1] >> s1) & 1;
  res = cc == (cc2 ^ mpfr_round_raw2 (tmp, bn, neg, rnd2, prec2));

 end:
  MPFR_TMP_FREE(marker);
  return res;
}