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 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 128 Linux kernel (what versions?) apparently uses untested code in its trap 129 handling routines, and gets the sign wrong. We don't use such a limb to 130 double cast, neither in the IEEE or generic code. */ 131 132 133 double 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 = GMP_LIMB_BITS - lshift; 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 == 64) 282 { 283 mlo = (mlo >> rshift) | (mhi << lshift); 284 mhi >>= rshift; 285 } 286 else 287 { 288 if (rshift >= 32) 289 { 290 mlo = mhi; 291 mhi = 0; 292 rshift -= 32; 293 } 294 lshift = GMP_LIMB_BITS - rshift; 295 mlo = (mlo >> rshift) | (rshift == 0 ? 0 : mhi << lshift); 296 mhi >>= rshift; 297 } 298 exp = -1023; 299 } 300 } 301 u.s.manh = mhi; 302 u.s.manl = mlo; 303 u.s.exp = exp + 1023; 304 u.s.sig = (sign < 0); 305 return u.d; 306 } 307 #else 308 309 310 #define ONE_LIMB (GMP_LIMB_BITS == 64 && 2*GMP_NUMB_BITS >= 53) 311 #define TWO_LIMBS (GMP_LIMB_BITS == 32 && 3*GMP_NUMB_BITS >= 53) 312 313 if (_GMP_IEEE_FLOATS && (ONE_LIMB || TWO_LIMBS)) 314 { 315 union ieee_double_extract u; 316 mp_limb_t m0, m1, m2, rmask; 317 int lshift, rshift; 318 319 m0 = up[size-1]; /* high limb */ 320 m1 = (size >= 2 ? up[size-2] : 0); /* second highest limb */ 321 count_leading_zeros (lshift, m0); 322 323 /* relative to just under high non-zero bit */ 324 exp -= (lshift - GMP_NAIL_BITS) + 1; 325 326 if (ONE_LIMB) 327 { 328 /* lshift to have high of m0 non-zero, and collapse nails */ 329 rshift = GMP_LIMB_BITS - lshift; 330 m1 <<= GMP_NAIL_BITS; 331 rmask = GMP_NAIL_BITS == 0 && lshift == 0 ? 0 : MP_LIMB_T_MAX; 332 m0 = (m0 << lshift) | ((m1 >> rshift) & rmask); 333 334 /* rshift back to have bit 53 of m0 the high non-zero */ 335 m0 >>= 11; 336 } 337 else /* TWO_LIMBS */ 338 { 339 m2 = (size >= 3 ? up[size-3] : 0); /* third highest limb */ 340 341 /* collapse nails from m1 and m2 */ 342 #if GMP_NAIL_BITS != 0 343 m1 = (m1 << GMP_NAIL_BITS) | (m2 >> (GMP_NUMB_BITS-GMP_NAIL_BITS)); 344 m2 <<= 2*GMP_NAIL_BITS; 345 #endif 346 347 /* lshift to have high of m0:m1 non-zero, collapse nails from m0 */ 348 rshift = GMP_LIMB_BITS - lshift; 349 rmask = (GMP_NAIL_BITS == 0 && lshift == 0 ? 0 : MP_LIMB_T_MAX); 350 m0 = (m0 << lshift) | ((m1 >> rshift) & rmask); 351 m1 = (m1 << lshift) | ((m2 >> rshift) & rmask); 352 353 /* rshift back to have bit 53 of m0:m1 the high non-zero */ 354 m1 = (m1 >> 11) | (m0 << (GMP_LIMB_BITS-11)); 355 m0 >>= 11; 356 } 357 358 if (UNLIKELY (exp >= CONST_1024)) 359 { 360 /* overflow, return infinity */ 361 ieee_infinity: 362 m0 = 0; 363 m1 = 0; 364 exp = 1024; 365 } 366 else if (UNLIKELY (exp <= CONST_NEG_1023)) 367 { 368 if (LIKELY (exp <= CONST_NEG_1022_SUB_53)) 369 return 0.0; /* denorm underflows to zero */ 370 371 rshift = -1022 - exp; 372 ASSERT (rshift > 0 && rshift < 53); 373 if (ONE_LIMB) 374 { 375 m0 >>= rshift; 376 } 377 else /* TWO_LIMBS */ 378 { 379 if (rshift >= 32) 380 { 381 m1 = m0; 382 m0 = 0; 383 rshift -= 32; 384 } 385 lshift = GMP_LIMB_BITS - rshift; 386 m1 = (m1 >> rshift) | (rshift == 0 ? 0 : m0 << lshift); 387 m0 >>= rshift; 388 } 389 exp = -1023; 390 } 391 392 if (ONE_LIMB) 393 { 394 #if GMP_LIMB_BITS > 32 /* avoid compiler warning about big shift */ 395 u.s.manh = m0 >> 32; 396 #endif 397 u.s.manl = m0; 398 } 399 else /* TWO_LIMBS */ 400 { 401 u.s.manh = m0; 402 u.s.manl = m1; 403 } 404 405 u.s.exp = exp + 1023; 406 u.s.sig = (sign < 0); 407 return u.d; 408 } 409 else 410 { 411 /* Non-IEEE or strange limb size, do something generic. */ 412 413 mp_size_t i; 414 mp_limb_t limb, bit; 415 int shift; 416 double base, factor, prev_factor, d, new_d, diff; 417 418 /* "limb" is "up[i]" the limb being examined, "bit" is a mask for the 419 bit being examined, initially the highest non-zero bit. */ 420 i = size-1; 421 limb = up[i]; 422 count_leading_zeros (shift, limb); 423 bit = GMP_LIMB_HIGHBIT >> shift; 424 425 /* relative to just under high non-zero bit */ 426 exp -= (shift - GMP_NAIL_BITS) + 1; 427 428 /* Power up "factor" to 2^exp, being the value of the "bit" in "limb" 429 being examined. */ 430 base = (exp >= 0 ? 2.0 : 0.5); 431 exp = ABS (exp); 432 factor = 1.0; 433 for (;;) 434 { 435 if (exp & 1) 436 { 437 prev_factor = factor; 438 factor *= base; 439 FORCE_DOUBLE (factor); 440 if (factor == 0.0) 441 return 0.0; /* underflow */ 442 if (factor == prev_factor) 443 { 444 d = factor; /* overflow, apparent infinity */ 445 goto generic_done; 446 } 447 } 448 exp >>= 1; 449 if (exp == 0) 450 break; 451 base *= base; 452 } 453 454 /* Add a "factor" for each non-zero bit, working from high to low. 455 Stop if any rounding occurs, hence implementing a truncation. 456 457 Note no attention is paid to DBL_MANT_DIG, since the effective 458 number of bits in the mantissa isn't constant when in denorm range. 459 We also encountered an ARM system with apparently somewhat doubtful 460 software floats where DBL_MANT_DIG claimed 53 bits but only 32 461 actually worked. */ 462 463 d = factor; /* high bit */ 464 for (;;) 465 { 466 factor *= 0.5; /* next bit */ 467 bit >>= 1; 468 if (bit == 0) 469 { 470 /* next limb, if any */ 471 i--; 472 if (i < 0) 473 break; 474 limb = up[i]; 475 bit = GMP_NUMB_HIGHBIT; 476 } 477 478 if (bit & limb) 479 { 480 new_d = d + factor; 481 FORCE_DOUBLE (new_d); 482 diff = new_d - d; 483 if (diff != factor) 484 break; /* rounding occured, stop now */ 485 d = new_d; 486 } 487 } 488 489 generic_done: 490 return (sign >= 0 ? d : -d); 491 } 492 #endif 493 } 494