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 } 439