1 /* mpfr_rint -- Round to an integer.
3 Copyright 1999-2015 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.LESSER. If not, see
20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23 #include "mpfr-impl.h"
25 /* Merge the following mpfr_rint code with mpfr_round_raw_generic? */
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.
33 * We also need to be careful with the flags.
37 mpfr_rint (mpfr_ptr r
, mpfr_srcptr u
, mpfr_rnd_t rnd_mode
)
43 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u
) ))
50 MPFR_SET_SAME_SIGN(r
, u
);
54 MPFR_RET(0); /* infinity is exact */
56 else /* now u is zero */
58 MPFR_ASSERTD(MPFR_IS_ZERO(u
));
60 MPFR_RET(0); /* zero is exact */
63 MPFR_SET_SAME_SIGN (r
, u
); /* Does nothing if r==u */
65 sign
= MPFR_INT_SIGN (u
);
66 exp
= MPFR_GET_EXP (u
);
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) */
76 1 if round away from zero,
78 -1 if not decided yet.
81 if (MPFR_UNLIKELY (exp
<= 0)) /* 0 < |u| < 1 ==> round |u| to 0 or 1 */
83 /* Note: in the MPFR_RNDN mode, 0.5 must be rounded to 0. */
86 (exp
== 0 && (rnd_mode
== MPFR_RNDNA
||
87 !mpfr_powerof2_raw (u
)))))
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);
96 MPFR_SET_ZERO(r
); /* r = 0 */
97 MPFR_RET(sign
> 0 ? -2 : 2);
100 else /* exp > 0, |u| >= 1 */
103 mp_size_t un
, rn
, ui
;
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.
117 un
= MPFR_LIMB_SIZE(u
);
118 rn
= MPFR_LIMB_SIZE(r
);
119 MPFR_UNSIGNED_MINUS_MODULO (sh
, MPFR_PREC (r
));
121 /* exp is in the current exponent range: obtained from the input. */
122 MPFR_SET_EXP (r
, exp
); /* Does nothing if r==u */
124 if ((exp
- 1) / GMP_NUMB_BITS
>= un
)
128 uflags
= 0; /* u is an integer, representable or not in r */
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 */
139 uflags
= idiff
== 0 || (up
[uj
] << idiff
) == 0 ? 0 : 2;
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 */
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 */
163 /* a: rounding bit and some of the following bits */
164 /* b: boundary for a (weight of the rounding bit in a) */
167 a
= rp
[0] & ((MPFR_LIMB_ONE
<< sh
) - 1);
168 b
= MPFR_LIMB_ONE
<< (sh
- 1);
173 b
= MPFR_LIMB_HIGHBIT
;
179 for (i
= un
- rn
- 1 - (sh
== 0); i
>= 0; i
--)
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));
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 */
199 for (i
= un
- rn
- 1; i
>= 0; i
--)
202 uflags
= 1; /* u is not representable in r */
213 uj
= un
- ui
; /* lowest limb of the integer part in u */
214 rj
= rn
- ui
; /* lowest limb of the integer part in r */
216 if (MPFR_LIKELY (rp
!= up
))
217 MPN_COPY(rp
+ rj
, up
+ uj
, ui
);
219 /* Ignore the lowest rj limbs, all equal to zero. */
223 /* number of fractional bits in whole rp[0] */
224 ush
= idiff
== 0 ? 0 : GMP_NUMB_BITS
- idiff
;
226 if (rj
== 0 && ush
< sh
)
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 */
235 else /* The integer part of u fits in r, we'll round to it. */
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 */
248 /* a: rounding bit and some of the following bits */
249 /* b: boundary for a (weight of the rounding bit in a) */
252 a
= rp
[0] & ((MPFR_LIMB_ONE
<< sh
) - 1);
253 b
= MPFR_LIMB_ONE
<< (sh
- 1);
257 MPFR_ASSERTD (uj
>= 1); /* see above */
259 b
= MPFR_LIMB_HIGHBIT
;
265 for (i
= uj
- 1 - (sh
== 0); i
>= 0; i
--)
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));
279 /* Now we can make the low rj limbs to 0 */
280 MPN_ZERO (rp
-rj
, rj
);
284 rp
[0] &= MP_LIMB_T_MAX
<< sh
;
286 /* If u is a representable integer, there is no rounding. */
290 MPFR_ASSERTD (rnd_away
>= 0); /* rounding direction is defined */
291 if (rnd_away
&& mpn_add_1(rp
, rp
, rn
, MPFR_LIMB_ONE
<< sh
))
293 if (exp
== __gmpfr_emax
)
294 return mpfr_overflow (r
, rnd_mode
, sign
) >= 0 ?
296 else /* no overflow */
298 MPFR_SET_EXP(r
, exp
+ 1);
299 rp
[rn
-1] = MPFR_LIMB_HIGHBIT
;
303 MPFR_RET (rnd_away
^ (sign
< 0) ? uflags
: -uflags
);
304 } /* exp > 0, |u| >= 1 */
310 mpfr_round (mpfr_ptr r
, mpfr_srcptr u
)
312 return mpfr_rint (r
, u
, MPFR_RNDNA
);
318 mpfr_trunc (mpfr_ptr r
, mpfr_srcptr u
)
320 return mpfr_rint (r
, u
, MPFR_RNDZ
);
326 mpfr_ceil (mpfr_ptr r
, mpfr_srcptr u
)
328 return mpfr_rint (r
, u
, MPFR_RNDU
);
334 mpfr_floor (mpfr_ptr r
, mpfr_srcptr u
)
336 return mpfr_rint (r
, u
, MPFR_RNDD
);
339 /* We need to save the flags and restore them after calling the mpfr_round,
340 * mpfr_trunc, mpfr_ceil, mpfr_floor functions because these functions set
341 * the inexact flag when the argument is not an integer.
344 #undef mpfr_rint_round
347 mpfr_rint_round (mpfr_ptr r
, mpfr_srcptr u
, mpfr_rnd_t rnd_mode
)
349 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u
) ) || mpfr_integer_p (u
))
350 return mpfr_set (r
, u
, rnd_mode
);
355 unsigned int saved_flags
= __gmpfr_flags
;
356 MPFR_BLOCK_DECL (flags
);
358 mpfr_init2 (tmp
, MPFR_PREC (u
));
359 /* round(u) is representable in tmp unless an overflow occurs */
360 MPFR_BLOCK (flags
, mpfr_round (tmp
, u
));
361 __gmpfr_flags
= saved_flags
;
362 inex
= (MPFR_OVERFLOW (flags
)
363 ? mpfr_overflow (r
, rnd_mode
, MPFR_SIGN (u
))
364 : mpfr_set (r
, tmp
, rnd_mode
));
370 #undef mpfr_rint_trunc
373 mpfr_rint_trunc (mpfr_ptr r
, mpfr_srcptr u
, mpfr_rnd_t rnd_mode
)
375 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u
) ) || mpfr_integer_p (u
))
376 return mpfr_set (r
, u
, rnd_mode
);
381 unsigned int saved_flags
= __gmpfr_flags
;
383 mpfr_init2 (tmp
, MPFR_PREC (u
));
384 /* trunc(u) is always representable in tmp */
386 __gmpfr_flags
= saved_flags
;
387 inex
= mpfr_set (r
, tmp
, rnd_mode
);
393 #undef mpfr_rint_ceil
396 mpfr_rint_ceil (mpfr_ptr r
, mpfr_srcptr u
, mpfr_rnd_t rnd_mode
)
398 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u
) ) || mpfr_integer_p (u
))
399 return mpfr_set (r
, u
, rnd_mode
);
404 unsigned int saved_flags
= __gmpfr_flags
;
405 MPFR_BLOCK_DECL (flags
);
407 mpfr_init2 (tmp
, MPFR_PREC (u
));
408 /* ceil(u) is representable in tmp unless an overflow occurs */
409 MPFR_BLOCK (flags
, mpfr_ceil (tmp
, u
));
410 __gmpfr_flags
= saved_flags
;
411 inex
= (MPFR_OVERFLOW (flags
)
412 ? mpfr_overflow (r
, rnd_mode
, MPFR_SIGN_POS
)
413 : mpfr_set (r
, tmp
, rnd_mode
));
419 #undef mpfr_rint_floor
422 mpfr_rint_floor (mpfr_ptr r
, mpfr_srcptr u
, mpfr_rnd_t rnd_mode
)
424 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u
) ) || mpfr_integer_p (u
))
425 return mpfr_set (r
, u
, rnd_mode
);
430 unsigned int saved_flags
= __gmpfr_flags
;
431 MPFR_BLOCK_DECL (flags
);
433 mpfr_init2 (tmp
, MPFR_PREC (u
));
434 /* floor(u) is representable in tmp unless an overflow occurs */
435 MPFR_BLOCK (flags
, mpfr_floor (tmp
, u
));
436 __gmpfr_flags
= saved_flags
;
437 inex
= (MPFR_OVERFLOW (flags
)
438 ? mpfr_overflow (r
, rnd_mode
, MPFR_SIGN_NEG
)
439 : mpfr_set (r
, tmp
, rnd_mode
));