contrib/mpfr/round_near_x.c

   1 /* mpfr_round_near_x -- Round a floating point number nears another one.
   2
   3 Copyright 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, and was contributed by Mathieu Dutour.
   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 /* Use MPFR_FAST_COMPUTE_IF_SMALL_INPUT instead (a simple wrapper) */
  26
  27 /* int mpfr_round_near_x (mpfr_ptr y, mpfr_srcptr v, mpfr_uexp_t err, int dir,
  28                           mp_rnd_t rnd)
  29
  30    TODO: fix this description.
  31    Assuming y = o(f(x)) = o(x + g(x)) with |g(x)| < 2^(EXP(v)-error)
  32    If x is small enough, y ~= v. This function checks and does this.
  33
  34    It assumes that f(x) is not representable exactly as a FP number.
  35    v must not be a singular value (NAN, INF or ZERO), usual values are
  36    v=1 or v=x.
  37
  38    y is the destination (a mpfr_t), v the value to set (a mpfr_t),
  39    err the error term (a mpfr_uexp_t) such that |g(x)| < 2^(EXP(x)-err),
  40    dir (an int) is the direction of the error (if dir = 0,
  41    it rounds towards 0, if dir=1, it rounds away from 0),
  42    rnd the rounding mode.
  43
  44    It returns 0 if it can't round.
  45    Otherwise it returns the ternary flag (It can't return an exact value).
  46 */
  47
  48 /* What "small enough" means?
  49
  50    We work with the positive values.
  51    Assuming err > Prec (y)+1
  52
  53    i = [ y = o(x)]   // i = inexact flag
  54    If i == 0
  55        Setting x in y is exact. We have:
  56        y = [XXXXXXXXX[...]]0[...] + error where [..] are optional zeros
  57       if dirError = ToInf,
  58         x < f(x) < x + 2^(EXP(x)-err)
  59         since x=y, and ulp (y)/2 > 2^(EXP(x)-err), we have:
  60         y < f(x) < y+ulp(y) and |y-f(x)| < ulp(y)/2
  61        if rnd = RNDN, nothing
  62        if rnd = RNDZ, nothing
  63        if rnd = RNDA, addoneulp
  64       elif dirError = ToZero
  65         x -2^(EXP(x)-err) < f(x) < x
  66         since x=y, and ulp (y)/2 > 2^(EXP(x)-err), we have:
  67         y-ulp(y) < f(x) < y and |y-f(x)| < ulp(y)/2
  68        if rnd = RNDN, nothing
  69        if rnd = RNDZ, nexttozero
  70        if rnd = RNDA, nothing
  71      NOTE: err > prec (y)+1 is needed only for RNDN.
  72    elif i > 0 and i = EVEN_ROUNDING
  73       So rnd = RNDN and we have y = x + ulp(y)/2
  74        if dirError = ToZero,
  75          we have x -2^(EXP(x)-err) < f(x) < x
  76          so y - ulp(y)/2 - 2^(EXP(x)-err) < f(x) < y-ulp(y)/2
  77          so y -ulp(y) < f(x) < y-ulp(y)/2
  78          => nexttozero(y)
  79        elif dirError = ToInf
  80          we have x < f(x) < x + 2^(EXP(x)-err)
  81          so y - ulp(y)/2 < f(x) < y+ulp(y)/2-ulp(y)/2
  82          so y - ulp(y)/2 < f(x) < y
  83          => do nothing
  84    elif i < 0 and i = -EVEN_ROUNDING
  85       So rnd = RNDN and we have y = x - ulp(y)/2
  86       if dirError = ToZero,
  87         y < f(x) < y + ulp(y)/2 => do nothing
  88       if dirError = ToInf
  89         y + ulp(y)/2 < f(x) < y + ulp(y) => AddOneUlp
  90    elif i > 0
  91      we can't have rnd = RNDZ, and prec(x) > prec(y), so ulp(x) < ulp(y)
  92      we have y - ulp (y) < x < y
  93      or more exactly y - ulp(y) + ulp(x)/2 <= x <= y - ulp(x)/2
  94      if rnd = RNDA,
  95       if dirError = ToInf,
  96        we have x < f(x) < x + 2^(EXP(x)-err)
  97        if err > prec (x),
  98          we have 2^(EXP(x)-err) < ulp(x), so 2^(EXP(x)-err) <= ulp(x)/2
  99          so f(x) <= y - ulp(x)/2+ulp(x)/2 <= y
 100          and y - ulp(y) < x < f(x)
 101          so we have y - ulp(y) < f(x) < y
 102          so do nothing.
 103        elif we can round, ie y - ulp(y) < x + 2^(EXP(x)-err) < y
 104          we have y - ulp(y) < x <  f(x) < x + 2^(EXP(x)-err) < y
 105          so do nothing
 106        otherwise
 107          Wrong. Example X=[0.11101]111111110000
 108                          +             1111111111111111111....
 109       elif dirError = ToZero
 110        we have x - 2^(EXP(x)-err) < f(x) < x
 111        so f(x) < x < y
 112        if err > prec (x)
 113          x-2^(EXP(x)-err) >= x-ulp(x)/2 >= y - ulp(y) + ulp(x)/2-ulp(x)/2
 114          so y - ulp(y) < f(x) < y
 115          so do nothing
 116        elif we can round, ie y - ulp(y) < x - 2^(EXP(x)-err) < y
 117          y - ulp(y) < x - 2^(EXP(x)-err) < f(x) < y
 118          so do nothing
 119        otherwise
 120         Wrong. Example: X=[1.111010]00000010
 121                          -             10000001000000000000100....
 122      elif rnd = RNDN,
 123       y - ulp(y)/2 < x < y and we can't have x = y-ulp(y)/2:
 124       so we have:
 125        y - ulp(y)/2 + ulp(x)/2 <= x <= y - ulp(x)/2
 126       if dirError = ToInf
 127         we have x < f(x) < x+2^(EXP(x)-err) and ulp(y) > 2^(EXP(x)-err)
 128         so y - ulp(y)/2 + ulp (x)/2 < f(x) < y + ulp (y)/2 - ulp (x)/2
 129         we can round but we can't compute inexact flag.
 130         if err > prec (x)
 131           y - ulp(y)/2 + ulp (x)/2 < f(x) < y + ulp(x)/2 - ulp(x)/2
 132           so y - ulp(y)/2 + ulp (x)/2 < f(x) < y
 133           we can round and compute inexact flag. do nothing
 134         elif we can round, ie y - ulp(y)/2 < x + 2^(EXP(x)-err) < y
 135           we have  y - ulp(y)/2 + ulp (x)/2 < f(x) < y
 136           so do nothing
 137         otherwise
 138           Wrong
 139       elif dirError = ToZero
 140         we have x -2^(EXP(x)-err) < f(x) < x and ulp(y)/2 > 2^(EXP(x)-err)
 141         so y-ulp(y)+ulp(x)/2 < f(x) < y - ulp(x)/2
 142         if err > prec (x)
 143            x- ulp(x)/2 < f(x) < x
 144            so y - ulp(y)/2+ulp(x)/2 - ulp(x)/2 < f(x) < x <= y - ulp(x)/2 < y
 145            do nothing
 146         elif we can round, ie y-ulp(y)/2 < x-2^(EXP(x)-err) < y
 147            we have y-ulp(y)/2 < x-2^(EXP(x)-err) < f(x) < x < y
 148            do nothing
 149         otherwise
 150           Wrong
 151    elif i < 0
 152      same thing?
 153  */
 154
 155 int
 156 mpfr_round_near_x (mpfr_ptr y, mpfr_srcptr v, mpfr_uexp_t err, int dir,
 157                    mp_rnd_t rnd)
 158 {
 159   int inexact, sign;
 160   unsigned int old_flags = __gmpfr_flags;
 161
 162   MPFR_ASSERTD (!MPFR_IS_SINGULAR (v));
 163   MPFR_ASSERTD (dir == 0 || dir == 1);
 164
 165   /* First check if we can round. The test is more restrictive than
 166      necessary. Note that if err is not representable in an mp_exp_t,
 167      then err > MPFR_PREC (v) and the conversion to mp_exp_t will not
 168      occur. */
 169   if (!(err > MPFR_PREC (y) + 1
 170         && (err > MPFR_PREC (v)
 171             || mpfr_round_p (MPFR_MANT (v), MPFR_LIMB_SIZE (v),
 172                              (mp_exp_t) err,
 173                              MPFR_PREC (y) + (rnd == GMP_RNDN)))))
 174     /* If we assume we can not round, return 0, and y is not modified */
 175     return 0;
 176
 177   /* First round v in y */
 178   sign = MPFR_SIGN (v);
 179   MPFR_SET_EXP (y, MPFR_GET_EXP (v));
 180   MPFR_SET_SIGN (y, sign);
 181   MPFR_RNDRAW_GEN (inexact, y, MPFR_MANT (v), MPFR_PREC (v), rnd, sign,
 182                    if (dir == 0)
 183                      {
 184                        inexact = -sign;
 185                        goto trunc_doit;
 186                      }
 187                    else
 188                      goto addoneulp;
 189                    , if (MPFR_UNLIKELY (++MPFR_EXP (y) > __gmpfr_emax))
 190                        mpfr_overflow (y, rnd, sign)
 191                   );
 192
 193   /* Fix it in some cases */
 194   MPFR_ASSERTD (!MPFR_IS_NAN (y) && !MPFR_IS_ZERO (y));
 195   /* If inexact == 0, setting y from v is exact but we haven't
 196      take into account yet the error term */
 197   if (inexact == 0)
 198     {
 199       if (dir == 0) /* The error term is negative for v positive */
 200         {
 201           inexact = sign;
 202           if (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG_SIGN (sign)))
 203             {
 204               /* case nexttozero */
 205               /* The underflow flag should be set if the result is zero */
 206               __gmpfr_flags = old_flags;
 207               inexact = -sign;
 208               mpfr_nexttozero (y);
 209               if (MPFR_UNLIKELY (MPFR_IS_ZERO (y)))
 210                 mpfr_set_underflow ();
 211             }
 212         }
 213       else /* The error term is positive for v positive */
 214         {
 215           inexact = -sign;
 216           /* Round Away */
 217           if (rnd != GMP_RNDN && rnd != GMP_RNDZ
 218               && MPFR_IS_RNDUTEST_OR_RNDDNOTTEST (rnd, MPFR_IS_POS_SIGN(sign)))
 219             {
 220               /* case nexttoinf */
 221               /* The overflow flag should be set if the result is infinity */
 222               inexact = sign;
 223               mpfr_nexttoinf (y);
 224               if (MPFR_UNLIKELY (MPFR_IS_INF (y)))
 225                 mpfr_set_overflow ();
 226             }
 227         }
 228     }
 229
 230   /* the inexact flag cannot be 0, since this would mean an exact value,
 231      and in this case we cannot round correctly */
 232   MPFR_ASSERTD(inexact != 0);
 233   MPFR_RET (inexact);
 234 }