1 /* mpfr_round_raw_generic -- Generic rounding function 2 3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4 Contributed by the AriC 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 = MPFR_PREC2LIMBS (xprec); 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 260