1 /* mpfr_round_near_x -- Round a floating point number nears another one.
3 Copyright 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
4 Contributed by the Arenaire and Cacao projects, INRIA.
6 This file is part of the GNU MPFR Library, and was contributed by Mathieu Dutour.
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 2.1 of the License, or (at your
11 option) any later version.
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.
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.LIB. If not, write to
20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
21 MA 02110-1301, USA. */
23 #include "mpfr-impl.h"
25 /* Use MPFR_FAST_COMPUTE_IF_SMALL_INPUT instead (a simple wrapper) */
27 /* int mpfr_round_near_x (mpfr_ptr y, mpfr_srcptr v, mpfr_uexp_t err, int dir,
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.
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
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 towards 0, if dir=1, it rounds away from 0),
42 rnd the rounding mode.
44 It returns 0 if it can't round.
45 Otherwise it returns the ternary flag (It can't return an exact value).
48 /* What "small enough" means?
50 We work with the positive values.
51 Assuming err > Prec (y)+1
53 i = [ y = o(x)] // i = inexact flag
55 Setting x in y is exact. We have:
56 y = [XXXXXXXXX[...]]0[...] + error where [..] are optional zeros
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
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
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
84 elif i < 0 and i = -EVEN_ROUNDING
85 So rnd = RNDN and we have y = x - ulp(y)/2
87 y < f(x) < y + ulp(y)/2 => do nothing
89 y + ulp(y)/2 < f(x) < y + ulp(y) => AddOneUlp
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
96 we have x < f(x) < x + 2^(EXP(x)-err)
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
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
107 Wrong. Example X=[0.11101]111111110000
108 + 1111111111111111111....
109 elif dirError = ToZero
110 we have x - 2^(EXP(x)-err) < f(x) < 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
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
120 Wrong. Example: X=[1.111010]00000010
121 - 10000001000000000000100....
123 y - ulp(y)/2 < x < y and we can't have x = y-ulp(y)/2:
125 y - ulp(y)/2 + ulp(x)/2 <= x <= y - ulp(x)/2
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.
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
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
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
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
156 mpfr_round_near_x (mpfr_ptr y, mpfr_srcptr v, mpfr_uexp_t err, int dir,
160 unsigned int old_flags = __gmpfr_flags;
162 MPFR_ASSERTD (!MPFR_IS_SINGULAR (v));
163 MPFR_ASSERTD (dir == 0 || dir == 1);
165 /* First check if we can round. The test is more restrictive than
166 necessary. Note that if err is not representable in an mp_exp_t,
167 then err > MPFR_PREC (v) and the conversion to mp_exp_t will not
169 if (!(err > MPFR_PREC (y) + 1
170 && (err > MPFR_PREC (v)
171 || mpfr_round_p (MPFR_MANT (v), MPFR_LIMB_SIZE (v),
173 MPFR_PREC (y) + (rnd == GMP_RNDN)))))
174 /* If we assume we can not round, return 0, and y is not modified */
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,
189 , if (MPFR_UNLIKELY (++MPFR_EXP (y) > __gmpfr_emax))
190 mpfr_overflow (y, rnd, sign)
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 */
199 if (dir == 0) /* The error term is negative for v positive */
202 if (MPFR_IS_LIKE_RNDZ (rnd, MPFR_IS_NEG_SIGN (sign)))
204 /* case nexttozero */
205 /* The underflow flag should be set if the result is zero */
206 __gmpfr_flags = old_flags;
209 if (MPFR_UNLIKELY (MPFR_IS_ZERO (y)))
210 mpfr_set_underflow ();
213 else /* The error term is positive for v positive */
217 if (rnd != GMP_RNDN && rnd != GMP_RNDZ
218 && MPFR_IS_RNDUTEST_OR_RNDDNOTTEST (rnd, MPFR_IS_POS_SIGN(sign)))
221 /* The overflow flag should be set if the result is infinity */
224 if (MPFR_UNLIKELY (MPFR_IS_INF (y)))
225 mpfr_set_overflow ();
230 /* the inexact flag cannot be 0, since this would mean an exact value,
231 and in this case we cannot round correctly */
232 MPFR_ASSERTD(inexact != 0);