contrib/mpfr/src/round_raw_generic.c

   1 /* mpfr_round_raw_generic -- Generic rounding function
   2
   3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
   4 Contributed by the Arenaire and Caramel projects, INRIA.
   5
   6 This file is part of the GNU MPFR Library.
   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 3 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.LESSER.  If not, see
  20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
  21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
  22
  23 #ifndef flag
  24 # error "ERROR: flag must be defined (0 / 1)"
  25 #endif
  26 #ifndef use_inexp
  27 # error "ERROR: use_enexp must be defined (0 / 1)"
  28 #endif
  29 #ifndef mpfr_round_raw_generic
  30 # error "ERROR: mpfr_round_raw_generic must be defined"
  31 #endif
  32
  33 /*
  34  * If flag = 0, puts in y the value of xp (with precision xprec and
  35  * sign 1 if negative=0, -1 otherwise) rounded to precision yprec and
  36  * direction rnd_mode. Supposes x is not zero nor NaN nor +/- Infinity
  37  * (i.e. *xp != 0). In that case, the return value is a possible carry
  38  * (0 or 1) that may happen during the rounding, in which case the result
  39  * is a power of two.
  40  *
  41  * If inexp != NULL, put in *inexp the inexact flag of the rounding (0, 1, -1).
  42  * In case of even rounding when rnd = MPFR_RNDN, put MPFR_EVEN_INEX (2) or
  43  * -MPFR_EVEN_INEX (-2) in *inexp.
  44  *
  45  * If flag = 1, just returns whether one should add 1 or not for rounding.
  46  *
  47  * Note: yprec may be < MPFR_PREC_MIN; in particular, it may be equal
  48  * to 1. In this case, the even rounding is done away from 0, which is
  49  * a natural generalization. Indeed, a number with 1-bit precision can
  50  * be seen as a denormalized number with more precision.
  51  */
  52
  53 int
  54 mpfr_round_raw_generic(
  55 #if flag == 0
  56                        mp_limb_t *yp,
  57 #endif
  58                        const mp_limb_t *xp, mpfr_prec_t xprec,
  59                        int neg, mpfr_prec_t yprec, mpfr_rnd_t rnd_mode
  60 #if use_inexp != 0
  61                        , int *inexp
  62 #endif
  63                        )
  64 {
  65   mp_size_t xsize, nw;
  66   mp_limb_t himask, lomask, sb;
  67   int rw;
  68 #if flag == 0
  69   int carry;
  70 #endif
  71 #if use_inexp == 0
  72   int *inexp;
  73 #endif
  74
  75   if (use_inexp)
  76     MPFR_ASSERTD(inexp != ((int*) 0));
  77   MPFR_ASSERTD(neg == 0 || neg == 1);
  78
  79   if (flag && !use_inexp &&
  80       (xprec <= yprec || MPFR_IS_LIKE_RNDZ (rnd_mode, neg)))
  81     return 0;
  82
  83   xsize = (xprec-1)/GMP_NUMB_BITS + 1;
  84   nw = yprec / GMP_NUMB_BITS;
  85   rw = yprec & (GMP_NUMB_BITS - 1);
  86
  87   if (MPFR_UNLIKELY(xprec <= yprec))
  88     { /* No rounding is necessary. */
  89       /* if yp=xp, maybe an overlap: MPN_COPY_DECR is ok when src <= dst */
  90       if (MPFR_LIKELY(rw))
  91         nw++;
  92       MPFR_ASSERTD(nw >= 1);
  93       MPFR_ASSERTD(nw >= xsize);
  94       if (use_inexp)
  95         *inexp = 0;
  96 #if flag == 0
  97       MPN_COPY_DECR(yp + (nw - xsize), xp, xsize);
  98       MPN_ZERO(yp, nw - xsize);
  99 #endif
 100       return 0;
 101     }
 102
 103   if (use_inexp || !MPFR_IS_LIKE_RNDZ(rnd_mode, neg))
 104     {
 105       mp_size_t k = xsize - nw - 1;
 106
 107       if (MPFR_LIKELY(rw))
 108         {
 109           nw++;
 110           lomask = MPFR_LIMB_MASK (GMP_NUMB_BITS - rw);
 111           himask = ~lomask;
 112         }
 113       else
 114         {
 115           lomask = ~(mp_limb_t) 0;
 116           himask = ~(mp_limb_t) 0;
 117         }
 118       MPFR_ASSERTD(k >= 0);
 119       sb = xp[k] & lomask;  /* First non-significant bits */
 120       /* Rounding to nearest ? */
 121       if (MPFR_LIKELY( rnd_mode == MPFR_RNDN) )
 122         {
 123           /* Rounding to nearest */
 124           mp_limb_t rbmask = MPFR_LIMB_ONE << (GMP_NUMB_BITS - 1 - rw);
 125           if (sb & rbmask) /* rounding bit */
 126             sb &= ~rbmask; /* it is 1, clear it */
 127           else
 128             {
 129               /* Rounding bit is 0, behave like rounding to 0 */
 130               goto rnd_RNDZ;
 131             }
 132           while (MPFR_UNLIKELY(sb == 0) && k > 0)
 133             sb = xp[--k];
 134           /* rounding to nearest, with rounding bit = 1 */
 135           if (MPFR_UNLIKELY(sb == 0)) /* Even rounding. */
 136             {
 137               /* sb == 0 && rnd_mode == MPFR_RNDN */
 138               sb = xp[xsize - nw] & (himask ^ (himask << 1));
 139               if (sb == 0)
 140                 {
 141                   if (use_inexp)
 142                     *inexp = 2*MPFR_EVEN_INEX*neg-MPFR_EVEN_INEX;
 143                   /* ((neg!=0)^(sb!=0)) ? MPFR_EVEN_INEX  : -MPFR_EVEN_INEX;*/
 144                   /* Since neg = 0 or 1 and sb=0*/
 145 #if flag == 1
 146                   return 0 /*sb != 0 && rnd_mode != MPFR_RNDZ */;
 147 #else
 148                   MPN_COPY_INCR(yp, xp + xsize - nw, nw);
 149                   yp[0] &= himask;
 150                   return 0;
 151 #endif
 152                 }
 153               else
 154                 {
 155                   /* sb != 0 && rnd_mode == MPFR_RNDN */
 156                   if (use_inexp)
 157                     *inexp = MPFR_EVEN_INEX-2*MPFR_EVEN_INEX*neg;
 158                   /*((neg!=0)^(sb!=0))? MPFR_EVEN_INEX  : -MPFR_EVEN_INEX; */
 159                   /*Since neg= 0 or 1 and sb != 0 */
 160                   goto rnd_RNDN_add_one_ulp;
 161                 }
 162             }
 163           else /* sb != 0  && rnd_mode == MPFR_RNDN*/
 164             {
 165               if (use_inexp)
 166                 /* *inexp = (neg == 0) ? 1 : -1; but since neg = 0 or 1 */
 167                 *inexp = 1-2*neg;
 168             rnd_RNDN_add_one_ulp:
 169 #if flag == 1
 170               return 1; /*sb != 0 && rnd_mode != MPFR_RNDZ;*/
 171 #else
 172               carry = mpn_add_1 (yp, xp + xsize - nw, nw,
 173                                  rw ?
 174                                  MPFR_LIMB_ONE << (GMP_NUMB_BITS - rw)
 175                                  : MPFR_LIMB_ONE);
 176               yp[0] &= himask;
 177               return carry;
 178 #endif
 179             }
 180         }
 181       /* Rounding to Zero ? */
 182       else if (MPFR_IS_LIKE_RNDZ(rnd_mode, neg))
 183         {
 184           /* rnd_mode == MPFR_RNDZ */
 185         rnd_RNDZ:
 186           while (MPFR_UNLIKELY(sb == 0) && k > 0)
 187             sb = xp[--k];
 188           if (use_inexp)
 189             /* rnd_mode == MPFR_RNDZ and neg = 0 or 1 */
 190             /* (neg != 0) ^ (rnd_mode != MPFR_RNDZ)) ? 1 : -1);*/
 191             *inexp = MPFR_UNLIKELY(sb == 0) ? 0 : (2*neg-1);
 192 #if flag == 1
 193           return 0; /*sb != 0 && rnd_mode != MPFR_RNDZ;*/
 194 #else
 195           MPN_COPY_INCR(yp, xp + xsize - nw, nw);
 196           yp[0] &= himask;
 197           return 0;
 198 #endif
 199         }
 200       else
 201         {
 202           /* rnd_mode = Away */
 203           while (MPFR_UNLIKELY(sb == 0) && k > 0)
 204             sb = xp[--k];
 205           if (MPFR_UNLIKELY(sb == 0))
 206             {
 207               /* sb = 0 && rnd_mode != MPFR_RNDZ */
 208               if (use_inexp)
 209                 /* (neg != 0) ^ (rnd_mode != MPFR_RNDZ)) ? 1 : -1);*/
 210                 *inexp = 0;
 211 #if flag == 1
 212               return 0;
 213 #else
 214               MPN_COPY_INCR(yp, xp + xsize - nw, nw);
 215               yp[0] &= himask;
 216               return 0;
 217 #endif
 218             }
 219           else
 220             {
 221               /* sb != 0 && rnd_mode != MPFR_RNDZ */
 222               if (use_inexp)
 223                 /* (neg != 0) ^ (rnd_mode != MPFR_RNDZ)) ? 1 : -1);*/
 224                 *inexp = 1-2*neg;
 225 #if flag == 1
 226               return 1;
 227 #else
 228               carry = mpn_add_1(yp, xp + xsize - nw, nw,
 229                                 rw ? MPFR_LIMB_ONE << (GMP_NUMB_BITS - rw)
 230                                 : 1);
 231               yp[0] &= himask;
 232               return carry;
 233 #endif
 234             }
 235         }
 236     }
 237   else
 238     {
 239       /* Roundind mode = Zero / No inexact flag */
 240 #if flag == 1
 241       return 0 /*sb != 0 && rnd_mode != MPFR_RNDZ*/;
 242 #else
 243       if (MPFR_LIKELY(rw))
 244         {
 245           nw++;
 246           himask = ~MPFR_LIMB_MASK (GMP_NUMB_BITS - rw);
 247         }
 248       else
 249         himask = ~(mp_limb_t) 0;
 250       MPN_COPY_INCR(yp, xp + xsize - nw, nw);
 251       yp[0] &= himask;
 252       return 0;
 253 #endif
 254     }
 255 }
 256
 257 #undef flag
 258 #undef use_inexp
 259 #undef mpfr_round_raw_generic