contrib/mpfr/src/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, 2010, 2011 Free Software Foundation, Inc.
   5 Contributed by the Arenaire and Caramel 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 3 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.LESSER.  If not, see
  21 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  22 51 Franklin St, Fifth Floor, Boston, 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, mpfr_prec_t prec, mpfr_rnd_t rnd_mode)
  50 {
  51   mp_limb_t *tmp, *xp;
  52   int carry, inexact;
  53   mpfr_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) / GMP_NUMB_BITS; /* needed allocated limbs */
  59
  60   /* check if x has enough allocated space for the significand */
  61   /* Get the number of limbs from the precision.
  62      (Compatible with all allocation methods) */
  63   ow = (MPFR_PREC (x) + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
  64   if (nw > ow)
  65     {
  66       /* FIXME: Variable can't be created using custom allocation,
  67          MPFR_DECL_INIT or GROUP_ALLOC: How to detect? */
  68       ow = MPFR_GET_ALLOC_SIZE(x);
  69       if (nw > ow)
  70        {
  71          /* Realloc significand */
  72          mpfr_limb_ptr tmpx = (mpfr_limb_ptr) (*__gmp_reallocate_func)
  73            (MPFR_GET_REAL_PTR(x), MPFR_MALLOC_SIZE(ow), MPFR_MALLOC_SIZE(nw));
  74          MPFR_SET_MANT_PTR(x, tmpx); /* mant ptr must be set
  75                                         before alloc size */
  76          MPFR_SET_ALLOC_SIZE(x, nw); /* new number of allocated limbs */
  77        }
  78     }
  79
  80   if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) ))
  81     {
  82       MPFR_PREC(x) = prec; /* Special value: need to set prec */
  83       if (MPFR_IS_NAN(x))
  84         MPFR_RET_NAN;
  85       MPFR_ASSERTD(MPFR_IS_INF(x) || MPFR_IS_ZERO(x));
  86       return 0; /* infinity and zero are exact */
  87     }
  88
  89   /* x is a non-zero real number */
  90
  91   MPFR_TMP_MARK(marker);
  92   tmp = MPFR_TMP_LIMBS_ALLOC (nw);
  93   xp = MPFR_MANT(x);
  94   carry = mpfr_round_raw (tmp, xp, MPFR_PREC(x), MPFR_IS_NEG(x),
  95                           prec, rnd_mode, &inexact);
  96   MPFR_PREC(x) = prec;
  97
  98   if (MPFR_UNLIKELY(carry))
  99     {
 100       mpfr_exp_t exp = MPFR_EXP (x);
 101
 102       if (MPFR_UNLIKELY(exp == __gmpfr_emax))
 103         (void) mpfr_overflow(x, rnd_mode, MPFR_SIGN(x));
 104       else
 105         {
 106           MPFR_ASSERTD (exp < __gmpfr_emax);
 107           MPFR_SET_EXP (x, exp + 1);
 108           xp[nw - 1] = MPFR_LIMB_HIGHBIT;
 109           if (nw - 1 > 0)
 110             MPN_ZERO(xp, nw - 1);
 111         }
 112     }
 113   else
 114     MPN_COPY(xp, tmp, nw);
 115
 116   MPFR_TMP_FREE(marker);
 117   return inexact;
 118 }
 119
 120 /* assumption: GMP_NUMB_BITS is a power of 2 */
 121
 122 /* assuming b is an approximation to x in direction rnd1 with error at
 123    most 2^(MPFR_EXP(b)-err), returns 1 if one is able to round exactly
 124    x to precision prec with direction rnd2, and 0 otherwise.
 125
 126    Side effects: none.
 127 */
 128
 129 int
 130 mpfr_can_round (mpfr_srcptr b, mpfr_exp_t err, mpfr_rnd_t rnd1,
 131                 mpfr_rnd_t rnd2, mpfr_prec_t prec)
 132 {
 133   if (MPFR_UNLIKELY(MPFR_IS_SINGULAR(b)))
 134     return 0; /* We cannot round if Zero, Nan or Inf */
 135   else
 136     return mpfr_can_round_raw (MPFR_MANT(b), MPFR_LIMB_SIZE(b),
 137                                MPFR_SIGN(b), err, rnd1, rnd2, prec);
 138 }
 139
 140 int
 141 mpfr_can_round_raw (const mp_limb_t *bp, mp_size_t bn, int neg, mpfr_exp_t err0,
 142                     mpfr_rnd_t rnd1, mpfr_rnd_t rnd2, mpfr_prec_t prec)
 143 {
 144   mpfr_prec_t err;
 145   mp_size_t k, k1, tn;
 146   int s, s1;
 147   mp_limb_t cc, cc2;
 148   mp_limb_t *tmp;
 149   MPFR_TMP_DECL(marker);
 150
 151   if (MPFR_UNLIKELY(err0 < 0 || (mpfr_uexp_t) err0 <= prec))
 152     return 0;  /* can't round */
 153   else if (MPFR_UNLIKELY (prec > (mpfr_prec_t) bn * GMP_NUMB_BITS))
 154     { /* then ulp(b) < precision < error */
 155       return rnd2 == MPFR_RNDN && (mpfr_uexp_t) err0 - 2 >= prec;
 156       /* can round only in rounding to the nearest and err0 >= prec + 2 */
 157     }
 158
 159   MPFR_ASSERT_SIGN(neg);
 160   neg = MPFR_IS_NEG_SIGN(neg);
 161
 162   /* if the error is smaller than ulp(b), then anyway it will propagate
 163      up to ulp(b) */
 164   err = ((mpfr_uexp_t) err0 > (mpfr_prec_t) bn * GMP_NUMB_BITS) ?
 165     (mpfr_prec_t) bn * GMP_NUMB_BITS : (mpfr_prec_t) err0;
 166
 167   /* warning: if k = m*GMP_NUMB_BITS, consider limb m-1 and not m */
 168   k = (err - 1) / GMP_NUMB_BITS;
 169   MPFR_UNSIGNED_MINUS_MODULO(s, err);
 170   /* the error corresponds to bit s in limb k, the most significant limb
 171      being limb 0 */
 172
 173   k1 = (prec - 1) / GMP_NUMB_BITS;
 174   MPFR_UNSIGNED_MINUS_MODULO(s1, prec);
 175   /* the last significant bit is bit s1 in limb k1 */
 176
 177   /* don't need to consider the k1 most significant limbs */
 178   k -= k1;
 179   bn -= k1;
 180   prec -= (mpfr_prec_t) k1 * GMP_NUMB_BITS;
 181
 182   /* if when adding or subtracting (1 << s) in bp[bn-1-k], it does not
 183      change bp[bn-1] >> s1, then we can round */
 184   MPFR_TMP_MARK(marker);
 185   tn = bn;
 186   k++; /* since we work with k+1 everywhere */
 187   tmp = MPFR_TMP_LIMBS_ALLOC (tn);
 188   if (bn > k)
 189     MPN_COPY (tmp, bp, bn - k);
 190
 191   MPFR_ASSERTD (k > 0);
 192
 193   /* Transform RNDD and RNDU to Zero / Away */
 194   MPFR_ASSERTD((neg == 0) || (neg ==1));
 195   if (MPFR_IS_RNDUTEST_OR_RNDDNOTTEST(rnd1, neg))
 196     rnd1 = MPFR_RNDZ;
 197
 198   switch (rnd1)
 199     {
 200     case MPFR_RNDZ:
 201       /* Round to Zero */
 202       cc = (bp[bn - 1] >> s1) & 1;
 203       /* mpfr_round_raw2 returns 1 if one should add 1 at ulp(b,prec),
 204          and 0 otherwise */
 205       cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec);
 206       /* cc is the new value of bit s1 in bp[bn-1] */
 207       /* now round b + 2^(MPFR_EXP(b)-err) */
 208       cc2 = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
 209       break;
 210     case MPFR_RNDN:
 211       /* Round to nearest */
 212        /* first round b+2^(MPFR_EXP(b)-err) */
 213       cc = mpn_add_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
 214       cc = (tmp[bn - 1] >> s1) & 1; /* gives 0 when cc=1 */
 215       cc ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec);
 216       /* now round b-2^(MPFR_EXP(b)-err) */
 217       cc2 = mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
 218       break;
 219     default:
 220       /* Round away */
 221       cc = (bp[bn - 1] >> s1) & 1;
 222       cc ^= mpfr_round_raw2 (bp, bn, neg, rnd2, prec);
 223       /* now round b +/- 2^(MPFR_EXP(b)-err) */
 224       cc2 = mpn_sub_1 (tmp + bn - k, bp + bn - k, k, MPFR_LIMB_ONE << s);
 225       break;
 226     }
 227
 228   /* if cc2 is 1, then a carry or borrow propagates to the next limb */
 229   if (cc2 && cc)
 230     {
 231       MPFR_TMP_FREE(marker);
 232       return 0;
 233     }
 234
 235   cc2 = (tmp[bn - 1] >> s1) & 1;
 236   cc2 ^= mpfr_round_raw2 (tmp, bn, neg, rnd2, prec);
 237
 238   MPFR_TMP_FREE(marker);
 239   return cc == cc2;
 240 }