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