1 /* mpfr_rint -- Round to an integer.
3 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 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? */
28 mpfr_rint (mpfr_ptr r
, mpfr_srcptr u
, mpfr_rnd_t rnd_mode
)
34 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u
) ))
41 MPFR_SET_SAME_SIGN(r
, u
);
45 MPFR_RET(0); /* infinity is exact */
47 else /* now u is zero */
49 MPFR_ASSERTD(MPFR_IS_ZERO(u
));
51 MPFR_RET(0); /* zero is exact */
54 MPFR_SET_SAME_SIGN (r
, u
); /* Does nothing if r==u */
56 sign
= MPFR_INT_SIGN (u
);
57 exp
= MPFR_GET_EXP (u
);
60 rnd_mode
== MPFR_RNDD
? sign
< 0 :
61 rnd_mode
== MPFR_RNDU
? sign
> 0 :
62 rnd_mode
== MPFR_RNDZ
? 0 :
63 rnd_mode
== MPFR_RNDA
? 1 :
64 -1; /* round to nearest-even (RNDN) or nearest-away (RNDNA) */
67 1 if round away from zero,
69 -1 if not decided yet.
72 if (MPFR_UNLIKELY (exp
<= 0)) /* 0 < |u| < 1 ==> round |u| to 0 or 1 */
74 /* Note: in the MPFR_RNDN mode, 0.5 must be rounded to 0. */
77 (exp
== 0 && (rnd_mode
== MPFR_RNDNA
||
78 !mpfr_powerof2_raw (u
)))))
84 rm
= (MPFR_PREC(r
) - 1) / GMP_NUMB_BITS
;
85 rp
[rm
] = MPFR_LIMB_HIGHBIT
;
87 MPFR_SET_EXP (r
, 1); /* |r| = 1 */
88 MPFR_RET(sign
> 0 ? 2 : -2);
92 MPFR_SET_ZERO(r
); /* r = 0 */
93 MPFR_RET(sign
> 0 ? -2 : 2);
96 else /* exp > 0, |u| >= 1 */
104 * uflags will contain:
105 * _ 0 if u is an integer representable in r,
106 * _ 1 if u is an integer not representable in r,
107 * _ 2 if u is not an integer.
113 un
= MPFR_LIMB_SIZE(u
);
114 rn
= MPFR_LIMB_SIZE(r
);
115 MPFR_UNSIGNED_MINUS_MODULO (sh
, MPFR_PREC (r
));
117 MPFR_SET_EXP (r
, exp
); /* Does nothing if r==u */
119 if ((exp
- 1) / GMP_NUMB_BITS
>= un
)
123 uflags
= 0; /* u is an integer, representable or not in r */
129 ui
= (exp
- 1) / GMP_NUMB_BITS
+ 1; /* #limbs of the int part */
130 MPFR_ASSERTD (un
>= ui
);
131 uj
= un
- ui
; /* lowest limb of the integer part */
132 idiff
= exp
% GMP_NUMB_BITS
; /* #int-part bits in up[uj] or 0 */
134 uflags
= idiff
== 0 || (up
[uj
] << idiff
) == 0 ? 0 : 2;
146 /* More limbs in the integer part of u than in r.
147 Just round u with the precision of r. */
148 MPFR_ASSERTD (rp
!= up
&& un
> rn
);
149 MPN_COPY (rp
, up
+ (un
- rn
), rn
); /* r != u */
152 /* This is a rounding to nearest mode (MPFR_RNDN or MPFR_RNDNA).
153 Decide the rounding direction here. */
154 if (rnd_mode
== MPFR_RNDN
&&
155 (rp
[0] & (MPFR_LIMB_ONE
<< sh
)) == 0)
156 { /* halfway cases rounded toward zero */
158 /* a: rounding bit and some of the following bits */
159 /* b: boundary for a (weight of the rounding bit in a) */
162 a
= rp
[0] & ((MPFR_LIMB_ONE
<< sh
) - 1);
163 b
= MPFR_LIMB_ONE
<< (sh
- 1);
168 b
= MPFR_LIMB_HIGHBIT
;
174 for (i
= un
- rn
- 1 - (sh
== 0); i
>= 0; i
--)
182 else /* halfway cases rounded away from zero */
183 rnd_away
= /* rounding bit */
184 ((sh
!= 0 && (rp
[0] & (MPFR_LIMB_ONE
<< (sh
- 1))) != 0) ||
185 (sh
== 0 && (up
[un
- rn
- 1] & MPFR_LIMB_HIGHBIT
) != 0));
188 { /* u is an integer; determine if it is representable in r */
189 if (sh
!= 0 && rp
[0] << (GMP_NUMB_BITS
- sh
) != 0)
190 uflags
= 1; /* u is not representable in r */
194 for (i
= un
- rn
- 1; i
>= 0; i
--)
197 uflags
= 1; /* u is not representable in r */
208 uj
= un
- ui
; /* lowest limb of the integer part in u */
209 rj
= rn
- ui
; /* lowest limb of the integer part in r */
211 if (MPFR_LIKELY (rp
!= up
))
212 MPN_COPY(rp
+ rj
, up
+ uj
, ui
);
214 /* Ignore the lowest rj limbs, all equal to zero. */
218 /* number of fractional bits in whole rp[0] */
219 ush
= idiff
== 0 ? 0 : GMP_NUMB_BITS
- idiff
;
221 if (rj
== 0 && ush
< sh
)
223 /* If u is an integer (uflags == 0), we need to determine
224 if it is representable in r, i.e. if its sh - ush bits
225 in the non-significant part of r are all 0. */
226 if (uflags
== 0 && (rp
[0] & ((MPFR_LIMB_ONE
<< sh
) -
227 (MPFR_LIMB_ONE
<< ush
))) != 0)
228 uflags
= 1; /* u is an integer not representable in r */
230 else /* The integer part of u fits in r, we'll round to it. */
235 /* This is a rounding to nearest mode.
236 Decide the rounding direction here. */
237 if (uj
== 0 && sh
== 0)
238 rnd_away
= 0; /* rounding bit = 0 (not represented in u) */
239 else if (rnd_mode
== MPFR_RNDN
&&
240 (rp
[0] & (MPFR_LIMB_ONE
<< sh
)) == 0)
241 { /* halfway cases rounded toward zero */
243 /* a: rounding bit and some of the following bits */
244 /* b: boundary for a (weight of the rounding bit in a) */
247 a
= rp
[0] & ((MPFR_LIMB_ONE
<< sh
) - 1);
248 b
= MPFR_LIMB_ONE
<< (sh
- 1);
252 MPFR_ASSERTD (uj
>= 1); /* see above */
254 b
= MPFR_LIMB_HIGHBIT
;
260 for (i
= uj
- 1 - (sh
== 0); i
>= 0; i
--)
268 else /* halfway cases rounded away from zero */
269 rnd_away
= /* rounding bit */
270 ((sh
!= 0 && (rp
[0] & (MPFR_LIMB_ONE
<< (sh
- 1))) != 0) ||
271 (sh
== 0 && (MPFR_ASSERTD (uj
>= 1),
272 up
[uj
- 1] & MPFR_LIMB_HIGHBIT
) != 0));
274 /* Now we can make the low rj limbs to 0 */
275 MPN_ZERO (rp
-rj
, rj
);
279 rp
[0] &= MP_LIMB_T_MAX
<< sh
;
281 /* If u is a representable integer, there is no rounding. */
285 MPFR_ASSERTD (rnd_away
>= 0); /* rounding direction is defined */
286 if (rnd_away
&& mpn_add_1(rp
, rp
, rn
, MPFR_LIMB_ONE
<< sh
))
288 if (exp
== __gmpfr_emax
)
289 return mpfr_overflow(r
, rnd_mode
, MPFR_SIGN(r
)) >= 0 ?
293 MPFR_SET_EXP(r
, exp
+ 1);
294 rp
[rn
-1] = MPFR_LIMB_HIGHBIT
;
298 MPFR_RET (rnd_away
^ (sign
< 0) ? uflags
: -uflags
);
299 } /* exp > 0, |u| >= 1 */
305 mpfr_round (mpfr_ptr r
, mpfr_srcptr u
)
307 return mpfr_rint (r
, u
, MPFR_RNDNA
);
313 mpfr_trunc (mpfr_ptr r
, mpfr_srcptr u
)
315 return mpfr_rint (r
, u
, MPFR_RNDZ
);
321 mpfr_ceil (mpfr_ptr r
, mpfr_srcptr u
)
323 return mpfr_rint (r
, u
, MPFR_RNDU
);
329 mpfr_floor (mpfr_ptr r
, mpfr_srcptr u
)
331 return mpfr_rint (r
, u
, MPFR_RNDD
);
334 #undef mpfr_rint_round
337 mpfr_rint_round (mpfr_ptr r
, mpfr_srcptr u
, mpfr_rnd_t rnd_mode
)
339 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u
) ) || mpfr_integer_p (u
))
340 return mpfr_set (r
, u
, rnd_mode
);
345 MPFR_SAVE_EXPO_DECL (expo
);
346 MPFR_BLOCK_DECL (flags
);
348 MPFR_SAVE_EXPO_MARK (expo
);
349 mpfr_init2 (tmp
, MPFR_PREC (u
));
350 /* round(u) is representable in tmp unless an overflow occurs */
351 MPFR_BLOCK (flags
, mpfr_round (tmp
, u
));
352 inex
= (MPFR_OVERFLOW (flags
)
353 ? mpfr_overflow (r
, rnd_mode
, MPFR_SIGN (u
))
354 : mpfr_set (r
, tmp
, rnd_mode
));
356 MPFR_SAVE_EXPO_FREE (expo
);
357 return mpfr_check_range (r
, inex
, rnd_mode
);
361 #undef mpfr_rint_trunc
364 mpfr_rint_trunc (mpfr_ptr r
, mpfr_srcptr u
, mpfr_rnd_t rnd_mode
)
366 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u
) ) || mpfr_integer_p (u
))
367 return mpfr_set (r
, u
, rnd_mode
);
372 MPFR_SAVE_EXPO_DECL (expo
);
374 MPFR_SAVE_EXPO_MARK (expo
);
375 mpfr_init2 (tmp
, MPFR_PREC (u
));
376 /* trunc(u) is always representable in tmp */
378 inex
= mpfr_set (r
, tmp
, rnd_mode
);
380 MPFR_SAVE_EXPO_FREE (expo
);
381 return mpfr_check_range (r
, inex
, rnd_mode
);
385 #undef mpfr_rint_ceil
388 mpfr_rint_ceil (mpfr_ptr r
, mpfr_srcptr u
, mpfr_rnd_t rnd_mode
)
390 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u
) ) || mpfr_integer_p (u
))
391 return mpfr_set (r
, u
, rnd_mode
);
396 MPFR_SAVE_EXPO_DECL (expo
);
397 MPFR_BLOCK_DECL (flags
);
399 MPFR_SAVE_EXPO_MARK (expo
);
400 mpfr_init2 (tmp
, MPFR_PREC (u
));
401 /* ceil(u) is representable in tmp unless an overflow occurs */
402 MPFR_BLOCK (flags
, mpfr_ceil (tmp
, u
));
403 inex
= (MPFR_OVERFLOW (flags
)
404 ? mpfr_overflow (r
, rnd_mode
, MPFR_SIGN_POS
)
405 : mpfr_set (r
, tmp
, rnd_mode
));
407 MPFR_SAVE_EXPO_FREE (expo
);
408 return mpfr_check_range (r
, inex
, rnd_mode
);
412 #undef mpfr_rint_floor
415 mpfr_rint_floor (mpfr_ptr r
, mpfr_srcptr u
, mpfr_rnd_t rnd_mode
)
417 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u
) ) || mpfr_integer_p (u
))
418 return mpfr_set (r
, u
, rnd_mode
);
423 MPFR_SAVE_EXPO_DECL (expo
);
424 MPFR_BLOCK_DECL (flags
);
426 MPFR_SAVE_EXPO_MARK (expo
);
427 mpfr_init2 (tmp
, MPFR_PREC (u
));
428 /* floor(u) is representable in tmp unless an overflow occurs */
429 MPFR_BLOCK (flags
, mpfr_floor (tmp
, u
));
430 inex
= (MPFR_OVERFLOW (flags
)
431 ? mpfr_overflow (r
, rnd_mode
, MPFR_SIGN_NEG
)
432 : mpfr_set (r
, tmp
, rnd_mode
));
434 MPFR_SAVE_EXPO_FREE (expo
);
435 return mpfr_check_range (r
, inex
, rnd_mode
);