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