1 /* mpfr_rint -- Round to an integer.
2
3 Copyright 1999-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 /* Merge the following mpfr_rint code with mpfr_round_raw_generic? */
26
27 /* For all the round-to-integer functions, we don't need to extend the
28 * exponent range. And it is better not to do so, so that we can test
29 * the flag setting for intermediate overflow in the test suite without
30 * involving huge non-integer numbers (thus in huge precision). This
31 * should also be faster.
32 *
33 * We also need to be careful with the flags.
34 */
35
36 int
37 mpfr_rint (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
38 {
39 int sign;
40 int rnd_away;
41 mpfr_exp_t exp;
42
43 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ))
44 {
45 if (MPFR_IS_NAN(u))
46 {
47 MPFR_SET_NAN(r);
48 MPFR_RET_NAN;
49 }
50 MPFR_SET_SAME_SIGN(r, u);
51 if (MPFR_IS_INF(u))
52 {
53 MPFR_SET_INF(r);
54 MPFR_RET(0); /* infinity is exact */
55 }
56 else /* now u is zero */
57 {
58 MPFR_ASSERTD(MPFR_IS_ZERO(u));
59 MPFR_SET_ZERO(r);
60 MPFR_RET(0); /* zero is exact */
61 }
62 }
63 MPFR_SET_SAME_SIGN (r, u); /* Does nothing if r==u */
64
65 sign = MPFR_INT_SIGN (u);
66 exp = MPFR_GET_EXP (u);
67
68 rnd_away =
69 rnd_mode == MPFR_RNDD ? sign < 0 :
70 rnd_mode == MPFR_RNDU ? sign > 0 :
71 rnd_mode == MPFR_RNDZ ? 0 :
72 rnd_mode == MPFR_RNDA ? 1 :
73 -1; /* round to nearest-even (RNDN) or nearest-away (RNDNA) */
74
75 /* rnd_away:
76 1 if round away from zero,
77 0 if round to zero,
78 -1 if not decided yet.
79 */
80
81 if (MPFR_UNLIKELY (exp <= 0)) /* 0 < |u| < 1 ==> round |u| to 0 or 1 */
82 {
83 /* Note: in the MPFR_RNDN mode, 0.5 must be rounded to 0. */
84 if (rnd_away != 0 &&
85 (rnd_away > 0 ||
86 (exp == 0 && (rnd_mode == MPFR_RNDNA ||
87 !mpfr_powerof2_raw (u)))))
88 {
89 /* The flags will correctly be set and overflow will correctly
90 be handled by mpfr_set_si. */
91 mpfr_set_si (r, sign, rnd_mode);
92 MPFR_RET(sign > 0 ? 2 : -2);
93 }
94 else
95 {
96 MPFR_SET_ZERO(r); /* r = 0 */
97 MPFR_RET(sign > 0 ? -2 : 2);
98 }
99 }
100 else /* exp > 0, |u| >= 1 */
101 {
102 mp_limb_t *up, *rp;
103 mp_size_t un, rn, ui;
104 int sh, idiff;
105 int uflags;
106
107 /*
108 * uflags will contain:
109 * _ 0 if u is an integer representable in r,
110 * _ 1 if u is an integer not representable in r,
111 * _ 2 if u is not an integer.
112 */
113
114 up = MPFR_MANT(u);
115 rp = MPFR_MANT(r);
116
117 un = MPFR_LIMB_SIZE(u);
118 rn = MPFR_LIMB_SIZE(r);
119 MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (r));
120
121 /* exp is in the current exponent range: obtained from the input. */
122 MPFR_SET_EXP (r, exp); /* Does nothing if r==u */
123
124 if ((exp - 1) / GMP_NUMB_BITS >= un)
125 {
126 ui = un;
127 idiff = 0;
128 uflags = 0; /* u is an integer, representable or not in r */
129 }
130 else
131 {
132 mp_size_t uj;
133
134 ui = (exp - 1) / GMP_NUMB_BITS + 1; /* #limbs of the int part */
135 MPFR_ASSERTD (un >= ui);
136 uj = un - ui; /* lowest limb of the integer part */
137 idiff = exp % GMP_NUMB_BITS; /* #int-part bits in up[uj] or 0 */
138
139 uflags = idiff == 0 || MPFR_LIMB_LSHIFT(up[uj],idiff) == 0 ? 0 : 2;
140 if (uflags == 0)
141 while (uj > 0)
142 if (up[--uj] != 0)
143 {
144 uflags = 2;
145 break;
146 }
147 }
148
149 if (ui > rn)
150 {
151 /* More limbs in the integer part of u than in r.
152 Just round u with the precision of r. */
153 MPFR_ASSERTD (rp != up && un > rn);
154 MPN_COPY (rp, up + (un - rn), rn); /* r != u */
155 if (rnd_away < 0)
156 {
157 /* This is a rounding to nearest mode (MPFR_RNDN or MPFR_RNDNA).
158 Decide the rounding direction here. */
159 if (rnd_mode == MPFR_RNDN &&
160 (rp[0] & (MPFR_LIMB_ONE << sh)) == 0)
161 { /* halfway cases rounded toward zero */
162 mp_limb_t a, b;
163 /* a: rounding bit and some of the following bits */
164 /* b: boundary for a (weight of the rounding bit in a) */
165 if (sh != 0)
166 {
167 a = rp[0] & ((MPFR_LIMB_ONE << sh) - 1);
168 b = MPFR_LIMB_ONE << (sh - 1);
169 }
170 else
171 {
172 a = up[un - rn - 1];
173 b = MPFR_LIMB_HIGHBIT;
174 }
175 rnd_away = a > b;
176 if (a == b)
177 {
178 mp_size_t i;
179 for (i = un - rn - 1 - (sh == 0); i >= 0; i--)
180 if (up[i] != 0)
181 {
182 rnd_away = 1;
183 break;
184 }
185 }
186 }
187 else /* halfway cases rounded away from zero */
188 rnd_away = /* rounding bit */
189 ((sh != 0 && (rp[0] & (MPFR_LIMB_ONE << (sh - 1))) != 0) ||
190 (sh == 0 && (up[un - rn - 1] & MPFR_LIMB_HIGHBIT) != 0));
191 }
192 if (uflags == 0)
193 { /* u is an integer; determine if it is representable in r */
194 if (sh != 0 && MPFR_LIMB_LSHIFT(rp[0], GMP_NUMB_BITS - sh) != 0)
195 uflags = 1; /* u is not representable in r */
196 else
197 {
198 mp_size_t i;
199 for (i = un - rn - 1; i >= 0; i--)
200 if (up[i] != 0)
201 {
202 uflags = 1; /* u is not representable in r */
203 break;
204 }
205 }
206 }
207 }
208 else /* ui <= rn */
209 {
210 mp_size_t uj, rj;
211 int ush;
212
213 uj = un - ui; /* lowest limb of the integer part in u */
214 rj = rn - ui; /* lowest limb of the integer part in r */
215
216 if (rp != up)
217 MPN_COPY(rp + rj, up + uj, ui);
218
219 /* Ignore the lowest rj limbs, all equal to zero. */
220 rp += rj;
221 rn = ui;
222
223 /* number of fractional bits in whole rp[0] */
224 ush = idiff == 0 ? 0 : GMP_NUMB_BITS - idiff;
225
226 if (rj == 0 && ush < sh)
227 {
228 /* If u is an integer (uflags == 0), we need to determine
229 if it is representable in r, i.e. if its sh - ush bits
230 in the non-significant part of r are all 0. */
231 if (uflags == 0 && (rp[0] & ((MPFR_LIMB_ONE << sh) -
232 (MPFR_LIMB_ONE << ush))) != 0)
233 uflags = 1; /* u is an integer not representable in r */
234 }
235 else /* The integer part of u fits in r, we'll round to it. */
236 sh = ush;
237
238 if (rnd_away < 0)
239 {
240 /* This is a rounding to nearest mode.
241 Decide the rounding direction here. */
242 if (uj == 0 && sh == 0)
243 rnd_away = 0; /* rounding bit = 0 (not represented in u) */
244 else if (rnd_mode == MPFR_RNDN &&
245 (rp[0] & (MPFR_LIMB_ONE << sh)) == 0)
246 { /* halfway cases rounded toward zero */
247 mp_limb_t a, b;
248 /* a: rounding bit and some of the following bits */
249 /* b: boundary for a (weight of the rounding bit in a) */
250 if (sh != 0)
251 {
252 a = rp[0] & ((MPFR_LIMB_ONE << sh) - 1);
253 b = MPFR_LIMB_ONE << (sh - 1);
254 }
255 else
256 {
257 MPFR_ASSERTD (uj >= 1); /* see above */
258 a = up[uj - 1];
259 b = MPFR_LIMB_HIGHBIT;
260 }
261 rnd_away = a > b;
262 if (a == b)
263 {
264 mp_size_t i;
265 for (i = uj - 1 - (sh == 0); i >= 0; i--)
266 if (up[i] != 0)
267 {
268 rnd_away = 1;
269 break;
270 }
271 }
272 }
273 else /* halfway cases rounded away from zero */
274 rnd_away = /* rounding bit */
275 ((sh != 0 && (rp[0] & (MPFR_LIMB_ONE << (sh - 1))) != 0) ||
276 (sh == 0 && (MPFR_ASSERTD (uj >= 1),
277 up[uj - 1] & MPFR_LIMB_HIGHBIT) != 0));
278 }
279 /* Now we can make the low rj limbs to 0 */
280 MPN_ZERO (rp-rj, rj);
281 }
282
283 if (sh != 0)
284 rp[0] &= MPFR_LIMB_MAX << sh;
285
286 /* If u is a representable integer, there is no rounding. */
287 if (uflags == 0)
288 MPFR_RET(0);
289
290 MPFR_ASSERTD (rnd_away >= 0); /* rounding direction is defined */
291 if (rnd_away && mpn_add_1 (rp, rp, rn, MPFR_LIMB_ONE << sh))
292 {
293 if (exp == __gmpfr_emax)
294 return mpfr_overflow (r, rnd_mode, sign) >= 0 ?
295 uflags : -uflags;
296 else /* no overflow */
297 {
298 MPFR_SET_EXP(r, exp + 1);
299 rp[rn-1] = MPFR_LIMB_HIGHBIT;
300 }
301 }
302
303 MPFR_RET (rnd_away ^ (sign < 0) ? uflags : -uflags);
304 } /* exp > 0, |u| >= 1 */
305 }
306
307 #undef mpfr_roundeven
308
309 int
310 mpfr_roundeven (mpfr_ptr r, mpfr_srcptr u)
311 {
312 return mpfr_rint (r, u, MPFR_RNDN);
313 }
314
315 #undef mpfr_round
316
317 int
318 mpfr_round (mpfr_ptr r, mpfr_srcptr u)
319 {
320 return mpfr_rint (r, u, MPFR_RNDNA);
321 }
322
323 #undef mpfr_trunc
324
325 int
326 mpfr_trunc (mpfr_ptr r, mpfr_srcptr u)
327 {
328 return mpfr_rint (r, u, MPFR_RNDZ);
329 }
330
331 #undef mpfr_ceil
332
333 int
334 mpfr_ceil (mpfr_ptr r, mpfr_srcptr u)
335 {
336 return mpfr_rint (r, u, MPFR_RNDU);
337 }
338
339 #undef mpfr_floor
340
341 int
342 mpfr_floor (mpfr_ptr r, mpfr_srcptr u)
343 {
344 return mpfr_rint (r, u, MPFR_RNDD);
345 }
346
347 /* We need to save the flags and restore them after calling the mpfr_round,
348 * mpfr_trunc, mpfr_ceil, mpfr_floor functions because these functions set
349 * the inexact flag when the argument is not an integer.
350 */
351
352 #undef mpfr_rint_roundeven
353
354 int
355 mpfr_rint_roundeven (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
356 {
357 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
358 return mpfr_set (r, u, rnd_mode);
359 else
360 {
361 mpfr_t tmp;
362 int inex;
363 mpfr_flags_t saved_flags = __gmpfr_flags;
364 MPFR_BLOCK_DECL (flags);
365
366 mpfr_init2 (tmp, MPFR_PREC (u));
367 /* round(u) is representable in tmp unless an overflow occurs */
368 MPFR_BLOCK (flags, mpfr_roundeven (tmp, u));
369 __gmpfr_flags = saved_flags;
370 inex = (MPFR_OVERFLOW (flags)
371 ? mpfr_overflow (r, rnd_mode, MPFR_SIGN (u))
372 : mpfr_set (r, tmp, rnd_mode));
373 mpfr_clear (tmp);
374 return inex;
375 }
376 }
377
378 #undef mpfr_rint_round
379
380 int
381 mpfr_rint_round (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
382 {
383 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
384 return mpfr_set (r, u, rnd_mode);
385 else
386 {
387 mpfr_t tmp;
388 int inex;
389 mpfr_flags_t saved_flags = __gmpfr_flags;
390 MPFR_BLOCK_DECL (flags);
391
392 mpfr_init2 (tmp, MPFR_PREC (u));
393 /* round(u) is representable in tmp unless an overflow occurs */
394 MPFR_BLOCK (flags, mpfr_round (tmp, u));
395 __gmpfr_flags = saved_flags;
396 inex = (MPFR_OVERFLOW (flags)
397 ? mpfr_overflow (r, rnd_mode, MPFR_SIGN (u))
398 : mpfr_set (r, tmp, rnd_mode));
399 mpfr_clear (tmp);
400 return inex;
401 }
402 }
403
404 #undef mpfr_rint_trunc
405
406 int
407 mpfr_rint_trunc (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
408 {
409 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
410 return mpfr_set (r, u, rnd_mode);
411 else
412 {
413 mpfr_t tmp;
414 int inex;
415 mpfr_flags_t saved_flags = __gmpfr_flags;
416
417 mpfr_init2 (tmp, MPFR_PREC (u));
418 /* trunc(u) is always representable in tmp */
419 mpfr_trunc (tmp, u);
420 __gmpfr_flags = saved_flags;
421 inex = mpfr_set (r, tmp, rnd_mode);
422 mpfr_clear (tmp);
423 return inex;
424 }
425 }
426
427 #undef mpfr_rint_ceil
428
429 int
430 mpfr_rint_ceil (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
431 {
432 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
433 return mpfr_set (r, u, rnd_mode);
434 else
435 {
436 mpfr_t tmp;
437 int inex;
438 mpfr_flags_t saved_flags = __gmpfr_flags;
439 MPFR_BLOCK_DECL (flags);
440
441 mpfr_init2 (tmp, MPFR_PREC (u));
442 /* ceil(u) is representable in tmp unless an overflow occurs */
443 MPFR_BLOCK (flags, mpfr_ceil (tmp, u));
444 __gmpfr_flags = saved_flags;
445 inex = (MPFR_OVERFLOW (flags)
446 ? mpfr_overflow (r, rnd_mode, MPFR_SIGN_POS)
447 : mpfr_set (r, tmp, rnd_mode));
448 mpfr_clear (tmp);
449 return inex;
450 }
451 }
452
453 #undef mpfr_rint_floor
454
455 int
456 mpfr_rint_floor (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
457 {
458 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
459 return mpfr_set (r, u, rnd_mode);
460 else
461 {
462 mpfr_t tmp;
463 int inex;
464 mpfr_flags_t saved_flags = __gmpfr_flags;
465 MPFR_BLOCK_DECL (flags);
466
467 mpfr_init2 (tmp, MPFR_PREC (u));
468 /* floor(u) is representable in tmp unless an overflow occurs */
469 MPFR_BLOCK (flags, mpfr_floor (tmp, u));
470 __gmpfr_flags = saved_flags;
471 inex = (MPFR_OVERFLOW (flags)
472 ? mpfr_overflow (r, rnd_mode, MPFR_SIGN_NEG)
473 : mpfr_set (r, tmp, rnd_mode));
474 mpfr_clear (tmp);
475 return inex;
476 }
477 }