1 /* mpfr_round_raw_generic -- Generic rounding function 2 3 Copyright 1999-2023 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 https://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 #if !(flag == 0 || flag == 1) 58 #error "flag must be 0 or 1" 59 #endif 60 61 int 62 mpfr_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
Indexes created Sun Mar 01 05:31:48 UTC 2026