1 //===-- APFloat.cpp - Implement APFloat class -----------------------------===// 2 // 3 // The LLVM Compiler Infrastructure 4 // 5 // This file is distributed under the University of Illinois Open Source 6 // License. See LICENSE.TXT for details. 7 // 8 //===----------------------------------------------------------------------===// 9 // 10 // This file implements a class to represent arbitrary precision floating 11 // point values and provide a variety of arithmetic operations on them. 12 // 13 //===----------------------------------------------------------------------===// 14 15 #include "llvm/ADT/APFloat.h" 16 #include "llvm/ADT/APSInt.h" 17 #include "llvm/ADT/ArrayRef.h" 18 #include "llvm/ADT/FoldingSet.h" 19 #include "llvm/ADT/Hashing.h" 20 #include "llvm/ADT/StringExtras.h" 21 #include "llvm/ADT/StringRef.h" 22 #include "llvm/Support/Debug.h" 23 #include "llvm/Support/ErrorHandling.h" 24 #include "llvm/Support/MathExtras.h" 25 #include "llvm/Support/raw_ostream.h" 26 #include <cstring> 27 #include <limits.h> 28 29 using namespace llvm; 30 31 /// A macro used to combine two fcCategory enums into one key which can be used 32 /// in a switch statement to classify how the interaction of two APFloat's 33 /// categories affects an operation. 34 /// 35 /// TODO: If clang source code is ever allowed to use constexpr in its own 36 /// codebase, change this into a static inline function. 37 #define PackCategoriesIntoKey(_lhs, _rhs) ((_lhs) * 4 + (_rhs)) 38 39 /* Assumed in hexadecimal significand parsing, and conversion to 40 hexadecimal strings. */ 41 static_assert(integerPartWidth % 4 == 0, "Part width must be divisible by 4!"); 42 43 namespace llvm { 44 /* Represents floating point arithmetic semantics. */ 45 struct fltSemantics { 46 /* The largest E such that 2^E is representable; this matches the 47 definition of IEEE 754. */ 48 APFloatBase::ExponentType maxExponent; 49 50 /* The smallest E such that 2^E is a normalized number; this 51 matches the definition of IEEE 754. */ 52 APFloatBase::ExponentType minExponent; 53 54 /* Number of bits in the significand. This includes the integer 55 bit. */ 56 unsigned int precision; 57 58 /* Number of bits actually used in the semantics. */ 59 unsigned int sizeInBits; 60 }; 61 62 static const fltSemantics semIEEEhalf = {15, -14, 11, 16}; 63 static const fltSemantics semIEEEsingle = {127, -126, 24, 32}; 64 static const fltSemantics semIEEEdouble = {1023, -1022, 53, 64}; 65 static const fltSemantics semIEEEquad = {16383, -16382, 113, 128}; 66 static const fltSemantics semX87DoubleExtended = {16383, -16382, 64, 80}; 67 static const fltSemantics semBogus = {0, 0, 0, 0}; 68 69 /* The PowerPC format consists of two doubles. It does not map cleanly 70 onto the usual format above. It is approximated using twice the 71 mantissa bits. Note that for exponents near the double minimum, 72 we no longer can represent the full 106 mantissa bits, so those 73 will be treated as denormal numbers. 74 75 FIXME: While this approximation is equivalent to what GCC uses for 76 compile-time arithmetic on PPC double-double numbers, it is not able 77 to represent all possible values held by a PPC double-double number, 78 for example: (long double) 1.0 + (long double) 0x1p-106 79 Should this be replaced by a full emulation of PPC double-double? */ 80 static const fltSemantics semPPCDoubleDouble = {0, 0, 0, 0}; 81 82 /* There are temporary semantics for the real PPCDoubleDouble implementation. 83 Currently, APFloat of PPCDoubleDouble holds one PPCDoubleDoubleImpl as the 84 high part of double double, and one IEEEdouble as the low part, so that 85 the old operations operate on PPCDoubleDoubleImpl, while the newly added 86 operations also populate the IEEEdouble. 87 88 TODO: Once all functions support DoubleAPFloat mode, we'll change all 89 PPCDoubleDoubleImpl to IEEEdouble and remove PPCDoubleDoubleImpl. */ 90 static const fltSemantics semPPCDoubleDoubleImpl = {1023, -1022 + 53, 53 + 53, 91 128}; 92 93 const fltSemantics &APFloatBase::IEEEhalf() { 94 return semIEEEhalf; 95 } 96 const fltSemantics &APFloatBase::IEEEsingle() { 97 return semIEEEsingle; 98 } 99 const fltSemantics &APFloatBase::IEEEdouble() { 100 return semIEEEdouble; 101 } 102 const fltSemantics &APFloatBase::IEEEquad() { 103 return semIEEEquad; 104 } 105 const fltSemantics &APFloatBase::x87DoubleExtended() { 106 return semX87DoubleExtended; 107 } 108 const fltSemantics &APFloatBase::Bogus() { 109 return semBogus; 110 } 111 const fltSemantics &APFloatBase::PPCDoubleDouble() { 112 return semPPCDoubleDouble; 113 } 114 115 /* A tight upper bound on number of parts required to hold the value 116 pow(5, power) is 117 118 power * 815 / (351 * integerPartWidth) + 1 119 120 However, whilst the result may require only this many parts, 121 because we are multiplying two values to get it, the 122 multiplication may require an extra part with the excess part 123 being zero (consider the trivial case of 1 * 1, tcFullMultiply 124 requires two parts to hold the single-part result). So we add an 125 extra one to guarantee enough space whilst multiplying. */ 126 const unsigned int maxExponent = 16383; 127 const unsigned int maxPrecision = 113; 128 const unsigned int maxPowerOfFiveExponent = maxExponent + maxPrecision - 1; 129 const unsigned int maxPowerOfFiveParts = 2 + ((maxPowerOfFiveExponent * 815) 130 / (351 * integerPartWidth)); 131 132 unsigned int APFloatBase::semanticsPrecision(const fltSemantics &semantics) { 133 return semantics.precision; 134 } 135 APFloatBase::ExponentType 136 APFloatBase::semanticsMaxExponent(const fltSemantics &semantics) { 137 return semantics.maxExponent; 138 } 139 APFloatBase::ExponentType 140 APFloatBase::semanticsMinExponent(const fltSemantics &semantics) { 141 return semantics.minExponent; 142 } 143 unsigned int APFloatBase::semanticsSizeInBits(const fltSemantics &semantics) { 144 return semantics.sizeInBits; 145 } 146 147 unsigned APFloatBase::getSizeInBits(const fltSemantics &Sem) { 148 return Sem.sizeInBits; 149 } 150 151 /* A bunch of private, handy routines. */ 152 153 static inline unsigned int 154 partCountForBits(unsigned int bits) 155 { 156 return ((bits) + integerPartWidth - 1) / integerPartWidth; 157 } 158 159 /* Returns 0U-9U. Return values >= 10U are not digits. */ 160 static inline unsigned int 161 decDigitValue(unsigned int c) 162 { 163 return c - '0'; 164 } 165 166 /* Return the value of a decimal exponent of the form 167 [+-]ddddddd. 168 169 If the exponent overflows, returns a large exponent with the 170 appropriate sign. */ 171 static int 172 readExponent(StringRef::iterator begin, StringRef::iterator end) 173 { 174 bool isNegative; 175 unsigned int absExponent; 176 const unsigned int overlargeExponent = 24000; /* FIXME. */ 177 StringRef::iterator p = begin; 178 179 assert(p != end && "Exponent has no digits"); 180 181 isNegative = (*p == '-'); 182 if (*p == '-' || *p == '+') { 183 p++; 184 assert(p != end && "Exponent has no digits"); 185 } 186 187 absExponent = decDigitValue(*p++); 188 assert(absExponent < 10U && "Invalid character in exponent"); 189 190 for (; p != end; ++p) { 191 unsigned int value; 192 193 value = decDigitValue(*p); 194 assert(value < 10U && "Invalid character in exponent"); 195 196 value += absExponent * 10; 197 if (absExponent >= overlargeExponent) { 198 absExponent = overlargeExponent; 199 p = end; /* outwit assert below */ 200 break; 201 } 202 absExponent = value; 203 } 204 205 assert(p == end && "Invalid exponent in exponent"); 206 207 if (isNegative) 208 return -(int) absExponent; 209 else 210 return (int) absExponent; 211 } 212 213 /* This is ugly and needs cleaning up, but I don't immediately see 214 how whilst remaining safe. */ 215 static int 216 totalExponent(StringRef::iterator p, StringRef::iterator end, 217 int exponentAdjustment) 218 { 219 int unsignedExponent; 220 bool negative, overflow; 221 int exponent = 0; 222 223 assert(p != end && "Exponent has no digits"); 224 225 negative = *p == '-'; 226 if (*p == '-' || *p == '+') { 227 p++; 228 assert(p != end && "Exponent has no digits"); 229 } 230 231 unsignedExponent = 0; 232 overflow = false; 233 for (; p != end; ++p) { 234 unsigned int value; 235 236 value = decDigitValue(*p); 237 assert(value < 10U && "Invalid character in exponent"); 238 239 unsignedExponent = unsignedExponent * 10 + value; 240 if (unsignedExponent > 32767) { 241 overflow = true; 242 break; 243 } 244 } 245 246 if (exponentAdjustment > 32767 || exponentAdjustment < -32768) 247 overflow = true; 248 249 if (!overflow) { 250 exponent = unsignedExponent; 251 if (negative) 252 exponent = -exponent; 253 exponent += exponentAdjustment; 254 if (exponent > 32767 || exponent < -32768) 255 overflow = true; 256 } 257 258 if (overflow) 259 exponent = negative ? -32768: 32767; 260 261 return exponent; 262 } 263 264 static StringRef::iterator 265 skipLeadingZeroesAndAnyDot(StringRef::iterator begin, StringRef::iterator end, 266 StringRef::iterator *dot) 267 { 268 StringRef::iterator p = begin; 269 *dot = end; 270 while (p != end && *p == '0') 271 p++; 272 273 if (p != end && *p == '.') { 274 *dot = p++; 275 276 assert(end - begin != 1 && "Significand has no digits"); 277 278 while (p != end && *p == '0') 279 p++; 280 } 281 282 return p; 283 } 284 285 /* Given a normal decimal floating point number of the form 286 287 dddd.dddd[eE][+-]ddd 288 289 where the decimal point and exponent are optional, fill out the 290 structure D. Exponent is appropriate if the significand is 291 treated as an integer, and normalizedExponent if the significand 292 is taken to have the decimal point after a single leading 293 non-zero digit. 294 295 If the value is zero, V->firstSigDigit points to a non-digit, and 296 the return exponent is zero. 297 */ 298 struct decimalInfo { 299 const char *firstSigDigit; 300 const char *lastSigDigit; 301 int exponent; 302 int normalizedExponent; 303 }; 304 305 static void 306 interpretDecimal(StringRef::iterator begin, StringRef::iterator end, 307 decimalInfo *D) 308 { 309 StringRef::iterator dot = end; 310 StringRef::iterator p = skipLeadingZeroesAndAnyDot (begin, end, &dot); 311 312 D->firstSigDigit = p; 313 D->exponent = 0; 314 D->normalizedExponent = 0; 315 316 for (; p != end; ++p) { 317 if (*p == '.') { 318 assert(dot == end && "String contains multiple dots"); 319 dot = p++; 320 if (p == end) 321 break; 322 } 323 if (decDigitValue(*p) >= 10U) 324 break; 325 } 326 327 if (p != end) { 328 assert((*p == 'e' || *p == 'E') && "Invalid character in significand"); 329 assert(p != begin && "Significand has no digits"); 330 assert((dot == end || p - begin != 1) && "Significand has no digits"); 331 332 /* p points to the first non-digit in the string */ 333 D->exponent = readExponent(p + 1, end); 334 335 /* Implied decimal point? */ 336 if (dot == end) 337 dot = p; 338 } 339 340 /* If number is all zeroes accept any exponent. */ 341 if (p != D->firstSigDigit) { 342 /* Drop insignificant trailing zeroes. */ 343 if (p != begin) { 344 do 345 do 346 p--; 347 while (p != begin && *p == '0'); 348 while (p != begin && *p == '.'); 349 } 350 351 /* Adjust the exponents for any decimal point. */ 352 D->exponent += static_cast<APFloat::ExponentType>((dot - p) - (dot > p)); 353 D->normalizedExponent = (D->exponent + 354 static_cast<APFloat::ExponentType>((p - D->firstSigDigit) 355 - (dot > D->firstSigDigit && dot < p))); 356 } 357 358 D->lastSigDigit = p; 359 } 360 361 /* Return the trailing fraction of a hexadecimal number. 362 DIGITVALUE is the first hex digit of the fraction, P points to 363 the next digit. */ 364 static lostFraction 365 trailingHexadecimalFraction(StringRef::iterator p, StringRef::iterator end, 366 unsigned int digitValue) 367 { 368 unsigned int hexDigit; 369 370 /* If the first trailing digit isn't 0 or 8 we can work out the 371 fraction immediately. */ 372 if (digitValue > 8) 373 return lfMoreThanHalf; 374 else if (digitValue < 8 && digitValue > 0) 375 return lfLessThanHalf; 376 377 // Otherwise we need to find the first non-zero digit. 378 while (p != end && (*p == '0' || *p == '.')) 379 p++; 380 381 assert(p != end && "Invalid trailing hexadecimal fraction!"); 382 383 hexDigit = hexDigitValue(*p); 384 385 /* If we ran off the end it is exactly zero or one-half, otherwise 386 a little more. */ 387 if (hexDigit == -1U) 388 return digitValue == 0 ? lfExactlyZero: lfExactlyHalf; 389 else 390 return digitValue == 0 ? lfLessThanHalf: lfMoreThanHalf; 391 } 392 393 /* Return the fraction lost were a bignum truncated losing the least 394 significant BITS bits. */ 395 static lostFraction 396 lostFractionThroughTruncation(const integerPart *parts, 397 unsigned int partCount, 398 unsigned int bits) 399 { 400 unsigned int lsb; 401 402 lsb = APInt::tcLSB(parts, partCount); 403 404 /* Note this is guaranteed true if bits == 0, or LSB == -1U. */ 405 if (bits <= lsb) 406 return lfExactlyZero; 407 if (bits == lsb + 1) 408 return lfExactlyHalf; 409 if (bits <= partCount * integerPartWidth && 410 APInt::tcExtractBit(parts, bits - 1)) 411 return lfMoreThanHalf; 412 413 return lfLessThanHalf; 414 } 415 416 /* Shift DST right BITS bits noting lost fraction. */ 417 static lostFraction 418 shiftRight(integerPart *dst, unsigned int parts, unsigned int bits) 419 { 420 lostFraction lost_fraction; 421 422 lost_fraction = lostFractionThroughTruncation(dst, parts, bits); 423 424 APInt::tcShiftRight(dst, parts, bits); 425 426 return lost_fraction; 427 } 428 429 /* Combine the effect of two lost fractions. */ 430 static lostFraction 431 combineLostFractions(lostFraction moreSignificant, 432 lostFraction lessSignificant) 433 { 434 if (lessSignificant != lfExactlyZero) { 435 if (moreSignificant == lfExactlyZero) 436 moreSignificant = lfLessThanHalf; 437 else if (moreSignificant == lfExactlyHalf) 438 moreSignificant = lfMoreThanHalf; 439 } 440 441 return moreSignificant; 442 } 443 444 /* The error from the true value, in half-ulps, on multiplying two 445 floating point numbers, which differ from the value they 446 approximate by at most HUE1 and HUE2 half-ulps, is strictly less 447 than the returned value. 448 449 See "How to Read Floating Point Numbers Accurately" by William D 450 Clinger. */ 451 static unsigned int 452 HUerrBound(bool inexactMultiply, unsigned int HUerr1, unsigned int HUerr2) 453 { 454 assert(HUerr1 < 2 || HUerr2 < 2 || (HUerr1 + HUerr2 < 8)); 455 456 if (HUerr1 + HUerr2 == 0) 457 return inexactMultiply * 2; /* <= inexactMultiply half-ulps. */ 458 else 459 return inexactMultiply + 2 * (HUerr1 + HUerr2); 460 } 461 462 /* The number of ulps from the boundary (zero, or half if ISNEAREST) 463 when the least significant BITS are truncated. BITS cannot be 464 zero. */ 465 static integerPart 466 ulpsFromBoundary(const integerPart *parts, unsigned int bits, bool isNearest) 467 { 468 unsigned int count, partBits; 469 integerPart part, boundary; 470 471 assert(bits != 0); 472 473 bits--; 474 count = bits / integerPartWidth; 475 partBits = bits % integerPartWidth + 1; 476 477 part = parts[count] & (~(integerPart) 0 >> (integerPartWidth - partBits)); 478 479 if (isNearest) 480 boundary = (integerPart) 1 << (partBits - 1); 481 else 482 boundary = 0; 483 484 if (count == 0) { 485 if (part - boundary <= boundary - part) 486 return part - boundary; 487 else 488 return boundary - part; 489 } 490 491 if (part == boundary) { 492 while (--count) 493 if (parts[count]) 494 return ~(integerPart) 0; /* A lot. */ 495 496 return parts[0]; 497 } else if (part == boundary - 1) { 498 while (--count) 499 if (~parts[count]) 500 return ~(integerPart) 0; /* A lot. */ 501 502 return -parts[0]; 503 } 504 505 return ~(integerPart) 0; /* A lot. */ 506 } 507 508 /* Place pow(5, power) in DST, and return the number of parts used. 509 DST must be at least one part larger than size of the answer. */ 510 static unsigned int 511 powerOf5(integerPart *dst, unsigned int power) 512 { 513 static const integerPart firstEightPowers[] = { 1, 5, 25, 125, 625, 3125, 514 15625, 78125 }; 515 integerPart pow5s[maxPowerOfFiveParts * 2 + 5]; 516 pow5s[0] = 78125 * 5; 517 518 unsigned int partsCount[16] = { 1 }; 519 integerPart scratch[maxPowerOfFiveParts], *p1, *p2, *pow5; 520 unsigned int result; 521 assert(power <= maxExponent); 522 523 p1 = dst; 524 p2 = scratch; 525 526 *p1 = firstEightPowers[power & 7]; 527 power >>= 3; 528 529 result = 1; 530 pow5 = pow5s; 531 532 for (unsigned int n = 0; power; power >>= 1, n++) { 533 unsigned int pc; 534 535 pc = partsCount[n]; 536 537 /* Calculate pow(5,pow(2,n+3)) if we haven't yet. */ 538 if (pc == 0) { 539 pc = partsCount[n - 1]; 540 APInt::tcFullMultiply(pow5, pow5 - pc, pow5 - pc, pc, pc); 541 pc *= 2; 542 if (pow5[pc - 1] == 0) 543 pc--; 544 partsCount[n] = pc; 545 } 546 547 if (power & 1) { 548 integerPart *tmp; 549 550 APInt::tcFullMultiply(p2, p1, pow5, result, pc); 551 result += pc; 552 if (p2[result - 1] == 0) 553 result--; 554 555 /* Now result is in p1 with partsCount parts and p2 is scratch 556 space. */ 557 tmp = p1; 558 p1 = p2; 559 p2 = tmp; 560 } 561 562 pow5 += pc; 563 } 564 565 if (p1 != dst) 566 APInt::tcAssign(dst, p1, result); 567 568 return result; 569 } 570 571 /* Zero at the end to avoid modular arithmetic when adding one; used 572 when rounding up during hexadecimal output. */ 573 static const char hexDigitsLower[] = "0123456789abcdef0"; 574 static const char hexDigitsUpper[] = "0123456789ABCDEF0"; 575 static const char infinityL[] = "infinity"; 576 static const char infinityU[] = "INFINITY"; 577 static const char NaNL[] = "nan"; 578 static const char NaNU[] = "NAN"; 579 580 /* Write out an integerPart in hexadecimal, starting with the most 581 significant nibble. Write out exactly COUNT hexdigits, return 582 COUNT. */ 583 static unsigned int 584 partAsHex (char *dst, integerPart part, unsigned int count, 585 const char *hexDigitChars) 586 { 587 unsigned int result = count; 588 589 assert(count != 0 && count <= integerPartWidth / 4); 590 591 part >>= (integerPartWidth - 4 * count); 592 while (count--) { 593 dst[count] = hexDigitChars[part & 0xf]; 594 part >>= 4; 595 } 596 597 return result; 598 } 599 600 /* Write out an unsigned decimal integer. */ 601 static char * 602 writeUnsignedDecimal (char *dst, unsigned int n) 603 { 604 char buff[40], *p; 605 606 p = buff; 607 do 608 *p++ = '0' + n % 10; 609 while (n /= 10); 610 611 do 612 *dst++ = *--p; 613 while (p != buff); 614 615 return dst; 616 } 617 618 /* Write out a signed decimal integer. */ 619 static char * 620 writeSignedDecimal (char *dst, int value) 621 { 622 if (value < 0) { 623 *dst++ = '-'; 624 dst = writeUnsignedDecimal(dst, -(unsigned) value); 625 } else 626 dst = writeUnsignedDecimal(dst, value); 627 628 return dst; 629 } 630 631 namespace detail { 632 /* Constructors. */ 633 void IEEEFloat::initialize(const fltSemantics *ourSemantics) { 634 unsigned int count; 635 636 semantics = ourSemantics; 637 count = partCount(); 638 if (count > 1) 639 significand.parts = new integerPart[count]; 640 } 641 642 void IEEEFloat::freeSignificand() { 643 if (needsCleanup()) 644 delete [] significand.parts; 645 } 646 647 void IEEEFloat::assign(const IEEEFloat &rhs) { 648 assert(semantics == rhs.semantics); 649 650 sign = rhs.sign; 651 category = rhs.category; 652 exponent = rhs.exponent; 653 if (isFiniteNonZero() || category == fcNaN) 654 copySignificand(rhs); 655 } 656 657 void IEEEFloat::copySignificand(const IEEEFloat &rhs) { 658 assert(isFiniteNonZero() || category == fcNaN); 659 assert(rhs.partCount() >= partCount()); 660 661 APInt::tcAssign(significandParts(), rhs.significandParts(), 662 partCount()); 663 } 664 665 /* Make this number a NaN, with an arbitrary but deterministic value 666 for the significand. If double or longer, this is a signalling NaN, 667 which may not be ideal. If float, this is QNaN(0). */ 668 void IEEEFloat::makeNaN(bool SNaN, bool Negative, const APInt *fill) { 669 category = fcNaN; 670 sign = Negative; 671 672 integerPart *significand = significandParts(); 673 unsigned numParts = partCount(); 674 675 // Set the significand bits to the fill. 676 if (!fill || fill->getNumWords() < numParts) 677 APInt::tcSet(significand, 0, numParts); 678 if (fill) { 679 APInt::tcAssign(significand, fill->getRawData(), 680 std::min(fill->getNumWords(), numParts)); 681 682 // Zero out the excess bits of the significand. 683 unsigned bitsToPreserve = semantics->precision - 1; 684 unsigned part = bitsToPreserve / 64; 685 bitsToPreserve %= 64; 686 significand[part] &= ((1ULL << bitsToPreserve) - 1); 687 for (part++; part != numParts; ++part) 688 significand[part] = 0; 689 } 690 691 unsigned QNaNBit = semantics->precision - 2; 692 693 if (SNaN) { 694 // We always have to clear the QNaN bit to make it an SNaN. 695 APInt::tcClearBit(significand, QNaNBit); 696 697 // If there are no bits set in the payload, we have to set 698 // *something* to make it a NaN instead of an infinity; 699 // conventionally, this is the next bit down from the QNaN bit. 700 if (APInt::tcIsZero(significand, numParts)) 701 APInt::tcSetBit(significand, QNaNBit - 1); 702 } else { 703 // We always have to set the QNaN bit to make it a QNaN. 704 APInt::tcSetBit(significand, QNaNBit); 705 } 706 707 // For x87 extended precision, we want to make a NaN, not a 708 // pseudo-NaN. Maybe we should expose the ability to make 709 // pseudo-NaNs? 710 if (semantics == &semX87DoubleExtended) 711 APInt::tcSetBit(significand, QNaNBit + 1); 712 } 713 714 IEEEFloat &IEEEFloat::operator=(const IEEEFloat &rhs) { 715 if (this != &rhs) { 716 if (semantics != rhs.semantics) { 717 freeSignificand(); 718 initialize(rhs.semantics); 719 } 720 assign(rhs); 721 } 722 723 return *this; 724 } 725 726 IEEEFloat &IEEEFloat::operator=(IEEEFloat &&rhs) { 727 freeSignificand(); 728 729 semantics = rhs.semantics; 730 significand = rhs.significand; 731 exponent = rhs.exponent; 732 category = rhs.category; 733 sign = rhs.sign; 734 735 rhs.semantics = &semBogus; 736 return *this; 737 } 738 739 bool IEEEFloat::isDenormal() const { 740 return isFiniteNonZero() && (exponent == semantics->minExponent) && 741 (APInt::tcExtractBit(significandParts(), 742 semantics->precision - 1) == 0); 743 } 744 745 bool IEEEFloat::isSmallest() const { 746 // The smallest number by magnitude in our format will be the smallest 747 // denormal, i.e. the floating point number with exponent being minimum 748 // exponent and significand bitwise equal to 1 (i.e. with MSB equal to 0). 749 return isFiniteNonZero() && exponent == semantics->minExponent && 750 significandMSB() == 0; 751 } 752 753 bool IEEEFloat::isSignificandAllOnes() const { 754 // Test if the significand excluding the integral bit is all ones. This allows 755 // us to test for binade boundaries. 756 const integerPart *Parts = significandParts(); 757 const unsigned PartCount = partCount(); 758 for (unsigned i = 0; i < PartCount - 1; i++) 759 if (~Parts[i]) 760 return false; 761 762 // Set the unused high bits to all ones when we compare. 763 const unsigned NumHighBits = 764 PartCount*integerPartWidth - semantics->precision + 1; 765 assert(NumHighBits <= integerPartWidth && "Can not have more high bits to " 766 "fill than integerPartWidth"); 767 const integerPart HighBitFill = 768 ~integerPart(0) << (integerPartWidth - NumHighBits); 769 if (~(Parts[PartCount - 1] | HighBitFill)) 770 return false; 771 772 return true; 773 } 774 775 bool IEEEFloat::isSignificandAllZeros() const { 776 // Test if the significand excluding the integral bit is all zeros. This 777 // allows us to test for binade boundaries. 778 const integerPart *Parts = significandParts(); 779 const unsigned PartCount = partCount(); 780 781 for (unsigned i = 0; i < PartCount - 1; i++) 782 if (Parts[i]) 783 return false; 784 785 const unsigned NumHighBits = 786 PartCount*integerPartWidth - semantics->precision + 1; 787 assert(NumHighBits <= integerPartWidth && "Can not have more high bits to " 788 "clear than integerPartWidth"); 789 const integerPart HighBitMask = ~integerPart(0) >> NumHighBits; 790 791 if (Parts[PartCount - 1] & HighBitMask) 792 return false; 793 794 return true; 795 } 796 797 bool IEEEFloat::isLargest() const { 798 // The largest number by magnitude in our format will be the floating point 799 // number with maximum exponent and with significand that is all ones. 800 return isFiniteNonZero() && exponent == semantics->maxExponent 801 && isSignificandAllOnes(); 802 } 803 804 bool IEEEFloat::isInteger() const { 805 // This could be made more efficient; I'm going for obviously correct. 806 if (!isFinite()) return false; 807 IEEEFloat truncated = *this; 808 truncated.roundToIntegral(rmTowardZero); 809 return compare(truncated) == cmpEqual; 810 } 811 812 bool IEEEFloat::bitwiseIsEqual(const IEEEFloat &rhs) const { 813 if (this == &rhs) 814 return true; 815 if (semantics != rhs.semantics || 816 category != rhs.category || 817 sign != rhs.sign) 818 return false; 819 if (category==fcZero || category==fcInfinity) 820 return true; 821 822 if (isFiniteNonZero() && exponent != rhs.exponent) 823 return false; 824 825 return std::equal(significandParts(), significandParts() + partCount(), 826 rhs.significandParts()); 827 } 828 829 IEEEFloat::IEEEFloat(const fltSemantics &ourSemantics, integerPart value) { 830 initialize(&ourSemantics); 831 sign = 0; 832 category = fcNormal; 833 zeroSignificand(); 834 exponent = ourSemantics.precision - 1; 835 significandParts()[0] = value; 836 normalize(rmNearestTiesToEven, lfExactlyZero); 837 } 838 839 IEEEFloat::IEEEFloat(const fltSemantics &ourSemantics) { 840 initialize(&ourSemantics); 841 category = fcZero; 842 sign = false; 843 } 844 845 // Delegate to the previous constructor, because later copy constructor may 846 // actually inspects category, which can't be garbage. 847 IEEEFloat::IEEEFloat(const fltSemantics &ourSemantics, uninitializedTag tag) 848 : IEEEFloat(ourSemantics) {} 849 850 IEEEFloat::IEEEFloat(const IEEEFloat &rhs) { 851 initialize(rhs.semantics); 852 assign(rhs); 853 } 854 855 IEEEFloat::IEEEFloat(IEEEFloat &&rhs) : semantics(&semBogus) { 856 *this = std::move(rhs); 857 } 858 859 IEEEFloat::~IEEEFloat() { freeSignificand(); } 860 861 // Profile - This method 'profiles' an APFloat for use with FoldingSet. 862 void IEEEFloat::Profile(FoldingSetNodeID &ID) const { 863 ID.Add(bitcastToAPInt()); 864 } 865 866 unsigned int IEEEFloat::partCount() const { 867 return partCountForBits(semantics->precision + 1); 868 } 869 870 const integerPart *IEEEFloat::significandParts() const { 871 return const_cast<IEEEFloat *>(this)->significandParts(); 872 } 873 874 integerPart *IEEEFloat::significandParts() { 875 if (partCount() > 1) 876 return significand.parts; 877 else 878 return &significand.part; 879 } 880 881 void IEEEFloat::zeroSignificand() { 882 APInt::tcSet(significandParts(), 0, partCount()); 883 } 884 885 /* Increment an fcNormal floating point number's significand. */ 886 void IEEEFloat::incrementSignificand() { 887 integerPart carry; 888 889 carry = APInt::tcIncrement(significandParts(), partCount()); 890 891 /* Our callers should never cause us to overflow. */ 892 assert(carry == 0); 893 (void)carry; 894 } 895 896 /* Add the significand of the RHS. Returns the carry flag. */ 897 integerPart IEEEFloat::addSignificand(const IEEEFloat &rhs) { 898 integerPart *parts; 899 900 parts = significandParts(); 901 902 assert(semantics == rhs.semantics); 903 assert(exponent == rhs.exponent); 904 905 return APInt::tcAdd(parts, rhs.significandParts(), 0, partCount()); 906 } 907 908 /* Subtract the significand of the RHS with a borrow flag. Returns 909 the borrow flag. */ 910 integerPart IEEEFloat::subtractSignificand(const IEEEFloat &rhs, 911 integerPart borrow) { 912 integerPart *parts; 913 914 parts = significandParts(); 915 916 assert(semantics == rhs.semantics); 917 assert(exponent == rhs.exponent); 918 919 return APInt::tcSubtract(parts, rhs.significandParts(), borrow, 920 partCount()); 921 } 922 923 /* Multiply the significand of the RHS. If ADDEND is non-NULL, add it 924 on to the full-precision result of the multiplication. Returns the 925 lost fraction. */ 926 lostFraction IEEEFloat::multiplySignificand(const IEEEFloat &rhs, 927 const IEEEFloat *addend) { 928 unsigned int omsb; // One, not zero, based MSB. 929 unsigned int partsCount, newPartsCount, precision; 930 integerPart *lhsSignificand; 931 integerPart scratch[4]; 932 integerPart *fullSignificand; 933 lostFraction lost_fraction; 934 bool ignored; 935 936 assert(semantics == rhs.semantics); 937 938 precision = semantics->precision; 939 940 // Allocate space for twice as many bits as the original significand, plus one 941 // extra bit for the addition to overflow into. 942 newPartsCount = partCountForBits(precision * 2 + 1); 943 944 if (newPartsCount > 4) 945 fullSignificand = new integerPart[newPartsCount]; 946 else 947 fullSignificand = scratch; 948 949 lhsSignificand = significandParts(); 950 partsCount = partCount(); 951 952 APInt::tcFullMultiply(fullSignificand, lhsSignificand, 953 rhs.significandParts(), partsCount, partsCount); 954 955 lost_fraction = lfExactlyZero; 956 omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1; 957 exponent += rhs.exponent; 958 959 // Assume the operands involved in the multiplication are single-precision 960 // FP, and the two multiplicants are: 961 // *this = a23 . a22 ... a0 * 2^e1 962 // rhs = b23 . b22 ... b0 * 2^e2 963 // the result of multiplication is: 964 // *this = c48 c47 c46 . c45 ... c0 * 2^(e1+e2) 965 // Note that there are three significant bits at the left-hand side of the 966 // radix point: two for the multiplication, and an overflow bit for the 967 // addition (that will always be zero at this point). Move the radix point 968 // toward left by two bits, and adjust exponent accordingly. 969 exponent += 2; 970 971 if (addend && addend->isNonZero()) { 972 // The intermediate result of the multiplication has "2 * precision" 973 // signicant bit; adjust the addend to be consistent with mul result. 974 // 975 Significand savedSignificand = significand; 976 const fltSemantics *savedSemantics = semantics; 977 fltSemantics extendedSemantics; 978 opStatus status; 979 unsigned int extendedPrecision; 980 981 // Normalize our MSB to one below the top bit to allow for overflow. 982 extendedPrecision = 2 * precision + 1; 983 if (omsb != extendedPrecision - 1) { 984 assert(extendedPrecision > omsb); 985 APInt::tcShiftLeft(fullSignificand, newPartsCount, 986 (extendedPrecision - 1) - omsb); 987 exponent -= (extendedPrecision - 1) - omsb; 988 } 989 990 /* Create new semantics. */ 991 extendedSemantics = *semantics; 992 extendedSemantics.precision = extendedPrecision; 993 994 if (newPartsCount == 1) 995 significand.part = fullSignificand[0]; 996 else 997 significand.parts = fullSignificand; 998 semantics = &extendedSemantics; 999 1000 IEEEFloat extendedAddend(*addend); 1001 status = extendedAddend.convert(extendedSemantics, rmTowardZero, &ignored); 1002 assert(status == opOK); 1003 (void)status; 1004 1005 // Shift the significand of the addend right by one bit. This guarantees 1006 // that the high bit of the significand is zero (same as fullSignificand), 1007 // so the addition will overflow (if it does overflow at all) into the top bit. 1008 lost_fraction = extendedAddend.shiftSignificandRight(1); 1009 assert(lost_fraction == lfExactlyZero && 1010 "Lost precision while shifting addend for fused-multiply-add."); 1011 1012 lost_fraction = addOrSubtractSignificand(extendedAddend, false); 1013 1014 /* Restore our state. */ 1015 if (newPartsCount == 1) 1016 fullSignificand[0] = significand.part; 1017 significand = savedSignificand; 1018 semantics = savedSemantics; 1019 1020 omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1; 1021 } 1022 1023 // Convert the result having "2 * precision" significant-bits back to the one 1024 // having "precision" significant-bits. First, move the radix point from 1025 // poision "2*precision - 1" to "precision - 1". The exponent need to be 1026 // adjusted by "2*precision - 1" - "precision - 1" = "precision". 1027 exponent -= precision + 1; 1028 1029 // In case MSB resides at the left-hand side of radix point, shift the 1030 // mantissa right by some amount to make sure the MSB reside right before 1031 // the radix point (i.e. "MSB . rest-significant-bits"). 1032 // 1033 // Note that the result is not normalized when "omsb < precision". So, the 1034 // caller needs to call IEEEFloat::normalize() if normalized value is 1035 // expected. 1036 if (omsb > precision) { 1037 unsigned int bits, significantParts; 1038 lostFraction lf; 1039 1040 bits = omsb - precision; 1041 significantParts = partCountForBits(omsb); 1042 lf = shiftRight(fullSignificand, significantParts, bits); 1043 lost_fraction = combineLostFractions(lf, lost_fraction); 1044 exponent += bits; 1045 } 1046 1047 APInt::tcAssign(lhsSignificand, fullSignificand, partsCount); 1048 1049 if (newPartsCount > 4) 1050 delete [] fullSignificand; 1051 1052 return lost_fraction; 1053 } 1054 1055 /* Multiply the significands of LHS and RHS to DST. */ 1056 lostFraction IEEEFloat::divideSignificand(const IEEEFloat &rhs) { 1057 unsigned int bit, i, partsCount; 1058 const integerPart *rhsSignificand; 1059 integerPart *lhsSignificand, *dividend, *divisor; 1060 integerPart scratch[4]; 1061 lostFraction lost_fraction; 1062 1063 assert(semantics == rhs.semantics); 1064 1065 lhsSignificand = significandParts(); 1066 rhsSignificand = rhs.significandParts(); 1067 partsCount = partCount(); 1068 1069 if (partsCount > 2) 1070 dividend = new integerPart[partsCount * 2]; 1071 else 1072 dividend = scratch; 1073 1074 divisor = dividend + partsCount; 1075 1076 /* Copy the dividend and divisor as they will be modified in-place. */ 1077 for (i = 0; i < partsCount; i++) { 1078 dividend[i] = lhsSignificand[i]; 1079 divisor[i] = rhsSignificand[i]; 1080 lhsSignificand[i] = 0; 1081 } 1082 1083 exponent -= rhs.exponent; 1084 1085 unsigned int precision = semantics->precision; 1086 1087 /* Normalize the divisor. */ 1088 bit = precision - APInt::tcMSB(divisor, partsCount) - 1; 1089 if (bit) { 1090 exponent += bit; 1091 APInt::tcShiftLeft(divisor, partsCount, bit); 1092 } 1093 1094 /* Normalize the dividend. */ 1095 bit = precision - APInt::tcMSB(dividend, partsCount) - 1; 1096 if (bit) { 1097 exponent -= bit; 1098 APInt::tcShiftLeft(dividend, partsCount, bit); 1099 } 1100 1101 /* Ensure the dividend >= divisor initially for the loop below. 1102 Incidentally, this means that the division loop below is 1103 guaranteed to set the integer bit to one. */ 1104 if (APInt::tcCompare(dividend, divisor, partsCount) < 0) { 1105 exponent--; 1106 APInt::tcShiftLeft(dividend, partsCount, 1); 1107 assert(APInt::tcCompare(dividend, divisor, partsCount) >= 0); 1108 } 1109 1110 /* Long division. */ 1111 for (bit = precision; bit; bit -= 1) { 1112 if (APInt::tcCompare(dividend, divisor, partsCount) >= 0) { 1113 APInt::tcSubtract(dividend, divisor, 0, partsCount); 1114 APInt::tcSetBit(lhsSignificand, bit - 1); 1115 } 1116 1117 APInt::tcShiftLeft(dividend, partsCount, 1); 1118 } 1119 1120 /* Figure out the lost fraction. */ 1121 int cmp = APInt::tcCompare(dividend, divisor, partsCount); 1122 1123 if (cmp > 0) 1124 lost_fraction = lfMoreThanHalf; 1125 else if (cmp == 0) 1126 lost_fraction = lfExactlyHalf; 1127 else if (APInt::tcIsZero(dividend, partsCount)) 1128 lost_fraction = lfExactlyZero; 1129 else 1130 lost_fraction = lfLessThanHalf; 1131 1132 if (partsCount > 2) 1133 delete [] dividend; 1134 1135 return lost_fraction; 1136 } 1137 1138 unsigned int IEEEFloat::significandMSB() const { 1139 return APInt::tcMSB(significandParts(), partCount()); 1140 } 1141 1142 unsigned int IEEEFloat::significandLSB() const { 1143 return APInt::tcLSB(significandParts(), partCount()); 1144 } 1145 1146 /* Note that a zero result is NOT normalized to fcZero. */ 1147 lostFraction IEEEFloat::shiftSignificandRight(unsigned int bits) { 1148 /* Our exponent should not overflow. */ 1149 assert((ExponentType) (exponent + bits) >= exponent); 1150 1151 exponent += bits; 1152 1153 return shiftRight(significandParts(), partCount(), bits); 1154 } 1155 1156 /* Shift the significand left BITS bits, subtract BITS from its exponent. */ 1157 void IEEEFloat::shiftSignificandLeft(unsigned int bits) { 1158 assert(bits < semantics->precision); 1159 1160 if (bits) { 1161 unsigned int partsCount = partCount(); 1162 1163 APInt::tcShiftLeft(significandParts(), partsCount, bits); 1164 exponent -= bits; 1165 1166 assert(!APInt::tcIsZero(significandParts(), partsCount)); 1167 } 1168 } 1169 1170 IEEEFloat::cmpResult 1171 IEEEFloat::compareAbsoluteValue(const IEEEFloat &rhs) const { 1172 int compare; 1173 1174 assert(semantics == rhs.semantics); 1175 assert(isFiniteNonZero()); 1176 assert(rhs.isFiniteNonZero()); 1177 1178 compare = exponent - rhs.exponent; 1179 1180 /* If exponents are equal, do an unsigned bignum comparison of the 1181 significands. */ 1182 if (compare == 0) 1183 compare = APInt::tcCompare(significandParts(), rhs.significandParts(), 1184 partCount()); 1185 1186 if (compare > 0) 1187 return cmpGreaterThan; 1188 else if (compare < 0) 1189 return cmpLessThan; 1190 else 1191 return cmpEqual; 1192 } 1193 1194 /* Handle overflow. Sign is preserved. We either become infinity or 1195 the largest finite number. */ 1196 IEEEFloat::opStatus IEEEFloat::handleOverflow(roundingMode rounding_mode) { 1197 /* Infinity? */ 1198 if (rounding_mode == rmNearestTiesToEven || 1199 rounding_mode == rmNearestTiesToAway || 1200 (rounding_mode == rmTowardPositive && !sign) || 1201 (rounding_mode == rmTowardNegative && sign)) { 1202 category = fcInfinity; 1203 return (opStatus) (opOverflow | opInexact); 1204 } 1205 1206 /* Otherwise we become the largest finite number. */ 1207 category = fcNormal; 1208 exponent = semantics->maxExponent; 1209 APInt::tcSetLeastSignificantBits(significandParts(), partCount(), 1210 semantics->precision); 1211 1212 return opInexact; 1213 } 1214 1215 /* Returns TRUE if, when truncating the current number, with BIT the 1216 new LSB, with the given lost fraction and rounding mode, the result 1217 would need to be rounded away from zero (i.e., by increasing the 1218 signficand). This routine must work for fcZero of both signs, and 1219 fcNormal numbers. */ 1220 bool IEEEFloat::roundAwayFromZero(roundingMode rounding_mode, 1221 lostFraction lost_fraction, 1222 unsigned int bit) const { 1223 /* NaNs and infinities should not have lost fractions. */ 1224 assert(isFiniteNonZero() || category == fcZero); 1225 1226 /* Current callers never pass this so we don't handle it. */ 1227 assert(lost_fraction != lfExactlyZero); 1228 1229 switch (rounding_mode) { 1230 case rmNearestTiesToAway: 1231 return lost_fraction == lfExactlyHalf || lost_fraction == lfMoreThanHalf; 1232 1233 case rmNearestTiesToEven: 1234 if (lost_fraction == lfMoreThanHalf) 1235 return true; 1236 1237 /* Our zeroes don't have a significand to test. */ 1238 if (lost_fraction == lfExactlyHalf && category != fcZero) 1239 return APInt::tcExtractBit(significandParts(), bit); 1240 1241 return false; 1242 1243 case rmTowardZero: 1244 return false; 1245 1246 case rmTowardPositive: 1247 return !sign; 1248 1249 case rmTowardNegative: 1250 return sign; 1251 } 1252 llvm_unreachable("Invalid rounding mode found"); 1253 } 1254 1255 IEEEFloat::opStatus IEEEFloat::normalize(roundingMode rounding_mode, 1256 lostFraction lost_fraction) { 1257 unsigned int omsb; /* One, not zero, based MSB. */ 1258 int exponentChange; 1259 1260 if (!isFiniteNonZero()) 1261 return opOK; 1262 1263 /* Before rounding normalize the exponent of fcNormal numbers. */ 1264 omsb = significandMSB() + 1; 1265 1266 if (omsb) { 1267 /* OMSB is numbered from 1. We want to place it in the integer 1268 bit numbered PRECISION if possible, with a compensating change in 1269 the exponent. */ 1270 exponentChange = omsb - semantics->precision; 1271 1272 /* If the resulting exponent is too high, overflow according to 1273 the rounding mode. */ 1274 if (exponent + exponentChange > semantics->maxExponent) 1275 return handleOverflow(rounding_mode); 1276 1277 /* Subnormal numbers have exponent minExponent, and their MSB 1278 is forced based on that. */ 1279 if (exponent + exponentChange < semantics->minExponent) 1280 exponentChange = semantics->minExponent - exponent; 1281 1282 /* Shifting left is easy as we don't lose precision. */ 1283 if (exponentChange < 0) { 1284 assert(lost_fraction == lfExactlyZero); 1285 1286 shiftSignificandLeft(-exponentChange); 1287 1288 return opOK; 1289 } 1290 1291 if (exponentChange > 0) { 1292 lostFraction lf; 1293 1294 /* Shift right and capture any new lost fraction. */ 1295 lf = shiftSignificandRight(exponentChange); 1296 1297 lost_fraction = combineLostFractions(lf, lost_fraction); 1298 1299 /* Keep OMSB up-to-date. */ 1300 if (omsb > (unsigned) exponentChange) 1301 omsb -= exponentChange; 1302 else 1303 omsb = 0; 1304 } 1305 } 1306 1307 /* Now round the number according to rounding_mode given the lost 1308 fraction. */ 1309 1310 /* As specified in IEEE 754, since we do not trap we do not report 1311 underflow for exact results. */ 1312 if (lost_fraction == lfExactlyZero) { 1313 /* Canonicalize zeroes. */ 1314 if (omsb == 0) 1315 category = fcZero; 1316 1317 return opOK; 1318 } 1319 1320 /* Increment the significand if we're rounding away from zero. */ 1321 if (roundAwayFromZero(rounding_mode, lost_fraction, 0)) { 1322 if (omsb == 0) 1323 exponent = semantics->minExponent; 1324 1325 incrementSignificand(); 1326 omsb = significandMSB() + 1; 1327 1328 /* Did the significand increment overflow? */ 1329 if (omsb == (unsigned) semantics->precision + 1) { 1330 /* Renormalize by incrementing the exponent and shifting our 1331 significand right one. However if we already have the 1332 maximum exponent we overflow to infinity. */ 1333 if (exponent == semantics->maxExponent) { 1334 category = fcInfinity; 1335 1336 return (opStatus) (opOverflow | opInexact); 1337 } 1338 1339 shiftSignificandRight(1); 1340 1341 return opInexact; 1342 } 1343 } 1344 1345 /* The normal case - we were and are not denormal, and any 1346 significand increment above didn't overflow. */ 1347 if (omsb == semantics->precision) 1348 return opInexact; 1349 1350 /* We have a non-zero denormal. */ 1351 assert(omsb < semantics->precision); 1352 1353 /* Canonicalize zeroes. */ 1354 if (omsb == 0) 1355 category = fcZero; 1356 1357 /* The fcZero case is a denormal that underflowed to zero. */ 1358 return (opStatus) (opUnderflow | opInexact); 1359 } 1360 1361 IEEEFloat::opStatus IEEEFloat::addOrSubtractSpecials(const IEEEFloat &rhs, 1362 bool subtract) { 1363 switch (PackCategoriesIntoKey(category, rhs.category)) { 1364 default: 1365 llvm_unreachable(nullptr); 1366 1367 case PackCategoriesIntoKey(fcNaN, fcZero): 1368 case PackCategoriesIntoKey(fcNaN, fcNormal): 1369 case PackCategoriesIntoKey(fcNaN, fcInfinity): 1370 case PackCategoriesIntoKey(fcNaN, fcNaN): 1371 case PackCategoriesIntoKey(fcNormal, fcZero): 1372 case PackCategoriesIntoKey(fcInfinity, fcNormal): 1373 case PackCategoriesIntoKey(fcInfinity, fcZero): 1374 return opOK; 1375 1376 case PackCategoriesIntoKey(fcZero, fcNaN): 1377 case PackCategoriesIntoKey(fcNormal, fcNaN): 1378 case PackCategoriesIntoKey(fcInfinity, fcNaN): 1379 // We need to be sure to flip the sign here for subtraction because we 1380 // don't have a separate negate operation so -NaN becomes 0 - NaN here. 1381 sign = rhs.sign ^ subtract; 1382 category = fcNaN; 1383 copySignificand(rhs); 1384 return opOK; 1385 1386 case PackCategoriesIntoKey(fcNormal, fcInfinity): 1387 case PackCategoriesIntoKey(fcZero, fcInfinity): 1388 category = fcInfinity; 1389 sign = rhs.sign ^ subtract; 1390 return opOK; 1391 1392 case PackCategoriesIntoKey(fcZero, fcNormal): 1393 assign(rhs); 1394 sign = rhs.sign ^ subtract; 1395 return opOK; 1396 1397 case PackCategoriesIntoKey(fcZero, fcZero): 1398 /* Sign depends on rounding mode; handled by caller. */ 1399 return opOK; 1400 1401 case PackCategoriesIntoKey(fcInfinity, fcInfinity): 1402 /* Differently signed infinities can only be validly 1403 subtracted. */ 1404 if (((sign ^ rhs.sign)!=0) != subtract) { 1405 makeNaN(); 1406 return opInvalidOp; 1407 } 1408 1409 return opOK; 1410 1411 case PackCategoriesIntoKey(fcNormal, fcNormal): 1412 return opDivByZero; 1413 } 1414 } 1415 1416 /* Add or subtract two normal numbers. */ 1417 lostFraction IEEEFloat::addOrSubtractSignificand(const IEEEFloat &rhs, 1418 bool subtract) { 1419 integerPart carry; 1420 lostFraction lost_fraction; 1421 int bits; 1422 1423 /* Determine if the operation on the absolute values is effectively 1424 an addition or subtraction. */ 1425 subtract ^= static_cast<bool>(sign ^ rhs.sign); 1426 1427 /* Are we bigger exponent-wise than the RHS? */ 1428 bits = exponent - rhs.exponent; 1429 1430 /* Subtraction is more subtle than one might naively expect. */ 1431 if (subtract) { 1432 IEEEFloat temp_rhs(rhs); 1433 bool reverse; 1434 1435 if (bits == 0) { 1436 reverse = compareAbsoluteValue(temp_rhs) == cmpLessThan; 1437 lost_fraction = lfExactlyZero; 1438 } else if (bits > 0) { 1439 lost_fraction = temp_rhs.shiftSignificandRight(bits - 1); 1440 shiftSignificandLeft(1); 1441 reverse = false; 1442 } else { 1443 lost_fraction = shiftSignificandRight(-bits - 1); 1444 temp_rhs.shiftSignificandLeft(1); 1445 reverse = true; 1446 } 1447 1448 if (reverse) { 1449 carry = temp_rhs.subtractSignificand 1450 (*this, lost_fraction != lfExactlyZero); 1451 copySignificand(temp_rhs); 1452 sign = !sign; 1453 } else { 1454 carry = subtractSignificand 1455 (temp_rhs, lost_fraction != lfExactlyZero); 1456 } 1457 1458 /* Invert the lost fraction - it was on the RHS and 1459 subtracted. */ 1460 if (lost_fraction == lfLessThanHalf) 1461 lost_fraction = lfMoreThanHalf; 1462 else if (lost_fraction == lfMoreThanHalf) 1463 lost_fraction = lfLessThanHalf; 1464 1465 /* The code above is intended to ensure that no borrow is 1466 necessary. */ 1467 assert(!carry); 1468 (void)carry; 1469 } else { 1470 if (bits > 0) { 1471 IEEEFloat temp_rhs(rhs); 1472 1473 lost_fraction = temp_rhs.shiftSignificandRight(bits); 1474 carry = addSignificand(temp_rhs); 1475 } else { 1476 lost_fraction = shiftSignificandRight(-bits); 1477 carry = addSignificand(rhs); 1478 } 1479 1480 /* We have a guard bit; generating a carry cannot happen. */ 1481 assert(!carry); 1482 (void)carry; 1483 } 1484 1485 return lost_fraction; 1486 } 1487 1488 IEEEFloat::opStatus IEEEFloat::multiplySpecials(const IEEEFloat &rhs) { 1489 switch (PackCategoriesIntoKey(category, rhs.category)) { 1490 default: 1491 llvm_unreachable(nullptr); 1492 1493 case PackCategoriesIntoKey(fcNaN, fcZero): 1494 case PackCategoriesIntoKey(fcNaN, fcNormal): 1495 case PackCategoriesIntoKey(fcNaN, fcInfinity): 1496 case PackCategoriesIntoKey(fcNaN, fcNaN): 1497 sign = false; 1498 return opOK; 1499 1500 case PackCategoriesIntoKey(fcZero, fcNaN): 1501 case PackCategoriesIntoKey(fcNormal, fcNaN): 1502 case PackCategoriesIntoKey(fcInfinity, fcNaN): 1503 sign = false; 1504 category = fcNaN; 1505 copySignificand(rhs); 1506 return opOK; 1507 1508 case PackCategoriesIntoKey(fcNormal, fcInfinity): 1509 case PackCategoriesIntoKey(fcInfinity, fcNormal): 1510 case PackCategoriesIntoKey(fcInfinity, fcInfinity): 1511 category = fcInfinity; 1512 return opOK; 1513 1514 case PackCategoriesIntoKey(fcZero, fcNormal): 1515 case PackCategoriesIntoKey(fcNormal, fcZero): 1516 case PackCategoriesIntoKey(fcZero, fcZero): 1517 category = fcZero; 1518 return opOK; 1519 1520 case PackCategoriesIntoKey(fcZero, fcInfinity): 1521 case PackCategoriesIntoKey(fcInfinity, fcZero): 1522 makeNaN(); 1523 return opInvalidOp; 1524 1525 case PackCategoriesIntoKey(fcNormal, fcNormal): 1526 return opOK; 1527 } 1528 } 1529 1530 IEEEFloat::opStatus IEEEFloat::divideSpecials(const IEEEFloat &rhs) { 1531 switch (PackCategoriesIntoKey(category, rhs.category)) { 1532 default: 1533 llvm_unreachable(nullptr); 1534 1535 case PackCategoriesIntoKey(fcZero, fcNaN): 1536 case PackCategoriesIntoKey(fcNormal, fcNaN): 1537 case PackCategoriesIntoKey(fcInfinity, fcNaN): 1538 category = fcNaN; 1539 copySignificand(rhs); 1540 case PackCategoriesIntoKey(fcNaN, fcZero): 1541 case PackCategoriesIntoKey(fcNaN, fcNormal): 1542 case PackCategoriesIntoKey(fcNaN, fcInfinity): 1543 case PackCategoriesIntoKey(fcNaN, fcNaN): 1544 sign = false; 1545 case PackCategoriesIntoKey(fcInfinity, fcZero): 1546 case PackCategoriesIntoKey(fcInfinity, fcNormal): 1547 case PackCategoriesIntoKey(fcZero, fcInfinity): 1548 case PackCategoriesIntoKey(fcZero, fcNormal): 1549 return opOK; 1550 1551 case PackCategoriesIntoKey(fcNormal, fcInfinity): 1552 category = fcZero; 1553 return opOK; 1554 1555 case PackCategoriesIntoKey(fcNormal, fcZero): 1556 category = fcInfinity; 1557 return opDivByZero; 1558 1559 case PackCategoriesIntoKey(fcInfinity, fcInfinity): 1560 case PackCategoriesIntoKey(fcZero, fcZero): 1561 makeNaN(); 1562 return opInvalidOp; 1563 1564 case PackCategoriesIntoKey(fcNormal, fcNormal): 1565 return opOK; 1566 } 1567 } 1568 1569 IEEEFloat::opStatus IEEEFloat::modSpecials(const IEEEFloat &rhs) { 1570 switch (PackCategoriesIntoKey(category, rhs.category)) { 1571 default: 1572 llvm_unreachable(nullptr); 1573 1574 case PackCategoriesIntoKey(fcNaN, fcZero): 1575 case PackCategoriesIntoKey(fcNaN, fcNormal): 1576 case PackCategoriesIntoKey(fcNaN, fcInfinity): 1577 case PackCategoriesIntoKey(fcNaN, fcNaN): 1578 case PackCategoriesIntoKey(fcZero, fcInfinity): 1579 case PackCategoriesIntoKey(fcZero, fcNormal): 1580 case PackCategoriesIntoKey(fcNormal, fcInfinity): 1581 return opOK; 1582 1583 case PackCategoriesIntoKey(fcZero, fcNaN): 1584 case PackCategoriesIntoKey(fcNormal, fcNaN): 1585 case PackCategoriesIntoKey(fcInfinity, fcNaN): 1586 sign = false; 1587 category = fcNaN; 1588 copySignificand(rhs); 1589 return opOK; 1590 1591 case PackCategoriesIntoKey(fcNormal, fcZero): 1592 case PackCategoriesIntoKey(fcInfinity, fcZero): 1593 case PackCategoriesIntoKey(fcInfinity, fcNormal): 1594 case PackCategoriesIntoKey(fcInfinity, fcInfinity): 1595 case PackCategoriesIntoKey(fcZero, fcZero): 1596 makeNaN(); 1597 return opInvalidOp; 1598 1599 case PackCategoriesIntoKey(fcNormal, fcNormal): 1600 return opOK; 1601 } 1602 } 1603 1604 /* Change sign. */ 1605 void IEEEFloat::changeSign() { 1606 /* Look mummy, this one's easy. */ 1607 sign = !sign; 1608 } 1609 1610 void IEEEFloat::clearSign() { 1611 /* So is this one. */ 1612 sign = 0; 1613 } 1614 1615 void IEEEFloat::copySign(const IEEEFloat &rhs) { 1616 /* And this one. */ 1617 sign = rhs.sign; 1618 } 1619 1620 /* Normalized addition or subtraction. */ 1621 IEEEFloat::opStatus IEEEFloat::addOrSubtract(const IEEEFloat &rhs, 1622 roundingMode rounding_mode, 1623 bool subtract) { 1624 opStatus fs; 1625 1626 fs = addOrSubtractSpecials(rhs, subtract); 1627 1628 /* This return code means it was not a simple case. */ 1629 if (fs == opDivByZero) { 1630 lostFraction lost_fraction; 1631 1632 lost_fraction = addOrSubtractSignificand(rhs, subtract); 1633 fs = normalize(rounding_mode, lost_fraction); 1634 1635 /* Can only be zero if we lost no fraction. */ 1636 assert(category != fcZero || lost_fraction == lfExactlyZero); 1637 } 1638 1639 /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a 1640 positive zero unless rounding to minus infinity, except that 1641 adding two like-signed zeroes gives that zero. */ 1642 if (category == fcZero) { 1643 if (rhs.category != fcZero || (sign == rhs.sign) == subtract) 1644 sign = (rounding_mode == rmTowardNegative); 1645 } 1646 1647 return fs; 1648 } 1649 1650 /* Normalized addition. */ 1651 IEEEFloat::opStatus IEEEFloat::add(const IEEEFloat &rhs, 1652 roundingMode rounding_mode) { 1653 return addOrSubtract(rhs, rounding_mode, false); 1654 } 1655 1656 /* Normalized subtraction. */ 1657 IEEEFloat::opStatus IEEEFloat::subtract(const IEEEFloat &rhs, 1658 roundingMode rounding_mode) { 1659 return addOrSubtract(rhs, rounding_mode, true); 1660 } 1661 1662 /* Normalized multiply. */ 1663 IEEEFloat::opStatus IEEEFloat::multiply(const IEEEFloat &rhs, 1664 roundingMode rounding_mode) { 1665 opStatus fs; 1666 1667 sign ^= rhs.sign; 1668 fs = multiplySpecials(rhs); 1669 1670 if (isFiniteNonZero()) { 1671 lostFraction lost_fraction = multiplySignificand(rhs, nullptr); 1672 fs = normalize(rounding_mode, lost_fraction); 1673 if (lost_fraction != lfExactlyZero) 1674 fs = (opStatus) (fs | opInexact); 1675 } 1676 1677 return fs; 1678 } 1679 1680 /* Normalized divide. */ 1681 IEEEFloat::opStatus IEEEFloat::divide(const IEEEFloat &rhs, 1682 roundingMode rounding_mode) { 1683 opStatus fs; 1684 1685 sign ^= rhs.sign; 1686 fs = divideSpecials(rhs); 1687 1688 if (isFiniteNonZero()) { 1689 lostFraction lost_fraction = divideSignificand(rhs); 1690 fs = normalize(rounding_mode, lost_fraction); 1691 if (lost_fraction != lfExactlyZero) 1692 fs = (opStatus) (fs | opInexact); 1693 } 1694 1695 return fs; 1696 } 1697 1698 /* Normalized remainder. This is not currently correct in all cases. */ 1699 IEEEFloat::opStatus IEEEFloat::remainder(const IEEEFloat &rhs) { 1700 opStatus fs; 1701 IEEEFloat V = *this; 1702 unsigned int origSign = sign; 1703 1704 fs = V.divide(rhs, rmNearestTiesToEven); 1705 if (fs == opDivByZero) 1706 return fs; 1707 1708 int parts = partCount(); 1709 integerPart *x = new integerPart[parts]; 1710 bool ignored; 1711 fs = V.convertToInteger(x, parts * integerPartWidth, true, 1712 rmNearestTiesToEven, &ignored); 1713 if (fs==opInvalidOp) { 1714 delete[] x; 1715 return fs; 1716 } 1717 1718 fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true, 1719 rmNearestTiesToEven); 1720 assert(fs==opOK); // should always work 1721 1722 fs = V.multiply(rhs, rmNearestTiesToEven); 1723 assert(fs==opOK || fs==opInexact); // should not overflow or underflow 1724 1725 fs = subtract(V, rmNearestTiesToEven); 1726 assert(fs==opOK || fs==opInexact); // likewise 1727 1728 if (isZero()) 1729 sign = origSign; // IEEE754 requires this 1730 delete[] x; 1731 return fs; 1732 } 1733 1734 /* Normalized llvm frem (C fmod). 1735 This is not currently correct in all cases. */ 1736 IEEEFloat::opStatus IEEEFloat::mod(const IEEEFloat &rhs) { 1737 opStatus fs; 1738 fs = modSpecials(rhs); 1739 1740 if (isFiniteNonZero() && rhs.isFiniteNonZero()) { 1741 IEEEFloat V = *this; 1742 unsigned int origSign = sign; 1743 1744 fs = V.divide(rhs, rmNearestTiesToEven); 1745 if (fs == opDivByZero) 1746 return fs; 1747 1748 int parts = partCount(); 1749 integerPart *x = new integerPart[parts]; 1750 bool ignored; 1751 fs = V.convertToInteger(x, parts * integerPartWidth, true, 1752 rmTowardZero, &ignored); 1753 if (fs==opInvalidOp) { 1754 delete[] x; 1755 return fs; 1756 } 1757 1758 fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true, 1759 rmNearestTiesToEven); 1760 assert(fs==opOK); // should always work 1761 1762 fs = V.multiply(rhs, rmNearestTiesToEven); 1763 assert(fs==opOK || fs==opInexact); // should not overflow or underflow 1764 1765 fs = subtract(V, rmNearestTiesToEven); 1766 assert(fs==opOK || fs==opInexact); // likewise 1767 1768 if (isZero()) 1769 sign = origSign; // IEEE754 requires this 1770 delete[] x; 1771 } 1772 return fs; 1773 } 1774 1775 /* Normalized fused-multiply-add. */ 1776 IEEEFloat::opStatus IEEEFloat::fusedMultiplyAdd(const IEEEFloat &multiplicand, 1777 const IEEEFloat &addend, 1778 roundingMode rounding_mode) { 1779 opStatus fs; 1780 1781 /* Post-multiplication sign, before addition. */ 1782 sign ^= multiplicand.sign; 1783 1784 /* If and only if all arguments are normal do we need to do an 1785 extended-precision calculation. */ 1786 if (isFiniteNonZero() && 1787 multiplicand.isFiniteNonZero() && 1788 addend.isFinite()) { 1789 lostFraction lost_fraction; 1790 1791 lost_fraction = multiplySignificand(multiplicand, &addend); 1792 fs = normalize(rounding_mode, lost_fraction); 1793 if (lost_fraction != lfExactlyZero) 1794 fs = (opStatus) (fs | opInexact); 1795 1796 /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a 1797 positive zero unless rounding to minus infinity, except that 1798 adding two like-signed zeroes gives that zero. */ 1799 if (category == fcZero && !(fs & opUnderflow) && sign != addend.sign) 1800 sign = (rounding_mode == rmTowardNegative); 1801 } else { 1802 fs = multiplySpecials(multiplicand); 1803 1804 /* FS can only be opOK or opInvalidOp. There is no more work 1805 to do in the latter case. The IEEE-754R standard says it is 1806 implementation-defined in this case whether, if ADDEND is a 1807 quiet NaN, we raise invalid op; this implementation does so. 1808 1809 If we need to do the addition we can do so with normal 1810 precision. */ 1811 if (fs == opOK) 1812 fs = addOrSubtract(addend, rounding_mode, false); 1813 } 1814 1815 return fs; 1816 } 1817 1818 /* Rounding-mode corrrect round to integral value. */ 1819 IEEEFloat::opStatus IEEEFloat::roundToIntegral(roundingMode rounding_mode) { 1820 opStatus fs; 1821 1822 // If the exponent is large enough, we know that this value is already 1823 // integral, and the arithmetic below would potentially cause it to saturate 1824 // to +/-Inf. Bail out early instead. 1825 if (isFiniteNonZero() && exponent+1 >= (int)semanticsPrecision(*semantics)) 1826 return opOK; 1827 1828 // The algorithm here is quite simple: we add 2^(p-1), where p is the 1829 // precision of our format, and then subtract it back off again. The choice 1830 // of rounding modes for the addition/subtraction determines the rounding mode 1831 // for our integral rounding as well. 1832 // NOTE: When the input value is negative, we do subtraction followed by 1833 // addition instead. 1834 APInt IntegerConstant(NextPowerOf2(semanticsPrecision(*semantics)), 1); 1835 IntegerConstant <<= semanticsPrecision(*semantics)-1; 1836 IEEEFloat MagicConstant(*semantics); 1837 fs = MagicConstant.convertFromAPInt(IntegerConstant, false, 1838 rmNearestTiesToEven); 1839 MagicConstant.copySign(*this); 1840 1841 if (fs != opOK) 1842 return fs; 1843 1844 // Preserve the input sign so that we can handle 0.0/-0.0 cases correctly. 1845 bool inputSign = isNegative(); 1846 1847 fs = add(MagicConstant, rounding_mode); 1848 if (fs != opOK && fs != opInexact) 1849 return fs; 1850 1851 fs = subtract(MagicConstant, rounding_mode); 1852 1853 // Restore the input sign. 1854 if (inputSign != isNegative()) 1855 changeSign(); 1856 1857 return fs; 1858 } 1859 1860 1861 /* Comparison requires normalized numbers. */ 1862 IEEEFloat::cmpResult IEEEFloat::compare(const IEEEFloat &rhs) const { 1863 cmpResult result; 1864 1865 assert(semantics == rhs.semantics); 1866 1867 switch (PackCategoriesIntoKey(category, rhs.category)) { 1868 default: 1869 llvm_unreachable(nullptr); 1870 1871 case PackCategoriesIntoKey(fcNaN, fcZero): 1872 case PackCategoriesIntoKey(fcNaN, fcNormal): 1873 case PackCategoriesIntoKey(fcNaN, fcInfinity): 1874 case PackCategoriesIntoKey(fcNaN, fcNaN): 1875 case PackCategoriesIntoKey(fcZero, fcNaN): 1876 case PackCategoriesIntoKey(fcNormal, fcNaN): 1877 case PackCategoriesIntoKey(fcInfinity, fcNaN): 1878 return cmpUnordered; 1879 1880 case PackCategoriesIntoKey(fcInfinity, fcNormal): 1881 case PackCategoriesIntoKey(fcInfinity, fcZero): 1882 case PackCategoriesIntoKey(fcNormal, fcZero): 1883 if (sign) 1884 return cmpLessThan; 1885 else 1886 return cmpGreaterThan; 1887 1888 case PackCategoriesIntoKey(fcNormal, fcInfinity): 1889 case PackCategoriesIntoKey(fcZero, fcInfinity): 1890 case PackCategoriesIntoKey(fcZero, fcNormal): 1891 if (rhs.sign) 1892 return cmpGreaterThan; 1893 else 1894 return cmpLessThan; 1895 1896 case PackCategoriesIntoKey(fcInfinity, fcInfinity): 1897 if (sign == rhs.sign) 1898 return cmpEqual; 1899 else if (sign) 1900 return cmpLessThan; 1901 else 1902 return cmpGreaterThan; 1903 1904 case PackCategoriesIntoKey(fcZero, fcZero): 1905 return cmpEqual; 1906 1907 case PackCategoriesIntoKey(fcNormal, fcNormal): 1908 break; 1909 } 1910 1911 /* Two normal numbers. Do they have the same sign? */ 1912 if (sign != rhs.sign) { 1913 if (sign) 1914 result = cmpLessThan; 1915 else 1916 result = cmpGreaterThan; 1917 } else { 1918 /* Compare absolute values; invert result if negative. */ 1919 result = compareAbsoluteValue(rhs); 1920 1921 if (sign) { 1922 if (result == cmpLessThan) 1923 result = cmpGreaterThan; 1924 else if (result == cmpGreaterThan) 1925 result = cmpLessThan; 1926 } 1927 } 1928 1929 return result; 1930 } 1931 1932 /// IEEEFloat::convert - convert a value of one floating point type to another. 1933 /// The return value corresponds to the IEEE754 exceptions. *losesInfo 1934 /// records whether the transformation lost information, i.e. whether 1935 /// converting the result back to the original type will produce the 1936 /// original value (this is almost the same as return value==fsOK, but there 1937 /// are edge cases where this is not so). 1938 1939 IEEEFloat::opStatus IEEEFloat::convert(const fltSemantics &toSemantics, 1940 roundingMode rounding_mode, 1941 bool *losesInfo) { 1942 lostFraction lostFraction; 1943 unsigned int newPartCount, oldPartCount; 1944 opStatus fs; 1945 int shift; 1946 const fltSemantics &fromSemantics = *semantics; 1947 1948 lostFraction = lfExactlyZero; 1949 newPartCount = partCountForBits(toSemantics.precision + 1); 1950 oldPartCount = partCount(); 1951 shift = toSemantics.precision - fromSemantics.precision; 1952 1953 bool X86SpecialNan = false; 1954 if (&fromSemantics == &semX87DoubleExtended && 1955 &toSemantics != &semX87DoubleExtended && category == fcNaN && 1956 (!(*significandParts() & 0x8000000000000000ULL) || 1957 !(*significandParts() & 0x4000000000000000ULL))) { 1958 // x86 has some unusual NaNs which cannot be represented in any other 1959 // format; note them here. 1960 X86SpecialNan = true; 1961 } 1962 1963 // If this is a truncation of a denormal number, and the target semantics 1964 // has larger exponent range than the source semantics (this can happen 1965 // when truncating from PowerPC double-double to double format), the 1966 // right shift could lose result mantissa bits. Adjust exponent instead 1967 // of performing excessive shift. 1968 if (shift < 0 && isFiniteNonZero()) { 1969 int exponentChange = significandMSB() + 1 - fromSemantics.precision; 1970 if (exponent + exponentChange < toSemantics.minExponent) 1971 exponentChange = toSemantics.minExponent - exponent; 1972 if (exponentChange < shift) 1973 exponentChange = shift; 1974 if (exponentChange < 0) { 1975 shift -= exponentChange; 1976 exponent += exponentChange; 1977 } 1978 } 1979 1980 // If this is a truncation, perform the shift before we narrow the storage. 1981 if (shift < 0 && (isFiniteNonZero() || category==fcNaN)) 1982 lostFraction = shiftRight(significandParts(), oldPartCount, -shift); 1983 1984 // Fix the storage so it can hold to new value. 1985 if (newPartCount > oldPartCount) { 1986 // The new type requires more storage; make it available. 1987 integerPart *newParts; 1988 newParts = new integerPart[newPartCount]; 1989 APInt::tcSet(newParts, 0, newPartCount); 1990 if (isFiniteNonZero() || category==fcNaN) 1991 APInt::tcAssign(newParts, significandParts(), oldPartCount); 1992 freeSignificand(); 1993 significand.parts = newParts; 1994 } else if (newPartCount == 1 && oldPartCount != 1) { 1995 // Switch to built-in storage for a single part. 1996 integerPart newPart = 0; 1997 if (isFiniteNonZero() || category==fcNaN) 1998 newPart = significandParts()[0]; 1999 freeSignificand(); 2000 significand.part = newPart; 2001 } 2002 2003 // Now that we have the right storage, switch the semantics. 2004 semantics = &toSemantics; 2005 2006 // If this is an extension, perform the shift now that the storage is 2007 // available. 2008 if (shift > 0 && (isFiniteNonZero() || category==fcNaN)) 2009 APInt::tcShiftLeft(significandParts(), newPartCount, shift); 2010 2011 if (isFiniteNonZero()) { 2012 fs = normalize(rounding_mode, lostFraction); 2013 *losesInfo = (fs != opOK); 2014 } else if (category == fcNaN) { 2015 *losesInfo = lostFraction != lfExactlyZero || X86SpecialNan; 2016 2017 // For x87 extended precision, we want to make a NaN, not a special NaN if 2018 // the input wasn't special either. 2019 if (!X86SpecialNan && semantics == &semX87DoubleExtended) 2020 APInt::tcSetBit(significandParts(), semantics->precision - 1); 2021 2022 // gcc forces the Quiet bit on, which means (float)(double)(float_sNan) 2023 // does not give you back the same bits. This is dubious, and we 2024 // don't currently do it. You're really supposed to get 2025 // an invalid operation signal at runtime, but nobody does that. 2026 fs = opOK; 2027 } else { 2028 *losesInfo = false; 2029 fs = opOK; 2030 } 2031 2032 return fs; 2033 } 2034 2035 /* Convert a floating point number to an integer according to the 2036 rounding mode. If the rounded integer value is out of range this 2037 returns an invalid operation exception and the contents of the 2038 destination parts are unspecified. If the rounded value is in 2039 range but the floating point number is not the exact integer, the C 2040 standard doesn't require an inexact exception to be raised. IEEE 2041 854 does require it so we do that. 2042 2043 Note that for conversions to integer type the C standard requires 2044 round-to-zero to always be used. */ 2045 IEEEFloat::opStatus IEEEFloat::convertToSignExtendedInteger( 2046 integerPart *parts, unsigned int width, bool isSigned, 2047 roundingMode rounding_mode, bool *isExact) const { 2048 lostFraction lost_fraction; 2049 const integerPart *src; 2050 unsigned int dstPartsCount, truncatedBits; 2051 2052 *isExact = false; 2053 2054 /* Handle the three special cases first. */ 2055 if (category == fcInfinity || category == fcNaN) 2056 return opInvalidOp; 2057 2058 dstPartsCount = partCountForBits(width); 2059 2060 if (category == fcZero) { 2061 APInt::tcSet(parts, 0, dstPartsCount); 2062 // Negative zero can't be represented as an int. 2063 *isExact = !sign; 2064 return opOK; 2065 } 2066 2067 src = significandParts(); 2068 2069 /* Step 1: place our absolute value, with any fraction truncated, in 2070 the destination. */ 2071 if (exponent < 0) { 2072 /* Our absolute value is less than one; truncate everything. */ 2073 APInt::tcSet(parts, 0, dstPartsCount); 2074 /* For exponent -1 the integer bit represents .5, look at that. 2075 For smaller exponents leftmost truncated bit is 0. */ 2076 truncatedBits = semantics->precision -1U - exponent; 2077 } else { 2078 /* We want the most significant (exponent + 1) bits; the rest are 2079 truncated. */ 2080 unsigned int bits = exponent + 1U; 2081 2082 /* Hopelessly large in magnitude? */ 2083 if (bits > width) 2084 return opInvalidOp; 2085 2086 if (bits < semantics->precision) { 2087 /* We truncate (semantics->precision - bits) bits. */ 2088 truncatedBits = semantics->precision - bits; 2089 APInt::tcExtract(parts, dstPartsCount, src, bits, truncatedBits); 2090 } else { 2091 /* We want at least as many bits as are available. */ 2092 APInt::tcExtract(parts, dstPartsCount, src, semantics->precision, 0); 2093 APInt::tcShiftLeft(parts, dstPartsCount, bits - semantics->precision); 2094 truncatedBits = 0; 2095 } 2096 } 2097 2098 /* Step 2: work out any lost fraction, and increment the absolute 2099 value if we would round away from zero. */ 2100 if (truncatedBits) { 2101 lost_fraction = lostFractionThroughTruncation(src, partCount(), 2102 truncatedBits); 2103 if (lost_fraction != lfExactlyZero && 2104 roundAwayFromZero(rounding_mode, lost_fraction, truncatedBits)) { 2105 if (APInt::tcIncrement(parts, dstPartsCount)) 2106 return opInvalidOp; /* Overflow. */ 2107 } 2108 } else { 2109 lost_fraction = lfExactlyZero; 2110 } 2111 2112 /* Step 3: check if we fit in the destination. */ 2113 unsigned int omsb = APInt::tcMSB(parts, dstPartsCount) + 1; 2114 2115 if (sign) { 2116 if (!isSigned) { 2117 /* Negative numbers cannot be represented as unsigned. */ 2118 if (omsb != 0) 2119 return opInvalidOp; 2120 } else { 2121 /* It takes omsb bits to represent the unsigned integer value. 2122 We lose a bit for the sign, but care is needed as the 2123 maximally negative integer is a special case. */ 2124 if (omsb == width && APInt::tcLSB(parts, dstPartsCount) + 1 != omsb) 2125 return opInvalidOp; 2126 2127 /* This case can happen because of rounding. */ 2128 if (omsb > width) 2129 return opInvalidOp; 2130 } 2131 2132 APInt::tcNegate (parts, dstPartsCount); 2133 } else { 2134 if (omsb >= width + !isSigned) 2135 return opInvalidOp; 2136 } 2137 2138 if (lost_fraction == lfExactlyZero) { 2139 *isExact = true; 2140 return opOK; 2141 } else 2142 return opInexact; 2143 } 2144 2145 /* Same as convertToSignExtendedInteger, except we provide 2146 deterministic values in case of an invalid operation exception, 2147 namely zero for NaNs and the minimal or maximal value respectively 2148 for underflow or overflow. 2149 The *isExact output tells whether the result is exact, in the sense 2150 that converting it back to the original floating point type produces 2151 the original value. This is almost equivalent to result==opOK, 2152 except for negative zeroes. 2153 */ 2154 IEEEFloat::opStatus IEEEFloat::convertToInteger(integerPart *parts, 2155 unsigned int width, 2156 bool isSigned, 2157 roundingMode rounding_mode, 2158 bool *isExact) const { 2159 opStatus fs; 2160 2161 fs = convertToSignExtendedInteger(parts, width, isSigned, rounding_mode, 2162 isExact); 2163 2164 if (fs == opInvalidOp) { 2165 unsigned int bits, dstPartsCount; 2166 2167 dstPartsCount = partCountForBits(width); 2168 2169 if (category == fcNaN) 2170 bits = 0; 2171 else if (sign) 2172 bits = isSigned; 2173 else 2174 bits = width - isSigned; 2175 2176 APInt::tcSetLeastSignificantBits(parts, dstPartsCount, bits); 2177 if (sign && isSigned) 2178 APInt::tcShiftLeft(parts, dstPartsCount, width - 1); 2179 } 2180 2181 return fs; 2182 } 2183 2184 /* Same as convertToInteger(integerPart*, ...), except the result is returned in 2185 an APSInt, whose initial bit-width and signed-ness are used to determine the 2186 precision of the conversion. 2187 */ 2188 IEEEFloat::opStatus IEEEFloat::convertToInteger(APSInt &result, 2189 roundingMode rounding_mode, 2190 bool *isExact) const { 2191 unsigned bitWidth = result.getBitWidth(); 2192 SmallVector<uint64_t, 4> parts(result.getNumWords()); 2193 opStatus status = convertToInteger( 2194 parts.data(), bitWidth, result.isSigned(), rounding_mode, isExact); 2195 // Keeps the original signed-ness. 2196 result = APInt(bitWidth, parts); 2197 return status; 2198 } 2199 2200 /* Convert an unsigned integer SRC to a floating point number, 2201 rounding according to ROUNDING_MODE. The sign of the floating 2202 point number is not modified. */ 2203 IEEEFloat::opStatus IEEEFloat::convertFromUnsignedParts( 2204 const integerPart *src, unsigned int srcCount, roundingMode rounding_mode) { 2205 unsigned int omsb, precision, dstCount; 2206 integerPart *dst; 2207 lostFraction lost_fraction; 2208 2209 category = fcNormal; 2210 omsb = APInt::tcMSB(src, srcCount) + 1; 2211 dst = significandParts(); 2212 dstCount = partCount(); 2213 precision = semantics->precision; 2214 2215 /* We want the most significant PRECISION bits of SRC. There may not 2216 be that many; extract what we can. */ 2217 if (precision <= omsb) { 2218 exponent = omsb - 1; 2219 lost_fraction = lostFractionThroughTruncation(src, srcCount, 2220 omsb - precision); 2221 APInt::tcExtract(dst, dstCount, src, precision, omsb - precision); 2222 } else { 2223 exponent = precision - 1; 2224 lost_fraction = lfExactlyZero; 2225 APInt::tcExtract(dst, dstCount, src, omsb, 0); 2226 } 2227 2228 return normalize(rounding_mode, lost_fraction); 2229 } 2230 2231 IEEEFloat::opStatus IEEEFloat::convertFromAPInt(const APInt &Val, bool isSigned, 2232 roundingMode rounding_mode) { 2233 unsigned int partCount = Val.getNumWords(); 2234 APInt api = Val; 2235 2236 sign = false; 2237 if (isSigned && api.isNegative()) { 2238 sign = true; 2239 api = -api; 2240 } 2241 2242 return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode); 2243 } 2244 2245 /* Convert a two's complement integer SRC to a floating point number, 2246 rounding according to ROUNDING_MODE. ISSIGNED is true if the 2247 integer is signed, in which case it must be sign-extended. */ 2248 IEEEFloat::opStatus 2249 IEEEFloat::convertFromSignExtendedInteger(const integerPart *src, 2250 unsigned int srcCount, bool isSigned, 2251 roundingMode rounding_mode) { 2252 opStatus status; 2253 2254 if (isSigned && 2255 APInt::tcExtractBit(src, srcCount * integerPartWidth - 1)) { 2256 integerPart *copy; 2257 2258 /* If we're signed and negative negate a copy. */ 2259 sign = true; 2260 copy = new integerPart[srcCount]; 2261 APInt::tcAssign(copy, src, srcCount); 2262 APInt::tcNegate(copy, srcCount); 2263 status = convertFromUnsignedParts(copy, srcCount, rounding_mode); 2264 delete [] copy; 2265 } else { 2266 sign = false; 2267 status = convertFromUnsignedParts(src, srcCount, rounding_mode); 2268 } 2269 2270 return status; 2271 } 2272 2273 /* FIXME: should this just take a const APInt reference? */ 2274 IEEEFloat::opStatus 2275 IEEEFloat::convertFromZeroExtendedInteger(const integerPart *parts, 2276 unsigned int width, bool isSigned, 2277 roundingMode rounding_mode) { 2278 unsigned int partCount = partCountForBits(width); 2279 APInt api = APInt(width, makeArrayRef(parts, partCount)); 2280 2281 sign = false; 2282 if (isSigned && APInt::tcExtractBit(parts, width - 1)) { 2283 sign = true; 2284 api = -api; 2285 } 2286 2287 return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode); 2288 } 2289 2290 IEEEFloat::opStatus 2291 IEEEFloat::convertFromHexadecimalString(StringRef s, 2292 roundingMode rounding_mode) { 2293 lostFraction lost_fraction = lfExactlyZero; 2294 2295 category = fcNormal; 2296 zeroSignificand(); 2297 exponent = 0; 2298 2299 integerPart *significand = significandParts(); 2300 unsigned partsCount = partCount(); 2301 unsigned bitPos = partsCount * integerPartWidth; 2302 bool computedTrailingFraction = false; 2303 2304 // Skip leading zeroes and any (hexa)decimal point. 2305 StringRef::iterator begin = s.begin(); 2306 StringRef::iterator end = s.end(); 2307 StringRef::iterator dot; 2308 StringRef::iterator p = skipLeadingZeroesAndAnyDot(begin, end, &dot); 2309 StringRef::iterator firstSignificantDigit = p; 2310 2311 while (p != end) { 2312 integerPart hex_value; 2313 2314 if (*p == '.') { 2315 assert(dot == end && "String contains multiple dots"); 2316 dot = p++; 2317 continue; 2318 } 2319 2320 hex_value = hexDigitValue(*p); 2321 if (hex_value == -1U) 2322 break; 2323 2324 p++; 2325 2326 // Store the number while we have space. 2327 if (bitPos) { 2328 bitPos -= 4; 2329 hex_value <<= bitPos % integerPartWidth; 2330 significand[bitPos / integerPartWidth] |= hex_value; 2331 } else if (!computedTrailingFraction) { 2332 lost_fraction = trailingHexadecimalFraction(p, end, hex_value); 2333 computedTrailingFraction = true; 2334 } 2335 } 2336 2337 /* Hex floats require an exponent but not a hexadecimal point. */ 2338 assert(p != end && "Hex strings require an exponent"); 2339 assert((*p == 'p' || *p == 'P') && "Invalid character in significand"); 2340 assert(p != begin && "Significand has no digits"); 2341 assert((dot == end || p - begin != 1) && "Significand has no digits"); 2342 2343 /* Ignore the exponent if we are zero. */ 2344 if (p != firstSignificantDigit) { 2345 int expAdjustment; 2346 2347 /* Implicit hexadecimal point? */ 2348 if (dot == end) 2349 dot = p; 2350 2351 /* Calculate the exponent adjustment implicit in the number of 2352 significant digits. */ 2353 expAdjustment = static_cast<int>(dot - firstSignificantDigit); 2354 if (expAdjustment < 0) 2355 expAdjustment++; 2356 expAdjustment = expAdjustment * 4 - 1; 2357 2358 /* Adjust for writing the significand starting at the most 2359 significant nibble. */ 2360 expAdjustment += semantics->precision; 2361 expAdjustment -= partsCount * integerPartWidth; 2362 2363 /* Adjust for the given exponent. */ 2364 exponent = totalExponent(p + 1, end, expAdjustment); 2365 } 2366 2367 return normalize(rounding_mode, lost_fraction); 2368 } 2369 2370 IEEEFloat::opStatus 2371 IEEEFloat::roundSignificandWithExponent(const integerPart *decSigParts, 2372 unsigned sigPartCount, int exp, 2373 roundingMode rounding_mode) { 2374 unsigned int parts, pow5PartCount; 2375 fltSemantics calcSemantics = { 32767, -32767, 0, 0 }; 2376 integerPart pow5Parts[maxPowerOfFiveParts]; 2377 bool isNearest; 2378 2379 isNearest = (rounding_mode == rmNearestTiesToEven || 2380 rounding_mode == rmNearestTiesToAway); 2381 2382 parts = partCountForBits(semantics->precision + 11); 2383 2384 /* Calculate pow(5, abs(exp)). */ 2385 pow5PartCount = powerOf5(pow5Parts, exp >= 0 ? exp: -exp); 2386 2387 for (;; parts *= 2) { 2388 opStatus sigStatus, powStatus; 2389 unsigned int excessPrecision, truncatedBits; 2390 2391 calcSemantics.precision = parts * integerPartWidth - 1; 2392 excessPrecision = calcSemantics.precision - semantics->precision; 2393 truncatedBits = excessPrecision; 2394 2395 IEEEFloat decSig(calcSemantics, uninitialized); 2396 decSig.makeZero(sign); 2397 IEEEFloat pow5(calcSemantics); 2398 2399 sigStatus = decSig.convertFromUnsignedParts(decSigParts, sigPartCount, 2400 rmNearestTiesToEven); 2401 powStatus = pow5.convertFromUnsignedParts(pow5Parts, pow5PartCount, 2402 rmNearestTiesToEven); 2403 /* Add exp, as 10^n = 5^n * 2^n. */ 2404 decSig.exponent += exp; 2405 2406 lostFraction calcLostFraction; 2407 integerPart HUerr, HUdistance; 2408 unsigned int powHUerr; 2409 2410 if (exp >= 0) { 2411 /* multiplySignificand leaves the precision-th bit set to 1. */ 2412 calcLostFraction = decSig.multiplySignificand(pow5, nullptr); 2413 powHUerr = powStatus != opOK; 2414 } else { 2415 calcLostFraction = decSig.divideSignificand(pow5); 2416 /* Denormal numbers have less precision. */ 2417 if (decSig.exponent < semantics->minExponent) { 2418 excessPrecision += (semantics->minExponent - decSig.exponent); 2419 truncatedBits = excessPrecision; 2420 if (excessPrecision > calcSemantics.precision) 2421 excessPrecision = calcSemantics.precision; 2422 } 2423 /* Extra half-ulp lost in reciprocal of exponent. */ 2424 powHUerr = (powStatus == opOK && calcLostFraction == lfExactlyZero) ? 0:2; 2425 } 2426 2427 /* Both multiplySignificand and divideSignificand return the 2428 result with the integer bit set. */ 2429 assert(APInt::tcExtractBit 2430 (decSig.significandParts(), calcSemantics.precision - 1) == 1); 2431 2432 HUerr = HUerrBound(calcLostFraction != lfExactlyZero, sigStatus != opOK, 2433 powHUerr); 2434 HUdistance = 2 * ulpsFromBoundary(decSig.significandParts(), 2435 excessPrecision, isNearest); 2436 2437 /* Are we guaranteed to round correctly if we truncate? */ 2438 if (HUdistance >= HUerr) { 2439 APInt::tcExtract(significandParts(), partCount(), decSig.significandParts(), 2440 calcSemantics.precision - excessPrecision, 2441 excessPrecision); 2442 /* Take the exponent of decSig. If we tcExtract-ed less bits 2443 above we must adjust our exponent to compensate for the 2444 implicit right shift. */ 2445 exponent = (decSig.exponent + semantics->precision 2446 - (calcSemantics.precision - excessPrecision)); 2447 calcLostFraction = lostFractionThroughTruncation(decSig.significandParts(), 2448 decSig.partCount(), 2449 truncatedBits); 2450 return normalize(rounding_mode, calcLostFraction); 2451 } 2452 } 2453 } 2454 2455 IEEEFloat::opStatus 2456 IEEEFloat::convertFromDecimalString(StringRef str, roundingMode rounding_mode) { 2457 decimalInfo D; 2458 opStatus fs; 2459 2460 /* Scan the text. */ 2461 StringRef::iterator p = str.begin(); 2462 interpretDecimal(p, str.end(), &D); 2463 2464 /* Handle the quick cases. First the case of no significant digits, 2465 i.e. zero, and then exponents that are obviously too large or too 2466 small. Writing L for log 10 / log 2, a number d.ddddd*10^exp 2467 definitely overflows if 2468 2469 (exp - 1) * L >= maxExponent 2470 2471 and definitely underflows to zero where 2472 2473 (exp + 1) * L <= minExponent - precision 2474 2475 With integer arithmetic the tightest bounds for L are 2476 2477 93/28 < L < 196/59 [ numerator <= 256 ] 2478 42039/12655 < L < 28738/8651 [ numerator <= 65536 ] 2479 */ 2480 2481 // Test if we have a zero number allowing for strings with no null terminators 2482 // and zero decimals with non-zero exponents. 2483 // 2484 // We computed firstSigDigit by ignoring all zeros and dots. Thus if 2485 // D->firstSigDigit equals str.end(), every digit must be a zero and there can 2486 // be at most one dot. On the other hand, if we have a zero with a non-zero 2487 // exponent, then we know that D.firstSigDigit will be non-numeric. 2488 if (D.firstSigDigit == str.end() || decDigitValue(*D.firstSigDigit) >= 10U) { 2489 category = fcZero; 2490 fs = opOK; 2491 2492 /* Check whether the normalized exponent is high enough to overflow 2493 max during the log-rebasing in the max-exponent check below. */ 2494 } else if (D.normalizedExponent - 1 > INT_MAX / 42039) { 2495 fs = handleOverflow(rounding_mode); 2496 2497 /* If it wasn't, then it also wasn't high enough to overflow max 2498 during the log-rebasing in the min-exponent check. Check that it 2499 won't overflow min in either check, then perform the min-exponent 2500 check. */ 2501 } else if (D.normalizedExponent - 1 < INT_MIN / 42039 || 2502 (D.normalizedExponent + 1) * 28738 <= 2503 8651 * (semantics->minExponent - (int) semantics->precision)) { 2504 /* Underflow to zero and round. */ 2505 category = fcNormal; 2506 zeroSignificand(); 2507 fs = normalize(rounding_mode, lfLessThanHalf); 2508 2509 /* We can finally safely perform the max-exponent check. */ 2510 } else if ((D.normalizedExponent - 1) * 42039 2511 >= 12655 * semantics->maxExponent) { 2512 /* Overflow and round. */ 2513 fs = handleOverflow(rounding_mode); 2514 } else { 2515 integerPart *decSignificand; 2516 unsigned int partCount; 2517 2518 /* A tight upper bound on number of bits required to hold an 2519 N-digit decimal integer is N * 196 / 59. Allocate enough space 2520 to hold the full significand, and an extra part required by 2521 tcMultiplyPart. */ 2522 partCount = static_cast<unsigned int>(D.lastSigDigit - D.firstSigDigit) + 1; 2523 partCount = partCountForBits(1 + 196 * partCount / 59); 2524 decSignificand = new integerPart[partCount + 1]; 2525 partCount = 0; 2526 2527 /* Convert to binary efficiently - we do almost all multiplication 2528 in an integerPart. When this would overflow do we do a single 2529 bignum multiplication, and then revert again to multiplication 2530 in an integerPart. */ 2531 do { 2532 integerPart decValue, val, multiplier; 2533 2534 val = 0; 2535 multiplier = 1; 2536 2537 do { 2538 if (*p == '.') { 2539 p++; 2540 if (p == str.end()) { 2541 break; 2542 } 2543 } 2544 decValue = decDigitValue(*p++); 2545 assert(decValue < 10U && "Invalid character in significand"); 2546 multiplier *= 10; 2547 val = val * 10 + decValue; 2548 /* The maximum number that can be multiplied by ten with any 2549 digit added without overflowing an integerPart. */ 2550 } while (p <= D.lastSigDigit && multiplier <= (~ (integerPart) 0 - 9) / 10); 2551 2552 /* Multiply out the current part. */ 2553 APInt::tcMultiplyPart(decSignificand, decSignificand, multiplier, val, 2554 partCount, partCount + 1, false); 2555 2556 /* If we used another part (likely but not guaranteed), increase 2557 the count. */ 2558 if (decSignificand[partCount]) 2559 partCount++; 2560 } while (p <= D.lastSigDigit); 2561 2562 category = fcNormal; 2563 fs = roundSignificandWithExponent(decSignificand, partCount, 2564 D.exponent, rounding_mode); 2565 2566 delete [] decSignificand; 2567 } 2568 2569 return fs; 2570 } 2571 2572 bool IEEEFloat::convertFromStringSpecials(StringRef str) { 2573 if (str.equals("inf") || str.equals("INFINITY")) { 2574 makeInf(false); 2575 return true; 2576 } 2577 2578 if (str.equals("-inf") || str.equals("-INFINITY")) { 2579 makeInf(true); 2580 return true; 2581 } 2582 2583 if (str.equals("nan") || str.equals("NaN")) { 2584 makeNaN(false, false); 2585 return true; 2586 } 2587 2588 if (str.equals("-nan") || str.equals("-NaN")) { 2589 makeNaN(false, true); 2590 return true; 2591 } 2592 2593 return false; 2594 } 2595 2596 IEEEFloat::opStatus IEEEFloat::convertFromString(StringRef str, 2597 roundingMode rounding_mode) { 2598 assert(!str.empty() && "Invalid string length"); 2599 2600 // Handle special cases. 2601 if (convertFromStringSpecials(str)) 2602 return opOK; 2603 2604 /* Handle a leading minus sign. */ 2605 StringRef::iterator p = str.begin(); 2606 size_t slen = str.size(); 2607 sign = *p == '-' ? 1 : 0; 2608 if (*p == '-' || *p == '+') { 2609 p++; 2610 slen--; 2611 assert(slen && "String has no digits"); 2612 } 2613 2614 if (slen >= 2 && p[0] == '0' && (p[1] == 'x' || p[1] == 'X')) { 2615 assert(slen - 2 && "Invalid string"); 2616 return convertFromHexadecimalString(StringRef(p + 2, slen - 2), 2617 rounding_mode); 2618 } 2619 2620 return convertFromDecimalString(StringRef(p, slen), rounding_mode); 2621 } 2622 2623 /* Write out a hexadecimal representation of the floating point value 2624 to DST, which must be of sufficient size, in the C99 form 2625 [-]0xh.hhhhp[+-]d. Return the number of characters written, 2626 excluding the terminating NUL. 2627 2628 If UPPERCASE, the output is in upper case, otherwise in lower case. 2629 2630 HEXDIGITS digits appear altogether, rounding the value if 2631 necessary. If HEXDIGITS is 0, the minimal precision to display the 2632 number precisely is used instead. If nothing would appear after 2633 the decimal point it is suppressed. 2634 2635 The decimal exponent is always printed and has at least one digit. 2636 Zero values display an exponent of zero. Infinities and NaNs 2637 appear as "infinity" or "nan" respectively. 2638 2639 The above rules are as specified by C99. There is ambiguity about 2640 what the leading hexadecimal digit should be. This implementation 2641 uses whatever is necessary so that the exponent is displayed as 2642 stored. This implies the exponent will fall within the IEEE format 2643 range, and the leading hexadecimal digit will be 0 (for denormals), 2644 1 (normal numbers) or 2 (normal numbers rounded-away-from-zero with 2645 any other digits zero). 2646 */ 2647 unsigned int IEEEFloat::convertToHexString(char *dst, unsigned int hexDigits, 2648 bool upperCase, 2649 roundingMode rounding_mode) const { 2650 char *p; 2651 2652 p = dst; 2653 if (sign) 2654 *dst++ = '-'; 2655 2656 switch (category) { 2657 case fcInfinity: 2658 memcpy (dst, upperCase ? infinityU: infinityL, sizeof infinityU - 1); 2659 dst += sizeof infinityL - 1; 2660 break; 2661 2662 case fcNaN: 2663 memcpy (dst, upperCase ? NaNU: NaNL, sizeof NaNU - 1); 2664 dst += sizeof NaNU - 1; 2665 break; 2666 2667 case fcZero: 2668 *dst++ = '0'; 2669 *dst++ = upperCase ? 'X': 'x'; 2670 *dst++ = '0'; 2671 if (hexDigits > 1) { 2672 *dst++ = '.'; 2673 memset (dst, '0', hexDigits - 1); 2674 dst += hexDigits - 1; 2675 } 2676 *dst++ = upperCase ? 'P': 'p'; 2677 *dst++ = '0'; 2678 break; 2679 2680 case fcNormal: 2681 dst = convertNormalToHexString (dst, hexDigits, upperCase, rounding_mode); 2682 break; 2683 } 2684 2685 *dst = 0; 2686 2687 return static_cast<unsigned int>(dst - p); 2688 } 2689 2690 /* Does the hard work of outputting the correctly rounded hexadecimal 2691 form of a normal floating point number with the specified number of 2692 hexadecimal digits. If HEXDIGITS is zero the minimum number of 2693 digits necessary to print the value precisely is output. */ 2694 char *IEEEFloat::convertNormalToHexString(char *dst, unsigned int hexDigits, 2695 bool upperCase, 2696 roundingMode rounding_mode) const { 2697 unsigned int count, valueBits, shift, partsCount, outputDigits; 2698 const char *hexDigitChars; 2699 const integerPart *significand; 2700 char *p; 2701 bool roundUp; 2702 2703 *dst++ = '0'; 2704 *dst++ = upperCase ? 'X': 'x'; 2705 2706 roundUp = false; 2707 hexDigitChars = upperCase ? hexDigitsUpper: hexDigitsLower; 2708 2709 significand = significandParts(); 2710 partsCount = partCount(); 2711 2712 /* +3 because the first digit only uses the single integer bit, so 2713 we have 3 virtual zero most-significant-bits. */ 2714 valueBits = semantics->precision + 3; 2715 shift = integerPartWidth - valueBits % integerPartWidth; 2716 2717 /* The natural number of digits required ignoring trailing 2718 insignificant zeroes. */ 2719 outputDigits = (valueBits - significandLSB () + 3) / 4; 2720 2721 /* hexDigits of zero means use the required number for the 2722 precision. Otherwise, see if we are truncating. If we are, 2723 find out if we need to round away from zero. */ 2724 if (hexDigits) { 2725 if (hexDigits < outputDigits) { 2726 /* We are dropping non-zero bits, so need to check how to round. 2727 "bits" is the number of dropped bits. */ 2728 unsigned int bits; 2729 lostFraction fraction; 2730 2731 bits = valueBits - hexDigits * 4; 2732 fraction = lostFractionThroughTruncation (significand, partsCount, bits); 2733 roundUp = roundAwayFromZero(rounding_mode, fraction, bits); 2734 } 2735 outputDigits = hexDigits; 2736 } 2737 2738 /* Write the digits consecutively, and start writing in the location 2739 of the hexadecimal point. We move the most significant digit 2740 left and add the hexadecimal point later. */ 2741 p = ++dst; 2742 2743 count = (valueBits + integerPartWidth - 1) / integerPartWidth; 2744 2745 while (outputDigits && count) { 2746 integerPart part; 2747 2748 /* Put the most significant integerPartWidth bits in "part". */ 2749 if (--count == partsCount) 2750 part = 0; /* An imaginary higher zero part. */ 2751 else 2752 part = significand[count] << shift; 2753 2754 if (count && shift) 2755 part |= significand[count - 1] >> (integerPartWidth - shift); 2756 2757 /* Convert as much of "part" to hexdigits as we can. */ 2758 unsigned int curDigits = integerPartWidth / 4; 2759 2760 if (curDigits > outputDigits) 2761 curDigits = outputDigits; 2762 dst += partAsHex (dst, part, curDigits, hexDigitChars); 2763 outputDigits -= curDigits; 2764 } 2765 2766 if (roundUp) { 2767 char *q = dst; 2768 2769 /* Note that hexDigitChars has a trailing '0'. */ 2770 do { 2771 q--; 2772 *q = hexDigitChars[hexDigitValue (*q) + 1]; 2773 } while (*q == '0'); 2774 assert(q >= p); 2775 } else { 2776 /* Add trailing zeroes. */ 2777 memset (dst, '0', outputDigits); 2778 dst += outputDigits; 2779 } 2780 2781 /* Move the most significant digit to before the point, and if there 2782 is something after the decimal point add it. This must come 2783 after rounding above. */ 2784 p[-1] = p[0]; 2785 if (dst -1 == p) 2786 dst--; 2787 else 2788 p[0] = '.'; 2789 2790 /* Finally output the exponent. */ 2791 *dst++ = upperCase ? 'P': 'p'; 2792 2793 return writeSignedDecimal (dst, exponent); 2794 } 2795 2796 hash_code hash_value(const IEEEFloat &Arg) { 2797 if (!Arg.isFiniteNonZero()) 2798 return hash_combine((uint8_t)Arg.category, 2799 // NaN has no sign, fix it at zero. 2800 Arg.isNaN() ? (uint8_t)0 : (uint8_t)Arg.sign, 2801 Arg.semantics->precision); 2802 2803 // Normal floats need their exponent and significand hashed. 2804 return hash_combine((uint8_t)Arg.category, (uint8_t)Arg.sign, 2805 Arg.semantics->precision, Arg.exponent, 2806 hash_combine_range( 2807 Arg.significandParts(), 2808 Arg.significandParts() + Arg.partCount())); 2809 } 2810 2811 // Conversion from APFloat to/from host float/double. It may eventually be 2812 // possible to eliminate these and have everybody deal with APFloats, but that 2813 // will take a while. This approach will not easily extend to long double. 2814 // Current implementation requires integerPartWidth==64, which is correct at 2815 // the moment but could be made more general. 2816 2817 // Denormals have exponent minExponent in APFloat, but minExponent-1 in 2818 // the actual IEEE respresentations. We compensate for that here. 2819 2820 APInt IEEEFloat::convertF80LongDoubleAPFloatToAPInt() const { 2821 assert(semantics == (const llvm::fltSemantics*)&semX87DoubleExtended); 2822 assert(partCount()==2); 2823 2824 uint64_t myexponent, mysignificand; 2825 2826 if (isFiniteNonZero()) { 2827 myexponent = exponent+16383; //bias 2828 mysignificand = significandParts()[0]; 2829 if (myexponent==1 && !(mysignificand & 0x8000000000000000ULL)) 2830 myexponent = 0; // denormal 2831 } else if (category==fcZero) { 2832 myexponent = 0; 2833 mysignificand = 0; 2834 } else if (category==fcInfinity) { 2835 myexponent = 0x7fff; 2836 mysignificand = 0x8000000000000000ULL; 2837 } else { 2838 assert(category == fcNaN && "Unknown category"); 2839 myexponent = 0x7fff; 2840 mysignificand = significandParts()[0]; 2841 } 2842 2843 uint64_t words[2]; 2844 words[0] = mysignificand; 2845 words[1] = ((uint64_t)(sign & 1) << 15) | 2846 (myexponent & 0x7fffLL); 2847 return APInt(80, words); 2848 } 2849 2850 APInt IEEEFloat::convertPPCDoubleDoubleAPFloatToAPInt() const { 2851 assert(semantics == (const llvm::fltSemantics *)&semPPCDoubleDoubleImpl); 2852 assert(partCount()==2); 2853 2854 uint64_t words[2]; 2855 opStatus fs; 2856 bool losesInfo; 2857 2858 // Convert number to double. To avoid spurious underflows, we re- 2859 // normalize against the "double" minExponent first, and only *then* 2860 // truncate the mantissa. The result of that second conversion 2861 // may be inexact, but should never underflow. 2862 // Declare fltSemantics before APFloat that uses it (and 2863 // saves pointer to it) to ensure correct destruction order. 2864 fltSemantics extendedSemantics = *semantics; 2865 extendedSemantics.minExponent = semIEEEdouble.minExponent; 2866 IEEEFloat extended(*this); 2867 fs = extended.convert(extendedSemantics, rmNearestTiesToEven, &losesInfo); 2868 assert(fs == opOK && !losesInfo); 2869 (void)fs; 2870 2871 IEEEFloat u(extended); 2872 fs = u.convert(semIEEEdouble, rmNearestTiesToEven, &losesInfo); 2873 assert(fs == opOK || fs == opInexact); 2874 (void)fs; 2875 words[0] = *u.convertDoubleAPFloatToAPInt().getRawData(); 2876 2877 // If conversion was exact or resulted in a special case, we're done; 2878 // just set the second double to zero. Otherwise, re-convert back to 2879 // the extended format and compute the difference. This now should 2880 // convert exactly to double. 2881 if (u.isFiniteNonZero() && losesInfo) { 2882 fs = u.convert(extendedSemantics, rmNearestTiesToEven, &losesInfo); 2883 assert(fs == opOK && !losesInfo); 2884 (void)fs; 2885 2886 IEEEFloat v(extended); 2887 v.subtract(u, rmNearestTiesToEven); 2888 fs = v.convert(semIEEEdouble, rmNearestTiesToEven, &losesInfo); 2889 assert(fs == opOK && !losesInfo); 2890 (void)fs; 2891 words[1] = *v.convertDoubleAPFloatToAPInt().getRawData(); 2892 } else { 2893 words[1] = 0; 2894 } 2895 2896 return APInt(128, words); 2897 } 2898 2899 APInt IEEEFloat::convertQuadrupleAPFloatToAPInt() const { 2900 assert(semantics == (const llvm::fltSemantics*)&semIEEEquad); 2901 assert(partCount()==2); 2902 2903 uint64_t myexponent, mysignificand, mysignificand2; 2904 2905 if (isFiniteNonZero()) { 2906 myexponent = exponent+16383; //bias 2907 mysignificand = significandParts()[0]; 2908 mysignificand2 = significandParts()[1]; 2909 if (myexponent==1 && !(mysignificand2 & 0x1000000000000LL)) 2910 myexponent = 0; // denormal 2911 } else if (category==fcZero) { 2912 myexponent = 0; 2913 mysignificand = mysignificand2 = 0; 2914 } else if (category==fcInfinity) { 2915 myexponent = 0x7fff; 2916 mysignificand = mysignificand2 = 0; 2917 } else { 2918 assert(category == fcNaN && "Unknown category!"); 2919 myexponent = 0x7fff; 2920 mysignificand = significandParts()[0]; 2921 mysignificand2 = significandParts()[1]; 2922 } 2923 2924 uint64_t words[2]; 2925 words[0] = mysignificand; 2926 words[1] = ((uint64_t)(sign & 1) << 63) | 2927 ((myexponent & 0x7fff) << 48) | 2928 (mysignificand2 & 0xffffffffffffLL); 2929 2930 return APInt(128, words); 2931 } 2932 2933 APInt IEEEFloat::convertDoubleAPFloatToAPInt() const { 2934 assert(semantics == (const llvm::fltSemantics*)&semIEEEdouble); 2935 assert(partCount()==1); 2936 2937 uint64_t myexponent, mysignificand; 2938 2939 if (isFiniteNonZero()) { 2940 myexponent = exponent+1023; //bias 2941 mysignificand = *significandParts(); 2942 if (myexponent==1 && !(mysignificand & 0x10000000000000LL)) 2943 myexponent = 0; // denormal 2944 } else if (category==fcZero) { 2945 myexponent = 0; 2946 mysignificand = 0; 2947 } else if (category==fcInfinity) { 2948 myexponent = 0x7ff; 2949 mysignificand = 0; 2950 } else { 2951 assert(category == fcNaN && "Unknown category!"); 2952 myexponent = 0x7ff; 2953 mysignificand = *significandParts(); 2954 } 2955 2956 return APInt(64, ((((uint64_t)(sign & 1) << 63) | 2957 ((myexponent & 0x7ff) << 52) | 2958 (mysignificand & 0xfffffffffffffLL)))); 2959 } 2960 2961 APInt IEEEFloat::convertFloatAPFloatToAPInt() const { 2962 assert(semantics == (const llvm::fltSemantics*)&semIEEEsingle); 2963 assert(partCount()==1); 2964 2965 uint32_t myexponent, mysignificand; 2966 2967 if (isFiniteNonZero()) { 2968 myexponent = exponent+127; //bias 2969 mysignificand = (uint32_t)*significandParts(); 2970 if (myexponent == 1 && !(mysignificand & 0x800000)) 2971 myexponent = 0; // denormal 2972 } else if (category==fcZero) { 2973 myexponent = 0; 2974 mysignificand = 0; 2975 } else if (category==fcInfinity) { 2976 myexponent = 0xff; 2977 mysignificand = 0; 2978 } else { 2979 assert(category == fcNaN && "Unknown category!"); 2980 myexponent = 0xff; 2981 mysignificand = (uint32_t)*significandParts(); 2982 } 2983 2984 return APInt(32, (((sign&1) << 31) | ((myexponent&0xff) << 23) | 2985 (mysignificand & 0x7fffff))); 2986 } 2987 2988 APInt IEEEFloat::convertHalfAPFloatToAPInt() const { 2989 assert(semantics == (const llvm::fltSemantics*)&semIEEEhalf); 2990 assert(partCount()==1); 2991 2992 uint32_t myexponent, mysignificand; 2993 2994 if (isFiniteNonZero()) { 2995 myexponent = exponent+15; //bias 2996 mysignificand = (uint32_t)*significandParts(); 2997 if (myexponent == 1 && !(mysignificand & 0x400)) 2998 myexponent = 0; // denormal 2999 } else if (category==fcZero) { 3000 myexponent = 0; 3001 mysignificand = 0; 3002 } else if (category==fcInfinity) { 3003 myexponent = 0x1f; 3004 mysignificand = 0; 3005 } else { 3006 assert(category == fcNaN && "Unknown category!"); 3007 myexponent = 0x1f; 3008 mysignificand = (uint32_t)*significandParts(); 3009 } 3010 3011 return APInt(16, (((sign&1) << 15) | ((myexponent&0x1f) << 10) | 3012 (mysignificand & 0x3ff))); 3013 } 3014 3015 // This function creates an APInt that is just a bit map of the floating 3016 // point constant as it would appear in memory. It is not a conversion, 3017 // and treating the result as a normal integer is unlikely to be useful. 3018 3019 APInt IEEEFloat::bitcastToAPInt() const { 3020 if (semantics == (const llvm::fltSemantics*)&semIEEEhalf) 3021 return convertHalfAPFloatToAPInt(); 3022 3023 if (semantics == (const llvm::fltSemantics*)&semIEEEsingle) 3024 return convertFloatAPFloatToAPInt(); 3025 3026 if (semantics == (const llvm::fltSemantics*)&semIEEEdouble) 3027 return convertDoubleAPFloatToAPInt(); 3028 3029 if (semantics == (const llvm::fltSemantics*)&semIEEEquad) 3030 return convertQuadrupleAPFloatToAPInt(); 3031 3032 if (semantics == (const llvm::fltSemantics *)&semPPCDoubleDoubleImpl) 3033 return convertPPCDoubleDoubleAPFloatToAPInt(); 3034 3035 assert(semantics == (const llvm::fltSemantics*)&semX87DoubleExtended && 3036 "unknown format!"); 3037 return convertF80LongDoubleAPFloatToAPInt(); 3038 } 3039 3040 float IEEEFloat::convertToFloat() const { 3041 assert(semantics == (const llvm::fltSemantics*)&semIEEEsingle && 3042 "Float semantics are not IEEEsingle"); 3043 APInt api = bitcastToAPInt(); 3044 return api.bitsToFloat(); 3045 } 3046 3047 double IEEEFloat::convertToDouble() const { 3048 assert(semantics == (const llvm::fltSemantics*)&semIEEEdouble && 3049 "Float semantics are not IEEEdouble"); 3050 APInt api = bitcastToAPInt(); 3051 return api.bitsToDouble(); 3052 } 3053 3054 /// Integer bit is explicit in this format. Intel hardware (387 and later) 3055 /// does not support these bit patterns: 3056 /// exponent = all 1's, integer bit 0, significand 0 ("pseudoinfinity") 3057 /// exponent = all 1's, integer bit 0, significand nonzero ("pseudoNaN") 3058 /// exponent = 0, integer bit 1 ("pseudodenormal") 3059 /// exponent!=0 nor all 1's, integer bit 0 ("unnormal") 3060 /// At the moment, the first two are treated as NaNs, the second two as Normal. 3061 void IEEEFloat::initFromF80LongDoubleAPInt(const APInt &api) { 3062 assert(api.getBitWidth()==80); 3063 uint64_t i1 = api.getRawData()[0]; 3064 uint64_t i2 = api.getRawData()[1]; 3065 uint64_t myexponent = (i2 & 0x7fff); 3066 uint64_t mysignificand = i1; 3067 3068 initialize(&semX87DoubleExtended); 3069 assert(partCount()==2); 3070 3071 sign = static_cast<unsigned int>(i2>>15); 3072 if (myexponent==0 && mysignificand==0) { 3073 // exponent, significand meaningless 3074 category = fcZero; 3075 } else if (myexponent==0x7fff && mysignificand==0x8000000000000000ULL) { 3076 // exponent, significand meaningless 3077 category = fcInfinity; 3078 } else if (myexponent==0x7fff && mysignificand!=0x8000000000000000ULL) { 3079 // exponent meaningless 3080 category = fcNaN; 3081 significandParts()[0] = mysignificand; 3082 significandParts()[1] = 0; 3083 } else { 3084 category = fcNormal; 3085 exponent = myexponent - 16383; 3086 significandParts()[0] = mysignificand; 3087 significandParts()[1] = 0; 3088 if (myexponent==0) // denormal 3089 exponent = -16382; 3090 } 3091 } 3092 3093 void IEEEFloat::initFromPPCDoubleDoubleAPInt(const APInt &api) { 3094 assert(api.getBitWidth()==128); 3095 uint64_t i1 = api.getRawData()[0]; 3096 uint64_t i2 = api.getRawData()[1]; 3097 opStatus fs; 3098 bool losesInfo; 3099 3100 // Get the first double and convert to our format. 3101 initFromDoubleAPInt(APInt(64, i1)); 3102 fs = convert(semPPCDoubleDoubleImpl, rmNearestTiesToEven, &losesInfo); 3103 assert(fs == opOK && !losesInfo); 3104 (void)fs; 3105 3106 // Unless we have a special case, add in second double. 3107 if (isFiniteNonZero()) { 3108 IEEEFloat v(semIEEEdouble, APInt(64, i2)); 3109 fs = v.convert(semPPCDoubleDoubleImpl, rmNearestTiesToEven, &losesInfo); 3110 assert(fs == opOK && !losesInfo); 3111 (void)fs; 3112 3113 add(v, rmNearestTiesToEven); 3114 } 3115 } 3116 3117 void IEEEFloat::initFromQuadrupleAPInt(const APInt &api) { 3118 assert(api.getBitWidth()==128); 3119 uint64_t i1 = api.getRawData()[0]; 3120 uint64_t i2 = api.getRawData()[1]; 3121 uint64_t myexponent = (i2 >> 48) & 0x7fff; 3122 uint64_t mysignificand = i1; 3123 uint64_t mysignificand2 = i2 & 0xffffffffffffLL; 3124 3125 initialize(&semIEEEquad); 3126 assert(partCount()==2); 3127 3128 sign = static_cast<unsigned int>(i2>>63); 3129 if (myexponent==0 && 3130 (mysignificand==0 && mysignificand2==0)) { 3131 // exponent, significand meaningless 3132 category = fcZero; 3133 } else if (myexponent==0x7fff && 3134 (mysignificand==0 && mysignificand2==0)) { 3135 // exponent, significand meaningless 3136 category = fcInfinity; 3137 } else if (myexponent==0x7fff && 3138 (mysignificand!=0 || mysignificand2 !=0)) { 3139 // exponent meaningless 3140 category = fcNaN; 3141 significandParts()[0] = mysignificand; 3142 significandParts()[1] = mysignificand2; 3143 } else { 3144 category = fcNormal; 3145 exponent = myexponent - 16383; 3146 significandParts()[0] = mysignificand; 3147 significandParts()[1] = mysignificand2; 3148 if (myexponent==0) // denormal 3149 exponent = -16382; 3150 else 3151 significandParts()[1] |= 0x1000000000000LL; // integer bit 3152 } 3153 } 3154 3155 void IEEEFloat::initFromDoubleAPInt(const APInt &api) { 3156 assert(api.getBitWidth()==64); 3157 uint64_t i = *api.getRawData(); 3158 uint64_t myexponent = (i >> 52) & 0x7ff; 3159 uint64_t mysignificand = i & 0xfffffffffffffLL; 3160 3161 initialize(&semIEEEdouble); 3162 assert(partCount()==1); 3163 3164 sign = static_cast<unsigned int>(i>>63); 3165 if (myexponent==0 && mysignificand==0) { 3166 // exponent, significand meaningless 3167 category = fcZero; 3168 } else if (myexponent==0x7ff && mysignificand==0) { 3169 // exponent, significand meaningless 3170 category = fcInfinity; 3171 } else if (myexponent==0x7ff && mysignificand!=0) { 3172 // exponent meaningless 3173 category = fcNaN; 3174 *significandParts() = mysignificand; 3175 } else { 3176 category = fcNormal; 3177 exponent = myexponent - 1023; 3178 *significandParts() = mysignificand; 3179 if (myexponent==0) // denormal 3180 exponent = -1022; 3181 else 3182 *significandParts() |= 0x10000000000000LL; // integer bit 3183 } 3184 } 3185 3186 void IEEEFloat::initFromFloatAPInt(const APInt &api) { 3187 assert(api.getBitWidth()==32); 3188 uint32_t i = (uint32_t)*api.getRawData(); 3189 uint32_t myexponent = (i >> 23) & 0xff; 3190 uint32_t mysignificand = i & 0x7fffff; 3191 3192 initialize(&semIEEEsingle); 3193 assert(partCount()==1); 3194 3195 sign = i >> 31; 3196 if (myexponent==0 && mysignificand==0) { 3197 // exponent, significand meaningless 3198 category = fcZero; 3199 } else if (myexponent==0xff && mysignificand==0) { 3200 // exponent, significand meaningless 3201 category = fcInfinity; 3202 } else if (myexponent==0xff && mysignificand!=0) { 3203 // sign, exponent, significand meaningless 3204 category = fcNaN; 3205 *significandParts() = mysignificand; 3206 } else { 3207 category = fcNormal; 3208 exponent = myexponent - 127; //bias 3209 *significandParts() = mysignificand; 3210 if (myexponent==0) // denormal 3211 exponent = -126; 3212 else 3213 *significandParts() |= 0x800000; // integer bit 3214 } 3215 } 3216 3217 void IEEEFloat::initFromHalfAPInt(const APInt &api) { 3218 assert(api.getBitWidth()==16); 3219 uint32_t i = (uint32_t)*api.getRawData(); 3220 uint32_t myexponent = (i >> 10) & 0x1f; 3221 uint32_t mysignificand = i & 0x3ff; 3222 3223 initialize(&semIEEEhalf); 3224 assert(partCount()==1); 3225 3226 sign = i >> 15; 3227 if (myexponent==0 && mysignificand==0) { 3228 // exponent, significand meaningless 3229 category = fcZero; 3230 } else if (myexponent==0x1f && mysignificand==0) { 3231 // exponent, significand meaningless 3232 category = fcInfinity; 3233 } else if (myexponent==0x1f && mysignificand!=0) { 3234 // sign, exponent, significand meaningless 3235 category = fcNaN; 3236 *significandParts() = mysignificand; 3237 } else { 3238 category = fcNormal; 3239 exponent = myexponent - 15; //bias 3240 *significandParts() = mysignificand; 3241 if (myexponent==0) // denormal 3242 exponent = -14; 3243 else 3244 *significandParts() |= 0x400; // integer bit 3245 } 3246 } 3247 3248 /// Treat api as containing the bits of a floating point number. Currently 3249 /// we infer the floating point type from the size of the APInt. The 3250 /// isIEEE argument distinguishes between PPC128 and IEEE128 (not meaningful 3251 /// when the size is anything else). 3252 void IEEEFloat::initFromAPInt(const fltSemantics *Sem, const APInt &api) { 3253 if (Sem == &semIEEEhalf) 3254 return initFromHalfAPInt(api); 3255 if (Sem == &semIEEEsingle) 3256 return initFromFloatAPInt(api); 3257 if (Sem == &semIEEEdouble) 3258 return initFromDoubleAPInt(api); 3259 if (Sem == &semX87DoubleExtended) 3260 return initFromF80LongDoubleAPInt(api); 3261 if (Sem == &semIEEEquad) 3262 return initFromQuadrupleAPInt(api); 3263 if (Sem == &semPPCDoubleDoubleImpl) 3264 return initFromPPCDoubleDoubleAPInt(api); 3265 3266 llvm_unreachable(nullptr); 3267 } 3268 3269 /// Make this number the largest magnitude normal number in the given 3270 /// semantics. 3271 void IEEEFloat::makeLargest(bool Negative) { 3272 // We want (in interchange format): 3273 // sign = {Negative} 3274 // exponent = 1..10 3275 // significand = 1..1 3276 category = fcNormal; 3277 sign = Negative; 3278 exponent = semantics->maxExponent; 3279 3280 // Use memset to set all but the highest integerPart to all ones. 3281 integerPart *significand = significandParts(); 3282 unsigned PartCount = partCount(); 3283 memset(significand, 0xFF, sizeof(integerPart)*(PartCount - 1)); 3284 3285 // Set the high integerPart especially setting all unused top bits for 3286 // internal consistency. 3287 const unsigned NumUnusedHighBits = 3288 PartCount*integerPartWidth - semantics->precision; 3289 significand[PartCount - 1] = (NumUnusedHighBits < integerPartWidth) 3290 ? (~integerPart(0) >> NumUnusedHighBits) 3291 : 0; 3292 } 3293 3294 /// Make this number the smallest magnitude denormal number in the given 3295 /// semantics. 3296 void IEEEFloat::makeSmallest(bool Negative) { 3297 // We want (in interchange format): 3298 // sign = {Negative} 3299 // exponent = 0..0 3300 // significand = 0..01 3301 category = fcNormal; 3302 sign = Negative; 3303 exponent = semantics->minExponent; 3304 APInt::tcSet(significandParts(), 1, partCount()); 3305 } 3306 3307 void IEEEFloat::makeSmallestNormalized(bool Negative) { 3308 // We want (in interchange format): 3309 // sign = {Negative} 3310 // exponent = 0..0 3311 // significand = 10..0 3312 3313 category = fcNormal; 3314 zeroSignificand(); 3315 sign = Negative; 3316 exponent = semantics->minExponent; 3317 significandParts()[partCountForBits(semantics->precision) - 1] |= 3318 (((integerPart)1) << ((semantics->precision - 1) % integerPartWidth)); 3319 } 3320 3321 IEEEFloat::IEEEFloat(const fltSemantics &Sem, const APInt &API) { 3322 initFromAPInt(&Sem, API); 3323 } 3324 3325 IEEEFloat::IEEEFloat(float f) { 3326 initFromAPInt(&semIEEEsingle, APInt::floatToBits(f)); 3327 } 3328 3329 IEEEFloat::IEEEFloat(double d) { 3330 initFromAPInt(&semIEEEdouble, APInt::doubleToBits(d)); 3331 } 3332 3333 namespace { 3334 void append(SmallVectorImpl<char> &Buffer, StringRef Str) { 3335 Buffer.append(Str.begin(), Str.end()); 3336 } 3337 3338 /// Removes data from the given significand until it is no more 3339 /// precise than is required for the desired precision. 3340 void AdjustToPrecision(APInt &significand, 3341 int &exp, unsigned FormatPrecision) { 3342 unsigned bits = significand.getActiveBits(); 3343 3344 // 196/59 is a very slight overestimate of lg_2(10). 3345 unsigned bitsRequired = (FormatPrecision * 196 + 58) / 59; 3346 3347 if (bits <= bitsRequired) return; 3348 3349 unsigned tensRemovable = (bits - bitsRequired) * 59 / 196; 3350 if (!tensRemovable) return; 3351 3352 exp += tensRemovable; 3353 3354 APInt divisor(significand.getBitWidth(), 1); 3355 APInt powten(significand.getBitWidth(), 10); 3356 while (true) { 3357 if (tensRemovable & 1) 3358 divisor *= powten; 3359 tensRemovable >>= 1; 3360 if (!tensRemovable) break; 3361 powten *= powten; 3362 } 3363 3364 significand = significand.udiv(divisor); 3365 3366 // Truncate the significand down to its active bit count. 3367 significand = significand.trunc(significand.getActiveBits()); 3368 } 3369 3370 3371 void AdjustToPrecision(SmallVectorImpl<char> &buffer, 3372 int &exp, unsigned FormatPrecision) { 3373 unsigned N = buffer.size(); 3374 if (N <= FormatPrecision) return; 3375 3376 // The most significant figures are the last ones in the buffer. 3377 unsigned FirstSignificant = N - FormatPrecision; 3378 3379 // Round. 3380 // FIXME: this probably shouldn't use 'round half up'. 3381 3382 // Rounding down is just a truncation, except we also want to drop 3383 // trailing zeros from the new result. 3384 if (buffer[FirstSignificant - 1] < '5') { 3385 while (FirstSignificant < N && buffer[FirstSignificant] == '0') 3386 FirstSignificant++; 3387 3388 exp += FirstSignificant; 3389 buffer.erase(&buffer[0], &buffer[FirstSignificant]); 3390 return; 3391 } 3392 3393 // Rounding up requires a decimal add-with-carry. If we continue 3394 // the carry, the newly-introduced zeros will just be truncated. 3395 for (unsigned I = FirstSignificant; I != N; ++I) { 3396 if (buffer[I] == '9') { 3397 FirstSignificant++; 3398 } else { 3399 buffer[I]++; 3400 break; 3401 } 3402 } 3403 3404 // If we carried through, we have exactly one digit of precision. 3405 if (FirstSignificant == N) { 3406 exp += FirstSignificant; 3407 buffer.clear(); 3408 buffer.push_back('1'); 3409 return; 3410 } 3411 3412 exp += FirstSignificant; 3413 buffer.erase(&buffer[0], &buffer[FirstSignificant]); 3414 } 3415 } 3416 3417 void IEEEFloat::toString(SmallVectorImpl<char> &Str, unsigned FormatPrecision, 3418 unsigned FormatMaxPadding) const { 3419 switch (category) { 3420 case fcInfinity: 3421 if (isNegative()) 3422 return append(Str, "-Inf"); 3423 else 3424 return append(Str, "+Inf"); 3425 3426 case fcNaN: return append(Str, "NaN"); 3427 3428 case fcZero: 3429 if (isNegative()) 3430 Str.push_back('-'); 3431 3432 if (!FormatMaxPadding) 3433 append(Str, "0.0E+0"); 3434 else 3435 Str.push_back('0'); 3436 return; 3437 3438 case fcNormal: 3439 break; 3440 } 3441 3442 if (isNegative()) 3443 Str.push_back('-'); 3444 3445 // Decompose the number into an APInt and an exponent. 3446 int exp = exponent - ((int) semantics->precision - 1); 3447 APInt significand(semantics->precision, 3448 makeArrayRef(significandParts(), 3449 partCountForBits(semantics->precision))); 3450 3451 // Set FormatPrecision if zero. We want to do this before we 3452 // truncate trailing zeros, as those are part of the precision. 3453 if (!FormatPrecision) { 3454 // We use enough digits so the number can be round-tripped back to an 3455 // APFloat. The formula comes from "How to Print Floating-Point Numbers 3456 // Accurately" by Steele and White. 3457 // FIXME: Using a formula based purely on the precision is conservative; 3458 // we can print fewer digits depending on the actual value being printed. 3459 3460 // FormatPrecision = 2 + floor(significandBits / lg_2(10)) 3461 FormatPrecision = 2 + semantics->precision * 59 / 196; 3462 } 3463 3464 // Ignore trailing binary zeros. 3465 int trailingZeros = significand.countTrailingZeros(); 3466 exp += trailingZeros; 3467 significand = significand.lshr(trailingZeros); 3468 3469 // Change the exponent from 2^e to 10^e. 3470 if (exp == 0) { 3471 // Nothing to do. 3472 } else if (exp > 0) { 3473 // Just shift left. 3474 significand = significand.zext(semantics->precision + exp); 3475 significand <<= exp; 3476 exp = 0; 3477 } else { /* exp < 0 */ 3478 int texp = -exp; 3479 3480 // We transform this using the identity: 3481 // (N)(2^-e) == (N)(5^e)(10^-e) 3482 // This means we have to multiply N (the significand) by 5^e. 3483 // To avoid overflow, we have to operate on numbers large 3484 // enough to store N * 5^e: 3485 // log2(N * 5^e) == log2(N) + e * log2(5) 3486 // <= semantics->precision + e * 137 / 59 3487 // (log_2(5) ~ 2.321928 < 2.322034 ~ 137/59) 3488 3489 unsigned precision = semantics->precision + (137 * texp + 136) / 59; 3490 3491 // Multiply significand by 5^e. 3492 // N * 5^0101 == N * 5^(1*1) * 5^(0*2) * 5^(1*4) * 5^(0*8) 3493 significand = significand.zext(precision); 3494 APInt five_to_the_i(precision, 5); 3495 while (true) { 3496 if (texp & 1) significand *= five_to_the_i; 3497 3498 texp >>= 1; 3499 if (!texp) break; 3500 five_to_the_i *= five_to_the_i; 3501 } 3502 } 3503 3504 AdjustToPrecision(significand, exp, FormatPrecision); 3505 3506 SmallVector<char, 256> buffer; 3507 3508 // Fill the buffer. 3509 unsigned precision = significand.getBitWidth(); 3510 APInt ten(precision, 10); 3511 APInt digit(precision, 0); 3512 3513 bool inTrail = true; 3514 while (significand != 0) { 3515 // digit <- significand % 10 3516 // significand <- significand / 10 3517 APInt::udivrem(significand, ten, significand, digit); 3518 3519 unsigned d = digit.getZExtValue(); 3520 3521 // Drop trailing zeros. 3522 if (inTrail && !d) exp++; 3523 else { 3524 buffer.push_back((char) ('0' + d)); 3525 inTrail = false; 3526 } 3527 } 3528 3529 assert(!buffer.empty() && "no characters in buffer!"); 3530 3531 // Drop down to FormatPrecision. 3532 // TODO: don't do more precise calculations above than are required. 3533 AdjustToPrecision(buffer, exp, FormatPrecision); 3534 3535 unsigned NDigits = buffer.size(); 3536 3537 // Check whether we should use scientific notation. 3538 bool FormatScientific; 3539 if (!FormatMaxPadding) 3540 FormatScientific = true; 3541 else { 3542 if (exp >= 0) { 3543 // 765e3 --> 765000 3544 // ^^^ 3545 // But we shouldn't make the number look more precise than it is. 3546 FormatScientific = ((unsigned) exp > FormatMaxPadding || 3547 NDigits + (unsigned) exp > FormatPrecision); 3548 } else { 3549 // Power of the most significant digit. 3550 int MSD = exp + (int) (NDigits - 1); 3551 if (MSD >= 0) { 3552 // 765e-2 == 7.65 3553 FormatScientific = false; 3554 } else { 3555 // 765e-5 == 0.00765 3556 // ^ ^^ 3557 FormatScientific = ((unsigned) -MSD) > FormatMaxPadding; 3558 } 3559 } 3560 } 3561 3562 // Scientific formatting is pretty straightforward. 3563 if (FormatScientific) { 3564 exp += (NDigits - 1); 3565 3566 Str.push_back(buffer[NDigits-1]); 3567 Str.push_back('.'); 3568 if (NDigits == 1) 3569 Str.push_back('0'); 3570 else 3571 for (unsigned I = 1; I != NDigits; ++I) 3572 Str.push_back(buffer[NDigits-1-I]); 3573 Str.push_back('E'); 3574 3575 Str.push_back(exp >= 0 ? '+' : '-'); 3576 if (exp < 0) exp = -exp; 3577 SmallVector<char, 6> expbuf; 3578 do { 3579 expbuf.push_back((char) ('0' + (exp % 10))); 3580 exp /= 10; 3581 } while (exp); 3582 for (unsigned I = 0, E = expbuf.size(); I != E; ++I) 3583 Str.push_back(expbuf[E-1-I]); 3584 return; 3585 } 3586 3587 // Non-scientific, positive exponents. 3588 if (exp >= 0) { 3589 for (unsigned I = 0; I != NDigits; ++I) 3590 Str.push_back(buffer[NDigits-1-I]); 3591 for (unsigned I = 0; I != (unsigned) exp; ++I) 3592 Str.push_back('0'); 3593 return; 3594 } 3595 3596 // Non-scientific, negative exponents. 3597 3598 // The number of digits to the left of the decimal point. 3599 int NWholeDigits = exp + (int) NDigits; 3600 3601 unsigned I = 0; 3602 if (NWholeDigits > 0) { 3603 for (; I != (unsigned) NWholeDigits; ++I) 3604 Str.push_back(buffer[NDigits-I-1]); 3605 Str.push_back('.'); 3606 } else { 3607 unsigned NZeros = 1 + (unsigned) -NWholeDigits; 3608 3609 Str.push_back('0'); 3610 Str.push_back('.'); 3611 for (unsigned Z = 1; Z != NZeros; ++Z) 3612 Str.push_back('0'); 3613 } 3614 3615 for (; I != NDigits; ++I) 3616 Str.push_back(buffer[NDigits-I-1]); 3617 } 3618 3619 bool IEEEFloat::getExactInverse(IEEEFloat *inv) const { 3620 // Special floats and denormals have no exact inverse. 3621 if (!isFiniteNonZero()) 3622 return false; 3623 3624 // Check that the number is a power of two by making sure that only the 3625 // integer bit is set in the significand. 3626 if (significandLSB() != semantics->precision - 1) 3627 return false; 3628 3629 // Get the inverse. 3630 IEEEFloat reciprocal(*semantics, 1ULL); 3631 if (reciprocal.divide(*this, rmNearestTiesToEven) != opOK) 3632 return false; 3633 3634 // Avoid multiplication with a denormal, it is not safe on all platforms and 3635 // may be slower than a normal division. 3636 if (reciprocal.isDenormal()) 3637 return false; 3638 3639 assert(reciprocal.isFiniteNonZero() && 3640 reciprocal.significandLSB() == reciprocal.semantics->precision - 1); 3641 3642 if (inv) 3643 *inv = reciprocal; 3644 3645 return true; 3646 } 3647 3648 bool IEEEFloat::isSignaling() const { 3649 if (!isNaN()) 3650 return false; 3651 3652 // IEEE-754R 2008 6.2.1: A signaling NaN bit string should be encoded with the 3653 // first bit of the trailing significand being 0. 3654 return !APInt::tcExtractBit(significandParts(), semantics->precision - 2); 3655 } 3656 3657 /// IEEE-754R 2008 5.3.1: nextUp/nextDown. 3658 /// 3659 /// *NOTE* since nextDown(x) = -nextUp(-x), we only implement nextUp with 3660 /// appropriate sign switching before/after the computation. 3661 IEEEFloat::opStatus IEEEFloat::next(bool nextDown) { 3662 // If we are performing nextDown, swap sign so we have -x. 3663 if (nextDown) 3664 changeSign(); 3665 3666 // Compute nextUp(x) 3667 opStatus result = opOK; 3668 3669 // Handle each float category separately. 3670 switch (category) { 3671 case fcInfinity: 3672 // nextUp(+inf) = +inf 3673 if (!isNegative()) 3674 break; 3675 // nextUp(-inf) = -getLargest() 3676 makeLargest(true); 3677 break; 3678 case fcNaN: 3679 // IEEE-754R 2008 6.2 Par 2: nextUp(sNaN) = qNaN. Set Invalid flag. 3680 // IEEE-754R 2008 6.2: nextUp(qNaN) = qNaN. Must be identity so we do not 3681 // change the payload. 3682 if (isSignaling()) { 3683 result = opInvalidOp; 3684 // For consistency, propagate the sign of the sNaN to the qNaN. 3685 makeNaN(false, isNegative(), nullptr); 3686 } 3687 break; 3688 case fcZero: 3689 // nextUp(pm 0) = +getSmallest() 3690 makeSmallest(false); 3691 break; 3692 case fcNormal: 3693 // nextUp(-getSmallest()) = -0 3694 if (isSmallest() && isNegative()) { 3695 APInt::tcSet(significandParts(), 0, partCount()); 3696 category = fcZero; 3697 exponent = 0; 3698 break; 3699 } 3700 3701 // nextUp(getLargest()) == INFINITY 3702 if (isLargest() && !isNegative()) { 3703 APInt::tcSet(significandParts(), 0, partCount()); 3704 category = fcInfinity; 3705 exponent = semantics->maxExponent + 1; 3706 break; 3707 } 3708 3709 // nextUp(normal) == normal + inc. 3710 if (isNegative()) { 3711 // If we are negative, we need to decrement the significand. 3712 3713 // We only cross a binade boundary that requires adjusting the exponent 3714 // if: 3715 // 1. exponent != semantics->minExponent. This implies we are not in the 3716 // smallest binade or are dealing with denormals. 3717 // 2. Our significand excluding the integral bit is all zeros. 3718 bool WillCrossBinadeBoundary = 3719 exponent != semantics->minExponent && isSignificandAllZeros(); 3720 3721 // Decrement the significand. 3722 // 3723 // We always do this since: 3724 // 1. If we are dealing with a non-binade decrement, by definition we 3725 // just decrement the significand. 3726 // 2. If we are dealing with a normal -> normal binade decrement, since 3727 // we have an explicit integral bit the fact that all bits but the 3728 // integral bit are zero implies that subtracting one will yield a 3729 // significand with 0 integral bit and 1 in all other spots. Thus we 3730 // must just adjust the exponent and set the integral bit to 1. 3731 // 3. If we are dealing with a normal -> denormal binade decrement, 3732 // since we set the integral bit to 0 when we represent denormals, we 3733 // just decrement the significand. 3734 integerPart *Parts = significandParts(); 3735 APInt::tcDecrement(Parts, partCount()); 3736 3737 if (WillCrossBinadeBoundary) { 3738 // Our result is a normal number. Do the following: 3739 // 1. Set the integral bit to 1. 3740 // 2. Decrement the exponent. 3741 APInt::tcSetBit(Parts, semantics->precision - 1); 3742 exponent--; 3743 } 3744 } else { 3745 // If we are positive, we need to increment the significand. 3746 3747 // We only cross a binade boundary that requires adjusting the exponent if 3748 // the input is not a denormal and all of said input's significand bits 3749 // are set. If all of said conditions are true: clear the significand, set 3750 // the integral bit to 1, and increment the exponent. If we have a 3751 // denormal always increment since moving denormals and the numbers in the 3752 // smallest normal binade have the same exponent in our representation. 3753 bool WillCrossBinadeBoundary = !isDenormal() && isSignificandAllOnes(); 3754 3755 if (WillCrossBinadeBoundary) { 3756 integerPart *Parts = significandParts(); 3757 APInt::tcSet(Parts, 0, partCount()); 3758 APInt::tcSetBit(Parts, semantics->precision - 1); 3759 assert(exponent != semantics->maxExponent && 3760 "We can not increment an exponent beyond the maxExponent allowed" 3761 " by the given floating point semantics."); 3762 exponent++; 3763 } else { 3764 incrementSignificand(); 3765 } 3766 } 3767 break; 3768 } 3769 3770 // If we are performing nextDown, swap sign so we have -nextUp(-x) 3771 if (nextDown) 3772 changeSign(); 3773 3774 return result; 3775 } 3776 3777 void IEEEFloat::makeInf(bool Negative) { 3778 category = fcInfinity; 3779 sign = Negative; 3780 exponent = semantics->maxExponent + 1; 3781 APInt::tcSet(significandParts(), 0, partCount()); 3782 } 3783 3784 void IEEEFloat::makeZero(bool Negative) { 3785 category = fcZero; 3786 sign = Negative; 3787 exponent = semantics->minExponent-1; 3788 APInt::tcSet(significandParts(), 0, partCount()); 3789 } 3790 3791 void IEEEFloat::makeQuiet() { 3792 assert(isNaN()); 3793 APInt::tcSetBit(significandParts(), semantics->precision - 2); 3794 } 3795 3796 int ilogb(const IEEEFloat &Arg) { 3797 if (Arg.isNaN()) 3798 return IEEEFloat::IEK_NaN; 3799 if (Arg.isZero()) 3800 return IEEEFloat::IEK_Zero; 3801 if (Arg.isInfinity()) 3802 return IEEEFloat::IEK_Inf; 3803 if (!Arg.isDenormal()) 3804 return Arg.exponent; 3805 3806 IEEEFloat Normalized(Arg); 3807 int SignificandBits = Arg.getSemantics().precision - 1; 3808 3809 Normalized.exponent += SignificandBits; 3810 Normalized.normalize(IEEEFloat::rmNearestTiesToEven, lfExactlyZero); 3811 return Normalized.exponent - SignificandBits; 3812 } 3813 3814 IEEEFloat scalbn(IEEEFloat X, int Exp, IEEEFloat::roundingMode RoundingMode) { 3815 auto MaxExp = X.getSemantics().maxExponent; 3816 auto MinExp = X.getSemantics().minExponent; 3817 3818 // If Exp is wildly out-of-scale, simply adding it to X.exponent will 3819 // overflow; clamp it to a safe range before adding, but ensure that the range 3820 // is large enough that the clamp does not change the result. The range we 3821 // need to support is the difference between the largest possible exponent and 3822 // the normalized exponent of half the smallest denormal. 3823 3824 int SignificandBits = X.getSemantics().precision - 1; 3825 int MaxIncrement = MaxExp - (MinExp - SignificandBits) + 1; 3826 3827 // Clamp to one past the range ends to let normalize handle overlflow. 3828 X.exponent += std::min(std::max(Exp, -MaxIncrement - 1), MaxIncrement); 3829 X.normalize(RoundingMode, lfExactlyZero); 3830 if (X.isNaN()) 3831 X.makeQuiet(); 3832 return X; 3833 } 3834 3835 IEEEFloat frexp(const IEEEFloat &Val, int &Exp, IEEEFloat::roundingMode RM) { 3836 Exp = ilogb(Val); 3837 3838 // Quiet signalling nans. 3839 if (Exp == IEEEFloat::IEK_NaN) { 3840 IEEEFloat Quiet(Val); 3841 Quiet.makeQuiet(); 3842 return Quiet; 3843 } 3844 3845 if (Exp == IEEEFloat::IEK_Inf) 3846 return Val; 3847 3848 // 1 is added because frexp is defined to return a normalized fraction in 3849 // +/-[0.5, 1.0), rather than the usual +/-[1.0, 2.0). 3850 Exp = Exp == IEEEFloat::IEK_Zero ? 0 : Exp + 1; 3851 return scalbn(Val, -Exp, RM); 3852 } 3853 3854 DoubleAPFloat::DoubleAPFloat(const fltSemantics &S) 3855 : Semantics(&S), Floats(new APFloat[2]{APFloat(semPPCDoubleDoubleImpl), 3856 APFloat(semIEEEdouble)}) { 3857 assert(Semantics == &semPPCDoubleDouble); 3858 } 3859 3860 DoubleAPFloat::DoubleAPFloat(const fltSemantics &S, uninitializedTag) 3861 : Semantics(&S), 3862 Floats(new APFloat[2]{APFloat(semPPCDoubleDoubleImpl, uninitialized), 3863 APFloat(semIEEEdouble, uninitialized)}) { 3864 assert(Semantics == &semPPCDoubleDouble); 3865 } 3866 3867 DoubleAPFloat::DoubleAPFloat(const fltSemantics &S, integerPart I) 3868 : Semantics(&S), Floats(new APFloat[2]{APFloat(semPPCDoubleDoubleImpl, I), 3869 APFloat(semIEEEdouble)}) { 3870 assert(Semantics == &semPPCDoubleDouble); 3871 } 3872 3873 DoubleAPFloat::DoubleAPFloat(const fltSemantics &S, const APInt &I) 3874 : Semantics(&S), Floats(new APFloat[2]{ 3875 APFloat(semPPCDoubleDoubleImpl, I), 3876 APFloat(semIEEEdouble, APInt(64, I.getRawData()[1]))}) { 3877 assert(Semantics == &semPPCDoubleDouble); 3878 } 3879 3880 DoubleAPFloat::DoubleAPFloat(const fltSemantics &S, APFloat &&First, 3881 APFloat &&Second) 3882 : Semantics(&S), 3883 Floats(new APFloat[2]{std::move(First), std::move(Second)}) { 3884 assert(Semantics == &semPPCDoubleDouble); 3885 // TODO Check for First == &IEEEdouble once the transition is done. 3886 assert(&Floats[0].getSemantics() == &semPPCDoubleDoubleImpl || 3887 &Floats[0].getSemantics() == &semIEEEdouble); 3888 assert(&Floats[1].getSemantics() == &semIEEEdouble); 3889 } 3890 3891 DoubleAPFloat::DoubleAPFloat(const DoubleAPFloat &RHS) 3892 : Semantics(RHS.Semantics), 3893 Floats(RHS.Floats ? new APFloat[2]{APFloat(RHS.Floats[0]), 3894 APFloat(RHS.Floats[1])} 3895 : nullptr) { 3896 assert(Semantics == &semPPCDoubleDouble); 3897 } 3898 3899 DoubleAPFloat::DoubleAPFloat(DoubleAPFloat &&RHS) 3900 : Semantics(RHS.Semantics), Floats(std::move(RHS.Floats)) { 3901 RHS.Semantics = &semBogus; 3902 assert(Semantics == &semPPCDoubleDouble); 3903 } 3904 3905 DoubleAPFloat &DoubleAPFloat::operator=(const DoubleAPFloat &RHS) { 3906 if (Semantics == RHS.Semantics && RHS.Floats) { 3907 Floats[0] = RHS.Floats[0]; 3908 Floats[1] = RHS.Floats[1]; 3909 } else if (this != &RHS) { 3910 this->~DoubleAPFloat(); 3911 new (this) DoubleAPFloat(RHS); 3912 } 3913 return *this; 3914 } 3915 3916 // "Software for Doubled-Precision Floating-Point Computations", 3917 // by Seppo Linnainmaa, ACM TOMS vol 7 no 3, September 1981, pages 272-283. 3918 APFloat::opStatus DoubleAPFloat::addImpl(const APFloat &a, const APFloat &aa, 3919 const APFloat &c, const APFloat &cc, 3920 roundingMode RM) { 3921 int Status = opOK; 3922 APFloat z = a; 3923 Status |= z.add(c, RM); 3924 if (!z.isFinite()) { 3925 if (!z.isInfinity()) { 3926 Floats[0] = std::move(z); 3927 Floats[1].makeZero(false); 3928 return (opStatus)Status; 3929 } 3930 Status = opOK; 3931 auto AComparedToC = a.compareAbsoluteValue(c); 3932 z = cc; 3933 Status |= z.add(aa, RM); 3934 if (AComparedToC == APFloat::cmpGreaterThan) { 3935 // z = cc + aa + c + a; 3936 Status |= z.add(c, RM); 3937 Status |= z.add(a, RM); 3938 } else { 3939 // z = cc + aa + a + c; 3940 Status |= z.add(a, RM); 3941 Status |= z.add(c, RM); 3942 } 3943 if (!z.isFinite()) { 3944 Floats[0] = std::move(z); 3945 Floats[1].makeZero(false); 3946 return (opStatus)Status; 3947 } 3948 Floats[0] = z; 3949 APFloat zz = aa; 3950 Status |= zz.add(cc, RM); 3951 if (AComparedToC == APFloat::cmpGreaterThan) { 3952 // Floats[1] = a - z + c + zz; 3953 Floats[1] = a; 3954 Status |= Floats[1].subtract(z, RM); 3955 Status |= Floats[1].add(c, RM); 3956 Status |= Floats[1].add(zz, RM); 3957 } else { 3958 // Floats[1] = c - z + a + zz; 3959 Floats[1] = c; 3960 Status |= Floats[1].subtract(z, RM); 3961 Status |= Floats[1].add(a, RM); 3962 Status |= Floats[1].add(zz, RM); 3963 } 3964 } else { 3965 // q = a - z; 3966 APFloat q = a; 3967 Status |= q.subtract(z, RM); 3968 3969 // zz = q + c + (a - (q + z)) + aa + cc; 3970 // Compute a - (q + z) as -((q + z) - a) to avoid temporary copies. 3971 auto zz = q; 3972 Status |= zz.add(c, RM); 3973 Status |= q.add(z, RM); 3974 Status |= q.subtract(a, RM); 3975 q.changeSign(); 3976 Status |= zz.add(q, RM); 3977 Status |= zz.add(aa, RM); 3978 Status |= zz.add(cc, RM); 3979 if (zz.isZero() && !zz.isNegative()) { 3980 Floats[0] = std::move(z); 3981 Floats[1].makeZero(false); 3982 return opOK; 3983 } 3984 Floats[0] = z; 3985 Status |= Floats[0].add(zz, RM); 3986 if (!Floats[0].isFinite()) { 3987 Floats[1].makeZero(false); 3988 return (opStatus)Status; 3989 } 3990 Floats[1] = std::move(z); 3991 Status |= Floats[1].subtract(Floats[0], RM); 3992 Status |= Floats[1].add(zz, RM); 3993 } 3994 return (opStatus)Status; 3995 } 3996 3997 APFloat::opStatus DoubleAPFloat::addWithSpecial(const DoubleAPFloat &LHS, 3998 const DoubleAPFloat &RHS, 3999 DoubleAPFloat &Out, 4000 roundingMode RM) { 4001 if (LHS.getCategory() == fcNaN) { 4002 Out = LHS; 4003 return opOK; 4004 } 4005 if (RHS.getCategory() == fcNaN) { 4006 Out = RHS; 4007 return opOK; 4008 } 4009 if (LHS.getCategory() == fcZero) { 4010 Out = RHS; 4011 return opOK; 4012 } 4013 if (RHS.getCategory() == fcZero) { 4014 Out = LHS; 4015 return opOK; 4016 } 4017 if (LHS.getCategory() == fcInfinity && RHS.getCategory() == fcInfinity && 4018 LHS.isNegative() != RHS.isNegative()) { 4019 Out.makeNaN(false, Out.isNegative(), nullptr); 4020 return opInvalidOp; 4021 } 4022 if (LHS.getCategory() == fcInfinity) { 4023 Out = LHS; 4024 return opOK; 4025 } 4026 if (RHS.getCategory() == fcInfinity) { 4027 Out = RHS; 4028 return opOK; 4029 } 4030 assert(LHS.getCategory() == fcNormal && RHS.getCategory() == fcNormal); 4031 4032 // These conversions will go away once PPCDoubleDoubleImpl goes away. 4033 // (PPCDoubleDoubleImpl, IEEEDouble) -> (IEEEDouble, IEEEDouble) 4034 APFloat A(semIEEEdouble, 4035 APInt(64, LHS.Floats[0].bitcastToAPInt().getRawData()[0])), 4036 AA(LHS.Floats[1]), 4037 C(semIEEEdouble, APInt(64, RHS.Floats[0].bitcastToAPInt().getRawData()[0])), 4038 CC(RHS.Floats[1]); 4039 assert(&AA.getSemantics() == &semIEEEdouble); 4040 assert(&CC.getSemantics() == &semIEEEdouble); 4041 Out.Floats[0] = APFloat(semIEEEdouble); 4042 assert(&Out.Floats[1].getSemantics() == &semIEEEdouble); 4043 4044 auto Ret = Out.addImpl(A, AA, C, CC, RM); 4045 4046 // (IEEEDouble, IEEEDouble) -> (PPCDoubleDoubleImpl, IEEEDouble) 4047 uint64_t Buffer[] = {Out.Floats[0].bitcastToAPInt().getRawData()[0], 4048 Out.Floats[1].bitcastToAPInt().getRawData()[0]}; 4049 Out.Floats[0] = APFloat(semPPCDoubleDoubleImpl, APInt(128, 2, Buffer)); 4050 return Ret; 4051 } 4052 4053 APFloat::opStatus DoubleAPFloat::add(const DoubleAPFloat &RHS, 4054 roundingMode RM) { 4055 return addWithSpecial(*this, RHS, *this, RM); 4056 } 4057 4058 APFloat::opStatus DoubleAPFloat::subtract(const DoubleAPFloat &RHS, 4059 roundingMode RM) { 4060 changeSign(); 4061 auto Ret = add(RHS, RM); 4062 changeSign(); 4063 return Ret; 4064 } 4065 4066 void DoubleAPFloat::changeSign() { 4067 Floats[0].changeSign(); 4068 Floats[1].changeSign(); 4069 } 4070 4071 APFloat::cmpResult 4072 DoubleAPFloat::compareAbsoluteValue(const DoubleAPFloat &RHS) const { 4073 auto Result = Floats[0].compareAbsoluteValue(RHS.Floats[0]); 4074 if (Result != cmpEqual) 4075 return Result; 4076 Result = Floats[1].compareAbsoluteValue(RHS.Floats[1]); 4077 if (Result == cmpLessThan || Result == cmpGreaterThan) { 4078 auto Against = Floats[0].isNegative() ^ Floats[1].isNegative(); 4079 auto RHSAgainst = RHS.Floats[0].isNegative() ^ RHS.Floats[1].isNegative(); 4080 if (Against && !RHSAgainst) 4081 return cmpLessThan; 4082 if (!Against && RHSAgainst) 4083 return cmpGreaterThan; 4084 if (!Against && !RHSAgainst) 4085 return Result; 4086 if (Against && RHSAgainst) 4087 return (cmpResult)(cmpLessThan + cmpGreaterThan - Result); 4088 } 4089 return Result; 4090 } 4091 4092 APFloat::fltCategory DoubleAPFloat::getCategory() const { 4093 return Floats[0].getCategory(); 4094 } 4095 4096 bool DoubleAPFloat::isNegative() const { return Floats[0].isNegative(); } 4097 4098 void DoubleAPFloat::makeInf(bool Neg) { 4099 Floats[0].makeInf(Neg); 4100 Floats[1].makeZero(false); 4101 } 4102 4103 void DoubleAPFloat::makeNaN(bool SNaN, bool Neg, const APInt *fill) { 4104 Floats[0].makeNaN(SNaN, Neg, fill); 4105 Floats[1].makeZero(false); 4106 } 4107 4108 } // End detail namespace 4109 4110 APFloat::Storage::Storage(IEEEFloat F, const fltSemantics &Semantics) { 4111 if (usesLayout<IEEEFloat>(Semantics)) { 4112 new (&IEEE) IEEEFloat(std::move(F)); 4113 } else if (usesLayout<DoubleAPFloat>(Semantics)) { 4114 new (&Double) 4115 DoubleAPFloat(Semantics, APFloat(std::move(F), F.getSemantics()), 4116 APFloat(semIEEEdouble)); 4117 } else { 4118 llvm_unreachable("Unexpected semantics"); 4119 } 4120 } 4121 4122 APFloat::opStatus APFloat::convertFromString(StringRef Str, roundingMode RM) { 4123 return getIEEE().convertFromString(Str, RM); 4124 } 4125 4126 hash_code hash_value(const APFloat &Arg) { return hash_value(Arg.getIEEE()); } 4127 4128 APFloat::APFloat(const fltSemantics &Semantics, StringRef S) 4129 : APFloat(Semantics) { 4130 convertFromString(S, rmNearestTiesToEven); 4131 } 4132 4133 APFloat::opStatus APFloat::convert(const fltSemantics &ToSemantics, 4134 roundingMode RM, bool *losesInfo) { 4135 if (&getSemantics() == &ToSemantics) 4136 return opOK; 4137 if (usesLayout<IEEEFloat>(getSemantics()) && 4138 usesLayout<IEEEFloat>(ToSemantics)) { 4139 return U.IEEE.convert(ToSemantics, RM, losesInfo); 4140 } else if (usesLayout<IEEEFloat>(getSemantics()) && 4141 usesLayout<DoubleAPFloat>(ToSemantics)) { 4142 assert(&ToSemantics == &semPPCDoubleDouble); 4143 auto Ret = U.IEEE.convert(semPPCDoubleDoubleImpl, RM, losesInfo); 4144 *this = APFloat( 4145 DoubleAPFloat(semPPCDoubleDouble, std::move(*this), APFloat(semIEEEdouble)), 4146 ToSemantics); 4147 return Ret; 4148 } else if (usesLayout<DoubleAPFloat>(getSemantics()) && 4149 usesLayout<IEEEFloat>(ToSemantics)) { 4150 auto Ret = getIEEE().convert(ToSemantics, RM, losesInfo); 4151 *this = APFloat(std::move(getIEEE()), ToSemantics); 4152 return Ret; 4153 } else { 4154 llvm_unreachable("Unexpected semantics"); 4155 } 4156 } 4157 4158 APFloat APFloat::getAllOnesValue(unsigned BitWidth, bool isIEEE) { 4159 if (isIEEE) { 4160 switch (BitWidth) { 4161 case 16: 4162 return APFloat(semIEEEhalf, APInt::getAllOnesValue(BitWidth)); 4163 case 32: 4164 return APFloat(semIEEEsingle, APInt::getAllOnesValue(BitWidth)); 4165 case 64: 4166 return APFloat(semIEEEdouble, APInt::getAllOnesValue(BitWidth)); 4167 case 80: 4168 return APFloat(semX87DoubleExtended, APInt::getAllOnesValue(BitWidth)); 4169 case 128: 4170 return APFloat(semIEEEquad, APInt::getAllOnesValue(BitWidth)); 4171 default: 4172 llvm_unreachable("Unknown floating bit width"); 4173 } 4174 } else { 4175 assert(BitWidth == 128); 4176 return APFloat(semPPCDoubleDouble, APInt::getAllOnesValue(BitWidth)); 4177 } 4178 } 4179 4180 void APFloat::print(raw_ostream &OS) const { 4181 SmallVector<char, 16> Buffer; 4182 toString(Buffer); 4183 OS << Buffer << "\n"; 4184 } 4185 4186 void APFloat::dump() const { print(dbgs()); } 4187 4188 } // End llvm namespace 4189