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