1 /* mpfr_get_ld, mpfr_get_ld_2exp -- convert a multiple precision floating-point
2 number to a machine long double
3
4 Copyright 2002-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> /* needed so that MPFR_LDBL_MANT_DIG is correctly defined */
25
26 #include "mpfr-impl.h"
27
28 #if defined(HAVE_LDOUBLE_IS_DOUBLE)
29
30 /* special code when "long double" is the same format as "double" */
31 long double
32 mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
33 {
34 return (long double) mpfr_get_d (x, rnd_mode);
35 }
36
37 #elif defined(HAVE_LDOUBLE_IEEE_EXT_LITTLE)
38
39 /* Note: The code will return a result with a 64-bit precision, even
40 if the rounding precision is only 53 bits like on FreeBSD and
41 NetBSD 6- (or with GCC's -mpc64 option to simulate this on other
42 platforms). This is consistent with how strtold behaves in these
43 cases, for instance. */
44
45 /* special code for IEEE 754 little-endian extended format */
46 long double
47 mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
48 {
49 mpfr_long_double_t ld;
50 mpfr_t tmp;
51 int inex;
52 MPFR_SAVE_EXPO_DECL (expo);
53
54 MPFR_SAVE_EXPO_MARK (expo);
55
56 mpfr_init2 (tmp, MPFR_LDBL_MANT_DIG);
57 inex = mpfr_set (tmp, x, rnd_mode);
58
59 mpfr_set_emin (-16381-63); /* emin=-16444, see below */
60 mpfr_set_emax (16384);
61 mpfr_subnormalize (tmp, mpfr_check_range (tmp, inex, rnd_mode), rnd_mode);
62 mpfr_prec_round (tmp, 64, MPFR_RNDZ); /* exact */
63 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (tmp)))
64 ld.ld = (long double) mpfr_get_d (tmp, rnd_mode);
65 else
66 {
67 mp_limb_t *tmpmant;
68 mpfr_exp_t e, denorm;
69
70 tmpmant = MPFR_MANT (tmp);
71 e = MPFR_GET_EXP (tmp);
72 /* The smallest positive normal number is 2^(-16382), which is
73 0.5*2^(-16381) in MPFR, thus any exponent <= -16382 corresponds to a
74 subnormal number. The smallest positive subnormal number is 2^(-16445)
75 which is 0.5*2^(-16444) in MPFR thus 0 <= denorm <= 63. */
76 denorm = MPFR_UNLIKELY (e <= -16382) ? - e - 16382 + 1 : 0;
77 MPFR_ASSERTD (0 <= denorm && denorm < 64);
78 #if GMP_NUMB_BITS >= 64
79 ld.s.manl = (tmpmant[0] >> denorm);
80 ld.s.manh = (tmpmant[0] >> denorm) >> 32;
81 #elif GMP_NUMB_BITS == 32
82 if (MPFR_LIKELY (denorm == 0))
83 {
84 ld.s.manl = tmpmant[0];
85 ld.s.manh = tmpmant[1];
86 }
87 else if (denorm < 32)
88 {
89 ld.s.manl = (tmpmant[0] >> denorm) | (tmpmant[1] << (32 - denorm));
90 ld.s.manh = tmpmant[1] >> denorm;
91 }
92 else /* 32 <= denorm < 64 */
93 {
94 ld.s.manl = tmpmant[1] >> (denorm - 32);
95 ld.s.manh = 0;
96 }
97 #elif GMP_NUMB_BITS == 16
98 if (MPFR_LIKELY (denorm == 0))
99 {
100 /* manl = tmpmant[1] | tmpmant[0]
101 manh = tmpmant[3] | tmpmant[2] */
102 ld.s.manl = tmpmant[0] | ((unsigned long) tmpmant[1] << 16);
103 ld.s.manh = tmpmant[2] | ((unsigned long) tmpmant[3] << 16);
104 }
105 else if (denorm < 16)
106 {
107 /* manl = low(mant[2],denorm) | mant[1] | high(mant[0],16-denorm)
108 manh = mant[3] | high(mant[2],16-denorm) */
109 ld.s.manl = (tmpmant[0] >> denorm)
110 | ((unsigned long) tmpmant[1] << (16 - denorm))
111 | ((unsigned long) tmpmant[2] << (32 - denorm));
112 ld.s.manh = (tmpmant[2] >> denorm)
113 | ((unsigned long) tmpmant[3] << (16 - denorm));
114 }
115 else if (denorm == 16)
116 {
117 /* manl = tmpmant[2] | tmpmant[1]
118 manh = 0000000000 | tmpmant[3] */
119 ld.s.manl = tmpmant[1] | ((unsigned long) tmpmant[2] << 16);
120 ld.s.manh = tmpmant[3];
121 }
122 else if (denorm < 32)
123 {
124 /* manl = low(mant[3],denorm-16) | mant[2] | high(mant[1],32-denorm)
125 manh = high(mant[3],32-denorm) */
126 ld.s.manl = (tmpmant[1] >> (denorm - 16))
127 | ((unsigned long) tmpmant[2] << (32 - denorm))
128 | ((unsigned long) tmpmant[3] << (48 - denorm));
129 ld.s.manh = tmpmant[3] >> (denorm - 16);
130 }
131 else if (denorm == 32)
132 {
133 /* manl = tmpmant[3] | tmpmant[2]
134 manh = 0 */
135 ld.s.manl = tmpmant[2] | ((unsigned long) tmpmant[3] << 16);
136 ld.s.manh = 0;
137 }
138 else if (denorm < 48)
139 {
140 /* manl = zero(denorm-32) | tmpmant[3] | high(tmpmant[2],48-denorm)
141 manh = 0 */
142 ld.s.manl = (tmpmant[2] >> (denorm - 32))
143 | ((unsigned long) tmpmant[3] << (48 - denorm));
144 ld.s.manh = 0;
145 }
146 else /* 48 <= denorm < 64 */
147 {
148 /* we assume a right shift of 0 is identity */
149 ld.s.manl = tmpmant[3] >> (denorm - 48);
150 ld.s.manh = 0;
151 }
152 #elif GMP_NUMB_BITS == 8
153 {
154 unsigned long long mant = 0;
155 int i;
156 for (i = 0; i < 8; i++)
157 mant |= (unsigned long long) tmpmant[i] << (8*i);
158 mant >>= denorm;
159 ld.s.manl = mant;
160 ld.s.manh = mant >> 32;
161 }
162 #else
163 # error "GMP_NUMB_BITS must be 16, 32 or >= 64"
164 /* Other values have never been supported anyway. */
165 #endif
166 if (MPFR_LIKELY (denorm == 0))
167 {
168 ld.s.exph = (e + 0x3FFE) >> 8;
169 ld.s.expl = (e + 0x3FFE);
170 }
171 else
172 ld.s.exph = ld.s.expl = 0;
173 ld.s.sign = MPFR_IS_NEG (x);
174 }
175
176 mpfr_clear (tmp);
177 MPFR_SAVE_EXPO_FREE (expo);
178 return ld.ld;
179 }
180
181 #else
182
183 /* generic code */
184 long double
185 mpfr_get_ld (mpfr_srcptr x, mpfr_rnd_t rnd_mode)
186 {
187 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
188 return (long double) mpfr_get_d (x, rnd_mode);
189 else /* now x is a normal non-zero number */
190 {
191 long double r; /* result */
192 double s; /* part of result */
193 MPFR_SAVE_EXPO_DECL (expo);
194
195 MPFR_SAVE_EXPO_MARK (expo);
196
197 #if defined(HAVE_LDOUBLE_MAYBE_DOUBLE_DOUBLE)
198 if (MPFR_LDBL_MANT_DIG == 106)
199 {
200 /* Assume double-double format (as found with the PowerPC ABI).
201 The generic code below isn't used because numbers with
202 precision > 106 would not be supported. */
203 s = mpfr_get_d (x, MPFR_RNDN); /* high part of x */
204 /* Let's first consider special cases separately. The test for
205 infinity is really needed to avoid a NaN result. The test
206 for NaN is mainly for optimization. The test for 0 is useful
207 to get the correct sign (assuming mpfr_get_d supports signed
208 zeros on the implementation). */
209 if (s == 0 || DOUBLE_ISNAN (s) || DOUBLE_ISINF (s))
210 {
211 /* we don't propagate the sign bit of NaN */
212 r = (long double) s;
213 }
214 else
215 {
216 mpfr_t y, z;
217
218 mpfr_init2 (y, mpfr_get_prec (x));
219 mpfr_init2 (z, IEEE_DBL_MANT_DIG); /* keep the precision small */
220 mpfr_set_d (z, s, MPFR_RNDN); /* exact */
221 mpfr_sub (y, x, z, MPFR_RNDN); /* exact */
222 /* Add the second part of y (in the correct rounding mode). */
223 r = (long double) s + (long double) mpfr_get_d (y, rnd_mode);
224 mpfr_clear (z);
225 mpfr_clear (y);
226 }
227 }
228 else
229 #endif
230 {
231 long double m;
232 mpfr_exp_t sh; /* exponent shift -> x/2^sh is in the double range */
233 mpfr_t y, z;
234 int sign;
235
236 /* First round x to the target long double precision, so that
237 all subsequent operations are exact (this avoids double rounding
238 problems). However, if the format contains numbers that have
239 more precision, MPFR won't be able to generate such numbers. */
240 mpfr_init2 (y, MPFR_LDBL_MANT_DIG);
241 mpfr_init2 (z, MPFR_LDBL_MANT_DIG);
242 /* Note about the precision of z: even though IEEE_DBL_MANT_DIG is
243 sufficient, z has been set to the same precision as y so that
244 the mpfr_sub below calls mpfr_sub1sp, which is faster than the
245 generic subtraction, even in this particular case (from tests
246 done by Patrick Pelissier on a 64-bit Core2 Duo against r7285).
247 But here there is an important cancellation in the subtraction.
248 TODO: get more information about what has been tested. */
249
250 mpfr_set (y, x, rnd_mode);
251 sh = MPFR_GET_EXP (y);
252 sign = MPFR_SIGN (y);
253 MPFR_SET_EXP (y, 0);
254 MPFR_SET_POS (y);
255
256 r = 0.0;
257 do
258 {
259 s = mpfr_get_d (y, MPFR_RNDN); /* high part of y */
260 r += (long double) s;
261 mpfr_set_d (z, s, MPFR_RNDN); /* exact */
262 mpfr_sub (y, y, z, MPFR_RNDN); /* exact */
263 }
264 while (!MPFR_IS_ZERO (y));
265
266 mpfr_clear (z);
267 mpfr_clear (y);
268
269 /* we now have to multiply back by 2^sh */
270 MPFR_ASSERTD (r > 0);
271 if (sh != 0)
272 {
273 /* An overflow may occur (example: 0.5*2^1024) */
274 while (r < 1.0)
275 {
276 r += r;
277 sh--;
278 }
279
280 if (sh > 0)
281 m = 2.0;
282 else
283 {
284 m = 0.5;
285 sh = -sh;
286 }
287
288 for (;;)
289 {
290 if (sh % 2)
291 r = r * m;
292 sh >>= 1;
293 if (sh == 0)
294 break;
295 m = m * m;
296 }
297 }
298 if (sign < 0)
299 r = -r;
300 }
301 MPFR_SAVE_EXPO_FREE (expo);
302 return r;
303 }
304 }
305
306 #endif
307
308 /* contributed by Damien Stehle */
309 long double
310 mpfr_get_ld_2exp (long *expptr, mpfr_srcptr src, mpfr_rnd_t rnd_mode)
311 {
312 long double ret;
313 mpfr_exp_t exp;
314 mpfr_t tmp;
315
316 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (src)))
317 return (long double) mpfr_get_d_2exp (expptr, src, rnd_mode);
318
319 MPFR_ALIAS (tmp, src, MPFR_SIGN (src), 0);
320 ret = mpfr_get_ld (tmp, rnd_mode);
321
322 exp = MPFR_GET_EXP (src);
323
324 /* rounding can give 1.0, adjust back to 0.5 <= abs(ret) < 1.0 */
325 if (ret == 1.0)
326 {
327 ret = 0.5;
328 exp ++;
329 }
330 else if (ret == -1.0)
331 {
332 ret = -0.5;
333 exp ++;
334 }
335
336 MPFR_ASSERTN ((ret >= 0.5 && ret < 1.0)
337 || (ret <= -0.5 && ret > -1.0));
338 MPFR_ASSERTN (exp >= LONG_MIN && exp <= LONG_MAX);
339
340 *expptr = exp;
341 return ret;
342 }