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