1/* mpfr_get_d, mpfr_get_d_2exp -- convert a multiple precision floating-point
2                                  number to a machine double precision float
3
4Copyright 1999, 2000, 2001, 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#include "ieee_floats.h"
30
31/* Assumes IEEE-754 double precision; otherwise, only an approximated
32   result will be returned, without any guaranty (and special cases
33   such as NaN must be avoided if not supported). */
34
35double
36mpfr_get_d (mpfr_srcptr src, mpfr_rnd_t rnd_mode)
37{
38  double d;
39  int negative;
40  mpfr_exp_t e;
41
42  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
43    {
44      if (MPFR_IS_NAN (src))
45        return MPFR_DBL_NAN;
46
47      negative = MPFR_IS_NEG (src);
48
49      if (MPFR_IS_INF (src))
50        return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
51
52      MPFR_ASSERTD (MPFR_IS_ZERO(src));
53      return negative ? DBL_NEG_ZERO : 0.0;
54    }
55
56  e = MPFR_GET_EXP (src);
57  negative = MPFR_IS_NEG (src);
58
59  if (MPFR_UNLIKELY(rnd_mode == MPFR_RNDA))
60    rnd_mode = negative ? MPFR_RNDD : MPFR_RNDU;
61
62  /* the smallest normalized number is 2^(-1022)=0.1e-1021, and the smallest
63     subnormal is 2^(-1074)=0.1e-1073 */
64  if (MPFR_UNLIKELY (e < -1073))
65    {
66      /* Note: Avoid using a constant expression DBL_MIN * DBL_EPSILON
67         as this gives 0 instead of the correct result with gcc on some
68         Alpha machines. */
69      d = negative ?
70        (rnd_mode == MPFR_RNDD ||
71         (rnd_mode == MPFR_RNDN && mpfr_cmp_si_2exp(src, -1, -1075) < 0)
72         ? -DBL_MIN : DBL_NEG_ZERO) :
73        (rnd_mode == MPFR_RNDU ||
74         (rnd_mode == MPFR_RNDN && mpfr_cmp_si_2exp(src, 1, -1075) > 0)
75         ? DBL_MIN : 0.0);
76      if (d != 0.0) /* we multiply DBL_MIN = 2^(-1022) by DBL_EPSILON = 2^(-52)
77                       to get +-2^(-1074) */
78        d *= DBL_EPSILON;
79    }
80  /* the largest normalized number is 2^1024*(1-2^(-53))=0.111...111e1024 */
81  else if (MPFR_UNLIKELY (e > 1024))
82    {
83      d = negative ?
84        (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDU ?
85         -DBL_MAX : MPFR_DBL_INFM) :
86        (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDD ?
87         DBL_MAX : MPFR_DBL_INFP);
88    }
89  else
90    {
91      int nbits;
92      mp_size_t np, i;
93      mp_limb_t tp[ MPFR_LIMBS_PER_DOUBLE ];
94      int carry;
95
96      nbits = IEEE_DBL_MANT_DIG; /* 53 */
97      if (MPFR_UNLIKELY (e < -1021))
98        /*In the subnormal case, compute the exact number of significant bits*/
99        {
100          nbits += (1021 + e);
101          MPFR_ASSERTD (nbits >= 1);
102        }
103      np = MPFR_PREC2LIMBS (nbits);
104      MPFR_ASSERTD ( np <= MPFR_LIMBS_PER_DOUBLE );
105      carry = mpfr_round_raw_4 (tp, MPFR_MANT(src), MPFR_PREC(src), negative,
106                                nbits, rnd_mode);
107      if (MPFR_UNLIKELY(carry))
108        d = 1.0;
109      else
110        {
111          /* The following computations are exact thanks to the previous
112             mpfr_round_raw. */
113          d = (double) tp[0] / MP_BASE_AS_DOUBLE;
114          for (i = 1 ; i < np ; i++)
115            d = (d + tp[i]) / MP_BASE_AS_DOUBLE;
116          /* d is the mantissa (between 1/2 and 1) of the argument rounded
117             to 53 bits */
118        }
119      d = mpfr_scale2 (d, e);
120      if (negative)
121        d = -d;
122    }
123
124  return d;
125}
126
127#undef mpfr_get_d1
128double
129mpfr_get_d1 (mpfr_srcptr src)
130{
131  return mpfr_get_d (src, __gmpfr_default_rounding_mode);
132}
133
134double
135mpfr_get_d_2exp (long *expptr, mpfr_srcptr src, mpfr_rnd_t rnd_mode)
136{
137  double ret;
138  mpfr_exp_t exp;
139  mpfr_t tmp;
140
141  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
142    {
143      int negative;
144      *expptr = 0;
145      if (MPFR_IS_NAN (src))
146        return MPFR_DBL_NAN;
147      negative = MPFR_IS_NEG (src);
148      if (MPFR_IS_INF (src))
149        return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
150      MPFR_ASSERTD (MPFR_IS_ZERO(src));
151      return negative ? DBL_NEG_ZERO : 0.0;
152    }
153
154  tmp[0] = *src;        /* Hack copy mpfr_t */
155  MPFR_SET_EXP (tmp, 0);
156  ret = mpfr_get_d (tmp, rnd_mode);
157
158  if (MPFR_IS_PURE_FP(src))
159    {
160      exp = MPFR_GET_EXP (src);
161
162      /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */
163      if (ret == 1.0)
164        {
165          ret = 0.5;
166          exp++;
167        }
168      else if (ret == -1.0)
169        {
170          ret = -0.5;
171          exp++;
172        }
173
174      MPFR_ASSERTN ((ret >= 0.5 && ret < 1.0)
175                    || (ret <= -0.5 && ret > -1.0));
176      MPFR_ASSERTN (exp >= LONG_MIN && exp <= LONG_MAX);
177    }
178  else
179    exp = 0;
180
181  *expptr = exp;
182  return ret;
183}
184