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-2023 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        {
46          /* we don't propagate the sign bit */
47          return MPFR_DBL_NAN;
48        }
49
50      negative = MPFR_IS_NEG (src);
51
52      if (MPFR_IS_INF (src))
53        return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
54
55      MPFR_ASSERTD (MPFR_IS_ZERO(src));
56      return negative ? DBL_NEG_ZERO : 0.0;
57    }
58
59  e = MPFR_GET_EXP (src);
60  negative = MPFR_IS_NEG (src);
61
62  if (MPFR_UNLIKELY(rnd_mode == MPFR_RNDA))
63    rnd_mode = negative ? MPFR_RNDD : MPFR_RNDU;
64
65  /* the smallest normalized number is 2^(-1022)=0.1e-1021, and the smallest
66     subnormal is 2^(-1074)=0.1e-1073 */
67  if (MPFR_UNLIKELY (e < -1073))
68    {
69      /* Note: Avoid using a constant expression DBL_MIN * DBL_EPSILON
70         as this gives 0 instead of the correct result with gcc on some
71         Alpha machines. */
72      d = negative ?
73        (rnd_mode == MPFR_RNDD ||
74         (rnd_mode == MPFR_RNDN && mpfr_cmp_si_2exp(src, -1, -1075) < 0)
75         ? -DBL_MIN : DBL_NEG_ZERO) :
76        (rnd_mode == MPFR_RNDU ||
77         (rnd_mode == MPFR_RNDN && mpfr_cmp_si_2exp(src, 1, -1075) > 0)
78         ? DBL_MIN : 0.0);
79      if (d != 0.0) /* we multiply DBL_MIN = 2^(-1022) by DBL_EPSILON = 2^(-52)
80                       to get +-2^(-1074) */
81        d *= DBL_EPSILON;
82    }
83  /* the largest normalized number is 2^1024*(1-2^(-53))=0.111...111e1024 */
84  else if (MPFR_UNLIKELY (e > 1024))
85    {
86      d = negative ?
87        (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDU ?
88         -DBL_MAX : MPFR_DBL_INFM) :
89        (rnd_mode == MPFR_RNDZ || rnd_mode == MPFR_RNDD ?
90         DBL_MAX : MPFR_DBL_INFP);
91    }
92  else
93    {
94      int nbits;
95      mp_limb_t tp[ MPFR_LIMBS_PER_DOUBLE ];
96      int carry;
97
98      nbits = IEEE_DBL_MANT_DIG; /* 53 */
99      if (MPFR_UNLIKELY (e < -1021))
100        /*In the subnormal case, compute the exact number of significant bits*/
101        {
102          nbits += 1021 + e;
103          MPFR_ASSERTD (1 <= nbits && nbits < IEEE_DBL_MANT_DIG);
104        }
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#if MPFR_LIMBS_PER_DOUBLE == 1
112          d = (double) tp[0] / MP_BASE_AS_DOUBLE;
113#else
114          mp_size_t np, i;
115          MPFR_ASSERTD (nbits <= IEEE_DBL_MANT_DIG);
116          np = MPFR_PREC2LIMBS (nbits);
117          MPFR_ASSERTD ( np <= MPFR_LIMBS_PER_DOUBLE );
118          /* The following computations are exact thanks to the previous
119             mpfr_round_raw. */
120          d = (double) tp[0] / MP_BASE_AS_DOUBLE;
121          for (i = 1 ; i < np ; i++)
122            d = (d + tp[i]) / MP_BASE_AS_DOUBLE;
123          /* d is the mantissa (between 1/2 and 1) of the argument rounded
124             to 53 bits */
125#endif
126        }
127      d = mpfr_scale2 (d, e);
128      if (negative)
129        d = -d;
130    }
131
132  return d;
133}
134
135#undef mpfr_get_d1
136double
137mpfr_get_d1 (mpfr_srcptr src)
138{
139  return mpfr_get_d (src, __gmpfr_default_rounding_mode);
140}
141
142double
143mpfr_get_d_2exp (long *expptr, mpfr_srcptr src, mpfr_rnd_t rnd_mode)
144{
145  double ret;
146  mpfr_exp_t exp;
147  mpfr_t tmp;
148
149  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
150    {
151      int negative;
152      *expptr = 0;
153      if (MPFR_IS_NAN (src))
154        {
155          /* we don't propagate the sign bit */
156          return MPFR_DBL_NAN;
157        }
158      negative = MPFR_IS_NEG (src);
159      if (MPFR_IS_INF (src))
160        return negative ? MPFR_DBL_INFM : MPFR_DBL_INFP;
161      MPFR_ASSERTD (MPFR_IS_ZERO(src));
162      return negative ? DBL_NEG_ZERO : 0.0;
163    }
164
165  MPFR_ALIAS (tmp, src, MPFR_SIGN (src), 0);
166  ret = mpfr_get_d (tmp, rnd_mode);
167
168  exp = MPFR_GET_EXP (src);
169
170  /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */
171  if (ret == 1.0)
172    {
173      ret = 0.5;
174      exp++;
175    }
176  else if (ret == -1.0)
177    {
178      ret = -0.5;
179      exp++;
180    }
181
182  MPFR_ASSERTN ((ret >= 0.5 && ret < 1.0)
183                || (ret <= -0.5 && ret > -1.0));
184  MPFR_ASSERTN (exp >= LONG_MIN && exp <= LONG_MAX);
185
186  *expptr = exp;
187  return ret;
188}
189