1 /* mpfr_get_d, mpfr_get_d_2exp -- convert a multiple precision floating-point
2 number to a machine double precision float
3
4 Copyright 1999-2023 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramba projects, INRIA.
6
7 This file is part of the GNU MPFR Library.
8
9 The GNU MPFR Library is free software; you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
13
14 The GNU MPFR Library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 License for more details.
18
19 You should have received a copy of the GNU Lesser General Public License
20 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
21 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
22 51 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
35 double
mpfr_get_d(mpfr_srcptr src,mpfr_rnd_t rnd_mode)36 mpfr_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
136 double
mpfr_get_d1(mpfr_srcptr src)137 mpfr_get_d1 (mpfr_srcptr src)
138 {
139 return mpfr_get_d (src, __gmpfr_default_rounding_mode);
140 }
141
142 double
mpfr_get_d_2exp(long * expptr,mpfr_srcptr src,mpfr_rnd_t rnd_mode)143 mpfr_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