1 /* mpfr_frac -- Fractional part of a floating-point number.
3 Copyright 2002, 2003, 2004, 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. */
24 #define MPFR_NEED_LONGLONG_H
25 #include "mpfr-impl.h"
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). */
31 mpfr_frac (mpfr_ptr r, mpfr_srcptr u, mpfr_rnd_t rnd_mode)
36 mp_limb_t *up, *tp, k;
41 MPFR_SAVE_EXPO_DECL (expo);
44 if (MPFR_UNLIKELY(MPFR_IS_NAN(u)))
49 else if (MPFR_UNLIKELY(MPFR_IS_INF(u) || mpfr_integer_p (u)))
51 MPFR_SET_SAME_SIGN(r, u);
53 MPFR_RET(0); /* zero is exact */
56 ue = MPFR_GET_EXP (u);
57 if (ue <= 0) /* |u| < 1 */
58 return mpfr_set (r, u, rnd_mode);
60 /* Now |u| >= 1, meaning that an overflow is not possible. */
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 */
68 sh = ue % GMP_NUMB_BITS;
70 /* the first bit of the fractional part is the MSB of k */
76 count_leading_zeros(cnt, k);
77 /* first bit 1 of the fractional part -> MSB of the number */
80 MPFR_ASSERTN (sh < GMP_NUMB_BITS);
85 re = sh - GMP_NUMB_BITS;
86 /* searching for the first bit 1 (exists since u isn't an integer) */
89 MPFR_ASSERTN(un >= 0);
91 count_leading_zeros(sh, k);
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 */
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);
105 /* Put the fractional part of u into t */
106 tn = (MPFR_PREC(t) - 1) / GMP_NUMB_BITS;
107 MPFR_ASSERTN(tn >= un);
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);
117 MPFR_SAVE_EXPO_MARK (expo);
121 MPFR_EXP (t) = 0; /* should be re, but not necessarily in the range */
122 inex = mpfr_set (r, t, rnd_mode); /* no underflow */
127 { /* There may be remaining non-significant bits in t (= r). */
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,
137 tp[tn] = MPFR_LIMB_HIGHBIT;
142 MPFR_SAVE_EXPO_FREE (expo);
143 return mpfr_check_range (r, inex, rnd_mode);