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