lib/Support/APFloat.cpp

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