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