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