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