1 /* mpfr_round_near_x -- Round a floating point number nears another one. 2 3 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramel 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 #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 unsigned int old_flags = __gmpfr_flags; 161 162 MPFR_ASSERTD (!MPFR_IS_SINGULAR (v)); 163 MPFR_ASSERTD (dir == 0 || dir == 1); 164 165 /* First check if we can round. The test is more restrictive than 166 necessary. Note that if err is not representable in an mpfr_exp_t, 167 then err > MPFR_PREC (v) and the conversion to mpfr_exp_t will not 168 occur. */ 169 if (!(err > MPFR_PREC (y) + 1 170 && (err > MPFR_PREC (v) 171 || mpfr_round_p (MPFR_MANT (v), MPFR_LIMB_SIZE (v), 172 (mpfr_exp_t) err, 173 MPFR_PREC (y) + (rnd == MPFR_RNDN))))) 174 /* If we assume we can not round, return 0, and y is not modified */ 175 return 0; 176 177 /* First round v in y */ 178 sign = MPFR_SIGN (v); 179 MPFR_SET_EXP (y, MPFR_GET_EXP (v)); 180 MPFR_SET_SIGN (y, sign); 181 MPFR_RNDRAW_GEN (inexact, y, MPFR_MANT (v), MPFR_PREC (v), rnd, sign, 182 if (dir == 0) 183 { 184 inexact = -sign; 185 goto trunc_doit; 186 } 187 else 188 goto addoneulp; 189 , if (MPFR_UNLIKELY (++MPFR_EXP (y) > __gmpfr_emax)) 190 mpfr_overflow (y, rnd, sign) 191 ); 192 193 /* Fix it in some cases */ 194 MPFR_ASSERTD (!MPFR_IS_NAN (y) && !MPFR_IS_ZERO (y)); 195 /* If inexact == 0, setting y from v is exact but we haven't 196 take into account yet the error term */ 197 if (inexact == 0) 198 { 199 if (dir == 0) /* The error term is negative for v positive */ 200 { 201 inexact = sign; 202 if (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG_SIGN (sign))) 203 { 204 /* case nexttozero */ 205 /* The underflow flag should be set if the result is zero */ 206 __gmpfr_flags = old_flags; 207 inexact = -sign; 208 mpfr_nexttozero (y); 209 if (MPFR_UNLIKELY (MPFR_IS_ZERO (y))) 210 mpfr_set_underflow (); 211 } 212 } 213 else /* The error term is positive for v positive */ 214 { 215 inexact = -sign; 216 /* Round Away */ 217 if (rnd != MPFR_RNDN && !MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG_SIGN(sign))) 218 { 219 /* case nexttoinf */ 220 /* The overflow flag should be set if the result is infinity */ 221 inexact = sign; 222 mpfr_nexttoinf (y); 223 if (MPFR_UNLIKELY (MPFR_IS_INF (y))) 224 mpfr_set_overflow (); 225 } 226 } 227 } 228 229 /* the inexact flag cannot be 0, since this would mean an exact value, 230 and in this case we cannot round correctly */ 231 MPFR_ASSERTD(inexact != 0); 232 MPFR_RET (inexact); 233 } 234