get_d.c revision 1.1.1.4 
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-2020 Free Software Foundation, Inc.
5Contributed by the AriC and Caramba 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
21https://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_limb_t tp[ MPFR_LIMBS_PER_DOUBLE ];
93      int carry;
94
95      nbits = IEEE_DBL_MANT_DIG; /* 53 */
96      if (MPFR_UNLIKELY (e < -1021))
97        /*In the subnormal case, compute the exact number of significant bits*/
98        {
99          nbits += 1021 + e;
100          MPFR_ASSERTD (1 <= nbits && nbits < IEEE_DBL_MANT_DIG);
101        }
102      carry = mpfr_round_raw_4 (tp, MPFR_MANT(src), MPFR_PREC(src), negative,
103                                nbits, rnd_mode);
104      if (MPFR_UNLIKELY(carry))
105        d = 1.0;
106      else
107        {
108#if MPFR_LIMBS_PER_DOUBLE == 1
109          d = (double) tp[0] / MP_BASE_AS_DOUBLE;
110#else
111          mp_size_t np, i;
112          MPFR_ASSERTD (nbits <= IEEE_DBL_MANT_DIG);
113          np = MPFR_PREC2LIMBS (nbits);
114          MPFR_ASSERTD ( np <= MPFR_LIMBS_PER_DOUBLE );
115          /* The following computations are exact thanks to the previous
116             mpfr_round_raw. */
117          d = (double) tp[0] / MP_BASE_AS_DOUBLE;
118          for (i = 1 ; i < np ; i++)
119            d = (d + tp[i]) / MP_BASE_AS_DOUBLE;
120          /* d is the mantissa (between 1/2 and 1) of the argument rounded
121             to 53 bits */
122#endif
123        }
124      d = mpfr_scale2 (d, e);
125      if (negative)
126        d = -d;
127    }
128
129  return d;
130}
131
132#undef mpfr_get_d1
133double
134mpfr_get_d1 (mpfr_srcptr src)
135{
136  return mpfr_get_d (src, __gmpfr_default_rounding_mode);
137}
138
139double
140mpfr_get_d_2exp (long *expptr, mpfr_srcptr src, mpfr_rnd_t rnd_mode)
141{
142  double ret;
143  mpfr_exp_t exp;
144  mpfr_t tmp;
145
146  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
147    {
148      int negative;
149      *expptr = 0;
150      if (MPFR_IS_NAN (src))
151        return MPFR_DBL_NAN;
152      negative = MPFR_IS_NEG (src);
153      if (MPFR_IS_INF (src))
154        return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
155      MPFR_ASSERTD (MPFR_IS_ZERO(src));
156      return negative ? DBL_NEG_ZERO : 0.0;
157    }
158
159  MPFR_ALIAS (tmp, src, MPFR_SIGN (src), 0);
160  ret = mpfr_get_d (tmp, rnd_mode);
161
162  exp = MPFR_GET_EXP (src);
163
164  /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */
165  if (ret == 1.0)
166    {
167      ret = 0.5;
168      exp++;
169    }
170  else if (ret == -1.0)
171    {
172      ret = -0.5;
173      exp++;
174    }
175
176  MPFR_ASSERTN ((ret >= 0.5 && ret < 1.0)
177                || (ret <= -0.5 && ret > -1.0));
178  MPFR_ASSERTN (exp >= LONG_MIN && exp <= LONG_MAX);
179
180  *expptr = exp;
181  return ret;
182}
183