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