MPFR v4.0.2
[mpfr.git] / src / rint.c
1 /* mpfr_rint -- Round to an integer.
2
3 Copyright 1999-2019 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 || (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 && 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 }