1 /* mpfr_round_p -- check if an approximation is roundable.
3 Copyright 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 /* Check against mpfr_can_round? */
28 int mpfr_round_p_2 (mp_limb_t
*, mp_size_t
, mpfr_exp_t
, mpfr_prec_t
);
30 mpfr_round_p (mp_limb_t
*bp
, mp_size_t bn
, mpfr_exp_t err0
, mpfr_prec_t prec
)
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
);
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
);
47 # define mpfr_round_p mpfr_round_p_2
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
58 mpfr_round_p (mp_limb_t
*bp
, mp_size_t bn
, mpfr_exp_t err0
, mpfr_prec_t prec
)
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
);
70 k
= prec
/ GMP_NUMB_BITS
;
71 s
= GMP_NUMB_BITS
- prec
%GMP_NUMB_BITS
;
72 n
= err
/ GMP_NUMB_BITS
- k
;
74 MPFR_ASSERTD (n
>= 0);
75 MPFR_ASSERTD (bn
> k
);
77 /* Check first limb */
80 mask
= s
== GMP_NUMB_BITS
? MP_LIMB_T_MAX
: MPFR_LIMB_MASK (s
);
83 if (MPFR_LIKELY (n
== 0))
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
);
90 return tmp
!= 0 && tmp
!= mask
;
92 else if (MPFR_UNLIKELY (tmp
== 0))
94 /* Check if all (n-1) limbs are 0 */
98 /* Check if final error limb is 0 */
99 s
= GMP_NUMB_BITS
- err
% GMP_NUMB_BITS
;
100 if (s
== GMP_NUMB_BITS
)
105 else if (MPFR_UNLIKELY (tmp
== mask
))
107 /* Check if all (n-1) limbs are 11111111111111111 */
109 if (*bp
-- != MP_LIMB_T_MAX
)
111 /* Check if final error limb is 0 */
112 s
= GMP_NUMB_BITS
- err
% GMP_NUMB_BITS
;
113 if (s
== GMP_NUMB_BITS
)
116 return tmp
!= (MP_LIMB_T_MAX
>> s
);
120 /* First limb is different from 000000 or 1111111 */