1 /* mpfr_rint -- Round to an integer.
2
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.
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 http://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 int
mpfr_rint(mpfr_ptr r,mpfr_srcptr u,mpfr_rnd_t rnd_mode)28 mpfr_rint (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
29 {
30 int sign;
31 int rnd_away;
32 mpfr_exp_t exp;
33
34 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ))
35 {
36 if (MPFR_IS_NAN(u))
37 {
38 MPFR_SET_NAN(r);
39 MPFR_RET_NAN;
40 }
41 MPFR_SET_SAME_SIGN(r, u);
42 if (MPFR_IS_INF(u))
43 {
44 MPFR_SET_INF(r);
45 MPFR_RET(0); /* infinity is exact */
46 }
47 else /* now u is zero */
48 {
49 MPFR_ASSERTD(MPFR_IS_ZERO(u));
50 MPFR_SET_ZERO(r);
51 MPFR_RET(0); /* zero is exact */
52 }
53 }
54 MPFR_SET_SAME_SIGN (r, u); /* Does nothing if r==u */
55
56 sign = MPFR_INT_SIGN (u);
57 exp = MPFR_GET_EXP (u);
58
59 rnd_away =
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) */
65
66 /* rnd_away:
67 1 if round away from zero,
68 0 if round to zero,
69 -1 if not decided yet.
70 */
71
72 if (MPFR_UNLIKELY (exp <= 0)) /* 0 < |u| < 1 ==> round |u| to 0 or 1 */
73 {
74 /* Note: in the MPFR_RNDN mode, 0.5 must be rounded to 0. */
75 if (rnd_away != 0 &&
76 (rnd_away > 0 ||
77 (exp == 0 && (rnd_mode == MPFR_RNDNA ||
78 !mpfr_powerof2_raw (u)))))
79 {
80 mp_limb_t *rp;
81 mp_size_t rm;
82
83 rp = MPFR_MANT(r);
84 rm = (MPFR_PREC(r) - 1) / GMP_NUMB_BITS;
85 rp[rm] = MPFR_LIMB_HIGHBIT;
86 MPN_ZERO(rp, rm);
87 MPFR_SET_EXP (r, 1); /* |r| = 1 */
88 MPFR_RET(sign > 0 ? 2 : -2);
89 }
90 else
91 {
92 MPFR_SET_ZERO(r); /* r = 0 */
93 MPFR_RET(sign > 0 ? -2 : 2);
94 }
95 }
96 else /* exp > 0, |u| >= 1 */
97 {
98 mp_limb_t *up, *rp;
99 mp_size_t un, rn, ui;
100 int sh, idiff;
101 int uflags;
102
103 /*
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.
108 */
109
110 up = MPFR_MANT(u);
111 rp = MPFR_MANT(r);
112
113 un = MPFR_LIMB_SIZE(u);
114 rn = MPFR_LIMB_SIZE(r);
115 MPFR_UNSIGNED_MINUS_MODULO (sh, MPFR_PREC (r));
116
117 MPFR_SET_EXP (r, exp); /* Does nothing if r==u */
118
119 if ((exp - 1) / GMP_NUMB_BITS >= un)
120 {
121 ui = un;
122 idiff = 0;
123 uflags = 0; /* u is an integer, representable or not in r */
124 }
125 else
126 {
127 mp_size_t uj;
128
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 */
133
134 uflags = idiff == 0 || (up[uj] << idiff) == 0 ? 0 : 2;
135 if (uflags == 0)
136 while (uj > 0)
137 if (up[--uj] != 0)
138 {
139 uflags = 2;
140 break;
141 }
142 }
143
144 if (ui > rn)
145 {
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 */
150 if (rnd_away < 0)
151 {
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 */
157 mp_limb_t a, b;
158 /* a: rounding bit and some of the following bits */
159 /* b: boundary for a (weight of the rounding bit in a) */
160 if (sh != 0)
161 {
162 a = rp[0] & ((MPFR_LIMB_ONE << sh) - 1);
163 b = MPFR_LIMB_ONE << (sh - 1);
164 }
165 else
166 {
167 a = up[un - rn - 1];
168 b = MPFR_LIMB_HIGHBIT;
169 }
170 rnd_away = a > b;
171 if (a == b)
172 {
173 mp_size_t i;
174 for (i = un - rn - 1 - (sh == 0); i >= 0; i--)
175 if (up[i] != 0)
176 {
177 rnd_away = 1;
178 break;
179 }
180 }
181 }
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));
186 }
187 if (uflags == 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 */
191 else
192 {
193 mp_size_t i;
194 for (i = un - rn - 1; i >= 0; i--)
195 if (up[i] != 0)
196 {
197 uflags = 1; /* u is not representable in r */
198 break;
199 }
200 }
201 }
202 }
203 else /* ui <= rn */
204 {
205 mp_size_t uj, rj;
206 int ush;
207
208 uj = un - ui; /* lowest limb of the integer part in u */
209 rj = rn - ui; /* lowest limb of the integer part in r */
210
211 if (MPFR_LIKELY (rp != up))
212 MPN_COPY(rp + rj, up + uj, ui);
213
214 /* Ignore the lowest rj limbs, all equal to zero. */
215 rp += rj;
216 rn = ui;
217
218 /* number of fractional bits in whole rp[0] */
219 ush = idiff == 0 ? 0 : GMP_NUMB_BITS - idiff;
220
221 if (rj == 0 && ush < sh)
222 {
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 */
229 }
230 else /* The integer part of u fits in r, we'll round to it. */
231 sh = ush;
232
233 if (rnd_away < 0)
234 {
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 */
242 mp_limb_t a, b;
243 /* a: rounding bit and some of the following bits */
244 /* b: boundary for a (weight of the rounding bit in a) */
245 if (sh != 0)
246 {
247 a = rp[0] & ((MPFR_LIMB_ONE << sh) - 1);
248 b = MPFR_LIMB_ONE << (sh - 1);
249 }
250 else
251 {
252 MPFR_ASSERTD (uj >= 1); /* see above */
253 a = up[uj - 1];
254 b = MPFR_LIMB_HIGHBIT;
255 }
256 rnd_away = a > b;
257 if (a == b)
258 {
259 mp_size_t i;
260 for (i = uj - 1 - (sh == 0); i >= 0; i--)
261 if (up[i] != 0)
262 {
263 rnd_away = 1;
264 break;
265 }
266 }
267 }
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));
273 }
274 /* Now we can make the low rj limbs to 0 */
275 MPN_ZERO (rp-rj, rj);
276 }
277
278 if (sh != 0)
279 rp[0] &= MP_LIMB_T_MAX << sh;
280
281 /* If u is a representable integer, there is no rounding. */
282 if (uflags == 0)
283 MPFR_RET(0);
284
285 MPFR_ASSERTD (rnd_away >= 0); /* rounding direction is defined */
286 if (rnd_away && mpn_add_1(rp, rp, rn, MPFR_LIMB_ONE << sh))
287 {
288 if (exp == __gmpfr_emax)
289 return mpfr_overflow(r, rnd_mode, MPFR_SIGN(r)) >= 0 ?
290 uflags : -uflags;
291 else
292 {
293 MPFR_SET_EXP(r, exp + 1);
294 rp[rn-1] = MPFR_LIMB_HIGHBIT;
295 }
296 }
297
298 MPFR_RET (rnd_away ^ (sign < 0) ? uflags : -uflags);
299 } /* exp > 0, |u| >= 1 */
300 }
301
302 #undef mpfr_round
303
304 int
mpfr_round(mpfr_ptr r,mpfr_srcptr u)305 mpfr_round (mpfr_ptr r, mpfr_srcptr u)
306 {
307 return mpfr_rint (r, u, MPFR_RNDNA);
308 }
309
310 #undef mpfr_trunc
311
312 int
mpfr_trunc(mpfr_ptr r,mpfr_srcptr u)313 mpfr_trunc (mpfr_ptr r, mpfr_srcptr u)
314 {
315 return mpfr_rint (r, u, MPFR_RNDZ);
316 }
317
318 #undef mpfr_ceil
319
320 int
mpfr_ceil(mpfr_ptr r,mpfr_srcptr u)321 mpfr_ceil (mpfr_ptr r, mpfr_srcptr u)
322 {
323 return mpfr_rint (r, u, MPFR_RNDU);
324 }
325
326 #undef mpfr_floor
327
328 int
mpfr_floor(mpfr_ptr r,mpfr_srcptr u)329 mpfr_floor (mpfr_ptr r, mpfr_srcptr u)
330 {
331 return mpfr_rint (r, u, MPFR_RNDD);
332 }
333
334 #undef mpfr_rint_round
335
336 int
mpfr_rint_round(mpfr_ptr r,mpfr_srcptr u,mpfr_rnd_t rnd_mode)337 mpfr_rint_round (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
338 {
339 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
340 return mpfr_set (r, u, rnd_mode);
341 else
342 {
343 mpfr_t tmp;
344 int inex;
345 MPFR_SAVE_EXPO_DECL (expo);
346 MPFR_BLOCK_DECL (flags);
347
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));
355 mpfr_clear (tmp);
356 MPFR_SAVE_EXPO_FREE (expo);
357 return mpfr_check_range (r, inex, rnd_mode);
358 }
359 }
360
361 #undef mpfr_rint_trunc
362
363 int
mpfr_rint_trunc(mpfr_ptr r,mpfr_srcptr u,mpfr_rnd_t rnd_mode)364 mpfr_rint_trunc (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
365 {
366 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
367 return mpfr_set (r, u, rnd_mode);
368 else
369 {
370 mpfr_t tmp;
371 int inex;
372 MPFR_SAVE_EXPO_DECL (expo);
373
374 MPFR_SAVE_EXPO_MARK (expo);
375 mpfr_init2 (tmp, MPFR_PREC (u));
376 /* trunc(u) is always representable in tmp */
377 mpfr_trunc (tmp, u);
378 inex = mpfr_set (r, tmp, rnd_mode);
379 mpfr_clear (tmp);
380 MPFR_SAVE_EXPO_FREE (expo);
381 return mpfr_check_range (r, inex, rnd_mode);
382 }
383 }
384
385 #undef mpfr_rint_ceil
386
387 int
mpfr_rint_ceil(mpfr_ptr r,mpfr_srcptr u,mpfr_rnd_t rnd_mode)388 mpfr_rint_ceil (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
389 {
390 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
391 return mpfr_set (r, u, rnd_mode);
392 else
393 {
394 mpfr_t tmp;
395 int inex;
396 MPFR_SAVE_EXPO_DECL (expo);
397 MPFR_BLOCK_DECL (flags);
398
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));
406 mpfr_clear (tmp);
407 MPFR_SAVE_EXPO_FREE (expo);
408 return mpfr_check_range (r, inex, rnd_mode);
409 }
410 }
411
412 #undef mpfr_rint_floor
413
414 int
mpfr_rint_floor(mpfr_ptr r,mpfr_srcptr u,mpfr_rnd_t rnd_mode)415 mpfr_rint_floor (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
416 {
417 if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(u) ) || mpfr_integer_p (u))
418 return mpfr_set (r, u, rnd_mode);
419 else
420 {
421 mpfr_t tmp;
422 int inex;
423 MPFR_SAVE_EXPO_DECL (expo);
424 MPFR_BLOCK_DECL (flags);
425
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));
433 mpfr_clear (tmp);
434 MPFR_SAVE_EXPO_FREE (expo);
435 return mpfr_check_range (r, inex, rnd_mode);
436 }
437 }
438