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