1 /* mpfr_set_d -- convert a machine double precision float to
2 a multiple precision floating-point number
3
4 Copyright 1999-2004, 2006-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> /* For DOUBLE_ISINF and DOUBLE_ISNAN */
25
26 #define MPFR_NEED_LONGLONG_H
27 #include "mpfr-impl.h"
28
29 /* Extracts the bits of |d| in rp[0..n-1] where n=ceil(53/GMP_NUMB_BITS).
30 Assumes d finite and <> 0.
31 Returns the corresponding exponent such that |d| = {rp, n} * 2^exp,
32 with the value of {rp, n} in [1/2, 1).
33 The int type should be sufficient for exp.
34 */
35 static int
36 extract_double (mpfr_limb_ptr rp, double d)
37 {
38 int exp;
39 mp_limb_t man[MPFR_LIMBS_PER_DOUBLE];
40
41 /* FIXME: Generalize to handle GMP_NUMB_BITS < 16. */
42
43 MPFR_ASSERTD(!DOUBLE_ISNAN(d));
44 MPFR_ASSERTD(!DOUBLE_ISINF(d));
45 MPFR_ASSERTD(d != 0.0);
46
47 #if _MPFR_IEEE_FLOATS
48
49 {
50 union mpfr_ieee_double_extract x;
51 x.d = d;
52
53 exp = x.s.exp;
54 if (exp)
55 {
56 /* x.s.manh has 20 bits (in its low bits), x.s.manl has 32 bits */
57 #if GMP_NUMB_BITS >= 64
58 man[0] = ((MPFR_LIMB_ONE << (GMP_NUMB_BITS - 1)) |
59 ((mp_limb_t) x.s.manh << (GMP_NUMB_BITS - 21)) |
60 ((mp_limb_t) x.s.manl << (GMP_NUMB_BITS - 53)));
61 #elif GMP_NUMB_BITS == 32
62 man[1] = (MPFR_LIMB_ONE << 31) | (x.s.manh << 11) | (x.s.manl >> 21);
63 man[0] = x.s.manl << 11;
64 #elif GMP_NUMB_BITS == 16
65 MPFR_STAT_STATIC_ASSERT (GMP_NUMB_BITS == 16);
66 man[3] = (MPFR_LIMB_ONE << 15) | (x.s.manh >> 5);
67 man[2] = (x.s.manh << 11) | (x.s.manl >> 21);
68 man[1] = x.s.manl >> 5;
69 man[0] = MPFR_LIMB_LSHIFT(x.s.manl,11);
70 #else
71 MPFR_STAT_STATIC_ASSERT (GMP_NUMB_BITS == 8);
72 man[6] = (MPFR_LIMB_ONE << 7) | (x.s.manh >> 13);
73 man[5] = (mp_limb_t) (x.s.manh >> 5);
74 man[4] = MPFR_LIMB_LSHIFT(x.s.manh, 3) | (mp_limb_t) (x.s.manl >> 29);
75 man[3] = (mp_limb_t) (x.s.manl >> 21);
76 man[2] = (mp_limb_t) (x.s.manl >> 13);
77 man[1] = (mp_limb_t) (x.s.manl >> 5);
78 man[0] = MPFR_LIMB_LSHIFT(x.s.manl,3);
79 #endif
80 exp -= 1022;
81 }
82 else /* subnormal number */
83 {
84 int cnt;
85 exp = -1021;
86 #if GMP_NUMB_BITS >= 64
87 man[0] = (((mp_limb_t) x.s.manh << (GMP_NUMB_BITS - 21)) |
88 ((mp_limb_t) x.s.manl << (GMP_NUMB_BITS - 53)));
89 count_leading_zeros (cnt, man[0]);
90 #elif GMP_NUMB_BITS == 32
91 man[1] = (x.s.manh << 11) /* high 21 bits */
92 | (x.s.manl >> 21); /* middle 11 bits */
93 man[0] = x.s.manl << 11; /* low 21 bits */
94 if (man[1] == 0)
95 {
96 man[1] = man[0];
97 man[0] = 0;
98 exp -= GMP_NUMB_BITS;
99 }
100 count_leading_zeros (cnt, man[1]);
101 man[1] = (man[1] << cnt) |
102 (cnt != 0 ? man[0] >> (GMP_NUMB_BITS - cnt) : 0);
103 #elif GMP_NUMB_BITS == 16
104 man[3] = x.s.manh >> 5;
105 man[2] = (x.s.manh << 11) | (x.s.manl >> 21);
106 man[1] = x.s.manl >> 5;
107 man[0] = x.s.manl << 11;
108 while (man[3] == 0) /* d is assumed <> 0 */
109 {
110 man[3] = man[2];
111 man[2] = man[1];
112 man[1] = man[0];
113 man[0] = 0;
114 exp -= GMP_NUMB_BITS;
115 }
116 count_leading_zeros (cnt, man[3]);
117 if (cnt)
118 {
119 man[3] = (man[3] << cnt) | (man[2] >> (GMP_NUMB_BITS - cnt));
120 man[2] = (man[2] << cnt) | (man[1] >> (GMP_NUMB_BITS - cnt));
121 man[1] = (man[1] << cnt) | (man[0] >> (GMP_NUMB_BITS - cnt));
122 }
123 #else
124 MPFR_STAT_STATIC_ASSERT (GMP_NUMB_BITS == 8);
125 man[6] = x.s.manh >> 13;
126 man[5] = x.s.manh >> 5;
127 man[4] = (x.s.manh << 3) | (x.s.manl >> 29);
128 man[3] = x.s.manl >> 21;
129 man[2] = x.s.manl >> 13;
130 man[1] = x.s.manl >> 5;
131 man[0] = x.s.manl << 3;
132 while (man[6] == 0) /* d is assumed <> 0 */
133 {
134 man[6] = man[5];
135 man[5] = man[4];
136 man[4] = man[3];
137 man[3] = man[2];
138 man[2] = man[1];
139 man[1] = man[0];
140 man[0] = 0;
141 exp -= GMP_NUMB_BITS;
142 }
143 count_leading_zeros (cnt, man[6]);
144 if (cnt)
145 {
146 int i;
147 for (i = 6; i > 0; i--)
148 man[i] = (man[i] << cnt) | (man[i-1] >> (GMP_NUMB_BITS - cnt));
149 }
150 #endif
151 man[0] <<= cnt;
152 exp -= cnt;
153 }
154 }
155
156 #else /* _MPFR_IEEE_FLOATS */
157
158 {
159 /* Unknown (or known to be non-IEEE) double format. */
160 exp = 0;
161 d = ABS (d);
162 if (d >= 1.0)
163 {
164 while (d >= 32768.0)
165 {
166 d *= (1.0 / 65536.0);
167 exp += 16;
168 }
169 while (d >= 1.0)
170 {
171 d *= 0.5;
172 exp += 1;
173 }
174 }
175 else if (d < 0.5)
176 {
177 while (d < (1.0 / 65536.0))
178 {
179 d *= 65536.0;
180 exp -= 16;
181 }
182 while (d < 0.5)
183 {
184 d *= 2.0;
185 exp -= 1;
186 }
187 }
188
189 d *= MP_BASE_AS_DOUBLE;
190 #if GMP_NUMB_BITS >= 64
191 #ifndef __clang__
192 man[0] = d;
193 #else
194 /* clang up to version 11 produces an invalid exception when d >= 2^63,
195 see <https://github.com/llvm/llvm-project/issues/18060>
196 (old URL: <https://bugs.llvm.org/show_bug.cgi?id=17686>).
197 Since this is always the case, here, we use the following patch. */
198 MPFR_STAT_STATIC_ASSERT (GMP_NUMB_BITS == 64);
199 man[0] = 0x8000000000000000 + (mp_limb_t) (d - 0x8000000000000000);
200 #endif /* __clang__ */
201 #elif GMP_NUMB_BITS == 32
202 man[1] = (mp_limb_t) d;
203 man[0] = (mp_limb_t) ((d - man[1]) * MP_BASE_AS_DOUBLE);
204 #else
205 MPFR_STAT_STATIC_ASSERT (GMP_NUMB_BITS == 16);
206 {
207 double r = d;
208 man[3] = (mp_limb_t) r;
209 r = (r - man[3]) * MP_BASE_AS_DOUBLE;
210 man[2] = (mp_limb_t) r;
211 r = (r - man[2]) * MP_BASE_AS_DOUBLE;
212 man[1] = (mp_limb_t) r;
213 r = (r - man[1]) * MP_BASE_AS_DOUBLE;
214 man[0] = (mp_limb_t) r;
215 }
216 #endif
217 }
218
219 #endif /* _MPFR_IEEE_FLOATS */
220
221 rp[0] = man[0];
222 #if GMP_NUMB_BITS <= 32
223 rp[1] = man[1];
224 #endif
225 #if GMP_NUMB_BITS <= 16
226 rp[2] = man[2];
227 rp[3] = man[3];
228 #endif
229 #if GMP_NUMB_BITS <= 8
230 rp[4] = man[4];
231 rp[5] = man[5];
232 rp[6] = man[6];
233 #endif
234
235 MPFR_ASSERTD((rp[MPFR_LIMBS_PER_DOUBLE - 1] & MPFR_LIMB_HIGHBIT) != 0);
236
237 return exp;
238 }
239
240 /* End of part included from gmp-2.0.2 */
241
242 int
243 mpfr_set_d (mpfr_ptr r, double d, mpfr_rnd_t rnd_mode)
244 {
245 int inexact;
246 mpfr_t tmp;
247 mp_limb_t tmpmant[MPFR_LIMBS_PER_DOUBLE];
248 MPFR_SAVE_EXPO_DECL (expo);
249
250 if (MPFR_UNLIKELY(DOUBLE_ISNAN(d)))
251 {
252 /* we don't propagate the sign bit */
253 MPFR_SET_NAN(r);
254 MPFR_RET_NAN;
255 }
256 else if (MPFR_UNLIKELY(d == 0))
257 {
258 #if _MPFR_IEEE_FLOATS
259 union mpfr_ieee_double_extract x;
260
261 MPFR_SET_ZERO(r);
262 /* set correct sign */
263 x.d = d;
264 if (x.s.sig == 1)
265 MPFR_SET_NEG(r);
266 else
267 MPFR_SET_POS(r);
268 #else /* _MPFR_IEEE_FLOATS */
269 MPFR_SET_ZERO(r);
270 {
271 /* This is to get the sign of zero on non-IEEE hardware
272 Some systems support +0.0, -0.0, and unsigned zero.
273 Some other systems may just have an unsigned zero.
274 We can't use d == +0.0 since it should be always true,
275 so we check that the memory representation of d is the
276 same as +0.0, etc.
277 Note: r is set to -0 only if d is detected as a negative zero
278 *and*, for the double type, -0 has a different representation
279 from +0. If -0.0 has several representations, the code below
280 may not work as expected, but this is hardly fixable in a
281 portable way (without depending on a math library) and only
282 the sign could be incorrect. Such systems should be taken
283 into account on a case-by-case basis. If the code is changed
284 here, set_d64.c code should be updated too. */
285 double poszero = +0.0, negzero = DBL_NEG_ZERO;
286 if (memcmp(&d, &poszero, sizeof(double)) == 0)
287 MPFR_SET_POS(r);
288 else if (memcmp(&d, &negzero, sizeof(double)) == 0)
289 MPFR_SET_NEG(r);
290 else
291 MPFR_SET_POS(r);
292 }
293 #endif /* _MPFR_IEEE_FLOATS */
294 return 0; /* 0 is exact */
295 }
296 else if (MPFR_UNLIKELY(DOUBLE_ISINF(d)))
297 {
298 MPFR_SET_INF(r);
299 if (d > 0)
300 MPFR_SET_POS(r);
301 else
302 MPFR_SET_NEG(r);
303 return 0; /* infinity is exact */
304 }
305
306 /* now d is neither 0, nor NaN nor Inf */
307
308 MPFR_SAVE_EXPO_MARK (expo);
309
310 /* warning: don't use tmp=r here, even if SIZE(r) >= MPFR_LIMBS_PER_DOUBLE,
311 since PREC(r) may be different from PREC(tmp), and then both variables
312 would have same precision in the mpfr_set4 call below. */
313 MPFR_MANT(tmp) = tmpmant;
314 MPFR_PREC(tmp) = IEEE_DBL_MANT_DIG;
315
316 /* don't use MPFR_SET_EXP here since the exponent may be out of range */
317 MPFR_EXP(tmp) = extract_double (tmpmant, d);
318
319 /* tmp is exact since PREC(tmp)=53 */
320 inexact = mpfr_set4 (r, tmp, rnd_mode,
321 (d < 0) ? MPFR_SIGN_NEG : MPFR_SIGN_POS);
322
323 MPFR_SAVE_EXPO_FREE (expo);
324 return mpfr_check_range (r, inexact, rnd_mode);
325 }
326
327
328