contrib/llvm/lib/Support/APFloat.cpp

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