1 /* mpfr_round_near_x -- Round a floating point number nears another one. 2 3 Copyright 2005-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 #include "mpfr-impl.h" 24 25 /* Use MPFR_FAST_COMPUTE_IF_SMALL_INPUT instead (a simple wrapper) */ 26 27 /* int mpfr_round_near_x (mpfr_ptr y, mpfr_srcptr v, mpfr_uexp_t err, int dir, 28 mpfr_rnd_t rnd) 29 30 TODO: fix this description. 31 Assuming y = o(f(x)) = o(x + g(x)) with |g(x)| < 2^(EXP(v)-error) 32 If x is small enough, y ~= v. This function checks and does this. 33 34 It assumes that f(x) is not representable exactly as a FP number. 35 v must not be a singular value (NAN, INF or ZERO), usual values are 36 v=1 or v=x. 37 38 y is the destination (a mpfr_t), v the value to set (a mpfr_t), 39 err the error term (a mpfr_uexp_t) such that |g(x)| < 2^(EXP(x)-err), 40 dir (an int) is the direction of the error (if dir = 0, 41 it rounds toward 0, if dir=1, it rounds away from 0), 42 rnd the rounding mode. 43 44 It returns 0 if it can't round. 45 Otherwise it returns the ternary flag (It can't return an exact value). 46 */ 47 48 /* What "small enough" means? 49 50 We work with the positive values. 51 Assuming err > Prec (y)+1 52 53 i = [ y = o(x)] // i = inexact flag 54 If i == 0 55 Setting x in y is exact. We have: 56 y = [XXXXXXXXX[...]]0[...] + error where [..] are optional zeros 57 if dirError = ToInf, 58 x < f(x) < x + 2^(EXP(x)-err) 59 since x=y, and ulp (y)/2 > 2^(EXP(x)-err), we have: 60 y < f(x) < y+ulp(y) and |y-f(x)| < ulp(y)/2 61 if rnd = RNDN, nothing 62 if rnd = RNDZ, nothing 63 if rnd = RNDA, addoneulp 64 elif dirError = ToZero 65 x -2^(EXP(x)-err) < f(x) < x 66 since x=y, and ulp (y)/2 > 2^(EXP(x)-err), we have: 67 y-ulp(y) < f(x) < y and |y-f(x)| < ulp(y)/2 68 if rnd = RNDN, nothing 69 if rnd = RNDZ, nexttozero 70 if rnd = RNDA, nothing 71 NOTE: err > prec (y)+1 is needed only for RNDN. 72 elif i > 0 and i = EVEN_ROUNDING 73 So rnd = RNDN and we have y = x + ulp(y)/2 74 if dirError = ToZero, 75 we have x -2^(EXP(x)-err) < f(x) < x 76 so y - ulp(y)/2 - 2^(EXP(x)-err) < f(x) < y-ulp(y)/2 77 so y -ulp(y) < f(x) < y-ulp(y)/2 78 => nexttozero(y) 79 elif dirError = ToInf 80 we have x < f(x) < x + 2^(EXP(x)-err) 81 so y - ulp(y)/2 < f(x) < y+ulp(y)/2-ulp(y)/2 82 so y - ulp(y)/2 < f(x) < y 83 => do nothing 84 elif i < 0 and i = -EVEN_ROUNDING 85 So rnd = RNDN and we have y = x - ulp(y)/2 86 if dirError = ToZero, 87 y < f(x) < y + ulp(y)/2 => do nothing 88 if dirError = ToInf 89 y + ulp(y)/2 < f(x) < y + ulp(y) => AddOneUlp 90 elif i > 0 91 we can't have rnd = RNDZ, and prec(x) > prec(y), so ulp(x) < ulp(y) 92 we have y - ulp (y) < x < y 93 or more exactly y - ulp(y) + ulp(x)/2 <= x <= y - ulp(x)/2 94 if rnd = RNDA, 95 if dirError = ToInf, 96 we have x < f(x) < x + 2^(EXP(x)-err) 97 if err > prec (x), 98 we have 2^(EXP(x)-err) < ulp(x), so 2^(EXP(x)-err) <= ulp(x)/2 99 so f(x) <= y - ulp(x)/2+ulp(x)/2 <= y 100 and y - ulp(y) < x < f(x) 101 so we have y - ulp(y) < f(x) < y 102 so do nothing. 103 elif we can round, ie y - ulp(y) < x + 2^(EXP(x)-err) < y 104 we have y - ulp(y) < x < f(x) < x + 2^(EXP(x)-err) < y 105 so do nothing 106 otherwise 107 Wrong. Example X=[0.11101]111111110000 108 + 1111111111111111111.... 109 elif dirError = ToZero 110 we have x - 2^(EXP(x)-err) < f(x) < x 111 so f(x) < x < y 112 if err > prec (x) 113 x-2^(EXP(x)-err) >= x-ulp(x)/2 >= y - ulp(y) + ulp(x)/2-ulp(x)/2 114 so y - ulp(y) < f(x) < y 115 so do nothing 116 elif we can round, ie y - ulp(y) < x - 2^(EXP(x)-err) < y 117 y - ulp(y) < x - 2^(EXP(x)-err) < f(x) < y 118 so do nothing 119 otherwise 120 Wrong. Example: X=[1.111010]00000010 121 - 10000001000000000000100.... 122 elif rnd = RNDN, 123 y - ulp(y)/2 < x < y and we can't have x = y-ulp(y)/2: 124 so we have: 125 y - ulp(y)/2 + ulp(x)/2 <= x <= y - ulp(x)/2 126 if dirError = ToInf 127 we have x < f(x) < x+2^(EXP(x)-err) and ulp(y) > 2^(EXP(x)-err) 128 so y - ulp(y)/2 + ulp (x)/2 < f(x) < y + ulp (y)/2 - ulp (x)/2 129 we can round but we can't compute inexact flag. 130 if err > prec (x) 131 y - ulp(y)/2 + ulp (x)/2 < f(x) < y + ulp(x)/2 - ulp(x)/2 132 so y - ulp(y)/2 + ulp (x)/2 < f(x) < y 133 we can round and compute inexact flag. do nothing 134 elif we can round, ie y - ulp(y)/2 < x + 2^(EXP(x)-err) < y 135 we have y - ulp(y)/2 + ulp (x)/2 < f(x) < y 136 so do nothing 137 otherwise 138 Wrong 139 elif dirError = ToZero 140 we have x -2^(EXP(x)-err) < f(x) < x and ulp(y)/2 > 2^(EXP(x)-err) 141 so y-ulp(y)+ulp(x)/2 < f(x) < y - ulp(x)/2 142 if err > prec (x) 143 x- ulp(x)/2 < f(x) < x 144 so y - ulp(y)/2+ulp(x)/2 - ulp(x)/2 < f(x) < x <= y - ulp(x)/2 < y 145 do nothing 146 elif we can round, ie y-ulp(y)/2 < x-2^(EXP(x)-err) < y 147 we have y-ulp(y)/2 < x-2^(EXP(x)-err) < f(x) < x < y 148 do nothing 149 otherwise 150 Wrong 151 elif i < 0 152 same thing? 153 */ 154 155 int 156 mpfr_round_near_x (mpfr_ptr y, mpfr_srcptr v, mpfr_uexp_t err, int dir, 157 mpfr_rnd_t rnd) 158 { 159 int inexact, sign; 160 mpfr_flags_t old_flags = __gmpfr_flags; 161 162 if (rnd == MPFR_RNDF) 163 rnd = MPFR_RNDZ; 164 165 MPFR_ASSERTD (!MPFR_IS_SINGULAR (v)); 166 MPFR_ASSERTD (dir == 0 || dir == 1); 167 168 /* First check if we can round. The test is more restrictive than 169 necessary. Note that if err is not representable in an mpfr_exp_t, 170 then err > MPFR_PREC (v) and the conversion to mpfr_exp_t will not 171 occur. */ 172 if (!(err > MPFR_PREC (y) + 1 173 && (err > MPFR_PREC (v) 174 || mpfr_round_p (MPFR_MANT (v), MPFR_LIMB_SIZE (v), 175 (mpfr_exp_t) err, 176 MPFR_PREC (y) + (rnd == MPFR_RNDN))))) 177 /* If we assume we can not round, return 0, and y is not modified */ 178 return 0; 179 180 /* First round v in y */ 181 sign = MPFR_SIGN (v); 182 MPFR_SET_EXP (y, MPFR_GET_EXP (v)); 183 MPFR_SET_SIGN (y, sign); 184 MPFR_RNDRAW_GEN (inexact, y, MPFR_MANT (v), MPFR_PREC (v), rnd, sign, 185 if (dir == 0) 186 { 187 inexact = -sign; 188 goto trunc_doit; 189 } 190 else 191 goto addoneulp; 192 , if (MPFR_UNLIKELY (++MPFR_EXP (y) > __gmpfr_emax)) 193 mpfr_overflow (y, rnd, sign) 194 ); 195 196 /* Fix it in some cases */ 197 MPFR_ASSERTD (!MPFR_IS_NAN (y) && !MPFR_IS_ZERO (y)); 198 /* If inexact == 0, setting y from v is exact but we haven't 199 take into account yet the error term */ 200 if (inexact == 0) 201 { 202 if (dir == 0) /* The error term is negative for v positive */ 203 { 204 inexact = sign; 205 if (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG_SIGN (sign))) 206 { 207 /* case nexttozero */ 208 /* The underflow flag should be set if the result is zero */ 209 __gmpfr_flags = old_flags; 210 inexact = -sign; 211 mpfr_nexttozero (y); 212 if (MPFR_UNLIKELY (MPFR_IS_ZERO (y))) 213 MPFR_SET_UNDERFLOW (); 214 } 215 } 216 else /* The error term is positive for v positive */ 217 { 218 inexact = -sign; 219 /* Round Away */ 220 if (MPFR_IS_LIKE_RNDA (rnd, MPFR_IS_NEG_SIGN(sign))) 221 { 222 /* case nexttoinf */ 223 /* The overflow flag should be set if the result is infinity */ 224 inexact = sign; 225 mpfr_nexttoinf (y); 226 if (MPFR_UNLIKELY (MPFR_IS_INF (y))) 227 MPFR_SET_OVERFLOW (); 228 } 229 } 230 } 231 232 /* the inexact flag cannot be 0, since this would mean an exact value, 233 and in this case we cannot round correctly */ 234 MPFR_ASSERTD(inexact != 0); 235 MPFR_RET (inexact); 236 } 237
Indexes created Fri Jun 05 00:26:10 UTC 2026