lib/Support/APFloat.cpp

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