contrib/llvm-project/llvm/lib/Support/APFloat.cpp

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