frac.c revision 1.1.1.1
1/* mpfr_frac -- Fractional part of a floating-point number. 2 3Copyright 2002, 2003, 2004, 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 24#define MPFR_NEED_LONGLONG_H 25#include "mpfr-impl.h" 26 27/* Optimization note: it is not a good idea to call mpfr_integer_p, 28 as some cases will take longer (the number may be parsed twice). */ 29 30int 31mpfr_frac (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode) 32{ 33 mpfr_exp_t re, ue; 34 mpfr_prec_t uq; 35 mp_size_t un, tn, t0; 36 mp_limb_t *up, *tp, k; 37 int sh; 38 mpfr_t tmp; 39 mpfr_ptr t; 40 int inex; 41 MPFR_SAVE_EXPO_DECL (expo); 42 43 /* Special cases */ 44 if (MPFR_UNLIKELY(MPFR_IS_NAN(u))) 45 { 46 MPFR_SET_NAN(r); 47 MPFR_RET_NAN; 48 } 49 else if (MPFR_UNLIKELY(MPFR_IS_INF(u) || mpfr_integer_p (u))) 50 { 51 MPFR_SET_SAME_SIGN(r, u); 52 MPFR_SET_ZERO(r); 53 MPFR_RET(0); /* zero is exact */ 54 } 55 56 ue = MPFR_GET_EXP (u); 57 if (ue <= 0) /* |u| < 1 */ 58 return mpfr_set (r, u, rnd_mode); 59 60 /* Now |u| >= 1, meaning that an overflow is not possible. */ 61 62 uq = MPFR_PREC(u); 63 un = (uq - 1) / GMP_NUMB_BITS; /* index of most significant limb */ 64 un -= (mp_size_t) (ue / GMP_NUMB_BITS); 65 /* now the index of the MSL containing bits of the fractional part */ 66 67 up = MPFR_MANT(u); 68 sh = ue % GMP_NUMB_BITS; 69 k = up[un] << sh; 70 /* the first bit of the fractional part is the MSB of k */ 71 72 if (k != 0) 73 { 74 int cnt; 75 76 count_leading_zeros(cnt, k); 77 /* first bit 1 of the fractional part -> MSB of the number */ 78 re = -cnt; 79 sh += cnt; 80 MPFR_ASSERTN (sh < GMP_NUMB_BITS); 81 k <<= cnt; 82 } 83 else 84 { 85 re = sh - GMP_NUMB_BITS; 86 /* searching for the first bit 1 (exists since u isn't an integer) */ 87 while (up[--un] == 0) 88 re -= GMP_NUMB_BITS; 89 MPFR_ASSERTN(un >= 0); 90 k = up[un]; 91 count_leading_zeros(sh, k); 92 re -= sh; 93 k <<= sh; 94 } 95 /* The exponent of r will be re */ 96 /* un: index of the limb of u that contains the first bit 1 of the FP */ 97 98 t = (mp_size_t) (MPFR_PREC(r) - 1) / GMP_NUMB_BITS < un ? 99 (mpfr_init2 (tmp, (un + 1) * GMP_NUMB_BITS), tmp) : r; 100 /* t has enough precision to contain the fractional part of u */ 101 /* If we use a temporary variable, we take the non-significant bits 102 of u into account, because of the mpn_lshift below. */ 103 MPFR_SET_SAME_SIGN(t, u); 104 105 /* Put the fractional part of u into t */ 106 tn = (MPFR_PREC(t) - 1) / GMP_NUMB_BITS; 107 MPFR_ASSERTN(tn >= un); 108 t0 = tn - un; 109 tp = MPFR_MANT(t); 110 if (sh == 0) 111 MPN_COPY_DECR(tp + t0, up, un + 1); 112 else /* warning: un may be 0 here */ 113 tp[tn] = k | ((un) ? mpn_lshift (tp + t0, up, un, sh) : (mp_limb_t) 0); 114 if (t0 > 0) 115 MPN_ZERO(tp, t0); 116 117 MPFR_SAVE_EXPO_MARK (expo); 118 119 if (t != r) 120 { /* t is tmp */ 121 MPFR_EXP (t) = 0; /* should be re, but not necessarily in the range */ 122 inex = mpfr_set (r, t, rnd_mode); /* no underflow */ 123 mpfr_clear (t); 124 MPFR_EXP (r) += re; 125 } 126 else 127 { /* There may be remaining non-significant bits in t (= r). */ 128 int carry; 129 130 MPFR_EXP (r) = re; 131 carry = mpfr_round_raw (tp, tp, 132 (mpfr_prec_t) (tn + 1) * GMP_NUMB_BITS, 133 MPFR_IS_NEG (r), MPFR_PREC (r), rnd_mode, 134 &inex); 135 if (carry) 136 { 137 tp[tn] = MPFR_LIMB_HIGHBIT; 138 MPFR_EXP (r) ++; 139 } 140 } 141 142 MPFR_SAVE_EXPO_FREE (expo); 143 return mpfr_check_range (r, inex, rnd_mode); 144} 145