1/* mpfr_set_d -- convert a machine double precision float to
2                 a multiple precision floating-point number
3
4Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
5Contributed by the AriC and Caramel projects, INRIA.
6
7This file is part of the GNU MPFR Library.
8
9The GNU MPFR Library is free software; you can redistribute it and/or modify
10it under the terms of the GNU Lesser General Public License as published by
11the Free Software Foundation; either version 3 of the License, or (at your
12option) any later version.
13
14The GNU MPFR Library is distributed in the hope that it will be useful, but
15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
17License for more details.
18
19You should have received a copy of the GNU Lesser General Public License
20along with the GNU MPFR Library; see the file COPYING.LESSER.  If not, see
21http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
2251 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23
24#include <float.h>  /* For DOUBLE_ISINF and DOUBLE_ISNAN */
25
26#define MPFR_NEED_LONGLONG_H
27#include "mpfr-impl.h"
28
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. */
31static long
32__gmpfr_extract_double (mpfr_limb_ptr rp, double d)
33     /* e=0 iff GMP_NUMB_BITS=32 and rp has only one limb */
34{
35  long exp;
36  mp_limb_t manl;
37#if GMP_NUMB_BITS == 32
38  mp_limb_t manh;
39#endif
40
41  /* BUGS
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.
46   */
47
48  MPFR_ASSERTD(!DOUBLE_ISNAN(d));
49  MPFR_ASSERTD(!DOUBLE_ISINF(d));
50  MPFR_ASSERTD(d != 0.0);
51
52#if _GMP_IEEE_FLOATS
53
54  {
55    union ieee_double_extract x;
56    x.d = d;
57
58    exp = x.s.exp;
59    if (exp)
60      {
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));
64#else
65        manh = (MPFR_LIMB_ONE << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
66        manl = x.s.manl << 11;
67#endif
68      }
69    else /* denormalized number */
70      {
71#if GMP_NUMB_BITS >= 64
72        manl = ((mp_limb_t) x.s.manh << 43) | ((mp_limb_t) x.s.manl << 11);
73#else
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 */
77#endif
78      }
79
80    if (exp)
81      exp -= 1022;
82    else
83      exp = -1021;
84  }
85
86#else /* _GMP_IEEE_FLOATS */
87
88  {
89    /* Unknown (or known to be non-IEEE) double format.  */
90    exp = 0;
91    if (d >= 1.0)
92      {
93        MPFR_ASSERTN (d * 0.5 != d);
94        while (d >= 32768.0)
95          {
96            d *= (1.0 / 65536.0);
97            exp += 16;
98          }
99        while (d >= 1.0)
100          {
101            d *= 0.5;
102            exp += 1;
103          }
104      }
105    else if (d < 0.5)
106      {
107        while (d < (1.0 / 65536.0))
108          {
109            d *=  65536.0;
110            exp -= 16;
111          }
112        while (d < 0.5)
113          {
114            d *= 2.0;
115            exp -= 1;
116          }
117      }
118
119    d *= MP_BASE_AS_DOUBLE;
120#if GMP_NUMB_BITS >= 64
121    manl = d;
122#else
123    manh = (mp_limb_t) d;
124    manl = (mp_limb_t) ((d - manh) * MP_BASE_AS_DOUBLE);
125#endif
126  }
127
128#endif /* _GMP_IEEE_FLOATS */
129
130#if GMP_NUMB_BITS >= 64
131  rp[0] = manl;
132#else
133  rp[1] = manh;
134  rp[0] = manl;
135#endif
136
137  return exp;
138}
139
140/* End of part included from gmp-2.0.2 */
141
142int
143mpfr_set_d (mpfr_ptr r, double d, mpfr_rnd_t rnd_mode)
144{
145  int signd, inexact;
146  unsigned int cnt;
147  mp_size_t i, k;
148  mpfr_t tmp;
149  mp_limb_t tmpmant[MPFR_LIMBS_PER_DOUBLE];
150  MPFR_SAVE_EXPO_DECL (expo);
151
152  if (MPFR_UNLIKELY(DOUBLE_ISNAN(d)))
153    {
154      MPFR_SET_NAN(r);
155      MPFR_RET_NAN;
156    }
157  else if (MPFR_UNLIKELY(d == 0))
158    {
159#if _GMP_IEEE_FLOATS
160      union ieee_double_extract x;
161
162      MPFR_SET_ZERO(r);
163      /* set correct sign */
164      x.d = d;
165      if (x.s.sig == 1)
166        MPFR_SET_NEG(r);
167      else
168        MPFR_SET_POS(r);
169#else /* _GMP_IEEE_FLOATS */
170      MPFR_SET_ZERO(r);
171      {
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
178           representations. */
179        double poszero = +0.0, negzero = DBL_NEG_ZERO;
180        if (memcmp(&d, &poszero, sizeof(double)) == 0)
181          MPFR_SET_POS(r);
182        else if (memcmp(&d, &negzero, sizeof(double)) == 0)
183          MPFR_SET_NEG(r);
184        else
185          MPFR_SET_POS(r);
186      }
187#endif
188      return 0; /* 0 is exact */
189    }
190  else if (MPFR_UNLIKELY(DOUBLE_ISINF(d)))
191    {
192      MPFR_SET_INF(r);
193      if (d > 0)
194        MPFR_SET_POS(r);
195      else
196        MPFR_SET_NEG(r);
197      return 0; /* infinity is exact */
198    }
199
200  /* now d is neither 0, nor NaN nor Inf */
201
202  MPFR_SAVE_EXPO_MARK (expo);
203
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;
209
210  signd = (d < 0) ? MPFR_SIGN_NEG : MPFR_SIGN_POS;
211  d = ABS (d);
212
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);
215
216#ifdef WANT_ASSERT
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. */
220  i = 0;
221  while (tmpmant[i] == 0)
222    {
223      i++;
224      MPFR_ASSERTN(i < MPFR_LIMBS_PER_DOUBLE);
225    }
226#endif
227
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;
233
234  count_leading_zeros (cnt, tmpmant[i - 1]);
235
236  if (MPFR_LIKELY(cnt != 0))
237    mpn_lshift (tmpmant + k, tmpmant, i, cnt);
238  else if (k != 0)
239    MPN_COPY (tmpmant + k, tmpmant, i);
240
241  if (MPFR_UNLIKELY(k != 0))
242    MPN_ZERO (tmpmant, k);
243
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);
246
247  /* tmp is exact since PREC(tmp)=53 */
248  inexact = mpfr_set4 (r, tmp, rnd_mode, signd);
249
250  MPFR_SAVE_EXPO_FREE (expo);
251  return mpfr_check_range (r, inexact, rnd_mode);
252}
253
254
255
256