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