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, 2010, 2012, 2018 Free Software
8 Foundation, Inc.
9
10 This file is part of the GNU MP Library.
11
12 The GNU MP Library is free software; you can redistribute it and/or modify
13 it under the terms of either:
14
15 * the GNU Lesser General Public License as published by the Free
16 Software Foundation; either version 3 of the License, or (at your
17 option) any later version.
18
19 or
20
21 * the GNU General Public License as published by the Free Software
22 Foundation; either version 2 of the License, or (at your option) any
23 later version.
24
25 or both in parallel, as here.
26
27 The GNU MP Library is distributed in the hope that it will be useful, but
28 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
29 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
30 for more details.
31
32 You should have received copies of the GNU General Public License and the
33 GNU Lesser General Public License along with the GNU MP Library. If not,
34 see https://www.gnu.org/licenses/. */
35
36 #include "config.h"
37
38 #if HAVE_FLOAT_H
39 #include <float.h> /* for DBL_MANT_DIG and FLT_RADIX */
40 #endif
41
42 #include "gmp-impl.h"
43 #include "longlong.h"
44
45 #ifndef _GMP_IEEE_FLOATS
46 #define _GMP_IEEE_FLOATS 0
47 #endif
48
49 /* To force use of the generic C code for testing, put
50 "#define _GMP_IEEE_FLOATS 0" at this point. */
51
52
53 /* In alpha gcc prior to 3.4, signed DI comparisons involving constants are
54 rearranged from "x < n" to "x+(-n) < 0", which is of course hopelessly
55 wrong if that addition overflows.
56
57 The workaround here avoids this bug by ensuring n is not a literal constant.
58 Note that this is alpha specific. The offending transformation is/was in
59 alpha.c alpha_emit_conditional_branch() under "We want to use cmpcc/bcc".
60
61 Bizarrely, this happens also with Cray cc on alphaev5-cray-unicosmk2.0.6.X,
62 and has the same solution. Don't know why or how. */
63
64 #if HAVE_HOST_CPU_FAMILY_alpha \
65 && ((defined (__GNUC__) && ! __GMP_GNUC_PREREQ(3,4)) \
66 || defined (_CRAY))
67 static volatile const long CONST_1024 = 1024;
68 static volatile const long CONST_NEG_1023 = -1023;
69 static volatile const long CONST_NEG_1022_SUB_53 = -1022 - 53;
70 #else
71 #define CONST_1024 (1024)
72 #define CONST_NEG_1023 (-1023)
73 #define CONST_NEG_1022_SUB_53 (-1022 - 53)
74 #endif
75
76
77 /* Return the value {ptr,size}*2^exp, and negative if sign<0. Must have
78 size>=1, and a non-zero high limb ptr[size-1].
79
80 When we know the fp format, the result is truncated towards zero. This is
81 consistent with other gmp conversions, like mpz_set_f or mpz_set_q, and is
82 easy to implement and test.
83
84 When we do not know the format, such truncation seems much harder. One
85 would need to defeat any rounding mode, including round-up.
86
87 It's felt that GMP is not primarily concerned with hardware floats, and
88 really isn't enhanced by getting involved with hardware rounding modes
89 (which could even be some weird unknown style), so something unambiguous and
90 straightforward is best.
91
92
93 The IEEE code below is the usual case, it knows either a 32-bit or 64-bit
94 limb and is done with shifts and masks. The 64-bit case in particular
95 should come out nice and compact.
96
97 The generic code used to work one bit at a time, which was not only slow,
98 but implicitly relied upon denorms for intermediates, since the lowest bits'
99 weight of a perfectly valid fp number underflows in non-denorm. Therefore,
100 the generic code now works limb-per-limb, initially creating a number x such
101 that 1 <= x <= BASE. (BASE is reached only as result of rounding.) Then
102 x's exponent is scaled with explicit code (not ldexp to avoid libm
103 dependency). It is a tap-dance to avoid underflow or overflow, beware!
104
105
106 Traps:
107
108 Hardware traps for overflow to infinity, underflow to zero, or unsupported
109 denorms may or may not be taken. The IEEE code works bitwise and so
110 probably won't trigger them, the generic code works by float operations and
111 so probably will. This difference might be thought less than ideal, but
112 again its felt straightforward code is better than trying to get intimate
113 with hardware exceptions (of perhaps unknown nature).
114
115
116 Not done:
117
118 mpz_get_d in the past handled size==1 with a cast limb->double. This might
119 still be worthwhile there (for up to the mantissa many bits), but for
120 mpn_get_d here, the cost of applying "exp" to the resulting exponent would
121 probably use up any benefit a cast may have over bit twiddling. Also, if
122 the exponent is pushed into denorm range then bit twiddling is the only
123 option, to ensure the desired truncation is obtained.
124
125
126 Other:
127
128 For reference, note that HPPA 8000, 8200, 8500 and 8600 trap FCNV,UDW,DBL
129 to the kernel for values >= 2^63. This makes it slow, and worse the kernel
130 Linux (what versions?) apparently uses untested code in its trap handling
131 routines, and gets the sign wrong. We don't use such a limb-to-double
132 cast, neither in the IEEE or generic code. */
133
134
135
136 #undef FORMAT_RECOGNIZED
137
138 double
139 mpn_get_d (mp_srcptr up, mp_size_t size, mp_size_t sign, long exp)
140 {
141 int lshift, nbits;
142 mp_limb_t x, mhi, mlo;
143
144 ASSERT (size >= 0);
145 ASSERT_MPN (up, size);
146 ASSERT (size == 0 || up[size-1] != 0);
147
148 if (size == 0)
149 return 0.0;
150
151 /* Adjust exp to a radix point just above {up,size}, guarding against
152 overflow. After this exp can of course be reduced to anywhere within
153 the {up,size} region without underflow. */
154 if (UNLIKELY ((unsigned long) (GMP_NUMB_BITS * size)
155 > ((unsigned long) LONG_MAX - exp)))
156 {
157 #if _GMP_IEEE_FLOATS
158 goto ieee_infinity;
159 #endif
160
161 /* generic */
162 exp = LONG_MAX;
163 }
164 else
165 {
166 exp += GMP_NUMB_BITS * size;
167 }
168
169 #if _GMP_IEEE_FLOATS
170 {
171 union ieee_double_extract u;
172
173 up += size;
174
175 #if GMP_LIMB_BITS == 64
176 mlo = up[-1];
177 count_leading_zeros (lshift, mlo);
178
179 exp -= (lshift - GMP_NAIL_BITS) + 1;
180 mlo <<= lshift;
181
182 nbits = GMP_LIMB_BITS - lshift;
183
184 if (nbits < 53 && size > 1)
185 {
186 x = up[-2];
187 x <<= GMP_NAIL_BITS;
188 x >>= nbits;
189 mlo |= x;
190 nbits += GMP_NUMB_BITS;
191
192 if (LIMBS_PER_DOUBLE >= 3 && nbits < 53 && size > 2)
193 {
194 x = up[-3];
195 x <<= GMP_NAIL_BITS;
196 x >>= nbits;
197 mlo |= x;
198 nbits += GMP_NUMB_BITS;
199 }
200 }
201 mhi = mlo >> (32 + 11);
202 mlo = mlo >> 11; /* later implicitly truncated to 32 bits */
203 #endif
204 #if GMP_LIMB_BITS == 32
205 x = *--up;
206 count_leading_zeros (lshift, x);
207
208 exp -= (lshift - GMP_NAIL_BITS) + 1;
209 x <<= lshift;
210 mhi = x >> 11;
211
212 if (lshift < 11) /* FIXME: never true if NUMB < 20 bits */
213 {
214 /* All 20 bits in mhi */
215 mlo = x << 21;
216 /* >= 1 bit in mlo */
217 nbits = GMP_LIMB_BITS - lshift - 21;
218 }
219 else
220 {
221 if (size > 1)
222 {
223 nbits = GMP_LIMB_BITS - lshift;
224
225 x = *--up, size--;
226 x <<= GMP_NAIL_BITS;
227 mhi |= x >> nbits >> 11;
228
229 mlo = x << (GMP_LIMB_BITS - nbits - 11);
230 nbits = nbits + 11 - GMP_NAIL_BITS;
231 }
232 else
233 {
234 mlo = 0;
235 goto done;
236 }
237 }
238
239 /* Now all needed bits in mhi have been accumulated. Add bits to mlo. */
240
241 if (LIMBS_PER_DOUBLE >= 2 && nbits < 32 && size > 1)
242 {
243 x = up[-1];
244 x <<= GMP_NAIL_BITS;
245 x >>= nbits;
246 mlo |= x;
247 nbits += GMP_NUMB_BITS;
248
249 if (LIMBS_PER_DOUBLE >= 3 && nbits < 32 && size > 2)
250 {
251 x = up[-2];
252 x <<= GMP_NAIL_BITS;
253 x >>= nbits;
254 mlo |= x;
255 nbits += GMP_NUMB_BITS;
256
257 if (LIMBS_PER_DOUBLE >= 4 && nbits < 32 && size > 3)
258 {
259 x = up[-3];
260 x <<= GMP_NAIL_BITS;
261 x >>= nbits;
262 mlo |= x;
263 nbits += GMP_NUMB_BITS;
264 }
265 }
266 }
267
268 done:;
269
270 #endif
271 if (UNLIKELY (exp >= CONST_1024))
272 {
273 /* overflow, return infinity */
274 ieee_infinity:
275 mhi = 0;
276 mlo = 0;
277 exp = 1024;
278 }
279 else if (UNLIKELY (exp <= CONST_NEG_1023))
280 {
281 int rshift;
282
283 if (LIKELY (exp <= CONST_NEG_1022_SUB_53))
284 return 0.0; /* denorm underflows to zero */
285
286 rshift = -1022 - exp;
287 ASSERT (rshift > 0 && rshift < 53);
288 #if GMP_LIMB_BITS > 53
289 mlo >>= rshift;
290 mhi = mlo >> 32;
291 #else
292 if (rshift >= 32)
293 {
294 mlo = mhi;
295 mhi = 0;
296 rshift -= 32;
297 }
298 lshift = GMP_LIMB_BITS - rshift;
299 mlo = (mlo >> rshift) | (rshift == 0 ? 0 : mhi << lshift);
300 mhi >>= rshift;
301 #endif
302 exp = -1023;
303 }
304 u.s.manh = mhi;
305 u.s.manl = mlo;
306 u.s.exp = exp + 1023;
307 u.s.sig = (sign < 0);
308 return u.d;
309 }
310 #define FORMAT_RECOGNIZED 1
311 #endif
312
313 #if HAVE_DOUBLE_VAX_D
314 {
315 union double_extract u;
316
317 up += size;
318
319 mhi = up[-1];
320
321 count_leading_zeros (lshift, mhi);
322 exp -= lshift;
323 mhi <<= lshift;
324
325 mlo = 0;
326 if (size > 1)
327 {
328 mlo = up[-2];
329 if (lshift != 0)
330 mhi += mlo >> (GMP_LIMB_BITS - lshift);
331 mlo <<= lshift;
332
333 if (size > 2 && lshift > 8)
334 {
335 x = up[-3];
336 mlo += x >> (GMP_LIMB_BITS - lshift);
337 }
338 }
339
340 if (UNLIKELY (exp >= 128))
341 {
342 /* overflow, return maximum number */
343 mhi = 0xffffffff;
344 mlo = 0xffffffff;
345 exp = 127;
346 }
347 else if (UNLIKELY (exp < -128))
348 {
349 return 0.0; /* underflows to zero */
350 }
351
352 u.s.man3 = mhi >> 24; /* drop msb, since implicit */
353 u.s.man2 = mhi >> 8;
354 u.s.man1 = (mhi << 8) + (mlo >> 24);
355 u.s.man0 = mlo >> 8;
356 u.s.exp = exp + 128;
357 u.s.sig = sign < 0;
358 return u.d;
359 }
360 #define FORMAT_RECOGNIZED 1
361 #endif
362
363 #if ! FORMAT_RECOGNIZED
364
365 #if !defined(GMP_DBL_MANT_BITS)
366 #if defined(DBL_MANT_DIG) && FLT_RADIX == 2
367 #define GMP_DBL_MANT_BITS DBL_MANT_DIG
368 #else
369 /* FIXME: Chose a smarter default value. */
370 #define GMP_DBL_MANT_BITS (16 * sizeof (double))
371 #endif
372 #endif
373
374 { /* Non-IEEE or strange limb size, generically convert
375 GMP_DBL_MANT_BITS bits. */
376 mp_limb_t l;
377 int m;
378 mp_size_t i;
379 double d, weight;
380 unsigned long uexp;
381
382 /* First generate an fp number disregarding exp, instead keeping things
383 within the numb base factor from 1, which should prevent overflow and
384 underflow even for the most exponent limited fp formats. */
385 i = size - 1;
386 l = up[i];
387 count_leading_zeros (m, l);
388 m = m + GMP_DBL_MANT_BITS - GMP_LIMB_BITS;
389 if (m < 0)
390 l &= GMP_NUMB_MAX << -m;
391 d = l;
392 for (weight = 1/MP_BASE_AS_DOUBLE; m > 0 && --i >= 0;)
393 {
394 l = up[i];
395 m -= GMP_NUMB_BITS;
396 if (m < 0)
397 l &= GMP_NUMB_MAX << -m;
398 d += l * weight;
399 weight /= MP_BASE_AS_DOUBLE;
400 if (weight == 0)
401 break;
402 }
403
404 /* Now apply exp. */
405 exp -= GMP_NUMB_BITS;
406 if (exp > 0)
407 {
408 weight = 2.0;
409 uexp = exp;
410 }
411 else
412 {
413 weight = 0.5;
414 uexp = NEG_CAST (unsigned long, exp);
415 }
416 #if 1
417 /* Square-and-multiply exponentiation. */
418 if (uexp & 1)
419 d *= weight;
420 while (uexp >>= 1)
421 {
422 weight *= weight;
423 if (uexp & 1)
424 d *= weight;
425 }
426 #else
427 /* Plain exponentiation. */
428 while (uexp > 0)
429 {
430 d *= weight;
431 uexp--;
432 }
433 #endif
434
435 return sign >= 0 ? d : -d;
436 }
437 #endif
438 }