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