1/* mpfr_round_p -- check if an approximation is roundable.
2
3Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4Contributed by the AriC and Caramel projects, INRIA.
5
6This file is part of the GNU MPFR Library.
7
8The GNU MPFR Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MPFR Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
20http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
22
23#include "mpfr-impl.h"
24
25/* Check against mpfr_can_round? */
26#ifdef WANT_ASSERT
27# if WANT_ASSERT >= 2
28int mpfr_round_p_2 (mp_limb_t *, mp_size_t, mpfr_exp_t, mpfr_prec_t);
29int
30mpfr_round_p (mp_limb_t *bp, mp_size_t bn, mpfr_exp_t err0, mpfr_prec_t prec)
31{
32  int i1, i2;
33
34  i1 = mpfr_round_p_2 (bp, bn, err0, prec);
35  i2 = mpfr_can_round_raw (bp, bn, MPFR_SIGN_POS, err0,
36                           MPFR_RNDN, MPFR_RNDZ, prec);
37  if (i1 != i2)
38    {
39      fprintf (stderr, "mpfr_round_p(%d) != mpfr_can_round(%d)!\n"
40               "bn = %lu, err0 = %ld, prec = %lu\nbp = ", i1, i2,
41               (unsigned long) bn, (long) err0, (unsigned long) prec);
42      gmp_fprintf (stderr, "%NX\n", bp, bn);
43      MPFR_ASSERTN (0);
44    }
45  return i1;
46}
47# define mpfr_round_p mpfr_round_p_2
48# endif
49#endif
50
51/*
52 * Assuming {bp, bn} is an approximation of a non-singular number
53 * with error at most equal to 2^(EXP(b)-err0) (`err0' bits of b are known)
54 * of direction unknown, check if we can round b toward zero with
55 * precision prec.
56 */
57int
58mpfr_round_p (mp_limb_t *bp, mp_size_t bn, mpfr_exp_t err0, mpfr_prec_t prec)
59{
60  mpfr_prec_t err;
61  mp_size_t k, n;
62  mp_limb_t tmp, mask;
63  int s;
64
65  err = (mpfr_prec_t) bn * GMP_NUMB_BITS;
66  if (MPFR_UNLIKELY (err0 <= 0 || (mpfr_uexp_t) err0 <= prec || prec >= err))
67    return 0;  /* can't round */
68  err = MIN (err, (mpfr_uexp_t) err0);
69
70  k = prec / GMP_NUMB_BITS;
71  s = GMP_NUMB_BITS - prec%GMP_NUMB_BITS;
72  n = err / GMP_NUMB_BITS - k;
73
74  MPFR_ASSERTD (n >= 0);
75  MPFR_ASSERTD (bn > k);
76
77  /* Check first limb */
78  bp += bn-1-k;
79  tmp = *bp--;
80  mask = s == GMP_NUMB_BITS ? MP_LIMB_T_MAX : MPFR_LIMB_MASK (s);
81  tmp &= mask;
82
83  if (MPFR_LIKELY (n == 0))
84    {
85      /* prec and error are in the same limb */
86      s = GMP_NUMB_BITS - err % GMP_NUMB_BITS;
87      MPFR_ASSERTD (s < GMP_NUMB_BITS);
88      tmp  >>= s;
89      mask >>= s;
90      return tmp != 0 && tmp != mask;
91    }
92  else if (MPFR_UNLIKELY (tmp == 0))
93    {
94      /* Check if all (n-1) limbs are 0 */
95      while (--n)
96        if (*bp-- != 0)
97          return 1;
98      /* Check if final error limb is 0 */
99      s = GMP_NUMB_BITS - err % GMP_NUMB_BITS;
100      if (s == GMP_NUMB_BITS)
101        return 0;
102      tmp = *bp >> s;
103      return tmp != 0;
104    }
105  else if (MPFR_UNLIKELY (tmp == mask))
106    {
107      /* Check if all (n-1) limbs are 11111111111111111 */
108      while (--n)
109        if (*bp-- != MP_LIMB_T_MAX)
110          return 1;
111      /* Check if final error limb is 0 */
112      s = GMP_NUMB_BITS - err % GMP_NUMB_BITS;
113      if (s == GMP_NUMB_BITS)
114        return 0;
115      tmp = *bp >> s;
116      return tmp != (MP_LIMB_T_MAX >> s);
117    }
118  else
119    {
120      /* First limb is different from 000000 or 1111111 */
121      return 1;
122    }
123}
124