lib/Support/APFloat.cpp

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