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 }