1/* mpfr_round_raw_generic -- Generic rounding function 2 3Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 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 53int 54mpfr_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