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 Free Software Foundation, Inc.
5 Contributed by the Arenaire and Cacao 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 2.1 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.LIB. If not, write to
21 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
22 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/BITS_PER_MP_LIMB).
30 Assumes d is neither 0 nor NaN nor Inf. */
32 __gmpfr_extract_double (mp_ptr rp, double d)
33 /* e=0 iff BITS_PER_MP_LIMB=32 and rp has only one limb */
37 #if BITS_PER_MP_LIMB == 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 BITS_PER_MP_LIMB.
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 BITS_PER_MP_LIMB >= 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 BITS_PER_MP_LIMB >= 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 BITS_PER_MP_LIMB >= 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 BITS_PER_MP_LIMB >= 64
140 /* End of part included from gmp-2.0.2 */
143 mpfr_set_d (mpfr_ptr r, double d, mp_rnd_t rnd_mode)
149 mp_limb_t tmpmant[MPFR_LIMBS_PER_DOUBLE];
150 MPFR_SAVE_EXPO_DECL (expo);
154 if (MPFR_UNLIKELY(DOUBLE_ISNAN(d)))
159 else if (MPFR_UNLIKELY(d == 0))
162 union ieee_double_extract x;
165 /* set correct sign */
171 #else /* _GMP_IEEE_FLOATS */
174 /* This is to get the sign of zero on non-IEEE hardware
175 Some systems support +0.0, -0.0 and unsigned zero.
176 We can't use d==+0.0 since it should be always true,
177 so we check that the memory representation of d is the
178 same than +0.0. etc */
179 /* FIXME: consider the case where +0.0 or -0.0 may have several
181 double poszero = +0.0, negzero = DBL_NEG_ZERO;
182 if (memcmp(&d, &poszero, sizeof(double)) == 0)
184 else if (memcmp(&d, &negzero, sizeof(double)) == 0)
190 return 0; /* 0 is exact */
192 else if (MPFR_UNLIKELY(DOUBLE_ISINF(d)))
199 return 0; /* infinity is exact */
202 /* now d is neither 0, nor NaN nor Inf */
204 MPFR_SAVE_EXPO_MARK (expo);
206 /* warning: don't use tmp=r here, even if SIZE(r) >= MPFR_LIMBS_PER_DOUBLE,
207 since PREC(r) may be different from PREC(tmp), and then both variables
208 would have same precision in the mpfr_set4 call below. */
209 MPFR_MANT(tmp) = tmpmant;
210 MPFR_PREC(tmp) = IEEE_DBL_MANT_DIG;
212 signd = (d < 0) ? MPFR_SIGN_NEG : MPFR_SIGN_POS;
215 /* don't use MPFR_SET_EXP here since the exponent may be out of range */
216 MPFR_EXP(tmp) = __gmpfr_extract_double (tmpmant, d);
219 /* Failed assertion if the stored value is 0 (e.g., if the exponent range
220 has been reduced at the wrong moment and an underflow to 0 occurred).
221 Probably a bug in the C implementation if this happens. */
223 while (tmpmant[i] == 0)
226 MPFR_ASSERTN(i < MPFR_LIMBS_PER_DOUBLE);
230 /* determine the index i-1 of the most significant non-zero limb
231 and the number k of zero high limbs */
232 i = MPFR_LIMBS_PER_DOUBLE;
233 MPN_NORMALIZE_NOT_ZERO(tmpmant, i);
234 k = MPFR_LIMBS_PER_DOUBLE - i;
236 count_leading_zeros (cnt, tmpmant[i - 1]);
238 if (MPFR_LIKELY(cnt != 0))
239 mpn_lshift (tmpmant + k, tmpmant, i, cnt);
241 MPN_COPY (tmpmant + k, tmpmant, i);
243 if (MPFR_UNLIKELY(k != 0))
244 MPN_ZERO (tmpmant, k);
246 /* don't use MPFR_SET_EXP here since the exponent may be out of range */
247 MPFR_EXP(tmp) -= (mp_exp_t) (cnt + k * BITS_PER_MP_LIMB);
249 /* tmp is exact since PREC(tmp)=53 */
250 inexact = mpfr_set4 (r, tmp, rnd_mode, signd);
252 MPFR_SAVE_EXPO_FREE (expo);
253 return mpfr_check_range (r, inexact, rnd_mode);