1/* mpfr_round_near_x -- Round a floating point number nears another one. 2 3Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc. 4Contributed by the AriC and Caramel 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 20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 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 155int 156mpfr_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