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