contrib/mpfr/round_prec.c

   1 /* mpfr_round_raw_generic, mpfr_round_raw2, mpfr_round_raw, mpfr_prec_round,
   2    mpfr_can_round, mpfr_can_round_raw -- various rounding functions
   3
   4 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
   5 Contributed by the Arenaire and Cacao projects, INRIA.
   6
   7 This file is part of the GNU MPFR Library.
   8
   9 The GNU MPFR Library is free software; you can redistribute it and/or modify
  10 it under the terms of the GNU Lesser General Public License as published by
  11 the Free Software Foundation; either version 2.1 of the License, or (at your
  12 option) any later version.
  13
  14 The GNU MPFR Library is distributed in the hope that it will be useful, but
  15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  16 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  17 License for more details.
  18
  19 You should have received a copy of the GNU Lesser General Public License
  20 along with the GNU MPFR Library; see the file COPYING.LIB.  If not, write to
  21 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
  22 MA 02110-1301, USA. */
  23
  24 #include "mpfr-impl.h"
  25
  26 #define mpfr_round_raw_generic mpfr_round_raw
  27 #define flag 0
  28 #define use_inexp 1
  29 #include "round_raw_generic.c"
  30
  31 #define mpfr_round_raw_generic mpfr_round_raw_2
  32 #define flag 1
  33 #define use_inexp 0
  34 #include "round_raw_generic.c"
  35
  36 /* Seems to be unused. Remove comment to implement it.
  37 #define mpfr_round_raw_generic mpfr_round_raw_3
  38 #define flag 1
  39 #define use_inexp 1
  40 #include "round_raw_generic.c"
  41 */
  42
  43 #define mpfr_round_raw_generic mpfr_round_raw_4
  44 #define flag 0
  45 #define use_inexp 0
  46 #include "round_raw_generic.c"
  47
  48 int
  49 mpfr_prec_round (mpfr_ptr x, mp_prec_t prec, mp_rnd_t rnd_mode)
  50 {
  51   mp_limb_t *tmp, *xp;
  52   int carry, inexact;
  53   mp_prec_t nw, ow;
  54   MPFR_TMP_DECL(marker);
  55
  56   MPFR_ASSERTN(prec >= MPFR_PREC_MIN && prec <= MPFR_PREC_MAX);
  57
  58   nw = 1 + (prec - 1) / BITS_PER_MP_LIMB; /* needed allocated limbs */
  59
  60   /* check if x has enough allocated space for the mantissa */
  61   ow = MPFR_GET_ALLOC_SIZE(x);
  62   if (nw > ow)
  63     {
  64       /* Realloc mantissa */
  65       mp_ptr tmpx = (mp_ptr) (*__gmp_reallocate_func)
  66         (MPFR_GET_REAL_PTR(x),  MPFR_MALLOC_SIZE(ow), MPFR_MALLOC_SIZE(nw));
  67       MPFR_SET_MANT_PTR(x, tmpx); /* mant ptr must be set before alloc size */
  68       MPFR_SET_ALLOC_SIZE(x, nw); /* new number of allocated limbs */
  69     }
  70
  71   if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) ))
  72     {
  73       MPFR_PREC(x) = prec; /* Special value: need to set prec */
  74       if (MPFR_IS_NAN(x))
  75         MPFR_RET_NAN;
  76       MPFR_ASSERTD(MPFR_IS_INF(x) || MPFR_IS_ZERO(x));
  77       return 0; /* infinity and zero are exact */
  78     }
  79
  80   /* x is a non-zero real number */
  81
  82   MPFR_TMP_MARK(marker);
  83   tmp = (mp_limb_t*) MPFR_TMP_ALLOC (nw * BYTES_PER_MP_LIMB);
  84   xp = MPFR_MANT(x);
  85   carry = mpfr_round_raw (tmp, xp, MPFR_PREC(x), MPFR_IS_NEG(x),
  86                           prec, rnd_mode, &inexact);
  87   MPFR_PREC(x) = prec;
  88
  89   if (MPFR_UNLIKELY(carry))
  90     {
  91       mp_exp_t exp = MPFR_EXP (x);
  92
  93       if (MPFR_UNLIKELY(exp == __gmpfr_emax))
  94         (void) mpfr_overflow(x, rnd_mode, MPFR_SIGN(x));
  95       else
  96         {
  97           MPFR_ASSERTD (exp < __gmpfr_emax);
  98           MPFR_SET_EXP (x, exp + 1);
  99           xp[nw - 1] = MPFR_LIMB_HIGHBIT;
 100           if (nw - 1 > 0)
 101             MPN_ZERO(xp, nw - 1);
 102         }
 103     }
 104   else
 105     MPN_COPY(xp, tmp, nw);
 106
 107   MPFR_TMP_FREE(marker);
 108   return inexact;
 109 }
 110
 111 /* assumption: BITS_PER_MP_LIMB is a power of 2 */
 112
 113 /* assuming b is an approximation to x in direction rnd1 with error at
 114    most 2^(MPFR_EXP(b)-err), returns 1 if one is able to round exactly
 115    x to precision prec with direction rnd2, and 0 otherwise.
 116
 117    Side effects: none.
 118 */
 119
 120 int
 121 mpfr_can_round (mpfr_srcptr b, mp_exp_t err, mp_rnd_t rnd1,
 122                 mp_rnd_t rnd2, mp_prec_t prec)
 123 {
 124   if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(b)))
 125     return 0; /* We cannot round if Zero, Nan or Inf */
 126   else
 127     return mpfr_can_round_raw(MPFR_MANT(b), MPFR_LIMB_SIZE(b),
 128                               MPFR_SIGN(b), err, rnd1, rnd2, prec);
 129 }
 130
 131 int
 132 mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mp_exp_t err0,
 133                     mp_rnd_t rnd1, mp_rnd_t rnd2, mp_prec_t prec)
 134 {
 135   mp_prec_t err;
 136   mp_size_t k, k1, tn;
 137   int s, s1;
 138   mp_limb_t cc, cc2;
 139   mp_limb_t *tmp;
 140   MPFR_TMP_DECL(marker);
 141
 142   if (MPFR_UNLIKELY(err0 < 0 || (mp_exp_unsigned_t) err0 <= prec))
 143     return 0;  /* can't round */
 144   else if (MPFR_UNLIKELY (prec > (mp_prec_t) bn * BITS_PER_MP_LIMB))
 145     { /* then ulp(b) < precision < error */
 146       return rnd2 == GMP_RNDN && (mp_exp_unsigned_t) err0 - 2 >= prec;
 147       /* can round only in rounding to the nearest and err0 >= prec + 2 */
 148     }
 149
 150   MPFR_ASSERT_SIGN(neg);
 151   neg = MPFR_IS_NEG_SIGN(neg);
 152
 153   /* if the error is smaller than ulp(b), then anyway it will propagate
 154      up to ulp(b) */
 155   err = ((mp_exp_unsigned_t) err0 > (mp_prec_t) bn * BITS_PER_MP_LIMB) ?
 156     (mp_prec_t) bn * BITS_PER_MP_LIMB : (mp_prec_t) err0;
 157
 158   /* warning: if k = m*BITS_PER_MP_LIMB, consider limb m-1 and not m */
 159   k = (err - 1) / BITS_PER_MP_LIMB;
 160   MPFR_UNSIGNED_MINUS_MODULO(s, err);
 161
 162   /* the error corresponds to bit s in limb k, the most significant limb
 163      being limb 0 */
 164   k1 = (prec - 1) / BITS_PER_MP_LIMB;
 165   MPFR_UNSIGNED_MINUS_MODULO(s1, prec);
 166
 167   /* the last significant bit is bit s1 in limb k1 */
 168
 169   /* don't need to consider the k1 most significant limbs */
 170   k -= k1;
 171   bn -= k1;
 172   prec -= (mp_prec_t) k1 * BITS_PER_MP_LIMB;
 173   /* if when adding or subtracting (1 << s) in bp[bn-1-k], it does not
 174      change bp[bn-1] >> s1, then we can round */
 175
 176   MPFR_TMP_MARK(marker);
 177   tn = bn;
 178   k++; /* since we work with k+1 everywhere */
 179   tmp = (mp_limb_t*) MPFR_TMP_ALLOC(tn * BYTES_PER_MP_LIMB);
 180   if (bn > k)
 181     MPN_COPY (tmp, bp, bn - k);
 182
 183   MPFR_ASSERTD (k > 0);
 184
 185   /* Transform RNDD and RNDU to Zero / Away */
 186   MPFR_ASSERTD((neg == 0) || (neg ==1));
 187   if (MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd1, neg))
 188     rnd1 = GMP_RNDZ;
 189
 190   switch (rnd1)
 191     {
 192     case GMP_RNDZ:
 193       /* Round to Zero */
 194       cc = (bp[bn - 1] >> s1) & 1;
 195       cc ^= mpfr_round_raw2(bp, bn, neg, rnd2, prec);
 196       /* now round b +/- 2^(MPFR_EXP(b)-err) */
 197       cc2 = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
 198       break;
 199     case GMP_RNDN:
 200       /* Round to nearest */
 201        /* first round b+2^(MPFR_EXP(b)-err) */
 202       cc = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
 203       cc = (tmp[bn - 1] >> s1) & 1; /* gives 0 when cc=1 */
 204       cc ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec);
 205       /* now round b-2^(MPFR_EXP(b)-err) */
 206       cc2 = mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
 207       break;
 208     default:
 209       /* Round away */
 210       cc = (bp[bn - 1] >> s1) & 1;
 211       cc ^= mpfr_round_raw2(bp, bn, neg, rnd2, prec);
 212       /* now round b +/- 2^(MPFR_EXP(b)-err) */
 213       cc2 = mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
 214       break;
 215     }
 216
 217   if (cc2 && cc)
 218     {
 219       MPFR_TMP_FREE(marker);
 220       return 0;
 221     }
 222
 223   cc2 = (tmp[bn - 1] >> s1) & 1;
 224   cc2 ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec);
 225
 226   MPFR_TMP_FREE(marker);
 227   return cc == cc2;
 228 }