lib/Support/APFloat.cpp

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