1/* mpfr_set_ld -- convert a machine long double to
2                  a multiple precision floating-point number
3
4Copyright 2002, 2003, 2004, 2005, 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>
25
26#define MPFR_NEED_LONGLONG_H
27#include "mpfr-impl.h"
28
29/* Various i386 systems have been seen with <float.h> LDBL constants equal
30   to the DBL ones, whereas they ought to be bigger, reflecting the 10-byte
31   IEEE extended format on that processor.  gcc 3.2.1 on FreeBSD and Solaris
32   has been seen with the problem, and gcc 2.95.4 on FreeBSD 4.7.  */
33
34#if HAVE_LDOUBLE_IEEE_EXT_LITTLE
35static const union {
36  char         bytes[10];
37  long double  d;
38} ldbl_max_struct = {
39  { '\377','\377','\377','\377',
40    '\377','\377','\377','\377',
41    '\376','\177' }
42};
43#define MPFR_LDBL_MAX   (ldbl_max_struct.d)
44#else
45#define MPFR_LDBL_MAX   LDBL_MAX
46#endif
47
48#ifndef HAVE_LDOUBLE_IEEE_EXT_LITTLE
49
50/* Generic code */
51int
52mpfr_set_ld (mpfr_ptr r, long double d, mpfr_rnd_t rnd_mode)
53{
54  mpfr_t t, u;
55  int inexact, shift_exp;
56  long double x;
57  MPFR_SAVE_EXPO_DECL (expo);
58
59  /* Check for NAN */
60  LONGDOUBLE_NAN_ACTION (d, goto nan);
61
62  /* Check for INF */
63  if (d > MPFR_LDBL_MAX)
64    {
65      mpfr_set_inf (r, 1);
66      return 0;
67    }
68  else if (d < -MPFR_LDBL_MAX)
69    {
70      mpfr_set_inf (r, -1);
71      return 0;
72    }
73  /* Check for ZERO */
74  else if (d == 0.0)
75    return mpfr_set_d (r, (double) d, rnd_mode);
76
77  mpfr_init2 (t, MPFR_LDBL_MANT_DIG);
78  mpfr_init2 (u, IEEE_DBL_MANT_DIG);
79
80  MPFR_SAVE_EXPO_MARK (expo);
81
82 convert:
83  x = d;
84  MPFR_SET_ZERO (t);  /* The sign doesn't matter. */
85  shift_exp = 0; /* invariant: remainder to deal with is d*2^shift_exp */
86  while (x != (long double) 0.0)
87    {
88      /* Check overflow of double */
89      if (x > (long double) DBL_MAX || (-x) > (long double) DBL_MAX)
90        {
91          long double div9, div10, div11, div12, div13;
92
93#define TWO_64 18446744073709551616.0 /* 2^64 */
94#define TWO_128 (TWO_64 * TWO_64)
95#define TWO_256 (TWO_128 * TWO_128)
96          div9 = (long double) (double) (TWO_256 * TWO_256); /* 2^(2^9) */
97          div10 = div9 * div9;
98          div11 = div10 * div10; /* 2^(2^11) */
99          div12 = div11 * div11; /* 2^(2^12) */
100          div13 = div12 * div12; /* 2^(2^13) */
101          if (ABS (x) >= div13)
102            {
103              x /= div13; /* exact */
104              shift_exp += 8192;
105              mpfr_div_2si (t, t, 8192, MPFR_RNDZ);
106            }
107          if (ABS (x) >= div12)
108            {
109              x /= div12; /* exact */
110              shift_exp += 4096;
111              mpfr_div_2si (t, t, 4096, MPFR_RNDZ);
112            }
113          if (ABS (x) >= div11)
114            {
115              x /= div11; /* exact */
116              shift_exp += 2048;
117              mpfr_div_2si (t, t, 2048, MPFR_RNDZ);
118            }
119          if (ABS (x) >= div10)
120            {
121              x /= div10; /* exact */
122              shift_exp += 1024;
123              mpfr_div_2si (t, t, 1024, MPFR_RNDZ);
124            }
125          /* warning: we may have DBL_MAX=2^1024*(1-2^(-53)) < x < 2^1024,
126             therefore we have one extra exponent reduction step */
127          if (ABS (x) >= div9)
128            {
129              x /= div9; /* exact */
130              shift_exp += 512;
131              mpfr_div_2si (t, t, 512, MPFR_RNDZ);
132            }
133        } /* Check overflow of double */
134      else /* no overflow on double */
135        {
136          long double div9, div10, div11;
137
138          div9 = (long double) (double) 7.4583407312002067432909653e-155;
139          /* div9 = 2^(-2^9) */
140          div10 = div9  * div9;  /* 2^(-2^10) */
141          div11 = div10 * div10; /* 2^(-2^11) if extended precision */
142          /* since -DBL_MAX <= x <= DBL_MAX, the cast to double should not
143             overflow here */
144          if (ABS(x) < div10 &&
145              div11 != (long double) 0.0 &&
146              div11 / div10 == div10) /* possible underflow */
147            {
148              long double div12, div13;
149              /* After the divisions, any bit of x must be >= div10,
150                 hence the possible division by div9. */
151              div12 = div11 * div11; /* 2^(-2^12) */
152              div13 = div12 * div12; /* 2^(-2^13) */
153              if (ABS (x) <= div13)
154                {
155                  x /= div13; /* exact */
156                  shift_exp -= 8192;
157                  mpfr_mul_2si (t, t, 8192, MPFR_RNDZ);
158                }
159              if (ABS (x) <= div12)
160                {
161                  x /= div12; /* exact */
162                  shift_exp -= 4096;
163                  mpfr_mul_2si (t, t, 4096, MPFR_RNDZ);
164                }
165              if (ABS (x) <= div11)
166                {
167                  x /= div11; /* exact */
168                  shift_exp -= 2048;
169                  mpfr_mul_2si (t, t, 2048, MPFR_RNDZ);
170                }
171              if (ABS (x) <= div10)
172                {
173                  x /= div10; /* exact */
174                  shift_exp -= 1024;
175                  mpfr_mul_2si (t, t, 1024, MPFR_RNDZ);
176                }
177              if (ABS(x) <= div9)
178                {
179                  x /= div9;  /* exact */
180                  shift_exp -= 512;
181                  mpfr_mul_2si (t, t, 512, MPFR_RNDZ);
182                }
183            }
184          else /* no underflow */
185            {
186              inexact = mpfr_set_d (u, (double) x, MPFR_RNDZ);
187              MPFR_ASSERTD (inexact == 0);
188              if (mpfr_add (t, t, u, MPFR_RNDZ) != 0)
189                {
190                  if (!mpfr_number_p (t))
191                    break;
192                  /* Inexact. This cannot happen unless the C implementation
193                     "lies" on the precision or when long doubles are
194                     implemented with FP expansions like under Mac OS X. */
195                  if (MPFR_PREC (t) != MPFR_PREC (r) + 1)
196                    {
197                      /* We assume that MPFR_PREC (r) < MPFR_PREC_MAX.
198                         The precision MPFR_PREC (r) + 1 allows us to
199                         deduce the rounding bit and the sticky bit. */
200                      mpfr_set_prec (t, MPFR_PREC (r) + 1);
201                      goto convert;
202                    }
203                  else
204                    {
205                      mp_limb_t *tp;
206                      int rb_mask;
207
208                      /* Since mpfr_add was inexact, the sticky bit is 1. */
209                      tp = MPFR_MANT (t);
210                      rb_mask = MPFR_LIMB_ONE <<
211                        (GMP_NUMB_BITS - 1 -
212                         (MPFR_PREC (r) & (GMP_NUMB_BITS - 1)));
213                      if (rnd_mode == MPFR_RNDN)
214                        rnd_mode = (*tp & rb_mask) ^ MPFR_IS_NEG (t) ?
215                          MPFR_RNDU : MPFR_RNDD;
216                      *tp |= rb_mask;
217                      break;
218                    }
219                }
220              x -= (long double) mpfr_get_d1 (u); /* exact */
221            }
222        }
223    }
224  inexact = mpfr_mul_2si (r, t, shift_exp, rnd_mode);
225  mpfr_clear (t);
226  mpfr_clear (u);
227
228  MPFR_SAVE_EXPO_FREE (expo);
229  return mpfr_check_range (r, inexact, rnd_mode);
230
231 nan:
232  MPFR_SET_NAN(r);
233  MPFR_RET_NAN;
234}
235
236#else /* IEEE Extended Little Endian Code */
237
238int
239mpfr_set_ld (mpfr_ptr r, long double d, mpfr_rnd_t rnd_mode)
240{
241  int inexact, i, k, cnt;
242  mpfr_t tmp;
243  mp_limb_t tmpmant[MPFR_LIMBS_PER_LONG_DOUBLE];
244  mpfr_long_double_t x;
245  mpfr_exp_t exp;
246  int signd;
247  MPFR_SAVE_EXPO_DECL (expo);
248
249  /* Check for NAN */
250  if (MPFR_UNLIKELY (d != d))
251    {
252      MPFR_SET_NAN (r);
253      MPFR_RET_NAN;
254    }
255  /* Check for INF */
256  else if (MPFR_UNLIKELY (d > MPFR_LDBL_MAX))
257    {
258      MPFR_SET_INF (r);
259      MPFR_SET_POS (r);
260      return 0;
261    }
262  else if (MPFR_UNLIKELY (d < -MPFR_LDBL_MAX))
263    {
264      MPFR_SET_INF (r);
265      MPFR_SET_NEG (r);
266      return 0;
267    }
268  /* Check for ZERO */
269  else if (MPFR_UNLIKELY (d == 0.0))
270    {
271      x.ld = d;
272      MPFR_SET_ZERO (r);
273      if (x.s.sign == 1)
274        MPFR_SET_NEG(r);
275      else
276        MPFR_SET_POS(r);
277      return 0;
278    }
279
280  /* now d is neither 0, nor NaN nor Inf */
281  MPFR_SAVE_EXPO_MARK (expo);
282
283  MPFR_MANT (tmp) = tmpmant;
284  MPFR_PREC (tmp) = 64;
285
286  /* Extract sign */
287  x.ld = d;
288  signd = MPFR_SIGN_POS;
289  if (x.ld < 0.0)
290    {
291      signd = MPFR_SIGN_NEG;
292      x.ld = -x.ld;
293    }
294
295  /* Extract mantissa */
296#if GMP_NUMB_BITS >= 64
297  tmpmant[0] = ((mp_limb_t) x.s.manh << 32) | ((mp_limb_t) x.s.manl);
298#else
299  tmpmant[0] = (mp_limb_t) x.s.manl;
300  tmpmant[1] = (mp_limb_t) x.s.manh;
301#endif
302
303  /* Normalize mantissa */
304  i = MPFR_LIMBS_PER_LONG_DOUBLE;
305  MPN_NORMALIZE_NOT_ZERO (tmpmant, i);
306  k = MPFR_LIMBS_PER_LONG_DOUBLE - i;
307  count_leading_zeros (cnt, tmpmant[i - 1]);
308  if (MPFR_LIKELY (cnt != 0))
309    mpn_lshift (tmpmant + k, tmpmant, i, cnt);
310  else if (k != 0)
311    MPN_COPY (tmpmant + k, tmpmant, i);
312  if (MPFR_UNLIKELY (k != 0))
313    MPN_ZERO (tmpmant, k);
314
315  /* Set exponent */
316  exp = (mpfr_exp_t) ((x.s.exph << 8) + x.s.expl);  /* 15-bit unsigned int */
317  if (MPFR_UNLIKELY (exp == 0))
318    exp -= 0x3FFD;
319  else
320    exp -= 0x3FFE;
321
322  MPFR_SET_EXP (tmp, exp - cnt - k * GMP_NUMB_BITS);
323
324  /* tmp is exact */
325  inexact = mpfr_set4 (r, tmp, rnd_mode, signd);
326
327  MPFR_SAVE_EXPO_FREE (expo);
328  return mpfr_check_range (r, inexact, rnd_mode);
329}
330
331#endif
332