1 /* mpfr_set_d -- convert a machine double precision float to
2 a multiple precision floating-point number
4 Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
5 Contributed by the Arenaire and Caramel projects, INRIA.
7 This file is part of the GNU MPFR Library.
9 The GNU MPFR Library is free software; you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
14 The GNU MPFR Library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 License for more details.
19 You should have received a copy of the GNU Lesser General Public License
20 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
21 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
22 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
24 #include <float.h> /* For DOUBLE_ISINF and DOUBLE_ISNAN */
26 #define MPFR_NEED_LONGLONG_H
27 #include "mpfr-impl.h"
29 /* extracts the bits of d in rp[0..n-1] where n=ceil(53/GMP_NUMB_BITS).
30 Assumes d is neither 0 nor NaN nor Inf. */
32 __gmpfr_extract_double (mpfr_limb_ptr rp, double d)
33 /* e=0 iff GMP_NUMB_BITS=32 and rp has only one limb */
37 #if GMP_NUMB_BITS == 32
42 1. Should handle Inf and NaN in IEEE specific code.
43 2. Handle Inf and NaN also in default code, to avoid hangs.
44 3. Generalize to handle all GMP_NUMB_BITS.
45 4. This lits is incomplete and misspelled.
48 MPFR_ASSERTD(!DOUBLE_ISNAN(d));
49 MPFR_ASSERTD(!DOUBLE_ISINF(d));
50 MPFR_ASSERTD(d != 0.0);
55 union ieee_double_extract x;
61 #if GMP_NUMB_BITS >= 64
62 manl = ((MPFR_LIMB_ONE << 63)
63 | ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11));
65 manh = (MPFR_LIMB_ONE << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
66 manl = x.s.manl << 11;
69 else /* denormalized number */
71 #if GMP_NUMB_BITS >= 64
72 manl = ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11);
74 manh = (x.s.manh << 11) /* high 21 bits */
75 | (x.s.manl >> 21); /* middle 11 bits */
76 manl = x.s.manl << 11; /* low 21 bits */
86 #else /* _GMP_IEEE_FLOATS */
89 /* Unknown (or known to be non-IEEE) double format. */
93 MPFR_ASSERTN (d * 0.5 != d);
107 while (d < (1.0 / 65536.0))
119 d *= MP_BASE_AS_DOUBLE;
120 #if GMP_NUMB_BITS >= 64
123 manh = (mp_limb_t) d;
124 manl = (mp_limb_t) ((d - manh) * MP_BASE_AS_DOUBLE);
128 #endif /* _GMP_IEEE_FLOATS */
130 #if GMP_NUMB_BITS >= 64
140 /* End of part included from gmp-2.0.2 */
143 mpfr_set_d (mpfr_ptr r, double d, mpfr_rnd_t rnd_mode)
149 mp_limb_t tmpmant[MPFR_LIMBS_PER_DOUBLE];
150 MPFR_SAVE_EXPO_DECL (expo);
152 if (MPFR_UNLIKELY(DOUBLE_ISNAN(d)))
157 else if (MPFR_UNLIKELY(d == 0))
160 union ieee_double_extract x;
163 /* set correct sign */
169 #else /* _GMP_IEEE_FLOATS */
172 /* This is to get the sign of zero on non-IEEE hardware
173 Some systems support +0.0, -0.0 and unsigned zero.
174 We can't use d==+0.0 since it should be always true,
175 so we check that the memory representation of d is the
176 same than +0.0. etc */
177 /* FIXME: consider the case where +0.0 or -0.0 may have several
179 double poszero = +0.0, negzero = DBL_NEG_ZERO;
180 if (memcmp(&d, &poszero, sizeof(double)) == 0)
182 else if (memcmp(&d, &negzero, sizeof(double)) == 0)
188 return 0; /* 0 is exact */
190 else if (MPFR_UNLIKELY(DOUBLE_ISINF(d)))
197 return 0; /* infinity is exact */
200 /* now d is neither 0, nor NaN nor Inf */
202 MPFR_SAVE_EXPO_MARK (expo);
204 /* warning: don't use tmp=r here, even if SIZE(r) >= MPFR_LIMBS_PER_DOUBLE,
205 since PREC(r) may be different from PREC(tmp), and then both variables
206 would have same precision in the mpfr_set4 call below. */
207 MPFR_MANT(tmp) = tmpmant;
208 MPFR_PREC(tmp) = IEEE_DBL_MANT_DIG;
210 signd = (d < 0) ? MPFR_SIGN_NEG : MPFR_SIGN_POS;
213 /* don't use MPFR_SET_EXP here since the exponent may be out of range */
214 MPFR_EXP(tmp) = __gmpfr_extract_double (tmpmant, d);
217 /* Failed assertion if the stored value is 0 (e.g., if the exponent range
218 has been reduced at the wrong moment and an underflow to 0 occurred).
219 Probably a bug in the C implementation if this happens. */
221 while (tmpmant[i] == 0)
224 MPFR_ASSERTN(i < MPFR_LIMBS_PER_DOUBLE);
228 /* determine the index i-1 of the most significant non-zero limb
229 and the number k of zero high limbs */
230 i = MPFR_LIMBS_PER_DOUBLE;
231 MPN_NORMALIZE_NOT_ZERO(tmpmant, i);
232 k = MPFR_LIMBS_PER_DOUBLE - i;
234 count_leading_zeros (cnt, tmpmant[i - 1]);
236 if (MPFR_LIKELY(cnt != 0))
237 mpn_lshift (tmpmant + k, tmpmant, i, cnt);
239 MPN_COPY (tmpmant + k, tmpmant, i);
241 if (MPFR_UNLIKELY(k != 0))
242 MPN_ZERO (tmpmant, k);
244 /* don't use MPFR_SET_EXP here since the exponent may be out of range */
245 MPFR_EXP(tmp) -= (mpfr_exp_t) (cnt + k * GMP_NUMB_BITS);
247 /* tmp is exact since PREC(tmp)=53 */
248 inexact = mpfr_set4 (r, tmp, rnd_mode, signd);
250 MPFR_SAVE_EXPO_FREE (expo);
251 return mpfr_check_range (r, inexact, rnd_mode);