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