1 /* mpfr_round_raw_generic -- Generic rounding function 2 3 Copyright 1999-2018 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramba 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 subnormal number with more precision. 51 * 52 * MPFR_RNDNA is now supported, but needs to be tested [TODO] and is 53 * still not part of the API. In particular, the MPFR_RNDNA value (-1) 54 * may change in the future without notice. 55 */ 56 57 int 58 mpfr_round_raw_generic( 59 #if flag == 0 60 mp_limb_t *yp, 61 #endif 62 const mp_limb_t *xp, mpfr_prec_t xprec, 63 int neg, mpfr_prec_t yprec, mpfr_rnd_t rnd_mode 64 #if use_inexp != 0 65 , int *inexp 66 #endif 67 ) 68 { 69 mp_size_t xsize, nw; 70 mp_limb_t himask, lomask, sb; 71 int rw, new_use_inexp; 72 #if flag == 0 73 int carry; 74 #endif 75 #if use_inexp == 0 76 int *inexp; 77 #endif 78 79 if (use_inexp) 80 MPFR_ASSERTD(inexp != ((int*) 0)); 81 MPFR_ASSERTD(neg == 0 || neg == 1); 82 83 if (rnd_mode == MPFR_RNDF) 84 { 85 if (use_inexp) 86 *inexp = 0; /* make sure it has a valid value */ 87 if (flag) 88 return 0; 89 rnd_mode = MPFR_RNDZ; /* faster */ 90 new_use_inexp = 0; 91 } 92 else 93 { 94 if (flag && !use_inexp && 95 (xprec <= yprec || MPFR_IS_LIKE_RNDZ (rnd_mode, neg))) 96 return 0; 97 new_use_inexp = use_inexp; 98 } 99 100 xsize = MPFR_PREC2LIMBS (xprec); 101 nw = yprec / GMP_NUMB_BITS; 102 rw = yprec & (GMP_NUMB_BITS - 1); 103 104 if (MPFR_UNLIKELY(xprec <= yprec)) 105 { /* No rounding is necessary. */ 106 /* if yp=xp, maybe an overlap: mpn_copyd is OK when src <= dst */ 107 if (MPFR_LIKELY(rw)) 108 nw++; 109 MPFR_ASSERTD(nw >= 1); 110 MPFR_ASSERTD(nw >= xsize); 111 if (use_inexp) 112 *inexp = 0; 113 #if flag == 0 114 mpn_copyd (yp + (nw - xsize), xp, xsize); 115 MPN_ZERO(yp, nw - xsize); 116 #endif 117 return 0; 118 } 119 120 if (new_use_inexp || !MPFR_IS_LIKE_RNDZ(rnd_mode, neg)) 121 { 122 mp_size_t k = xsize - nw - 1; 123 124 if (MPFR_LIKELY(rw)) 125 { 126 nw++; 127 lomask = MPFR_LIMB_MASK (GMP_NUMB_BITS - rw); 128 himask = ~lomask; 129 } 130 else 131 { 132 lomask = MPFR_LIMB_MAX; 133 himask = MPFR_LIMB_MAX; 134 } 135 MPFR_ASSERTD(k >= 0); 136 sb = xp[k] & lomask; /* First non-significant bits */ 137 /* Rounding to nearest? */ 138 if (rnd_mode == MPFR_RNDN || rnd_mode == MPFR_RNDNA) 139 { 140 /* Rounding to nearest */ 141 mp_limb_t rbmask = MPFR_LIMB_ONE << (GMP_NUMB_BITS - 1 - rw); 142 143 if ((sb & rbmask) == 0) /* rounding bit = 0 ? */ 144 goto rnd_RNDZ; /* yes, behave like rounding toward zero */ 145 /* Rounding to nearest with rounding bit = 1 */ 146 if (MPFR_UNLIKELY (rnd_mode == MPFR_RNDNA)) 147 goto away_addone_ulp; /* like rounding away from zero */ 148 sb &= ~rbmask; /* first bits after the rounding bit */ 149 while (MPFR_UNLIKELY(sb == 0) && k > 0) 150 sb = xp[--k]; 151 if (MPFR_UNLIKELY(sb == 0)) /* Even rounding. */ 152 { 153 /* sb == 0 && rnd_mode == MPFR_RNDN */ 154 sb = xp[xsize - nw] & (himask ^ (himask << 1)); 155 if (sb == 0) 156 { 157 if (use_inexp) 158 *inexp = 2*MPFR_EVEN_INEX*neg-MPFR_EVEN_INEX; 159 /* ((neg!=0)^(sb!=0)) ? MPFR_EVEN_INEX : -MPFR_EVEN_INEX */ 160 /* since neg = 0 or 1 and sb = 0 */ 161 #if flag == 0 162 mpn_copyi (yp, xp + xsize - nw, nw); 163 yp[0] &= himask; 164 #endif 165 return 0; /* sb != 0 && rnd_mode != MPFR_RNDZ */ 166 } 167 else 168 { 169 away_addone_ulp: 170 /* sb != 0 && rnd_mode == MPFR_RNDN */ 171 if (use_inexp) 172 *inexp = MPFR_EVEN_INEX-2*MPFR_EVEN_INEX*neg; 173 /* ((neg!=0)^(sb!=0)) ? MPFR_EVEN_INEX : -MPFR_EVEN_INEX */ 174 /* since neg = 0 or 1 and sb != 0 */ 175 goto rnd_RNDN_add_one_ulp; 176 } 177 } 178 else /* sb != 0 && rnd_mode == MPFR_RNDN */ 179 { 180 if (use_inexp) 181 *inexp = 1-2*neg; /* neg == 0 ? 1 : -1 */ 182 rnd_RNDN_add_one_ulp: 183 #if flag == 1 184 return 1; /* sb != 0 && rnd_mode != MPFR_RNDZ */ 185 #else 186 carry = mpn_add_1 (yp, xp + xsize - nw, nw, 187 rw ? 188 MPFR_LIMB_ONE << (GMP_NUMB_BITS - rw) 189 : MPFR_LIMB_ONE); 190 yp[0] &= himask; 191 return carry; 192 #endif 193 } 194 } 195 /* Rounding toward zero? */ 196 else if (MPFR_IS_LIKE_RNDZ(rnd_mode, neg)) 197 { 198 /* rnd_mode == MPFR_RNDZ */ 199 rnd_RNDZ: 200 while (MPFR_UNLIKELY(sb == 0) && k > 0) 201 sb = xp[--k]; 202 if (use_inexp) 203 /* rnd_mode == MPFR_RNDZ and neg = 0 or 1 */ 204 /* ((neg != 0) ^ (rnd_mode != MPFR_RNDZ)) ? 1 : -1 */ 205 *inexp = MPFR_UNLIKELY(sb == 0) ? 0 : (2*neg-1); 206 #if flag == 0 207 mpn_copyi (yp, xp + xsize - nw, nw); 208 yp[0] &= himask; 209 #endif 210 return 0; /* sb != 0 && rnd_mode != MPFR_RNDZ */ 211 } 212 else 213 { 214 /* Rounding away from zero */ 215 while (MPFR_UNLIKELY(sb == 0) && k > 0) 216 sb = xp[--k]; 217 if (MPFR_UNLIKELY(sb == 0)) 218 { 219 /* sb = 0 && rnd_mode != MPFR_RNDZ */ 220 if (use_inexp) 221 /* ((neg != 0) ^ (rnd_mode != MPFR_RNDZ)) ? 1 : -1 */ 222 *inexp = 0; 223 #if flag == 0 224 mpn_copyi (yp, xp + xsize - nw, nw); 225 yp[0] &= himask; 226 #endif 227 return 0; 228 } 229 else 230 { 231 /* sb != 0 && rnd_mode != MPFR_RNDZ */ 232 if (use_inexp) 233 *inexp = 1-2*neg; /* neg == 0 ? 1 : -1 */ 234 #if flag == 1 235 return 1; 236 #else 237 carry = mpn_add_1(yp, xp + xsize - nw, nw, 238 rw ? MPFR_LIMB_ONE << (GMP_NUMB_BITS - rw) 239 : 1); 240 yp[0] &= himask; 241 return carry; 242 #endif 243 } 244 } 245 } 246 else 247 { 248 /* Rounding toward zero / no inexact flag */ 249 #if flag == 0 250 if (MPFR_LIKELY(rw)) 251 { 252 nw++; 253 himask = ~MPFR_LIMB_MASK (GMP_NUMB_BITS - rw); 254 } 255 else 256 himask = MPFR_LIMB_MAX; 257 mpn_copyi (yp, xp + xsize - nw, nw); 258 yp[0] &= himask; 259 #endif 260 return 0; 261 } 262 } 263 264 #undef flag 265 #undef use_inexp 266 #undef mpfr_round_raw_generic 267