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       p = end;  /* outwit assert below */
 157       break;
 158     }
 159     absExponent = value;
 160   }
 161
 162   assert(p == end && "Invalid exponent in exponent");
 163
 164   if (isNegative)
 165     return -(int) absExponent;
 166   else
 167     return (int) absExponent;
 168 }
 169
 170 /* This is ugly and needs cleaning up, but I don't immediately see
 171    how whilst remaining safe.  */
 172 static int
 173 totalExponent(StringRef::iterator p, StringRef::iterator end,
 174               int exponentAdjustment)
 175 {
 176   int unsignedExponent;
 177   bool negative, overflow;
 178   int exponent = 0;
 179
 180   assert(p != end && "Exponent has no digits");
 181
 182   negative = *p == '-';
 183   if (*p == '-' || *p == '+') {
 184     p++;
 185     assert(p != end && "Exponent has no digits");
 186   }
 187
 188   unsignedExponent = 0;
 189   overflow = false;
 190   for (; p != end; ++p) {
 191     unsigned int value;
 192
 193     value = decDigitValue(*p);
 194     assert(value < 10U && "Invalid character in exponent");
 195
 196     unsignedExponent = unsignedExponent * 10 + value;
 197     if (unsignedExponent > 32767)
 198       overflow = true;
 199   }
 200
 201   if (exponentAdjustment > 32767 || exponentAdjustment < -32768)
 202     overflow = true;
 203
 204   if (!overflow) {
 205     exponent = unsignedExponent;
 206     if (negative)
 207       exponent = -exponent;
 208     exponent += exponentAdjustment;
 209     if (exponent > 32767 || exponent < -32768)
 210       overflow = true;
 211   }
 212
 213   if (overflow)
 214     exponent = negative ? -32768: 32767;
 215
 216   return exponent;
 217 }
 218
 219 static StringRef::iterator
 220 skipLeadingZeroesAndAnyDot(StringRef::iterator begin, StringRef::iterator end,
 221                            StringRef::iterator *dot)
 222 {
 223   StringRef::iterator p = begin;
 224   *dot = end;
 225   while (*p == '0' && p != end)
 226     p++;
 227
 228   if (*p == '.') {
 229     *dot = p++;
 230
 231     assert(end - begin != 1 && "Significand has no digits");
 232
 233     while (*p == '0' && p != end)
 234       p++;
 235   }
 236
 237   return p;
 238 }
 239
 240 /* Given a normal decimal floating point number of the form
 241
 242      dddd.dddd[eE][+-]ddd
 243
 244    where the decimal point and exponent are optional, fill out the
 245    structure D.  Exponent is appropriate if the significand is
 246    treated as an integer, and normalizedExponent if the significand
 247    is taken to have the decimal point after a single leading
 248    non-zero digit.
 249
 250    If the value is zero, V->firstSigDigit points to a non-digit, and
 251    the return exponent is zero.
 252 */
 253 struct decimalInfo {
 254   const char *firstSigDigit;
 255   const char *lastSigDigit;
 256   int exponent;
 257   int normalizedExponent;
 258 };
 259
 260 static void
 261 interpretDecimal(StringRef::iterator begin, StringRef::iterator end,
 262                  decimalInfo *D)
 263 {
 264   StringRef::iterator dot = end;
 265   StringRef::iterator p = skipLeadingZeroesAndAnyDot (begin, end, &dot);
 266
 267   D->firstSigDigit = p;
 268   D->exponent = 0;
 269   D->normalizedExponent = 0;
 270
 271   for (; p != end; ++p) {
 272     if (*p == '.') {
 273       assert(dot == end && "String contains multiple dots");
 274       dot = p++;
 275       if (p == end)
 276         break;
 277     }
 278     if (decDigitValue(*p) >= 10U)
 279       break;
 280   }
 281
 282   if (p != end) {
 283     assert((*p == 'e' || *p == 'E') && "Invalid character in significand");
 284     assert(p != begin && "Significand has no digits");
 285     assert((dot == end || p - begin != 1) && "Significand has no digits");
 286
 287     /* p points to the first non-digit in the string */
 288     D->exponent = readExponent(p + 1, end);
 289
 290     /* Implied decimal point?  */
 291     if (dot == end)
 292       dot = p;
 293   }
 294
 295   /* If number is all zeroes accept any exponent.  */
 296   if (p != D->firstSigDigit) {
 297     /* Drop insignificant trailing zeroes.  */
 298     if (p != begin) {
 299       do
 300         do
 301           p--;
 302         while (p != begin && *p == '0');
 303       while (p != begin && *p == '.');
 304     }
 305
 306     /* Adjust the exponents for any decimal point.  */
 307     D->exponent += static_cast<exponent_t>((dot - p) - (dot > p));
 308     D->normalizedExponent = (D->exponent +
 309               static_cast<exponent_t>((p - D->firstSigDigit)
 310                                       - (dot > D->firstSigDigit && dot < p)));
 311   }
 312
 313   D->lastSigDigit = p;
 314 }
 315
 316 /* Return the trailing fraction of a hexadecimal number.
 317    DIGITVALUE is the first hex digit of the fraction, P points to
 318    the next digit.  */
 319 static lostFraction
 320 trailingHexadecimalFraction(StringRef::iterator p, StringRef::iterator end,
 321                             unsigned int digitValue)
 322 {
 323   unsigned int hexDigit;
 324
 325   /* If the first trailing digit isn't 0 or 8 we can work out the
 326      fraction immediately.  */
 327   if (digitValue > 8)
 328     return lfMoreThanHalf;
 329   else if (digitValue < 8 && digitValue > 0)
 330     return lfLessThanHalf;
 331
 332   /* Otherwise we need to find the first non-zero digit.  */
 333   while (*p == '0')
 334     p++;
 335
 336   assert(p != end && "Invalid trailing hexadecimal fraction!");
 337
 338   hexDigit = hexDigitValue(*p);
 339
 340   /* If we ran off the end it is exactly zero or one-half, otherwise
 341      a little more.  */
 342   if (hexDigit == -1U)
 343     return digitValue == 0 ? lfExactlyZero: lfExactlyHalf;
 344   else
 345     return digitValue == 0 ? lfLessThanHalf: lfMoreThanHalf;
 346 }
 347
 348 /* Return the fraction lost were a bignum truncated losing the least
 349    significant BITS bits.  */
 350 static lostFraction
 351 lostFractionThroughTruncation(const integerPart *parts,
 352                               unsigned int partCount,
 353                               unsigned int bits)
 354 {
 355   unsigned int lsb;
 356
 357   lsb = APInt::tcLSB(parts, partCount);
 358
 359   /* Note this is guaranteed true if bits == 0, or LSB == -1U.  */
 360   if (bits <= lsb)
 361     return lfExactlyZero;
 362   if (bits == lsb + 1)
 363     return lfExactlyHalf;
 364   if (bits <= partCount * integerPartWidth &&
 365       APInt::tcExtractBit(parts, bits - 1))
 366     return lfMoreThanHalf;
 367
 368   return lfLessThanHalf;
 369 }
 370
 371 /* Shift DST right BITS bits noting lost fraction.  */
 372 static lostFraction
 373 shiftRight(integerPart *dst, unsigned int parts, unsigned int bits)
 374 {
 375   lostFraction lost_fraction;
 376
 377   lost_fraction = lostFractionThroughTruncation(dst, parts, bits);
 378
 379   APInt::tcShiftRight(dst, parts, bits);
 380
 381   return lost_fraction;
 382 }
 383
 384 /* Combine the effect of two lost fractions.  */
 385 static lostFraction
 386 combineLostFractions(lostFraction moreSignificant,
 387                      lostFraction lessSignificant)
 388 {
 389   if (lessSignificant != lfExactlyZero) {
 390     if (moreSignificant == lfExactlyZero)
 391       moreSignificant = lfLessThanHalf;
 392     else if (moreSignificant == lfExactlyHalf)
 393       moreSignificant = lfMoreThanHalf;
 394   }
 395
 396   return moreSignificant;
 397 }
 398
 399 /* The error from the true value, in half-ulps, on multiplying two
 400    floating point numbers, which differ from the value they
 401    approximate by at most HUE1 and HUE2 half-ulps, is strictly less
 402    than the returned value.
 403
 404    See "How to Read Floating Point Numbers Accurately" by William D
 405    Clinger.  */
 406 static unsigned int
 407 HUerrBound(bool inexactMultiply, unsigned int HUerr1, unsigned int HUerr2)
 408 {
 409   assert(HUerr1 < 2 || HUerr2 < 2 || (HUerr1 + HUerr2 < 8));
 410
 411   if (HUerr1 + HUerr2 == 0)
 412     return inexactMultiply * 2;  /* <= inexactMultiply half-ulps.  */
 413   else
 414     return inexactMultiply + 2 * (HUerr1 + HUerr2);
 415 }
 416
 417 /* The number of ulps from the boundary (zero, or half if ISNEAREST)
 418    when the least significant BITS are truncated.  BITS cannot be
 419    zero.  */
 420 static integerPart
 421 ulpsFromBoundary(const integerPart *parts, unsigned int bits, bool isNearest)
 422 {
 423   unsigned int count, partBits;
 424   integerPart part, boundary;
 425
 426   assert(bits != 0);
 427
 428   bits--;
 429   count = bits / integerPartWidth;
 430   partBits = bits % integerPartWidth + 1;
 431
 432   part = parts[count] & (~(integerPart) 0 >> (integerPartWidth - partBits));
 433
 434   if (isNearest)
 435     boundary = (integerPart) 1 << (partBits - 1);
 436   else
 437     boundary = 0;
 438
 439   if (count == 0) {
 440     if (part - boundary <= boundary - part)
 441       return part - boundary;
 442     else
 443       return boundary - part;
 444   }
 445
 446   if (part == boundary) {
 447     while (--count)
 448       if (parts[count])
 449         return ~(integerPart) 0; /* A lot.  */
 450
 451     return parts[0];
 452   } else if (part == boundary - 1) {
 453     while (--count)
 454       if (~parts[count])
 455         return ~(integerPart) 0; /* A lot.  */
 456
 457     return -parts[0];
 458   }
 459
 460   return ~(integerPart) 0; /* A lot.  */
 461 }
 462
 463 /* Place pow(5, power) in DST, and return the number of parts used.
 464    DST must be at least one part larger than size of the answer.  */
 465 static unsigned int
 466 powerOf5(integerPart *dst, unsigned int power)
 467 {
 468   static const integerPart firstEightPowers[] = { 1, 5, 25, 125, 625, 3125,
 469                                                   15625, 78125 };
 470   integerPart pow5s[maxPowerOfFiveParts * 2 + 5];
 471   pow5s[0] = 78125 * 5;
 472
 473   unsigned int partsCount[16] = { 1 };
 474   integerPart scratch[maxPowerOfFiveParts], *p1, *p2, *pow5;
 475   unsigned int result;
 476   assert(power <= maxExponent);
 477
 478   p1 = dst;
 479   p2 = scratch;
 480
 481   *p1 = firstEightPowers[power & 7];
 482   power >>= 3;
 483
 484   result = 1;
 485   pow5 = pow5s;
 486
 487   for (unsigned int n = 0; power; power >>= 1, n++) {
 488     unsigned int pc;
 489
 490     pc = partsCount[n];
 491
 492     /* Calculate pow(5,pow(2,n+3)) if we haven't yet.  */
 493     if (pc == 0) {
 494       pc = partsCount[n - 1];
 495       APInt::tcFullMultiply(pow5, pow5 - pc, pow5 - pc, pc, pc);
 496       pc *= 2;
 497       if (pow5[pc - 1] == 0)
 498         pc--;
 499       partsCount[n] = pc;
 500     }
 501
 502     if (power & 1) {
 503       integerPart *tmp;
 504
 505       APInt::tcFullMultiply(p2, p1, pow5, result, pc);
 506       result += pc;
 507       if (p2[result - 1] == 0)
 508         result--;
 509
 510       /* Now result is in p1 with partsCount parts and p2 is scratch
 511          space.  */
 512       tmp = p1, p1 = p2, p2 = tmp;
 513     }
 514
 515     pow5 += pc;
 516   }
 517
 518   if (p1 != dst)
 519     APInt::tcAssign(dst, p1, result);
 520
 521   return result;
 522 }
 523
 524 /* Zero at the end to avoid modular arithmetic when adding one; used
 525    when rounding up during hexadecimal output.  */
 526 static const char hexDigitsLower[] = "0123456789abcdef0";
 527 static const char hexDigitsUpper[] = "0123456789ABCDEF0";
 528 static const char infinityL[] = "infinity";
 529 static const char infinityU[] = "INFINITY";
 530 static const char NaNL[] = "nan";
 531 static const char NaNU[] = "NAN";
 532
 533 /* Write out an integerPart in hexadecimal, starting with the most
 534    significant nibble.  Write out exactly COUNT hexdigits, return
 535    COUNT.  */
 536 static unsigned int
 537 partAsHex (char *dst, integerPart part, unsigned int count,
 538            const char *hexDigitChars)
 539 {
 540   unsigned int result = count;
 541
 542   assert(count != 0 && count <= integerPartWidth / 4);
 543
 544   part >>= (integerPartWidth - 4 * count);
 545   while (count--) {
 546     dst[count] = hexDigitChars[part & 0xf];
 547     part >>= 4;
 548   }
 549
 550   return result;
 551 }
 552
 553 /* Write out an unsigned decimal integer.  */
 554 static char *
 555 writeUnsignedDecimal (char *dst, unsigned int n)
 556 {
 557   char buff[40], *p;
 558
 559   p = buff;
 560   do
 561     *p++ = '0' + n % 10;
 562   while (n /= 10);
 563
 564   do
 565     *dst++ = *--p;
 566   while (p != buff);
 567
 568   return dst;
 569 }
 570
 571 /* Write out a signed decimal integer.  */
 572 static char *
 573 writeSignedDecimal (char *dst, int value)
 574 {
 575   if (value < 0) {
 576     *dst++ = '-';
 577     dst = writeUnsignedDecimal(dst, -(unsigned) value);
 578   } else
 579     dst = writeUnsignedDecimal(dst, value);
 580
 581   return dst;
 582 }
 583
 584 /* Constructors.  */
 585 void
 586 APFloat::initialize(const fltSemantics *ourSemantics)
 587 {
 588   unsigned int count;
 589
 590   semantics = ourSemantics;
 591   count = partCount();
 592   if (count > 1)
 593     significand.parts = new integerPart[count];
 594 }
 595
 596 void
 597 APFloat::freeSignificand()
 598 {
 599   if (partCount() > 1)
 600     delete [] significand.parts;
 601 }
 602
 603 void
 604 APFloat::assign(const APFloat &rhs)
 605 {
 606   assert(semantics == rhs.semantics);
 607
 608   sign = rhs.sign;
 609   category = rhs.category;
 610   exponent = rhs.exponent;
 611   sign2 = rhs.sign2;
 612   exponent2 = rhs.exponent2;
 613   if (category == fcNormal || category == fcNaN)
 614     copySignificand(rhs);
 615 }
 616
 617 void
 618 APFloat::copySignificand(const APFloat &rhs)
 619 {
 620   assert(category == fcNormal || category == fcNaN);
 621   assert(rhs.partCount() >= partCount());
 622
 623   APInt::tcAssign(significandParts(), rhs.significandParts(),
 624                   partCount());
 625 }
 626
 627 /* Make this number a NaN, with an arbitrary but deterministic value
 628    for the significand.  If double or longer, this is a signalling NaN,
 629    which may not be ideal.  If float, this is QNaN(0).  */
 630 void APFloat::makeNaN(bool SNaN, bool Negative, const APInt *fill)
 631 {
 632   category = fcNaN;
 633   sign = Negative;
 634
 635   integerPart *significand = significandParts();
 636   unsigned numParts = partCount();
 637
 638   // Set the significand bits to the fill.
 639   if (!fill || fill->getNumWords() < numParts)
 640     APInt::tcSet(significand, 0, numParts);
 641   if (fill) {
 642     APInt::tcAssign(significand, fill->getRawData(),
 643                     std::min(fill->getNumWords(), numParts));
 644
 645     // Zero out the excess bits of the significand.
 646     unsigned bitsToPreserve = semantics->precision - 1;
 647     unsigned part = bitsToPreserve / 64;
 648     bitsToPreserve %= 64;
 649     significand[part] &= ((1ULL << bitsToPreserve) - 1);
 650     for (part++; part != numParts; ++part)
 651       significand[part] = 0;
 652   }
 653
 654   unsigned QNaNBit = semantics->precision - 2;
 655
 656   if (SNaN) {
 657     // We always have to clear the QNaN bit to make it an SNaN.
 658     APInt::tcClearBit(significand, QNaNBit);
 659
 660     // If there are no bits set in the payload, we have to set
 661     // *something* to make it a NaN instead of an infinity;
 662     // conventionally, this is the next bit down from the QNaN bit.
 663     if (APInt::tcIsZero(significand, numParts))
 664       APInt::tcSetBit(significand, QNaNBit - 1);
 665   } else {
 666     // We always have to set the QNaN bit to make it a QNaN.
 667     APInt::tcSetBit(significand, QNaNBit);
 668   }
 669
 670   // For x87 extended precision, we want to make a NaN, not a
 671   // pseudo-NaN.  Maybe we should expose the ability to make
 672   // pseudo-NaNs?
 673   if (semantics == &APFloat::x87DoubleExtended)
 674     APInt::tcSetBit(significand, QNaNBit + 1);
 675 }
 676
 677 APFloat APFloat::makeNaN(const fltSemantics &Sem, bool SNaN, bool Negative,
 678                          const APInt *fill) {
 679   APFloat value(Sem, uninitialized);
 680   value.makeNaN(SNaN, Negative, fill);
 681   return value;
 682 }
 683
 684 APFloat &
 685 APFloat::operator=(const APFloat &rhs)
 686 {
 687   if (this != &rhs) {
 688     if (semantics != rhs.semantics) {
 689       freeSignificand();
 690       initialize(rhs.semantics);
 691     }
 692     assign(rhs);
 693   }
 694
 695   return *this;
 696 }
 697
 698 bool
 699 APFloat::bitwiseIsEqual(const APFloat &rhs) const {
 700   if (this == &rhs)
 701     return true;
 702   if (semantics != rhs.semantics ||
 703       category != rhs.category ||
 704       sign != rhs.sign)
 705     return false;
 706   if (semantics==(const llvm::fltSemantics*)&PPCDoubleDouble &&
 707       sign2 != rhs.sign2)
 708     return false;
 709   if (category==fcZero || category==fcInfinity)
 710     return true;
 711   else if (category==fcNormal && exponent!=rhs.exponent)
 712     return false;
 713   else if (semantics==(const llvm::fltSemantics*)&PPCDoubleDouble &&
 714            exponent2!=rhs.exponent2)
 715     return false;
 716   else {
 717     int i= partCount();
 718     const integerPart* p=significandParts();
 719     const integerPart* q=rhs.significandParts();
 720     for (; i>0; i--, p++, q++) {
 721       if (*p != *q)
 722         return false;
 723     }
 724     return true;
 725   }
 726 }
 727
 728 APFloat::APFloat(const fltSemantics &ourSemantics, integerPart value)
 729   : exponent2(0), sign2(0) {
 730   assertArithmeticOK(ourSemantics);
 731   initialize(&ourSemantics);
 732   sign = 0;
 733   zeroSignificand();
 734   exponent = ourSemantics.precision - 1;
 735   significandParts()[0] = value;
 736   normalize(rmNearestTiesToEven, lfExactlyZero);
 737 }
 738
 739 APFloat::APFloat(const fltSemantics &ourSemantics) : exponent2(0), sign2(0) {
 740   assertArithmeticOK(ourSemantics);
 741   initialize(&ourSemantics);
 742   category = fcZero;
 743   sign = false;
 744 }
 745
 746 APFloat::APFloat(const fltSemantics &ourSemantics, uninitializedTag tag)
 747   : exponent2(0), sign2(0) {
 748   assertArithmeticOK(ourSemantics);
 749   // Allocates storage if necessary but does not initialize it.
 750   initialize(&ourSemantics);
 751 }
 752
 753 APFloat::APFloat(const fltSemantics &ourSemantics,
 754                  fltCategory ourCategory, bool negative)
 755   : exponent2(0), sign2(0) {
 756   assertArithmeticOK(ourSemantics);
 757   initialize(&ourSemantics);
 758   category = ourCategory;
 759   sign = negative;
 760   if (category == fcNormal)
 761     category = fcZero;
 762   else if (ourCategory == fcNaN)
 763     makeNaN();
 764 }
 765
 766 APFloat::APFloat(const fltSemantics &ourSemantics, StringRef text)
 767   : exponent2(0), sign2(0) {
 768   assertArithmeticOK(ourSemantics);
 769   initialize(&ourSemantics);
 770   convertFromString(text, rmNearestTiesToEven);
 771 }
 772
 773 APFloat::APFloat(const APFloat &rhs) : exponent2(0), sign2(0) {
 774   initialize(rhs.semantics);
 775   assign(rhs);
 776 }
 777
 778 APFloat::~APFloat()
 779 {
 780   freeSignificand();
 781 }
 782
 783 // Profile - This method 'profiles' an APFloat for use with FoldingSet.
 784 void APFloat::Profile(FoldingSetNodeID& ID) const {
 785   ID.Add(bitcastToAPInt());
 786 }
 787
 788 unsigned int
 789 APFloat::partCount() const
 790 {
 791   return partCountForBits(semantics->precision + 1);
 792 }
 793
 794 unsigned int
 795 APFloat::semanticsPrecision(const fltSemantics &semantics)
 796 {
 797   return semantics.precision;
 798 }
 799
 800 const integerPart *
 801 APFloat::significandParts() const
 802 {
 803   return const_cast<APFloat *>(this)->significandParts();
 804 }
 805
 806 integerPart *
 807 APFloat::significandParts()
 808 {
 809   assert(category == fcNormal || category == fcNaN);
 810
 811   if (partCount() > 1)
 812     return significand.parts;
 813   else
 814     return &significand.part;
 815 }
 816
 817 void
 818 APFloat::zeroSignificand()
 819 {
 820   category = fcNormal;
 821   APInt::tcSet(significandParts(), 0, partCount());
 822 }
 823
 824 /* Increment an fcNormal floating point number's significand.  */
 825 void
 826 APFloat::incrementSignificand()
 827 {
 828   integerPart carry;
 829
 830   carry = APInt::tcIncrement(significandParts(), partCount());
 831
 832   /* Our callers should never cause us to overflow.  */
 833   assert(carry == 0);
 834 }
 835
 836 /* Add the significand of the RHS.  Returns the carry flag.  */
 837 integerPart
 838 APFloat::addSignificand(const APFloat &rhs)
 839 {
 840   integerPart *parts;
 841
 842   parts = significandParts();
 843
 844   assert(semantics == rhs.semantics);
 845   assert(exponent == rhs.exponent);
 846
 847   return APInt::tcAdd(parts, rhs.significandParts(), 0, partCount());
 848 }
 849
 850 /* Subtract the significand of the RHS with a borrow flag.  Returns
 851    the borrow flag.  */
 852 integerPart
 853 APFloat::subtractSignificand(const APFloat &rhs, integerPart borrow)
 854 {
 855   integerPart *parts;
 856
 857   parts = significandParts();
 858
 859   assert(semantics == rhs.semantics);
 860   assert(exponent == rhs.exponent);
 861
 862   return APInt::tcSubtract(parts, rhs.significandParts(), borrow,
 863                            partCount());
 864 }
 865
 866 /* Multiply the significand of the RHS.  If ADDEND is non-NULL, add it
 867    on to the full-precision result of the multiplication.  Returns the
 868    lost fraction.  */
 869 lostFraction
 870 APFloat::multiplySignificand(const APFloat &rhs, const APFloat *addend)
 871 {
 872   unsigned int omsb;        // One, not zero, based MSB.
 873   unsigned int partsCount, newPartsCount, precision;
 874   integerPart *lhsSignificand;
 875   integerPart scratch[4];
 876   integerPart *fullSignificand;
 877   lostFraction lost_fraction;
 878   bool ignored;
 879
 880   assert(semantics == rhs.semantics);
 881
 882   precision = semantics->precision;
 883   newPartsCount = partCountForBits(precision * 2);
 884
 885   if (newPartsCount > 4)
 886     fullSignificand = new integerPart[newPartsCount];
 887   else
 888     fullSignificand = scratch;
 889
 890   lhsSignificand = significandParts();
 891   partsCount = partCount();
 892
 893   APInt::tcFullMultiply(fullSignificand, lhsSignificand,
 894                         rhs.significandParts(), partsCount, partsCount);
 895
 896   lost_fraction = lfExactlyZero;
 897   omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
 898   exponent += rhs.exponent;
 899
 900   if (addend) {
 901     Significand savedSignificand = significand;
 902     const fltSemantics *savedSemantics = semantics;
 903     fltSemantics extendedSemantics;
 904     opStatus status;
 905     unsigned int extendedPrecision;
 906
 907     /* Normalize our MSB.  */
 908     extendedPrecision = precision + precision - 1;
 909     if (omsb != extendedPrecision) {
 910       APInt::tcShiftLeft(fullSignificand, newPartsCount,
 911                          extendedPrecision - omsb);
 912       exponent -= extendedPrecision - omsb;
 913     }
 914
 915     /* Create new semantics.  */
 916     extendedSemantics = *semantics;
 917     extendedSemantics.precision = extendedPrecision;
 918
 919     if (newPartsCount == 1)
 920       significand.part = fullSignificand[0];
 921     else
 922       significand.parts = fullSignificand;
 923     semantics = &extendedSemantics;
 924
 925     APFloat extendedAddend(*addend);
 926     status = extendedAddend.convert(extendedSemantics, rmTowardZero, &ignored);
 927     assert(status == opOK);
 928     lost_fraction = addOrSubtractSignificand(extendedAddend, false);
 929
 930     /* Restore our state.  */
 931     if (newPartsCount == 1)
 932       fullSignificand[0] = significand.part;
 933     significand = savedSignificand;
 934     semantics = savedSemantics;
 935
 936     omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
 937   }
 938
 939   exponent -= (precision - 1);
 940
 941   if (omsb > precision) {
 942     unsigned int bits, significantParts;
 943     lostFraction lf;
 944
 945     bits = omsb - precision;
 946     significantParts = partCountForBits(omsb);
 947     lf = shiftRight(fullSignificand, significantParts, bits);
 948     lost_fraction = combineLostFractions(lf, lost_fraction);
 949     exponent += bits;
 950   }
 951
 952   APInt::tcAssign(lhsSignificand, fullSignificand, partsCount);
 953
 954   if (newPartsCount > 4)
 955     delete [] fullSignificand;
 956
 957   return lost_fraction;
 958 }
 959
 960 /* Multiply the significands of LHS and RHS to DST.  */
 961 lostFraction
 962 APFloat::divideSignificand(const APFloat &rhs)
 963 {
 964   unsigned int bit, i, partsCount;
 965   const integerPart *rhsSignificand;
 966   integerPart *lhsSignificand, *dividend, *divisor;
 967   integerPart scratch[4];
 968   lostFraction lost_fraction;
 969
 970   assert(semantics == rhs.semantics);
 971
 972   lhsSignificand = significandParts();
 973   rhsSignificand = rhs.significandParts();
 974   partsCount = partCount();
 975
 976   if (partsCount > 2)
 977     dividend = new integerPart[partsCount * 2];
 978   else
 979     dividend = scratch;
 980
 981   divisor = dividend + partsCount;
 982
 983   /* Copy the dividend and divisor as they will be modified in-place.  */
 984   for (i = 0; i < partsCount; i++) {
 985     dividend[i] = lhsSignificand[i];
 986     divisor[i] = rhsSignificand[i];
 987     lhsSignificand[i] = 0;
 988   }
 989
 990   exponent -= rhs.exponent;
 991
 992   unsigned int precision = semantics->precision;
 993
 994   /* Normalize the divisor.  */
 995   bit = precision - APInt::tcMSB(divisor, partsCount) - 1;
 996   if (bit) {
 997     exponent += bit;
 998     APInt::tcShiftLeft(divisor, partsCount, bit);
 999   }
1000
1001   /* Normalize the dividend.  */
1002   bit = precision - APInt::tcMSB(dividend, partsCount) - 1;
1003   if (bit) {
1004     exponent -= bit;
1005     APInt::tcShiftLeft(dividend, partsCount, bit);
1006   }
1007
1008   /* Ensure the dividend >= divisor initially for the loop below.
1009      Incidentally, this means that the division loop below is
1010      guaranteed to set the integer bit to one.  */
1011   if (APInt::tcCompare(dividend, divisor, partsCount) < 0) {
1012     exponent--;
1013     APInt::tcShiftLeft(dividend, partsCount, 1);
1014     assert(APInt::tcCompare(dividend, divisor, partsCount) >= 0);
1015   }
1016
1017   /* Long division.  */
1018   for (bit = precision; bit; bit -= 1) {
1019     if (APInt::tcCompare(dividend, divisor, partsCount) >= 0) {
1020       APInt::tcSubtract(dividend, divisor, 0, partsCount);
1021       APInt::tcSetBit(lhsSignificand, bit - 1);
1022     }
1023
1024     APInt::tcShiftLeft(dividend, partsCount, 1);
1025   }
1026
1027   /* Figure out the lost fraction.  */
1028   int cmp = APInt::tcCompare(dividend, divisor, partsCount);
1029
1030   if (cmp > 0)
1031     lost_fraction = lfMoreThanHalf;
1032   else if (cmp == 0)
1033     lost_fraction = lfExactlyHalf;
1034   else if (APInt::tcIsZero(dividend, partsCount))
1035     lost_fraction = lfExactlyZero;
1036   else
1037     lost_fraction = lfLessThanHalf;
1038
1039   if (partsCount > 2)
1040     delete [] dividend;
1041
1042   return lost_fraction;
1043 }
1044
1045 unsigned int
1046 APFloat::significandMSB() const
1047 {
1048   return APInt::tcMSB(significandParts(), partCount());
1049 }
1050
1051 unsigned int
1052 APFloat::significandLSB() const
1053 {
1054   return APInt::tcLSB(significandParts(), partCount());
1055 }
1056
1057 /* Note that a zero result is NOT normalized to fcZero.  */
1058 lostFraction
1059 APFloat::shiftSignificandRight(unsigned int bits)
1060 {
1061   /* Our exponent should not overflow.  */
1062   assert((exponent_t) (exponent + bits) >= exponent);
1063
1064   exponent += bits;
1065
1066   return shiftRight(significandParts(), partCount(), bits);
1067 }
1068
1069 /* Shift the significand left BITS bits, subtract BITS from its exponent.  */
1070 void
1071 APFloat::shiftSignificandLeft(unsigned int bits)
1072 {
1073   assert(bits < semantics->precision);
1074
1075   if (bits) {
1076     unsigned int partsCount = partCount();
1077
1078     APInt::tcShiftLeft(significandParts(), partsCount, bits);
1079     exponent -= bits;
1080
1081     assert(!APInt::tcIsZero(significandParts(), partsCount));
1082   }
1083 }
1084
1085 APFloat::cmpResult
1086 APFloat::compareAbsoluteValue(const APFloat &rhs) const
1087 {
1088   int compare;
1089
1090   assert(semantics == rhs.semantics);
1091   assert(category == fcNormal);
1092   assert(rhs.category == fcNormal);
1093
1094   compare = exponent - rhs.exponent;
1095
1096   /* If exponents are equal, do an unsigned bignum comparison of the
1097      significands.  */
1098   if (compare == 0)
1099     compare = APInt::tcCompare(significandParts(), rhs.significandParts(),
1100                                partCount());
1101
1102   if (compare > 0)
1103     return cmpGreaterThan;
1104   else if (compare < 0)
1105     return cmpLessThan;
1106   else
1107     return cmpEqual;
1108 }
1109
1110 /* Handle overflow.  Sign is preserved.  We either become infinity or
1111    the largest finite number.  */
1112 APFloat::opStatus
1113 APFloat::handleOverflow(roundingMode rounding_mode)
1114 {
1115   /* Infinity?  */
1116   if (rounding_mode == rmNearestTiesToEven ||
1117       rounding_mode == rmNearestTiesToAway ||
1118       (rounding_mode == rmTowardPositive && !sign) ||
1119       (rounding_mode == rmTowardNegative && sign)) {
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(StringRef s, roundingMode rounding_mode)
2190 {
2191   lostFraction lost_fraction = lfExactlyZero;
2192   integerPart *significand;
2193   unsigned int bitPos, partsCount;
2194   StringRef::iterator dot, firstSignificantDigit;
2195
2196   zeroSignificand();
2197   exponent = 0;
2198   category = fcNormal;
2199
2200   significand = significandParts();
2201   partsCount = partCount();
2202   bitPos = partsCount * integerPartWidth;
2203
2204   /* Skip leading zeroes and any (hexa)decimal point.  */
2205   StringRef::iterator begin = s.begin();
2206   StringRef::iterator end = s.end();
2207   StringRef::iterator p = skipLeadingZeroesAndAnyDot(begin, end, &dot);
2208   firstSignificantDigit = p;
2209
2210   for (; p != end;) {
2211     integerPart hex_value;
2212
2213     if (*p == '.') {
2214       assert(dot == end && "String contains multiple dots");
2215       dot = p++;
2216       if (p == end) {
2217         break;
2218       }
2219     }
2220
2221     hex_value = hexDigitValue(*p);
2222     if (hex_value == -1U) {
2223       break;
2224     }
2225
2226     p++;
2227
2228     if (p == end) {
2229       break;
2230     } else {
2231       /* Store the number whilst 4-bit nibbles remain.  */
2232       if (bitPos) {
2233         bitPos -= 4;
2234         hex_value <<= bitPos % integerPartWidth;
2235         significand[bitPos / integerPartWidth] |= hex_value;
2236       } else {
2237         lost_fraction = trailingHexadecimalFraction(p, end, hex_value);
2238         while (p != end && hexDigitValue(*p) != -1U)
2239           p++;
2240         break;
2241       }
2242     }
2243   }
2244
2245   /* Hex floats require an exponent but not a hexadecimal point.  */
2246   assert(p != end && "Hex strings require an exponent");
2247   assert((*p == 'p' || *p == 'P') && "Invalid character in significand");
2248   assert(p != begin && "Significand has no digits");
2249   assert((dot == end || p - begin != 1) && "Significand has no digits");
2250
2251   /* Ignore the exponent if we are zero.  */
2252   if (p != firstSignificantDigit) {
2253     int expAdjustment;
2254
2255     /* Implicit hexadecimal point?  */
2256     if (dot == end)
2257       dot = p;
2258
2259     /* Calculate the exponent adjustment implicit in the number of
2260        significant digits.  */
2261     expAdjustment = static_cast<int>(dot - firstSignificantDigit);
2262     if (expAdjustment < 0)
2263       expAdjustment++;
2264     expAdjustment = expAdjustment * 4 - 1;
2265
2266     /* Adjust for writing the significand starting at the most
2267        significant nibble.  */
2268     expAdjustment += semantics->precision;
2269     expAdjustment -= partsCount * integerPartWidth;
2270
2271     /* Adjust for the given exponent.  */
2272     exponent = totalExponent(p + 1, end, expAdjustment);
2273   }
2274
2275   return normalize(rounding_mode, lost_fraction);
2276 }
2277
2278 APFloat::opStatus
2279 APFloat::roundSignificandWithExponent(const integerPart *decSigParts,
2280                                       unsigned sigPartCount, int exp,
2281                                       roundingMode rounding_mode)
2282 {
2283   unsigned int parts, pow5PartCount;
2284   fltSemantics calcSemantics = { 32767, -32767, 0, true };
2285   integerPart pow5Parts[maxPowerOfFiveParts];
2286   bool isNearest;
2287
2288   isNearest = (rounding_mode == rmNearestTiesToEven ||
2289                rounding_mode == rmNearestTiesToAway);
2290
2291   parts = partCountForBits(semantics->precision + 11);
2292
2293   /* Calculate pow(5, abs(exp)).  */
2294   pow5PartCount = powerOf5(pow5Parts, exp >= 0 ? exp: -exp);
2295
2296   for (;; parts *= 2) {
2297     opStatus sigStatus, powStatus;
2298     unsigned int excessPrecision, truncatedBits;
2299
2300     calcSemantics.precision = parts * integerPartWidth - 1;
2301     excessPrecision = calcSemantics.precision - semantics->precision;
2302     truncatedBits = excessPrecision;
2303
2304     APFloat decSig(calcSemantics, fcZero, sign);
2305     APFloat pow5(calcSemantics, fcZero, false);
2306
2307     sigStatus = decSig.convertFromUnsignedParts(decSigParts, sigPartCount,
2308                                                 rmNearestTiesToEven);
2309     powStatus = pow5.convertFromUnsignedParts(pow5Parts, pow5PartCount,
2310                                               rmNearestTiesToEven);
2311     /* Add exp, as 10^n = 5^n * 2^n.  */
2312     decSig.exponent += exp;
2313
2314     lostFraction calcLostFraction;
2315     integerPart HUerr, HUdistance;
2316     unsigned int powHUerr;
2317
2318     if (exp >= 0) {
2319       /* multiplySignificand leaves the precision-th bit set to 1.  */
2320       calcLostFraction = decSig.multiplySignificand(pow5, NULL);
2321       powHUerr = powStatus != opOK;
2322     } else {
2323       calcLostFraction = decSig.divideSignificand(pow5);
2324       /* Denormal numbers have less precision.  */
2325       if (decSig.exponent < semantics->minExponent) {
2326         excessPrecision += (semantics->minExponent - decSig.exponent);
2327         truncatedBits = excessPrecision;
2328         if (excessPrecision > calcSemantics.precision)
2329           excessPrecision = calcSemantics.precision;
2330       }
2331       /* Extra half-ulp lost in reciprocal of exponent.  */
2332       powHUerr = (powStatus == opOK && calcLostFraction == lfExactlyZero) ? 0:2;
2333     }
2334
2335     /* Both multiplySignificand and divideSignificand return the
2336        result with the integer bit set.  */
2337     assert(APInt::tcExtractBit
2338            (decSig.significandParts(), calcSemantics.precision - 1) == 1);
2339
2340     HUerr = HUerrBound(calcLostFraction != lfExactlyZero, sigStatus != opOK,
2341                        powHUerr);
2342     HUdistance = 2 * ulpsFromBoundary(decSig.significandParts(),
2343                                       excessPrecision, isNearest);
2344
2345     /* Are we guaranteed to round correctly if we truncate?  */
2346     if (HUdistance >= HUerr) {
2347       APInt::tcExtract(significandParts(), partCount(), decSig.significandParts(),
2348                        calcSemantics.precision - excessPrecision,
2349                        excessPrecision);
2350       /* Take the exponent of decSig.  If we tcExtract-ed less bits
2351          above we must adjust our exponent to compensate for the
2352          implicit right shift.  */
2353       exponent = (decSig.exponent + semantics->precision
2354                   - (calcSemantics.precision - excessPrecision));
2355       calcLostFraction = lostFractionThroughTruncation(decSig.significandParts(),
2356                                                        decSig.partCount(),
2357                                                        truncatedBits);
2358       return normalize(rounding_mode, calcLostFraction);
2359     }
2360   }
2361 }
2362
2363 APFloat::opStatus
2364 APFloat::convertFromDecimalString(StringRef str, roundingMode rounding_mode)
2365 {
2366   decimalInfo D;
2367   opStatus fs;
2368
2369   /* Scan the text.  */
2370   StringRef::iterator p = str.begin();
2371   interpretDecimal(p, str.end(), &D);
2372
2373   /* Handle the quick cases.  First the case of no significant digits,
2374      i.e. zero, and then exponents that are obviously too large or too
2375      small.  Writing L for log 10 / log 2, a number d.ddddd*10^exp
2376      definitely overflows if
2377
2378            (exp - 1) * L >= maxExponent
2379
2380      and definitely underflows to zero where
2381
2382            (exp + 1) * L <= minExponent - precision
2383
2384      With integer arithmetic the tightest bounds for L are
2385
2386            93/28 < L < 196/59            [ numerator <= 256 ]
2387            42039/12655 < L < 28738/8651  [ numerator <= 65536 ]
2388   */
2389
2390   if (decDigitValue(*D.firstSigDigit) >= 10U) {
2391     category = fcZero;
2392     fs = opOK;
2393
2394   /* Check whether the normalized exponent is high enough to overflow
2395      max during the log-rebasing in the max-exponent check below. */
2396   } else if (D.normalizedExponent - 1 > INT_MAX / 42039) {
2397     fs = handleOverflow(rounding_mode);
2398
2399   /* If it wasn't, then it also wasn't high enough to overflow max
2400      during the log-rebasing in the min-exponent check.  Check that it
2401      won't overflow min in either check, then perform the min-exponent
2402      check. */
2403   } else if (D.normalizedExponent - 1 < INT_MIN / 42039 ||
2404              (D.normalizedExponent + 1) * 28738 <=
2405                8651 * (semantics->minExponent - (int) semantics->precision)) {
2406     /* Underflow to zero and round.  */
2407     zeroSignificand();
2408     fs = normalize(rounding_mode, lfLessThanHalf);
2409
2410   /* We can finally safely perform the max-exponent check. */
2411   } else if ((D.normalizedExponent - 1) * 42039
2412              >= 12655 * semantics->maxExponent) {
2413     /* Overflow and round.  */
2414     fs = handleOverflow(rounding_mode);
2415   } else {
2416     integerPart *decSignificand;
2417     unsigned int partCount;
2418
2419     /* A tight upper bound on number of bits required to hold an
2420        N-digit decimal integer is N * 196 / 59.  Allocate enough space
2421        to hold the full significand, and an extra part required by
2422        tcMultiplyPart.  */
2423     partCount = static_cast<unsigned int>(D.lastSigDigit - D.firstSigDigit) + 1;
2424     partCount = partCountForBits(1 + 196 * partCount / 59);
2425     decSignificand = new integerPart[partCount + 1];
2426     partCount = 0;
2427
2428     /* Convert to binary efficiently - we do almost all multiplication
2429        in an integerPart.  When this would overflow do we do a single
2430        bignum multiplication, and then revert again to multiplication
2431        in an integerPart.  */
2432     do {
2433       integerPart decValue, val, multiplier;
2434
2435       val = 0;
2436       multiplier = 1;
2437
2438       do {
2439         if (*p == '.') {
2440           p++;
2441           if (p == str.end()) {
2442             break;
2443           }
2444         }
2445         decValue = decDigitValue(*p++);
2446         assert(decValue < 10U && "Invalid character in significand");
2447         multiplier *= 10;
2448         val = val * 10 + decValue;
2449         /* The maximum number that can be multiplied by ten with any
2450            digit added without overflowing an integerPart.  */
2451       } while (p <= D.lastSigDigit && multiplier <= (~ (integerPart) 0 - 9) / 10);
2452
2453       /* Multiply out the current part.  */
2454       APInt::tcMultiplyPart(decSignificand, decSignificand, multiplier, val,
2455                             partCount, partCount + 1, false);
2456
2457       /* If we used another part (likely but not guaranteed), increase
2458          the count.  */
2459       if (decSignificand[partCount])
2460         partCount++;
2461     } while (p <= D.lastSigDigit);
2462
2463     category = fcNormal;
2464     fs = roundSignificandWithExponent(decSignificand, partCount,
2465                                       D.exponent, rounding_mode);
2466
2467     delete [] decSignificand;
2468   }
2469
2470   return fs;
2471 }
2472
2473 APFloat::opStatus
2474 APFloat::convertFromString(StringRef str, roundingMode rounding_mode)
2475 {
2476   assertArithmeticOK(*semantics);
2477   assert(!str.empty() && "Invalid string length");
2478
2479   /* Handle a leading minus sign.  */
2480   StringRef::iterator p = str.begin();
2481   size_t slen = str.size();
2482   sign = *p == '-' ? 1 : 0;
2483   if (*p == '-' || *p == '+') {
2484     p++;
2485     slen--;
2486     assert(slen && "String has no digits");
2487   }
2488
2489   if (slen >= 2 && p[0] == '0' && (p[1] == 'x' || p[1] == 'X')) {
2490     assert(slen - 2 && "Invalid string");
2491     return convertFromHexadecimalString(StringRef(p + 2, slen - 2),
2492                                         rounding_mode);
2493   }
2494
2495   return convertFromDecimalString(StringRef(p, slen), rounding_mode);
2496 }
2497
2498 /* Write out a hexadecimal representation of the floating point value
2499    to DST, which must be of sufficient size, in the C99 form
2500    [-]0xh.hhhhp[+-]d.  Return the number of characters written,
2501    excluding the terminating NUL.
2502
2503    If UPPERCASE, the output is in upper case, otherwise in lower case.
2504
2505    HEXDIGITS digits appear altogether, rounding the value if
2506    necessary.  If HEXDIGITS is 0, the minimal precision to display the
2507    number precisely is used instead.  If nothing would appear after
2508    the decimal point it is suppressed.
2509
2510    The decimal exponent is always printed and has at least one digit.
2511    Zero values display an exponent of zero.  Infinities and NaNs
2512    appear as "infinity" or "nan" respectively.
2513
2514    The above rules are as specified by C99.  There is ambiguity about
2515    what the leading hexadecimal digit should be.  This implementation
2516    uses whatever is necessary so that the exponent is displayed as
2517    stored.  This implies the exponent will fall within the IEEE format
2518    range, and the leading hexadecimal digit will be 0 (for denormals),
2519    1 (normal numbers) or 2 (normal numbers rounded-away-from-zero with
2520    any other digits zero).
2521 */
2522 unsigned int
2523 APFloat::convertToHexString(char *dst, unsigned int hexDigits,
2524                             bool upperCase, roundingMode rounding_mode) const
2525 {
2526   char *p;
2527
2528   assertArithmeticOK(*semantics);
2529
2530   p = dst;
2531   if (sign)
2532     *dst++ = '-';
2533
2534   switch (category) {
2535   case fcInfinity:
2536     memcpy (dst, upperCase ? infinityU: infinityL, sizeof infinityU - 1);
2537     dst += sizeof infinityL - 1;
2538     break;
2539
2540   case fcNaN:
2541     memcpy (dst, upperCase ? NaNU: NaNL, sizeof NaNU - 1);
2542     dst += sizeof NaNU - 1;
2543     break;
2544
2545   case fcZero:
2546     *dst++ = '0';
2547     *dst++ = upperCase ? 'X': 'x';
2548     *dst++ = '0';
2549     if (hexDigits > 1) {
2550       *dst++ = '.';
2551       memset (dst, '0', hexDigits - 1);
2552       dst += hexDigits - 1;
2553     }
2554     *dst++ = upperCase ? 'P': 'p';
2555     *dst++ = '0';
2556     break;
2557
2558   case fcNormal:
2559     dst = convertNormalToHexString (dst, hexDigits, upperCase, rounding_mode);
2560     break;
2561   }
2562
2563   *dst = 0;
2564
2565   return static_cast<unsigned int>(dst - p);
2566 }
2567
2568 /* Does the hard work of outputting the correctly rounded hexadecimal
2569    form of a normal floating point number with the specified number of
2570    hexadecimal digits.  If HEXDIGITS is zero the minimum number of
2571    digits necessary to print the value precisely is output.  */
2572 char *
2573 APFloat::convertNormalToHexString(char *dst, unsigned int hexDigits,
2574                                   bool upperCase,
2575                                   roundingMode rounding_mode) const
2576 {
2577   unsigned int count, valueBits, shift, partsCount, outputDigits;
2578   const char *hexDigitChars;
2579   const integerPart *significand;
2580   char *p;
2581   bool roundUp;
2582
2583   *dst++ = '0';
2584   *dst++ = upperCase ? 'X': 'x';
2585
2586   roundUp = false;
2587   hexDigitChars = upperCase ? hexDigitsUpper: hexDigitsLower;
2588
2589   significand = significandParts();
2590   partsCount = partCount();
2591
2592   /* +3 because the first digit only uses the single integer bit, so
2593      we have 3 virtual zero most-significant-bits.  */
2594   valueBits = semantics->precision + 3;
2595   shift = integerPartWidth - valueBits % integerPartWidth;
2596
2597   /* The natural number of digits required ignoring trailing
2598      insignificant zeroes.  */
2599   outputDigits = (valueBits - significandLSB () + 3) / 4;
2600
2601   /* hexDigits of zero means use the required number for the
2602      precision.  Otherwise, see if we are truncating.  If we are,
2603      find out if we need to round away from zero.  */
2604   if (hexDigits) {
2605     if (hexDigits < outputDigits) {
2606       /* We are dropping non-zero bits, so need to check how to round.
2607          "bits" is the number of dropped bits.  */
2608       unsigned int bits;
2609       lostFraction fraction;
2610
2611       bits = valueBits - hexDigits * 4;
2612       fraction = lostFractionThroughTruncation (significand, partsCount, bits);
2613       roundUp = roundAwayFromZero(rounding_mode, fraction, bits);
2614     }
2615     outputDigits = hexDigits;
2616   }
2617
2618   /* Write the digits consecutively, and start writing in the location
2619      of the hexadecimal point.  We move the most significant digit
2620      left and add the hexadecimal point later.  */
2621   p = ++dst;
2622
2623   count = (valueBits + integerPartWidth - 1) / integerPartWidth;
2624
2625   while (outputDigits && count) {
2626     integerPart part;
2627
2628     /* Put the most significant integerPartWidth bits in "part".  */
2629     if (--count == partsCount)
2630       part = 0;  /* An imaginary higher zero part.  */
2631     else
2632       part = significand[count] << shift;
2633
2634     if (count && shift)
2635       part |= significand[count - 1] >> (integerPartWidth - shift);
2636
2637     /* Convert as much of "part" to hexdigits as we can.  */
2638     unsigned int curDigits = integerPartWidth / 4;
2639
2640     if (curDigits > outputDigits)
2641       curDigits = outputDigits;
2642     dst += partAsHex (dst, part, curDigits, hexDigitChars);
2643     outputDigits -= curDigits;
2644   }
2645
2646   if (roundUp) {
2647     char *q = dst;
2648
2649     /* Note that hexDigitChars has a trailing '0'.  */
2650     do {
2651       q--;
2652       *q = hexDigitChars[hexDigitValue (*q) + 1];
2653     } while (*q == '0');
2654     assert(q >= p);
2655   } else {
2656     /* Add trailing zeroes.  */
2657     memset (dst, '0', outputDigits);
2658     dst += outputDigits;
2659   }
2660
2661   /* Move the most significant digit to before the point, and if there
2662      is something after the decimal point add it.  This must come
2663      after rounding above.  */
2664   p[-1] = p[0];
2665   if (dst -1 == p)
2666     dst--;
2667   else
2668     p[0] = '.';
2669
2670   /* Finally output the exponent.  */
2671   *dst++ = upperCase ? 'P': 'p';
2672
2673   return writeSignedDecimal (dst, exponent);
2674 }
2675
2676 // For good performance it is desirable for different APFloats
2677 // to produce different integers.
2678 uint32_t
2679 APFloat::getHashValue() const
2680 {
2681   if (category==fcZero) return sign<<8 | semantics->precision ;
2682   else if (category==fcInfinity) return sign<<9 | semantics->precision;
2683   else if (category==fcNaN) return 1<<10 | semantics->precision;
2684   else {
2685     uint32_t hash = sign<<11 | semantics->precision | exponent<<12;
2686     const integerPart* p = significandParts();
2687     for (int i=partCount(); i>0; i--, p++)
2688       hash ^= ((uint32_t)*p) ^ (uint32_t)((*p)>>32);
2689     return hash;
2690   }
2691 }
2692
2693 // Conversion from APFloat to/from host float/double.  It may eventually be
2694 // possible to eliminate these and have everybody deal with APFloats, but that
2695 // will take a while.  This approach will not easily extend to long double.
2696 // Current implementation requires integerPartWidth==64, which is correct at
2697 // the moment but could be made more general.
2698
2699 // Denormals have exponent minExponent in APFloat, but minExponent-1 in
2700 // the actual IEEE respresentations.  We compensate for that here.
2701
2702 APInt
2703 APFloat::convertF80LongDoubleAPFloatToAPInt() const
2704 {
2705   assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended);
2706   assert(partCount()==2);
2707
2708   uint64_t myexponent, mysignificand;
2709
2710   if (category==fcNormal) {
2711     myexponent = exponent+16383; //bias
2712     mysignificand = significandParts()[0];
2713     if (myexponent==1 && !(mysignificand & 0x8000000000000000ULL))
2714       myexponent = 0;   // denormal
2715   } else if (category==fcZero) {
2716     myexponent = 0;
2717     mysignificand = 0;
2718   } else if (category==fcInfinity) {
2719     myexponent = 0x7fff;
2720     mysignificand = 0x8000000000000000ULL;
2721   } else {
2722     assert(category == fcNaN && "Unknown category");
2723     myexponent = 0x7fff;
2724     mysignificand = significandParts()[0];
2725   }
2726
2727   uint64_t words[2];
2728   words[0] = mysignificand;
2729   words[1] =  ((uint64_t)(sign & 1) << 15) |
2730               (myexponent & 0x7fffLL);
2731   return APInt(80, 2, words);
2732 }
2733
2734 APInt
2735 APFloat::convertPPCDoubleDoubleAPFloatToAPInt() const
2736 {
2737   assert(semantics == (const llvm::fltSemantics*)&PPCDoubleDouble);
2738   assert(partCount()==2);
2739
2740   uint64_t myexponent, mysignificand, myexponent2, mysignificand2;
2741
2742   if (category==fcNormal) {
2743     myexponent = exponent + 1023; //bias
2744     myexponent2 = exponent2 + 1023;
2745     mysignificand = significandParts()[0];
2746     mysignificand2 = significandParts()[1];
2747     if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2748       myexponent = 0;   // denormal
2749     if (myexponent2==1 && !(mysignificand2 & 0x10000000000000LL))
2750       myexponent2 = 0;   // denormal
2751   } else if (category==fcZero) {
2752     myexponent = 0;
2753     mysignificand = 0;
2754     myexponent2 = 0;
2755     mysignificand2 = 0;
2756   } else if (category==fcInfinity) {
2757     myexponent = 0x7ff;
2758     myexponent2 = 0;
2759     mysignificand = 0;
2760     mysignificand2 = 0;
2761   } else {
2762     assert(category == fcNaN && "Unknown category");
2763     myexponent = 0x7ff;
2764     mysignificand = significandParts()[0];
2765     myexponent2 = exponent2;
2766     mysignificand2 = significandParts()[1];
2767   }
2768
2769   uint64_t words[2];
2770   words[0] =  ((uint64_t)(sign & 1) << 63) |
2771               ((myexponent & 0x7ff) <<  52) |
2772               (mysignificand & 0xfffffffffffffLL);
2773   words[1] =  ((uint64_t)(sign2 & 1) << 63) |
2774               ((myexponent2 & 0x7ff) <<  52) |
2775               (mysignificand2 & 0xfffffffffffffLL);
2776   return APInt(128, 2, words);
2777 }
2778
2779 APInt
2780 APFloat::convertQuadrupleAPFloatToAPInt() const
2781 {
2782   assert(semantics == (const llvm::fltSemantics*)&IEEEquad);
2783   assert(partCount()==2);
2784
2785   uint64_t myexponent, mysignificand, mysignificand2;
2786
2787   if (category==fcNormal) {
2788     myexponent = exponent+16383; //bias
2789     mysignificand = significandParts()[0];
2790     mysignificand2 = significandParts()[1];
2791     if (myexponent==1 && !(mysignificand2 & 0x1000000000000LL))
2792       myexponent = 0;   // denormal
2793   } else if (category==fcZero) {
2794     myexponent = 0;
2795     mysignificand = mysignificand2 = 0;
2796   } else if (category==fcInfinity) {
2797     myexponent = 0x7fff;
2798     mysignificand = mysignificand2 = 0;
2799   } else {
2800     assert(category == fcNaN && "Unknown category!");
2801     myexponent = 0x7fff;
2802     mysignificand = significandParts()[0];
2803     mysignificand2 = significandParts()[1];
2804   }
2805
2806   uint64_t words[2];
2807   words[0] = mysignificand;
2808   words[1] = ((uint64_t)(sign & 1) << 63) |
2809              ((myexponent & 0x7fff) << 48) |
2810              (mysignificand2 & 0xffffffffffffLL);
2811
2812   return APInt(128, 2, words);
2813 }
2814
2815 APInt
2816 APFloat::convertDoubleAPFloatToAPInt() const
2817 {
2818   assert(semantics == (const llvm::fltSemantics*)&IEEEdouble);
2819   assert(partCount()==1);
2820
2821   uint64_t myexponent, mysignificand;
2822
2823   if (category==fcNormal) {
2824     myexponent = exponent+1023; //bias
2825     mysignificand = *significandParts();
2826     if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2827       myexponent = 0;   // denormal
2828   } else if (category==fcZero) {
2829     myexponent = 0;
2830     mysignificand = 0;
2831   } else if (category==fcInfinity) {
2832     myexponent = 0x7ff;
2833     mysignificand = 0;
2834   } else {
2835     assert(category == fcNaN && "Unknown category!");
2836     myexponent = 0x7ff;
2837     mysignificand = *significandParts();
2838   }
2839
2840   return APInt(64, ((((uint64_t)(sign & 1) << 63) |
2841                      ((myexponent & 0x7ff) <<  52) |
2842                      (mysignificand & 0xfffffffffffffLL))));
2843 }
2844
2845 APInt
2846 APFloat::convertFloatAPFloatToAPInt() const
2847 {
2848   assert(semantics == (const llvm::fltSemantics*)&IEEEsingle);
2849   assert(partCount()==1);
2850
2851   uint32_t myexponent, mysignificand;
2852
2853   if (category==fcNormal) {
2854     myexponent = exponent+127; //bias
2855     mysignificand = (uint32_t)*significandParts();
2856     if (myexponent == 1 && !(mysignificand & 0x800000))
2857       myexponent = 0;   // denormal
2858   } else if (category==fcZero) {
2859     myexponent = 0;
2860     mysignificand = 0;
2861   } else if (category==fcInfinity) {
2862     myexponent = 0xff;
2863     mysignificand = 0;
2864   } else {
2865     assert(category == fcNaN && "Unknown category!");
2866     myexponent = 0xff;
2867     mysignificand = (uint32_t)*significandParts();
2868   }
2869
2870   return APInt(32, (((sign&1) << 31) | ((myexponent&0xff) << 23) |
2871                     (mysignificand & 0x7fffff)));
2872 }
2873
2874 APInt
2875 APFloat::convertHalfAPFloatToAPInt() const
2876 {
2877   assert(semantics == (const llvm::fltSemantics*)&IEEEhalf);
2878   assert(partCount()==1);
2879
2880   uint32_t myexponent, mysignificand;
2881
2882   if (category==fcNormal) {
2883     myexponent = exponent+15; //bias
2884     mysignificand = (uint32_t)*significandParts();
2885     if (myexponent == 1 && !(mysignificand & 0x400))
2886       myexponent = 0;   // denormal
2887   } else if (category==fcZero) {
2888     myexponent = 0;
2889     mysignificand = 0;
2890   } else if (category==fcInfinity) {
2891     myexponent = 0x1f;
2892     mysignificand = 0;
2893   } else {
2894     assert(category == fcNaN && "Unknown category!");
2895     myexponent = 0x1f;
2896     mysignificand = (uint32_t)*significandParts();
2897   }
2898
2899   return APInt(16, (((sign&1) << 15) | ((myexponent&0x1f) << 10) |
2900                     (mysignificand & 0x3ff)));
2901 }
2902
2903 // This function creates an APInt that is just a bit map of the floating
2904 // point constant as it would appear in memory.  It is not a conversion,
2905 // and treating the result as a normal integer is unlikely to be useful.
2906
2907 APInt
2908 APFloat::bitcastToAPInt() const
2909 {
2910   if (semantics == (const llvm::fltSemantics*)&IEEEhalf)
2911     return convertHalfAPFloatToAPInt();
2912
2913   if (semantics == (const llvm::fltSemantics*)&IEEEsingle)
2914     return convertFloatAPFloatToAPInt();
2915
2916   if (semantics == (const llvm::fltSemantics*)&IEEEdouble)
2917     return convertDoubleAPFloatToAPInt();
2918
2919   if (semantics == (const llvm::fltSemantics*)&IEEEquad)
2920     return convertQuadrupleAPFloatToAPInt();
2921
2922   if (semantics == (const llvm::fltSemantics*)&PPCDoubleDouble)
2923     return convertPPCDoubleDoubleAPFloatToAPInt();
2924
2925   assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended &&
2926          "unknown format!");
2927   return convertF80LongDoubleAPFloatToAPInt();
2928 }
2929
2930 float
2931 APFloat::convertToFloat() const
2932 {
2933   assert(semantics == (const llvm::fltSemantics*)&IEEEsingle &&
2934          "Float semantics are not IEEEsingle");
2935   APInt api = bitcastToAPInt();
2936   return api.bitsToFloat();
2937 }
2938
2939 double
2940 APFloat::convertToDouble() const
2941 {
2942   assert(semantics == (const llvm::fltSemantics*)&IEEEdouble &&
2943          "Float semantics are not IEEEdouble");
2944   APInt api = bitcastToAPInt();
2945   return api.bitsToDouble();
2946 }
2947
2948 /// Integer bit is explicit in this format.  Intel hardware (387 and later)
2949 /// does not support these bit patterns:
2950 ///  exponent = all 1's, integer bit 0, significand 0 ("pseudoinfinity")
2951 ///  exponent = all 1's, integer bit 0, significand nonzero ("pseudoNaN")
2952 ///  exponent = 0, integer bit 1 ("pseudodenormal")
2953 ///  exponent!=0 nor all 1's, integer bit 0 ("unnormal")
2954 /// At the moment, the first two are treated as NaNs, the second two as Normal.
2955 void
2956 APFloat::initFromF80LongDoubleAPInt(const APInt &api)
2957 {
2958   assert(api.getBitWidth()==80);
2959   uint64_t i1 = api.getRawData()[0];
2960   uint64_t i2 = api.getRawData()[1];
2961   uint64_t myexponent = (i2 & 0x7fff);
2962   uint64_t mysignificand = i1;
2963
2964   initialize(&APFloat::x87DoubleExtended);
2965   assert(partCount()==2);
2966
2967   sign = static_cast<unsigned int>(i2>>15);
2968   if (myexponent==0 && mysignificand==0) {
2969     // exponent, significand meaningless
2970     category = fcZero;
2971   } else if (myexponent==0x7fff && mysignificand==0x8000000000000000ULL) {
2972     // exponent, significand meaningless
2973     category = fcInfinity;
2974   } else if (myexponent==0x7fff && mysignificand!=0x8000000000000000ULL) {
2975     // exponent meaningless
2976     category = fcNaN;
2977     significandParts()[0] = mysignificand;
2978     significandParts()[1] = 0;
2979   } else {
2980     category = fcNormal;
2981     exponent = myexponent - 16383;
2982     significandParts()[0] = mysignificand;
2983     significandParts()[1] = 0;
2984     if (myexponent==0)          // denormal
2985       exponent = -16382;
2986   }
2987 }
2988
2989 void
2990 APFloat::initFromPPCDoubleDoubleAPInt(const APInt &api)
2991 {
2992   assert(api.getBitWidth()==128);
2993   uint64_t i1 = api.getRawData()[0];
2994   uint64_t i2 = api.getRawData()[1];
2995   uint64_t myexponent = (i1 >> 52) & 0x7ff;
2996   uint64_t mysignificand = i1 & 0xfffffffffffffLL;
2997   uint64_t myexponent2 = (i2 >> 52) & 0x7ff;
2998   uint64_t mysignificand2 = i2 & 0xfffffffffffffLL;
2999
3000   initialize(&APFloat::PPCDoubleDouble);
3001   assert(partCount()==2);
3002
3003   sign = static_cast<unsigned int>(i1>>63);
3004   sign2 = static_cast<unsigned int>(i2>>63);
3005   if (myexponent==0 && mysignificand==0) {
3006     // exponent, significand meaningless
3007     // exponent2 and significand2 are required to be 0; we don't check
3008     category = fcZero;
3009   } else if (myexponent==0x7ff && mysignificand==0) {
3010     // exponent, significand meaningless
3011     // exponent2 and significand2 are required to be 0; we don't check
3012     category = fcInfinity;
3013   } else if (myexponent==0x7ff && mysignificand!=0) {
3014     // exponent meaningless.  So is the whole second word, but keep it
3015     // for determinism.
3016     category = fcNaN;
3017     exponent2 = myexponent2;
3018     significandParts()[0] = mysignificand;
3019     significandParts()[1] = mysignificand2;
3020   } else {
3021     category = fcNormal;
3022     // Note there is no category2; the second word is treated as if it is
3023     // fcNormal, although it might be something else considered by itself.
3024     exponent = myexponent - 1023;
3025     exponent2 = myexponent2 - 1023;
3026     significandParts()[0] = mysignificand;
3027     significandParts()[1] = mysignificand2;
3028     if (myexponent==0)          // denormal
3029       exponent = -1022;
3030     else
3031       significandParts()[0] |= 0x10000000000000LL;  // integer bit
3032     if (myexponent2==0)
3033       exponent2 = -1022;
3034     else
3035       significandParts()[1] |= 0x10000000000000LL;  // integer bit
3036   }
3037 }
3038
3039 void
3040 APFloat::initFromQuadrupleAPInt(const APInt &api)
3041 {
3042   assert(api.getBitWidth()==128);
3043   uint64_t i1 = api.getRawData()[0];
3044   uint64_t i2 = api.getRawData()[1];
3045   uint64_t myexponent = (i2 >> 48) & 0x7fff;
3046   uint64_t mysignificand  = i1;
3047   uint64_t mysignificand2 = i2 & 0xffffffffffffLL;
3048
3049   initialize(&APFloat::IEEEquad);
3050   assert(partCount()==2);
3051
3052   sign = static_cast<unsigned int>(i2>>63);
3053   if (myexponent==0 &&
3054       (mysignificand==0 && mysignificand2==0)) {
3055     // exponent, significand meaningless
3056     category = fcZero;
3057   } else if (myexponent==0x7fff &&
3058              (mysignificand==0 && mysignificand2==0)) {
3059     // exponent, significand meaningless
3060     category = fcInfinity;
3061   } else if (myexponent==0x7fff &&
3062              (mysignificand!=0 || mysignificand2 !=0)) {
3063     // exponent meaningless
3064     category = fcNaN;
3065     significandParts()[0] = mysignificand;
3066     significandParts()[1] = mysignificand2;
3067   } else {
3068     category = fcNormal;
3069     exponent = myexponent - 16383;
3070     significandParts()[0] = mysignificand;
3071     significandParts()[1] = mysignificand2;
3072     if (myexponent==0)          // denormal
3073       exponent = -16382;
3074     else
3075       significandParts()[1] |= 0x1000000000000LL;  // integer bit
3076   }
3077 }
3078
3079 void
3080 APFloat::initFromDoubleAPInt(const APInt &api)
3081 {
3082   assert(api.getBitWidth()==64);
3083   uint64_t i = *api.getRawData();
3084   uint64_t myexponent = (i >> 52) & 0x7ff;
3085   uint64_t mysignificand = i & 0xfffffffffffffLL;
3086
3087   initialize(&APFloat::IEEEdouble);
3088   assert(partCount()==1);
3089
3090   sign = static_cast<unsigned int>(i>>63);
3091   if (myexponent==0 && mysignificand==0) {
3092     // exponent, significand meaningless
3093     category = fcZero;
3094   } else if (myexponent==0x7ff && mysignificand==0) {
3095     // exponent, significand meaningless
3096     category = fcInfinity;
3097   } else if (myexponent==0x7ff && mysignificand!=0) {
3098     // exponent meaningless
3099     category = fcNaN;
3100     *significandParts() = mysignificand;
3101   } else {
3102     category = fcNormal;
3103     exponent = myexponent - 1023;
3104     *significandParts() = mysignificand;
3105     if (myexponent==0)          // denormal
3106       exponent = -1022;
3107     else
3108       *significandParts() |= 0x10000000000000LL;  // integer bit
3109   }
3110 }
3111
3112 void
3113 APFloat::initFromFloatAPInt(const APInt & api)
3114 {
3115   assert(api.getBitWidth()==32);
3116   uint32_t i = (uint32_t)*api.getRawData();
3117   uint32_t myexponent = (i >> 23) & 0xff;
3118   uint32_t mysignificand = i & 0x7fffff;
3119
3120   initialize(&APFloat::IEEEsingle);
3121   assert(partCount()==1);
3122
3123   sign = i >> 31;
3124   if (myexponent==0 && mysignificand==0) {
3125     // exponent, significand meaningless
3126     category = fcZero;
3127   } else if (myexponent==0xff && mysignificand==0) {
3128     // exponent, significand meaningless
3129     category = fcInfinity;
3130   } else if (myexponent==0xff && mysignificand!=0) {
3131     // sign, exponent, significand meaningless
3132     category = fcNaN;
3133     *significandParts() = mysignificand;
3134   } else {
3135     category = fcNormal;
3136     exponent = myexponent - 127;  //bias
3137     *significandParts() = mysignificand;
3138     if (myexponent==0)    // denormal
3139       exponent = -126;
3140     else
3141       *significandParts() |= 0x800000; // integer bit
3142   }
3143 }
3144
3145 void
3146 APFloat::initFromHalfAPInt(const APInt & api)
3147 {
3148   assert(api.getBitWidth()==16);
3149   uint32_t i = (uint32_t)*api.getRawData();
3150   uint32_t myexponent = (i >> 10) & 0x1f;
3151   uint32_t mysignificand = i & 0x3ff;
3152
3153   initialize(&APFloat::IEEEhalf);
3154   assert(partCount()==1);
3155
3156   sign = i >> 15;
3157   if (myexponent==0 && mysignificand==0) {
3158     // exponent, significand meaningless
3159     category = fcZero;
3160   } else if (myexponent==0x1f && mysignificand==0) {
3161     // exponent, significand meaningless
3162     category = fcInfinity;
3163   } else if (myexponent==0x1f && mysignificand!=0) {
3164     // sign, exponent, significand meaningless
3165     category = fcNaN;
3166     *significandParts() = mysignificand;
3167   } else {
3168     category = fcNormal;
3169     exponent = myexponent - 15;  //bias
3170     *significandParts() = mysignificand;
3171     if (myexponent==0)    // denormal
3172       exponent = -14;
3173     else
3174       *significandParts() |= 0x400; // integer bit
3175   }
3176 }
3177
3178 /// Treat api as containing the bits of a floating point number.  Currently
3179 /// we infer the floating point type from the size of the APInt.  The
3180 /// isIEEE argument distinguishes between PPC128 and IEEE128 (not meaningful
3181 /// when the size is anything else).
3182 void
3183 APFloat::initFromAPInt(const APInt& api, bool isIEEE)
3184 {
3185   if (api.getBitWidth() == 16)
3186     return initFromHalfAPInt(api);
3187   else if (api.getBitWidth() == 32)
3188     return initFromFloatAPInt(api);
3189   else if (api.getBitWidth()==64)
3190     return initFromDoubleAPInt(api);
3191   else if (api.getBitWidth()==80)
3192     return initFromF80LongDoubleAPInt(api);
3193   else if (api.getBitWidth()==128)
3194     return (isIEEE ?
3195             initFromQuadrupleAPInt(api) : initFromPPCDoubleDoubleAPInt(api));
3196   else
3197     llvm_unreachable(0);
3198 }
3199
3200 APFloat
3201 APFloat::getAllOnesValue(unsigned BitWidth, bool isIEEE)
3202 {
3203   return APFloat(APInt::getAllOnesValue(BitWidth), isIEEE);
3204 }
3205
3206 APFloat APFloat::getLargest(const fltSemantics &Sem, bool Negative) {
3207   APFloat Val(Sem, fcNormal, Negative);
3208
3209   // We want (in interchange format):
3210   //   sign = {Negative}
3211   //   exponent = 1..10
3212   //   significand = 1..1
3213
3214   Val.exponent = Sem.maxExponent; // unbiased
3215
3216   // 1-initialize all bits....
3217   Val.zeroSignificand();
3218   integerPart *significand = Val.significandParts();
3219   unsigned N = partCountForBits(Sem.precision);
3220   for (unsigned i = 0; i != N; ++i)
3221     significand[i] = ~((integerPart) 0);
3222
3223   // ...and then clear the top bits for internal consistency.
3224   significand[N-1] &=
3225     (((integerPart) 1) << ((Sem.precision % integerPartWidth) - 1)) - 1;
3226
3227   return Val;
3228 }
3229
3230 APFloat APFloat::getSmallest(const fltSemantics &Sem, bool Negative) {
3231   APFloat Val(Sem, fcNormal, Negative);
3232
3233   // We want (in interchange format):
3234   //   sign = {Negative}
3235   //   exponent = 0..0
3236   //   significand = 0..01
3237
3238   Val.exponent = Sem.minExponent; // unbiased
3239   Val.zeroSignificand();
3240   Val.significandParts()[0] = 1;
3241   return Val;
3242 }
3243
3244 APFloat APFloat::getSmallestNormalized(const fltSemantics &Sem, bool Negative) {
3245   APFloat Val(Sem, fcNormal, Negative);
3246
3247   // We want (in interchange format):
3248   //   sign = {Negative}
3249   //   exponent = 0..0
3250   //   significand = 10..0
3251
3252   Val.exponent = Sem.minExponent;
3253   Val.zeroSignificand();
3254   Val.significandParts()[partCountForBits(Sem.precision)-1] |=
3255     (((integerPart) 1) << ((Sem.precision % integerPartWidth) - 1));
3256
3257   return Val;
3258 }
3259
3260 APFloat::APFloat(const APInt& api, bool isIEEE) : exponent2(0), sign2(0) {
3261   initFromAPInt(api, isIEEE);
3262 }
3263
3264 APFloat::APFloat(float f) : exponent2(0), sign2(0) {
3265   initFromAPInt(APInt::floatToBits(f));
3266 }
3267
3268 APFloat::APFloat(double d) : exponent2(0), sign2(0) {
3269   initFromAPInt(APInt::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 = 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) const {
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 = 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 = 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 }
3565
3566 bool APFloat::getExactInverse(APFloat *inv) const {
3567   // We can only guarantee the existence of an exact inverse for IEEE floats.
3568   if (semantics != &IEEEhalf && semantics != &IEEEsingle &&
3569       semantics != &IEEEdouble && semantics != &IEEEquad)
3570     return false;
3571
3572   // Special floats and denormals have no exact inverse.
3573   if (category != fcNormal)
3574     return false;
3575
3576   // Check that the number is a power of two by making sure that only the
3577   // integer bit is set in the significand.
3578   if (significandLSB() != semantics->precision - 1)
3579     return false;
3580
3581   // Get the inverse.
3582   APFloat reciprocal(*semantics, 1ULL);
3583   if (reciprocal.divide(*this, rmNearestTiesToEven) != opOK)
3584     return false;
3585
3586   // Avoid multiplication with a denormal, it is not safe on all platforms and
3587   // may be slower than a normal division.
3588   if (reciprocal.significandMSB() + 1 < reciprocal.semantics->precision)
3589     return false;
3590
3591   assert(reciprocal.category == fcNormal &&
3592          reciprocal.significandLSB() == reciprocal.semantics->precision - 1);
3593
3594   if (inv)
3595     *inv = reciprocal;
3596
3597   return true;
3598 }