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