1 /* mpfr_set_ld -- convert a machine long double to
2 a multiple precision floating-point number
3
4 Copyright 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
5 Contributed by the AriC and Caramel 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 http://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 /* Various i386 systems have been seen with <float.h> LDBL constants equal
30 to the DBL ones, whereas they ought to be bigger, reflecting the 10-byte
31 IEEE extended format on that processor. gcc 3.2.1 on FreeBSD and Solaris
32 has been seen with the problem, and gcc 2.95.4 on FreeBSD 4.7. */
33
34 #if HAVE_LDOUBLE_IEEE_EXT_LITTLE
35 static const union {
36 char bytes[10];
37 long double d;
38 } ldbl_max_struct = {
39 { '\377','\377','\377','\377',
40 '\377','\377','\377','\377',
41 '\376','\177' }
42 };
43 #define MPFR_LDBL_MAX (ldbl_max_struct.d)
44 #else
45 #define MPFR_LDBL_MAX LDBL_MAX
46 #endif
47
48 #ifndef HAVE_LDOUBLE_IEEE_EXT_LITTLE
49
50 /* Generic code */
51 int
mpfr_set_ld(mpfr_ptr r,long double d,mpfr_rnd_t rnd_mode)52 mpfr_set_ld (mpfr_ptr r, long double d, mpfr_rnd_t rnd_mode)
53 {
54 mpfr_t t, u;
55 int inexact, shift_exp;
56 long double x;
57 MPFR_SAVE_EXPO_DECL (expo);
58
59 /* Check for NAN */
60 LONGDOUBLE_NAN_ACTION (d, goto nan);
61
62 /* Check for INF */
63 if (d > MPFR_LDBL_MAX)
64 {
65 mpfr_set_inf (r, 1);
66 return 0;
67 }
68 else if (d < -MPFR_LDBL_MAX)
69 {
70 mpfr_set_inf (r, -1);
71 return 0;
72 }
73 /* Check for ZERO */
74 else if (d == 0.0)
75 return mpfr_set_d (r, (double) d, rnd_mode);
76
77 mpfr_init2 (t, MPFR_LDBL_MANT_DIG);
78 mpfr_init2 (u, IEEE_DBL_MANT_DIG);
79
80 MPFR_SAVE_EXPO_MARK (expo);
81
82 convert:
83 x = d;
84 MPFR_SET_ZERO (t); /* The sign doesn't matter. */
85 shift_exp = 0; /* invariant: remainder to deal with is d*2^shift_exp */
86 while (x != (long double) 0.0)
87 {
88 /* Check overflow of double */
89 if (x > (long double) DBL_MAX || (-x) > (long double) DBL_MAX)
90 {
91 long double div9, div10, div11, div12, div13;
92
93 #define TWO_64 18446744073709551616.0 /* 2^64 */
94 #define TWO_128 (TWO_64 * TWO_64)
95 #define TWO_256 (TWO_128 * TWO_128)
96 div9 = (long double) (double) (TWO_256 * TWO_256); /* 2^(2^9) */
97 div10 = div9 * div9;
98 div11 = div10 * div10; /* 2^(2^11) */
99 div12 = div11 * div11; /* 2^(2^12) */
100 div13 = div12 * div12; /* 2^(2^13) */
101 if (ABS (x) >= div13)
102 {
103 x /= div13; /* exact */
104 shift_exp += 8192;
105 mpfr_div_2si (t, t, 8192, MPFR_RNDZ);
106 }
107 if (ABS (x) >= div12)
108 {
109 x /= div12; /* exact */
110 shift_exp += 4096;
111 mpfr_div_2si (t, t, 4096, MPFR_RNDZ);
112 }
113 if (ABS (x) >= div11)
114 {
115 x /= div11; /* exact */
116 shift_exp += 2048;
117 mpfr_div_2si (t, t, 2048, MPFR_RNDZ);
118 }
119 if (ABS (x) >= div10)
120 {
121 x /= div10; /* exact */
122 shift_exp += 1024;
123 mpfr_div_2si (t, t, 1024, MPFR_RNDZ);
124 }
125 /* warning: we may have DBL_MAX=2^1024*(1-2^(-53)) < x < 2^1024,
126 therefore we have one extra exponent reduction step */
127 if (ABS (x) >= div9)
128 {
129 x /= div9; /* exact */
130 shift_exp += 512;
131 mpfr_div_2si (t, t, 512, MPFR_RNDZ);
132 }
133 } /* Check overflow of double */
134 else /* no overflow on double */
135 {
136 long double div9, div10, div11;
137
138 div9 = (long double) (double) 7.4583407312002067432909653e-155;
139 /* div9 = 2^(-2^9) */
140 div10 = div9 * div9; /* 2^(-2^10) */
141 div11 = div10 * div10; /* 2^(-2^11) if extended precision */
142 /* since -DBL_MAX <= x <= DBL_MAX, the cast to double should not
143 overflow here */
144 if (ABS(x) < div10 &&
145 div11 != (long double) 0.0 &&
146 div11 / div10 == div10) /* possible underflow */
147 {
148 long double div12, div13;
149 /* After the divisions, any bit of x must be >= div10,
150 hence the possible division by div9. */
151 div12 = div11 * div11; /* 2^(-2^12) */
152 div13 = div12 * div12; /* 2^(-2^13) */
153 if (ABS (x) <= div13)
154 {
155 x /= div13; /* exact */
156 shift_exp -= 8192;
157 mpfr_mul_2si (t, t, 8192, MPFR_RNDZ);
158 }
159 if (ABS (x) <= div12)
160 {
161 x /= div12; /* exact */
162 shift_exp -= 4096;
163 mpfr_mul_2si (t, t, 4096, MPFR_RNDZ);
164 }
165 if (ABS (x) <= div11)
166 {
167 x /= div11; /* exact */
168 shift_exp -= 2048;
169 mpfr_mul_2si (t, t, 2048, MPFR_RNDZ);
170 }
171 if (ABS (x) <= div10)
172 {
173 x /= div10; /* exact */
174 shift_exp -= 1024;
175 mpfr_mul_2si (t, t, 1024, MPFR_RNDZ);
176 }
177 if (ABS(x) <= div9)
178 {
179 x /= div9; /* exact */
180 shift_exp -= 512;
181 mpfr_mul_2si (t, t, 512, MPFR_RNDZ);
182 }
183 }
184 else /* no underflow */
185 {
186 inexact = mpfr_set_d (u, (double) x, MPFR_RNDZ);
187 MPFR_ASSERTD (inexact == 0);
188 if (mpfr_add (t, t, u, MPFR_RNDZ) != 0)
189 {
190 if (!mpfr_number_p (t))
191 break;
192 /* Inexact. This cannot happen unless the C implementation
193 "lies" on the precision or when long doubles are
194 implemented with FP expansions like under Mac OS X. */
195 if (MPFR_PREC (t) != MPFR_PREC (r) + 1)
196 {
197 /* We assume that MPFR_PREC (r) < MPFR_PREC_MAX.
198 The precision MPFR_PREC (r) + 1 allows us to
199 deduce the rounding bit and the sticky bit. */
200 mpfr_set_prec (t, MPFR_PREC (r) + 1);
201 goto convert;
202 }
203 else
204 {
205 mp_limb_t *tp;
206 int rb_mask;
207
208 /* Since mpfr_add was inexact, the sticky bit is 1. */
209 tp = MPFR_MANT (t);
210 rb_mask = MPFR_LIMB_ONE <<
211 (GMP_NUMB_BITS - 1 -
212 (MPFR_PREC (r) & (GMP_NUMB_BITS - 1)));
213 if (rnd_mode == MPFR_RNDN)
214 rnd_mode = (*tp & rb_mask) ^ MPFR_IS_NEG (t) ?
215 MPFR_RNDU : MPFR_RNDD;
216 *tp |= rb_mask;
217 break;
218 }
219 }
220 x -= (long double) mpfr_get_d1 (u); /* exact */
221 }
222 }
223 }
224 inexact = mpfr_mul_2si (r, t, shift_exp, rnd_mode);
225 mpfr_clear (t);
226 mpfr_clear (u);
227
228 MPFR_SAVE_EXPO_FREE (expo);
229 return mpfr_check_range (r, inexact, rnd_mode);
230
231 nan:
232 MPFR_SET_NAN(r);
233 MPFR_RET_NAN;
234 }
235
236 #else /* IEEE Extended Little Endian Code */
237
238 int
mpfr_set_ld(mpfr_ptr r,long double d,mpfr_rnd_t rnd_mode)239 mpfr_set_ld (mpfr_ptr r, long double d, mpfr_rnd_t rnd_mode)
240 {
241 int inexact, i, k, cnt;
242 mpfr_t tmp;
243 mp_limb_t tmpmant[MPFR_LIMBS_PER_LONG_DOUBLE];
244 mpfr_long_double_t x;
245 mpfr_exp_t exp;
246 int signd;
247 MPFR_SAVE_EXPO_DECL (expo);
248
249 /* Check for NAN */
250 if (MPFR_UNLIKELY (d != d))
251 {
252 MPFR_SET_NAN (r);
253 MPFR_RET_NAN;
254 }
255 /* Check for INF */
256 else if (MPFR_UNLIKELY (d > MPFR_LDBL_MAX))
257 {
258 MPFR_SET_INF (r);
259 MPFR_SET_POS (r);
260 return 0;
261 }
262 else if (MPFR_UNLIKELY (d < -MPFR_LDBL_MAX))
263 {
264 MPFR_SET_INF (r);
265 MPFR_SET_NEG (r);
266 return 0;
267 }
268 /* Check for ZERO */
269 else if (MPFR_UNLIKELY (d == 0.0))
270 {
271 x.ld = d;
272 MPFR_SET_ZERO (r);
273 if (x.s.sign == 1)
274 MPFR_SET_NEG(r);
275 else
276 MPFR_SET_POS(r);
277 return 0;
278 }
279
280 /* now d is neither 0, nor NaN nor Inf */
281 MPFR_SAVE_EXPO_MARK (expo);
282
283 MPFR_MANT (tmp) = tmpmant;
284 MPFR_PREC (tmp) = 64;
285
286 /* Extract sign */
287 x.ld = d;
288 signd = MPFR_SIGN_POS;
289 if (x.ld < 0.0)
290 {
291 signd = MPFR_SIGN_NEG;
292 x.ld = -x.ld;
293 }
294
295 /* Extract mantissa */
296 #if GMP_NUMB_BITS >= 64
297 tmpmant[0] = ((mp_limb_t) x.s.manh << 32) | ((mp_limb_t) x.s.manl);
298 #else
299 tmpmant[0] = (mp_limb_t) x.s.manl;
300 tmpmant[1] = (mp_limb_t) x.s.manh;
301 #endif
302
303 /* Normalize mantissa */
304 i = MPFR_LIMBS_PER_LONG_DOUBLE;
305 MPN_NORMALIZE_NOT_ZERO (tmpmant, i);
306 k = MPFR_LIMBS_PER_LONG_DOUBLE - i;
307 count_leading_zeros (cnt, tmpmant[i - 1]);
308 if (MPFR_LIKELY (cnt != 0))
309 mpn_lshift (tmpmant + k, tmpmant, i, cnt);
310 else if (k != 0)
311 MPN_COPY (tmpmant + k, tmpmant, i);
312 if (MPFR_UNLIKELY (k != 0))
313 MPN_ZERO (tmpmant, k);
314
315 /* Set exponent */
316 exp = (mpfr_exp_t) ((x.s.exph << 8) + x.s.expl); /* 15-bit unsigned int */
317 if (MPFR_UNLIKELY (exp == 0))
318 exp -= 0x3FFD;
319 else
320 exp -= 0x3FFE;
321
322 MPFR_SET_EXP (tmp, exp - cnt - k * GMP_NUMB_BITS);
323
324 /* tmp is exact */
325 inexact = mpfr_set4 (r, tmp, rnd_mode, signd);
326
327 MPFR_SAVE_EXPO_FREE (expo);
328 return mpfr_check_range (r, inexact, rnd_mode);
329 }
330
331 #endif
332