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
);