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