1 /* mpn_get_d -- limbs to double conversion.
2
3 THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST
4 CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
5 FUTURE GNU MP RELEASES.
6
7 Copyright 2003, 2004, 2007, 2009 Free Software Foundation, Inc.
8
9 This file is part of the GNU MP Library.
10
11 The GNU MP Library is free software; you can redistribute it and/or modify
12 it under the terms of the GNU Lesser General Public License as published by
13 the Free Software Foundation; either version 3 of the License, or (at your
14 option) any later version.
15
16 The GNU MP Library is distributed in the hope that it will be useful, but
17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19 License for more details.
20
21 You should have received a copy of the GNU Lesser General Public License
22 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
23
24 #include "gmp.h"
25 #include "gmp-impl.h"
26 #include "longlong.h"
27
28 #ifndef _GMP_IEEE_FLOATS
29 #define _GMP_IEEE_FLOATS 0
30 #endif
31
32 #if ! _GMP_IEEE_FLOATS
33 /* dummy definition, just to let dead code compile */
34 union ieee_double_extract {
35 struct {
36 int manh, manl, sig, exp;
37 } s;
38 double d;
39 };
40 #endif
41
42 /* To force use of the generic C code for testing, put
43 "#define _GMP_IEEE_FLOATS 0" at this point. */
44
45
46
47 /* In alpha gcc prior to 3.4, signed DI comparisons involving constants are
48 rearranged from "x < n" to "x+(-n) < 0", which is of course hopelessly
49 wrong if that addition overflows.
50
51 The workaround here avoids this bug by ensuring n is not a literal
52 constant. Note that this is alpha specific. The offending transformation
53 is/was in alpha.c alpha_emit_conditional_branch() under "We want to use
54 cmpcc/bcc".
55
56 Bizarrely, it turns out this happens also with Cray cc on
57 alphaev5-cray-unicosmk2.0.6.X, and has the same solution. Don't know why
58 or how. */
59
60 #if HAVE_HOST_CPU_FAMILY_alpha \
61 && ((defined (__GNUC__) && ! __GMP_GNUC_PREREQ(3,4)) \
62 || defined (_CRAY))
63 static volatile const long CONST_1024 = 1024;
64 static volatile const long CONST_NEG_1023 = -1023;
65 static volatile const long CONST_NEG_1022_SUB_53 = -1022 - 53;
66 #else
67 #define CONST_1024 (1024)
68 #define CONST_NEG_1023 (-1023)
69 #define CONST_NEG_1022_SUB_53 (-1022 - 53)
70 #endif
71
72
73
74 /* Return the value {ptr,size}*2^exp, and negative if sign<0.
75 Must have size>=1, and a non-zero high limb ptr[size-1].
76
77 {ptr,size} is truncated towards zero. This is consistent with other gmp
78 conversions, like mpz_set_f or mpz_set_q, and is easy to implement and
79 test.
80
81 In the past conversions had attempted (imperfectly) to let the hardware
82 float rounding mode take effect, but that gets tricky since multiple
83 roundings need to be avoided, or taken into account, and denorms mean the
84 effective precision of the mantissa is not constant. (For reference,
85 mpz_get_d on IEEE systems was ok, except it operated on the absolute
86 value. mpf_get_d and mpq_get_d suffered from multiple roundings and from
87 not always using enough bits to get the rounding right.)
88
89 It's felt that GMP is not primarily concerned with hardware floats, and
90 really isn't enhanced by getting involved with hardware rounding modes
91 (which could even be some weird unknown style), so something unambiguous
92 and straightforward is best.
93
94
95 The IEEE code below is the usual case, it knows either a 32-bit or 64-bit
96 limb and is done with shifts and masks. The 64-bit case in particular
97 should come out nice and compact.
98
99 The generic code works one bit at a time, which will be quite slow, but
100 should support any binary-based "double" and be safe against any rounding
101 mode. Note in particular it works on IEEE systems too.
102
103
104 Traps:
105
106 Hardware traps for overflow to infinity, underflow to zero, or
107 unsupported denorms may or may not be taken. The IEEE code works bitwise
108 and so probably won't trigger them, the generic code works by float
109 operations and so probably will. This difference might be thought less
110 than ideal, but again its felt straightforward code is better than trying
111 to get intimate with hardware exceptions (of perhaps unknown nature).
112
113
114 Not done:
115
116 mpz_get_d in the past handled size==1 with a cast limb->double. This
117 might still be worthwhile there (for up to the mantissa many bits), but
118 for mpn_get_d here, the cost of applying "exp" to the resulting exponent
119 would probably use up any benefit a cast may have over bit twiddling.
120 Also, if the exponent is pushed into denorm range then bit twiddling is
121 the only option, to ensure the desired truncation is obtained.
122
123
124 Other:
125
126 For reference, note that HPPA 8000, 8200, 8500 and 8600 trap FCNV,UDW,DBL
127 to the kernel for values >= 2^63. This makes it slow, and worse the kernel
128 Linux (what versions?) apparently uses untested code in its trap handling
129 routines, and gets the sign wrong. We don't use such a limb-to-double
130 cast, neither in the IEEE or generic code. */
131
132
133 double
mpn_get_d(mp_srcptr up,mp_size_t size,mp_size_t sign,long exp)134 mpn_get_d (mp_srcptr up, mp_size_t size, mp_size_t sign, long exp)
135 {
136 ASSERT (size >= 0);
137 ASSERT_MPN (up, size);
138 ASSERT (size == 0 || up[size-1] != 0);
139
140 if (size == 0)
141 return 0.0;
142
143 /* Adjust exp to a radix point just above {up,size}, guarding against
144 overflow. After this exp can of course be reduced to anywhere within
145 the {up,size} region without underflow. */
146 if (UNLIKELY ((unsigned long) (GMP_NUMB_BITS * size)
147 > (unsigned long) (LONG_MAX - exp)))
148 {
149 if (_GMP_IEEE_FLOATS)
150 goto ieee_infinity;
151
152 /* generic */
153 exp = LONG_MAX;
154 }
155 else
156 {
157 exp += GMP_NUMB_BITS * size;
158 }
159
160
161 #if 1
162 {
163 int lshift, nbits;
164 union ieee_double_extract u;
165 mp_limb_t x, mhi, mlo;
166 #if GMP_LIMB_BITS == 64
167 mp_limb_t m;
168 up += size;
169 m = *--up;
170 count_leading_zeros (lshift, m);
171
172 exp -= (lshift - GMP_NAIL_BITS) + 1;
173 m <<= lshift;
174
175 nbits = GMP_LIMB_BITS - lshift;
176
177 if (nbits < 53 && size > 1)
178 {
179 x = *--up;
180 x <<= GMP_NAIL_BITS;
181 x >>= nbits;
182 m |= x;
183 nbits += GMP_NUMB_BITS;
184
185 if (LIMBS_PER_DOUBLE >= 3 && nbits < 53 && size > 2)
186 {
187 x = *--up;
188 x <<= GMP_NAIL_BITS;
189 x >>= nbits;
190 m |= x;
191 nbits += GMP_NUMB_BITS;
192 }
193 }
194 mhi = m >> (32 + 11);
195 mlo = m >> 11;
196 #endif
197 #if GMP_LIMB_BITS == 32
198 up += size;
199 x = *--up, size--;
200 count_leading_zeros (lshift, x);
201
202 exp -= (lshift - GMP_NAIL_BITS) + 1;
203 x <<= lshift;
204 mhi = x >> 11;
205
206 if (lshift < 11) /* FIXME: never true if NUMB < 20 bits */
207 {
208 /* All 20 bits in mhi */
209 mlo = x << 21;
210 /* >= 1 bit in mlo */
211 nbits = GMP_LIMB_BITS - lshift - 21;
212 }
213 else
214 {
215 if (size != 0)
216 {
217 nbits = GMP_LIMB_BITS - lshift;
218
219 x = *--up, size--;
220 x <<= GMP_NAIL_BITS;
221 mhi |= x >> nbits >> 11;
222
223 mlo = x << GMP_LIMB_BITS - nbits - 11;
224 nbits = nbits + 11 - GMP_NAIL_BITS;
225 }
226 else
227 {
228 mlo = 0;
229 goto done;
230 }
231 }
232
233 if (LIMBS_PER_DOUBLE >= 2 && nbits < 32 && size != 0)
234 {
235 x = *--up, size--;
236 x <<= GMP_NAIL_BITS;
237 x >>= nbits;
238 mlo |= x;
239 nbits += GMP_NUMB_BITS;
240
241 if (LIMBS_PER_DOUBLE >= 3 && nbits < 32 && size != 0)
242 {
243 x = *--up, size--;
244 x <<= GMP_NAIL_BITS;
245 x >>= nbits;
246 mlo |= x;
247 nbits += GMP_NUMB_BITS;
248
249 if (LIMBS_PER_DOUBLE >= 4 && nbits < 32 && size != 0)
250 {
251 x = *--up;
252 x <<= GMP_NAIL_BITS;
253 x >>= nbits;
254 mlo |= x;
255 nbits += GMP_NUMB_BITS;
256 }
257 }
258 }
259
260 done:;
261
262 #endif
263 {
264 if (UNLIKELY (exp >= CONST_1024))
265 {
266 /* overflow, return infinity */
267 ieee_infinity:
268 mhi = 0;
269 mlo = 0;
270 exp = 1024;
271 }
272 else if (UNLIKELY (exp <= CONST_NEG_1023))
273 {
274 int rshift;
275
276 if (LIKELY (exp <= CONST_NEG_1022_SUB_53))
277 return 0.0; /* denorm underflows to zero */
278
279 rshift = -1022 - exp;
280 ASSERT (rshift > 0 && rshift < 53);
281 #if GMP_LIMB_BITS > 53
282 mlo >>= rshift;
283 mhi = mlo >> 32;
284 #else
285 if (rshift >= 32)
286 {
287 mlo = mhi;
288 mhi = 0;
289 rshift -= 32;
290 }
291 lshift = GMP_LIMB_BITS - rshift;
292 mlo = (mlo >> rshift) | (rshift == 0 ? 0 : mhi << lshift);
293 mhi >>= rshift;
294 #endif
295 exp = -1023;
296 }
297 }
298 u.s.manh = mhi;
299 u.s.manl = mlo;
300 u.s.exp = exp + 1023;
301 u.s.sig = (sign < 0);
302 return u.d;
303 }
304 #else
305
306
307 #define ONE_LIMB (GMP_LIMB_BITS == 64 && 2*GMP_NUMB_BITS >= 53)
308 #define TWO_LIMBS (GMP_LIMB_BITS == 32 && 3*GMP_NUMB_BITS >= 53)
309
310 if (_GMP_IEEE_FLOATS && (ONE_LIMB || TWO_LIMBS))
311 {
312 union ieee_double_extract u;
313 mp_limb_t m0, m1, m2, rmask;
314 int lshift, rshift;
315
316 m0 = up[size-1]; /* high limb */
317 m1 = (size >= 2 ? up[size-2] : 0); /* second highest limb */
318 count_leading_zeros (lshift, m0);
319
320 /* relative to just under high non-zero bit */
321 exp -= (lshift - GMP_NAIL_BITS) + 1;
322
323 if (ONE_LIMB)
324 {
325 /* lshift to have high of m0 non-zero, and collapse nails */
326 rshift = GMP_LIMB_BITS - lshift;
327 m1 <<= GMP_NAIL_BITS;
328 rmask = GMP_NAIL_BITS == 0 && lshift == 0 ? 0 : MP_LIMB_T_MAX;
329 m0 = (m0 << lshift) | ((m1 >> rshift) & rmask);
330
331 /* rshift back to have bit 53 of m0 the high non-zero */
332 m0 >>= 11;
333 }
334 else /* TWO_LIMBS */
335 {
336 m2 = (size >= 3 ? up[size-3] : 0); /* third highest limb */
337
338 /* collapse nails from m1 and m2 */
339 #if GMP_NAIL_BITS != 0
340 m1 = (m1 << GMP_NAIL_BITS) | (m2 >> (GMP_NUMB_BITS-GMP_NAIL_BITS));
341 m2 <<= 2*GMP_NAIL_BITS;
342 #endif
343
344 /* lshift to have high of m0:m1 non-zero, collapse nails from m0 */
345 rshift = GMP_LIMB_BITS - lshift;
346 rmask = (GMP_NAIL_BITS == 0 && lshift == 0 ? 0 : MP_LIMB_T_MAX);
347 m0 = (m0 << lshift) | ((m1 >> rshift) & rmask);
348 m1 = (m1 << lshift) | ((m2 >> rshift) & rmask);
349
350 /* rshift back to have bit 53 of m0:m1 the high non-zero */
351 m1 = (m1 >> 11) | (m0 << (GMP_LIMB_BITS-11));
352 m0 >>= 11;
353 }
354
355 if (UNLIKELY (exp >= CONST_1024))
356 {
357 /* overflow, return infinity */
358 ieee_infinity:
359 m0 = 0;
360 m1 = 0;
361 exp = 1024;
362 }
363 else if (UNLIKELY (exp <= CONST_NEG_1023))
364 {
365 if (LIKELY (exp <= CONST_NEG_1022_SUB_53))
366 return 0.0; /* denorm underflows to zero */
367
368 rshift = -1022 - exp;
369 ASSERT (rshift > 0 && rshift < 53);
370 if (ONE_LIMB)
371 {
372 m0 >>= rshift;
373 }
374 else /* TWO_LIMBS */
375 {
376 if (rshift >= 32)
377 {
378 m1 = m0;
379 m0 = 0;
380 rshift -= 32;
381 }
382 lshift = GMP_LIMB_BITS - rshift;
383 m1 = (m1 >> rshift) | (rshift == 0 ? 0 : m0 << lshift);
384 m0 >>= rshift;
385 }
386 exp = -1023;
387 }
388
389 if (ONE_LIMB)
390 {
391 #if GMP_LIMB_BITS > 32 /* avoid compiler warning about big shift */
392 u.s.manh = m0 >> 32;
393 #endif
394 u.s.manl = m0;
395 }
396 else /* TWO_LIMBS */
397 {
398 u.s.manh = m0;
399 u.s.manl = m1;
400 }
401
402 u.s.exp = exp + 1023;
403 u.s.sig = (sign < 0);
404 return u.d;
405 }
406 else
407 {
408 /* Non-IEEE or strange limb size, do something generic. */
409
410 mp_size_t i;
411 mp_limb_t limb, bit;
412 int shift;
413 double base, factor, prev_factor, d, new_d, diff;
414
415 /* "limb" is "up[i]" the limb being examined, "bit" is a mask for the
416 bit being examined, initially the highest non-zero bit. */
417 i = size-1;
418 limb = up[i];
419 count_leading_zeros (shift, limb);
420 bit = GMP_LIMB_HIGHBIT >> shift;
421
422 /* relative to just under high non-zero bit */
423 exp -= (shift - GMP_NAIL_BITS) + 1;
424
425 /* Power up "factor" to 2^exp, being the value of the "bit" in "limb"
426 being examined. */
427 base = (exp >= 0 ? 2.0 : 0.5);
428 exp = ABS (exp);
429 factor = 1.0;
430 for (;;)
431 {
432 if (exp & 1)
433 {
434 prev_factor = factor;
435 factor *= base;
436 FORCE_DOUBLE (factor);
437 if (factor == 0.0)
438 return 0.0; /* underflow */
439 if (factor == prev_factor)
440 {
441 d = factor; /* overflow, apparent infinity */
442 goto generic_done;
443 }
444 }
445 exp >>= 1;
446 if (exp == 0)
447 break;
448 base *= base;
449 }
450
451 /* Add a "factor" for each non-zero bit, working from high to low.
452 Stop if any rounding occurs, hence implementing a truncation.
453
454 Note no attention is paid to DBL_MANT_DIG, since the effective
455 number of bits in the mantissa isn't constant when in denorm range.
456 We also encountered an ARM system with apparently somewhat doubtful
457 software floats where DBL_MANT_DIG claimed 53 bits but only 32
458 actually worked. */
459
460 d = factor; /* high bit */
461 for (;;)
462 {
463 factor *= 0.5; /* next bit */
464 bit >>= 1;
465 if (bit == 0)
466 {
467 /* next limb, if any */
468 i--;
469 if (i < 0)
470 break;
471 limb = up[i];
472 bit = GMP_NUMB_HIGHBIT;
473 }
474
475 if (bit & limb)
476 {
477 new_d = d + factor;
478 FORCE_DOUBLE (new_d);
479 diff = new_d - d;
480 if (diff != factor)
481 break; /* rounding occured, stop now */
482 d = new_d;
483 }
484 }
485
486 generic_done:
487 return (sign >= 0 ? d : -d);
488 }
489 #endif
490 }
491