initial checkin of Neil's APFloat work.

git-svn-id: https://llvm.org/svn/llvm-project/llvm/trunk@41203 91177308-0d34-0410-b5e6-96231b3b80d8

4 files changed, 1780 insertions, 13 deletions

diff --git a/include/llvm/ADT/APFloat.h b/include/llvm/ADT/APFloat.h
new file mode 100644
index 0000000000..89614d5d36
--- /dev/null
+++ b/include/llvm/ADT/APFloat.h
@@ -0,0 +1,258 @@
+//== llvm/Support/APFloat.h - Arbitrary Precision Floating Point -*- C++ -*-==//
+//
+//                     The LLVM Compiler Infrastructure
+//
+// This file was developed by Neil Booth and is distributed under the
+// University of Illinois Open Source License. See LICENSE.TXT for details.
+//
+//===----------------------------------------------------------------------===//
+//
+// This file declares a class to represent arbitrary precision floating
+// point values and provide a variety of arithmetic operations on them.
+//
+//===----------------------------------------------------------------------===//
+
+/*  A self-contained host- and target-independent arbitrary-precision
+    floating-point software implementation using bignum integer
+    arithmetic, as provided by static functions in the APInt class.
+    The library will work with bignum integers whose parts are any
+    unsigned type at least 16 bits wide.  64 bits is recommended.
+
+    Written for clarity rather than speed, in particular with a view
+    to use in the front-end of a cross compiler so that target
+    arithmetic can be correctly performed on the host.  Performance
+    should nonetheless be reasonable, particularly for its intended
+    use.  It may be useful as a base implementation for a run-time
+    library during development of a faster target-specific one.
+
+    All 5 rounding modes in the IEEE-754R draft are handled correctly
+    for all implemented operations.  Currently implemented operations
+    are add, subtract, multiply, divide, fused-multiply-add,
+    conversion-to-float, conversion-to-integer and
+    conversion-from-integer.  New rounding modes (e.g. away from zero)
+    can be added with three or four lines of code.  The library reads
+    and correctly rounds hexadecimal floating point numbers as per
+    C99; syntax is required to have been validated by the caller.
+    Conversion from decimal is not currently implemented.
+
+    Four formats are built-in: IEEE single precision, double
+    precision, quadruple precision, and x87 80-bit extended double
+    (when operating with full extended precision).  Adding a new
+    format that obeys IEEE semantics only requires adding two lines of
+    code: a declaration and definition of the format.
+
+    All operations return the status of that operation as an exception
+    bit-mask, so multiple operations can be done consecutively with
+    their results or-ed together.  The returned status can be useful
+    for compiler diagnostics; e.g., inexact, underflow and overflow
+    can be easily diagnosed on constant folding, and compiler
+    optimizers can determine what exceptions would be raised by
+    folding operations and optimize, or perhaps not optimize,
+    accordingly.
+
+    At present, underflow tininess is detected after rounding; it
+    should be straight forward to add support for the before-rounding
+    case too.
+
+    Non-zero finite numbers are represented internally as a sign bit,
+    a 16-bit signed exponent, and the significand as an array of
+    integer parts.  After normalization of a number of precision P the
+    exponent is within the range of the format, and if the number is
+    not denormal the P-th bit of the significand is set as an explicit
+    integer bit.  For denormals the most significant bit is shifted
+    right so that the exponent is maintained at the format's minimum,
+    so that the smallest denormal has just the least significant bit
+    of the significand set.  The sign of zeroes and infinities is
+    significant; the exponent and significand of such numbers is
+    indeterminate and meaningless.  For QNaNs the sign bit, as well as
+    the exponent and significand are indeterminate and meaningless.
+
+    TODO
+    ====
+
+    Some features that may or may not be worth adding:
+
+    Conversions to and from decimal strings (hard).
+
+    Conversions to hexadecimal string.
+
+    Read and write IEEE-format in-memory representations.
+
+    Optional ability to detect underflow tininess before rounding.
+
+    New formats: x87 in single and double precision mode (IEEE apart
+    from extended exponent range) and IBM two-double extended
+    precision (hard).
+
+    New operations: sqrt, copysign, nextafter, nexttoward.
+*/
+
+#ifndef LLVM_FLOAT_H
+#define LLVM_FLOAT_H
+
+// APInt contains static functions implementing bignum arithmetic.
+#include "llvm/ADT/APInt.h"
+
+namespace llvm {
+
+  /* Exponents are stored as signed numbers.  */
+  typedef signed short exponent_t;
+
+  struct fltSemantics;
+
+  /* When bits of a floating point number are truncated, this enum is
+     used to indicate what fraction of the LSB those bits represented.
+     It essentially combines the roles of guard and sticky bits.  */
+  enum lostFraction {		// Example of truncated bits:
+    lfExactlyZero,		// 000000
+    lfLessThanHalf,		// 0xxxxx  x's not all zero
+    lfExactlyHalf,		// 100000
+    lfMoreThanHalf		// 1xxxxx  x's not all zero
+  };
+
+  class APFloat {
+  public:
+
+    /* We support the following floating point semantics.  */
+    static const fltSemantics IEEEsingle;
+    static const fltSemantics IEEEdouble;
+    static const fltSemantics IEEEquad;
+    static const fltSemantics x87DoubleExtended;
+
+    static unsigned int semanticsPrecision(const fltSemantics &);
+
+    /* Floating point numbers have a four-state comparison relation.  */
+    enum cmpResult {
+      cmpLessThan,
+      cmpEqual,
+      cmpGreaterThan,
+      cmpUnordered
+    };
+
+    /* IEEE-754R gives five rounding modes.  */
+    enum roundingMode {
+      rmNearestTiesToEven,
+      rmTowardPositive,
+      rmTowardNegative,
+      rmTowardZero,
+      rmNearestTiesToAway
+    };
+
+    /* Operation status.  opUnderflow or opOverflow are always returned
+       or-ed with opInexact.  */
+    enum opStatus {
+      opOK          = 0x00,
+      opInvalidOp   = 0x01,
+      opDivByZero   = 0x02,
+      opOverflow    = 0x04,
+      opUnderflow   = 0x08,
+      opInexact     = 0x10
+    };
+
+    /* Category of internally-represented number.  */
+    enum fltCategory {
+      fcInfinity,
+      fcQNaN,
+      fcNormal,
+      fcZero
+    };
+
+    /* Constructors.  */
+    APFloat(const fltSemantics &, const char *);
+    APFloat(const fltSemantics &, integerPart);
+    APFloat(const fltSemantics &, fltCategory, bool negative);
+    APFloat(const APFloat &);
+    ~APFloat();
+
+    /* Arithmetic.  */
+    opStatus add(const APFloat &, roundingMode);
+    opStatus subtract(const APFloat &, roundingMode);
+    opStatus multiply(const APFloat &, roundingMode);
+    opStatus divide(const APFloat &, roundingMode);
+    opStatus fusedMultiplyAdd(const APFloat &, const APFloat &, roundingMode);
+    void changeSign();
+
+    /* Conversions.  */
+    opStatus convert(const fltSemantics &, roundingMode);
+    opStatus convertToInteger(integerPart *, unsigned int, bool,
+			      roundingMode) const;
+    opStatus convertFromInteger(const integerPart *, unsigned int, bool,
+				roundingMode);
+    opStatus convertFromString(const char *, roundingMode);
+
+    /* Comparison with another floating point number.  */
+    cmpResult compare(const APFloat &) const;
+
+    /* Simple queries.  */
+    fltCategory getCategory() const { return category; }
+    const fltSemantics &getSemantics() const { return *semantics; }
+    bool isZero() const { return category == fcZero; }
+    bool isNonZero() const { return category != fcZero; }
+    bool isNegative() const { return sign; }
+
+    APFloat& operator=(const APFloat &);
+
+  private:
+
+    /* Trivial queries.  */
+    integerPart *significandParts();
+    const integerPart *significandParts() const;
+    unsigned int partCount() const;
+
+    /* Significand operations.  */
+    integerPart addSignificand(const APFloat &);
+    integerPart subtractSignificand(const APFloat &, integerPart);
+    lostFraction addOrSubtractSignificand(const APFloat &, bool subtract);
+    lostFraction multiplySignificand(const APFloat &, const APFloat *);
+    lostFraction divideSignificand(const APFloat &);
+    void incrementSignificand();
+    void initialize(const fltSemantics *);
+    void shiftSignificandLeft(unsigned int);
+    lostFraction shiftSignificandRight(unsigned int);
+    unsigned int significandLSB() const;
+    unsigned int significandMSB() const;
+    void zeroSignificand();
+
+    /* Arithmetic on special values.  */
+    opStatus addOrSubtractSpecials(const APFloat &, bool subtract);
+    opStatus divideSpecials(const APFloat &);
+    opStatus multiplySpecials(const APFloat &);
+
+    /* Miscellany.  */
+    opStatus normalize(roundingMode, lostFraction);
+    opStatus addOrSubtract(const APFloat &, roundingMode, bool subtract);
+    cmpResult compareAbsoluteValue(const APFloat &) const;
+    opStatus handleOverflow(roundingMode);
+    bool roundAwayFromZero(roundingMode, lostFraction);
+    opStatus convertFromUnsignedInteger(integerPart *, unsigned int,
+					roundingMode);
+    lostFraction combineLostFractions(lostFraction, lostFraction);
+    opStatus convertFromHexadecimalString(const char *, roundingMode);
+
+    void assign(const APFloat &);
+    void copySignificand(const APFloat &);
+    void freeSignificand();
+
+    /* What kind of semantics does this value obey?  */
+    const fltSemantics *semantics;
+
+    /* Significand - the fraction with an explicit integer bit.  Must be
+       at least one bit wider than the target precision.  */
+    union Significand
+    {
+      integerPart part;
+      integerPart *parts;
+    } significand;
+
+    /* The exponent - a signed number.  */
+    exponent_t exponent;
+
+    /* What kind of floating point number this is.  */
+    fltCategory category: 2;
+
+    /* The sign bit of this number.  */
+    unsigned int sign: 1;
+  };
+} /* namespace llvm */
+
+#endif /* LLVM_FLOAT_H */
diff --git a/include/llvm/ADT/APInt.h b/include/llvm/ADT/APInt.h
index ce3696fc1e..971822f1bd 100644
--- a/include/llvm/ADT/APInt.h
+++ b/include/llvm/ADT/APInt.h
@@ -19,9 +19,7 @@
 #include <cassert>
 #include <string>
 
-#define HOST_CHAR_BIT 8
-#define compileTimeAssert(cond) extern int CTAssert[(cond) ? 1 : -1]
-#define integerPartWidth (HOST_CHAR_BIT * sizeof(llvm::integerPart))
+#define COMPILE_TIME_ASSERT(cond) extern int CTAssert[(cond) ? 1 : -1]
 
 namespace llvm {
 
@@ -29,6 +27,9 @@ namespace llvm {
      bignum.  */
   typedef uint64_t integerPart;
 
+  const unsigned int host_char_bit = 8;
+  const unsigned int integerPartWidth = host_char_bit * sizeof(integerPart);
+
 //===----------------------------------------------------------------------===//
 //                              APInt Class
 //===----------------------------------------------------------------------===//
diff --git a/lib/Support/APFloat.cpp b/lib/Support/APFloat.cpp
new file mode 100644
index 0000000000..8ac92475db
--- /dev/null
+++ b/lib/Support/APFloat.cpp
@@ -0,0 +1,1488 @@
+//===-- APFloat.cpp - Implement APFloat class -----------------------------===//
+//
+//                     The LLVM Compiler Infrastructure
+//
+// This file was developed by Neil Booth and is distributed under the
+// University of Illinois Open Source License. See LICENSE.TXT for details.
+//
+//===----------------------------------------------------------------------===//
+//
+// This file implements a class to represent arbitrary precision floating
+// point values and provide a variety of arithmetic operations on them.
+//
+//===----------------------------------------------------------------------===//
+
+#include <cassert>
+#include "llvm/ADT/APFloat.h"
+
+using namespace llvm;
+
+#define convolve(lhs, rhs) ((lhs) * 4 + (rhs))
+
+/* Assumed in hexadecimal significand parsing.  */
+COMPILE_TIME_ASSERT(integerPartWidth % 4 == 0);
+
+namespace llvm {
+
+  /* Represents floating point arithmetic semantics.  */
+  struct fltSemantics {
+    /* The largest E such that 2^E is representable; this matches the
+       definition of IEEE 754.  */
+    exponent_t maxExponent;
+
+    /* The smallest E such that 2^E is a normalized number; this
+       matches the definition of IEEE 754.  */
+    exponent_t minExponent;
+
+    /* Number of bits in the significand.  This includes the integer
+       bit.  */
+    unsigned char precision;
+
+    /* If the target format has an implicit integer bit.  */
+    bool implicitIntegerBit;
+  };
+
+  const fltSemantics APFloat::IEEEsingle = { 127, -126, 24, true };
+  const fltSemantics APFloat::IEEEdouble = { 1023, -1022, 53, true };
+  const fltSemantics APFloat::IEEEquad = { 16383, -16382, 113, true };
+  const fltSemantics APFloat::x87DoubleExtended = { 16383, -16382, 64, false };
+}
+
+/* Put a bunch of private, handy routines in an anonymous namespace.  */
+namespace {
+
+  inline unsigned int
+  partCountForBits(unsigned int bits)
+  {
+    return ((bits) + integerPartWidth - 1) / integerPartWidth;
+  }
+
+  unsigned int
+  digitValue(unsigned int c)
+  {
+    unsigned int r;
+
+    r = c - '0';
+    if(r <= 9)
+      return r;
+
+    return -1U;
+  }
+
+  unsigned int
+  hexDigitValue (unsigned int c)
+  {
+    unsigned int r;
+
+    r = c - '0';
+    if(r <= 9)
+      return r;
+
+    r = c - 'A';
+    if(r <= 5)
+      return r + 10;
+
+    r = c - 'a';
+    if(r <= 5)
+      return r + 10;
+
+    return -1U;
+  }
+
+  /* This is ugly and needs cleaning up, but I don't immediately see
+     how whilst remaining safe.  */
+  static int
+  totalExponent(const char *p, int exponentAdjustment)
+  {
+    integerPart unsignedExponent;
+    bool negative, overflow;
+    long exponent;
+
+    /* Move past the exponent letter and sign to the digits.  */
+    p++;
+    negative = *p == '-';
+    if(*p == '-' || *p == '+')
+      p++;
+
+    unsignedExponent = 0;
+    overflow = false;
+    for(;;) {
+      unsigned int value;
+
+      value = digitValue(*p);
+      if(value == -1U)
+	break;
+
+      p++;
+      unsignedExponent = unsignedExponent * 10 + value;
+      if(unsignedExponent > 65535)
+	overflow = true;
+    }
+
+    if(exponentAdjustment > 65535 || exponentAdjustment < -65536)
+      overflow = true;
+
+    if(!overflow) {
+      exponent = unsignedExponent;
+      if(negative)
+	exponent = -exponent;
+      exponent += exponentAdjustment;
+      if(exponent > 65535 || exponent < -65536)
+	overflow = true;
+    }
+
+    if(overflow)
+      exponent = negative ? -65536: 65535;
+
+    return exponent;
+  }
+
+  const char *
+  skipLeadingZeroesAndAnyDot(const char *p, const char **dot)
+  {
+    *dot = 0;
+    while(*p == '0')
+      p++;
+
+    if(*p == '.') {
+      *dot = p++;
+      while(*p == '0')
+	p++;
+    }
+
+    return p;
+  }
+
+  /* Return the trailing fraction of a hexadecimal number.
+     DIGITVALUE is the first hex digit of the fraction, P points to
+     the next digit.  */
+  lostFraction
+  trailingHexadecimalFraction(const char *p, unsigned int digitValue)
+  {
+    unsigned int hexDigit;
+
+    /* If the first trailing digit isn't 0 or 8 we can work out the
+       fraction immediately.  */
+    if(digitValue > 8)
+      return lfMoreThanHalf;
+    else if(digitValue < 8 && digitValue > 0)
+      return lfLessThanHalf;
+
+    /* Otherwise we need to find the first non-zero digit.  */
+    while(*p == '0')
+      p++;
+
+    hexDigit = hexDigitValue(*p);
+
+    /* If we ran off the end it is exactly zero or one-half, otherwise
+       a little more.  */
+    if(hexDigit == -1U)
+      return digitValue == 0 ? lfExactlyZero: lfExactlyHalf;
+    else
+      return digitValue == 0 ? lfLessThanHalf: lfMoreThanHalf;
+  }
+
+  /* Return the fraction lost were a bignum truncated.  */
+  lostFraction
+  lostFractionThroughTruncation(integerPart *parts,
+				unsigned int partCount,
+				unsigned int bits)
+  {
+    unsigned int lsb;
+
+    lsb = APInt::tcLSB(parts, partCount);
+
+    /* Note this is guaranteed true if bits == 0, or LSB == -1U.  */
+    if(bits <= lsb)
+      return lfExactlyZero;
+    if(bits == lsb + 1)
+      return lfExactlyHalf;
+    if(bits <= partCount * integerPartWidth
+       && APInt::tcExtractBit(parts, bits - 1))
+      return lfMoreThanHalf;
+
+    return lfLessThanHalf;
+  }
+
+  /* Shift DST right BITS bits noting lost fraction.  */
+  lostFraction
+  shiftRight(integerPart *dst, unsigned int parts, unsigned int bits)
+  {
+    lostFraction lost_fraction;
+
+    lost_fraction = lostFractionThroughTruncation(dst, parts, bits);
+
+    APInt::tcShiftRight(dst, parts, bits);
+
+    return lost_fraction;
+  }
+}
+
+/* Constructors.  */
+void
+APFloat::initialize(const fltSemantics *ourSemantics)
+{
+  unsigned int count;
+
+  semantics = ourSemantics;
+  count = partCount();
+  if(count > 1)
+    significand.parts = new integerPart[count];
+}
+
+void
+APFloat::freeSignificand()
+{
+  if(partCount() > 1)
+    delete [] significand.parts;
+}
+
+void
+APFloat::assign(const APFloat &rhs)
+{
+  assert(semantics == rhs.semantics);
+
+  sign = rhs.sign;
+  category = rhs.category;
+  exponent = rhs.exponent;
+  if(category == fcNormal)
+    copySignificand(rhs);
+}
+
+void
+APFloat::copySignificand(const APFloat &rhs)
+{
+  assert(category == fcNormal);
+  assert(rhs.partCount() >= partCount());
+
+  APInt::tcAssign(significandParts(), rhs.significandParts(),
+		  partCount());
+}
+
+APFloat &
+APFloat::operator=(const APFloat &rhs)
+{
+  if(this != &rhs) {
+    if(semantics != rhs.semantics) {
+      freeSignificand();
+      initialize(rhs.semantics);
+    }
+    assign(rhs);
+  }
+
+  return *this;
+}
+
+APFloat::APFloat(const fltSemantics &ourSemantics, integerPart value)
+{
+  initialize(&ourSemantics);
+  sign = 0;
+  zeroSignificand();
+  exponent = ourSemantics.precision - 1;
+  significandParts()[0] = value;
+  normalize(rmNearestTiesToEven, lfExactlyZero);
+}
+
+APFloat::APFloat(const fltSemantics &ourSemantics,
+		 fltCategory ourCategory, bool negative)
+{
+  initialize(&ourSemantics);
+  category = ourCategory;
+  sign = negative;
+  if(category == fcNormal)
+    category = fcZero;
+}
+
+APFloat::APFloat(const fltSemantics &ourSemantics, const char *text)
+{
+  initialize(&ourSemantics);
+  convertFromString(text, rmNearestTiesToEven);
+}
+
+APFloat::APFloat(const APFloat &rhs)
+{
+  initialize(rhs.semantics);
+  assign(rhs);
+}
+
+APFloat::~APFloat()
+{
+  freeSignificand();
+}
+
+unsigned int
+APFloat::partCount() const
+{
+  return partCountForBits(semantics->precision + 1);
+}
+
+unsigned int
+APFloat::semanticsPrecision(const fltSemantics &semantics)
+{
+  return semantics.precision;
+}
+
+const integerPart *
+APFloat::significandParts() const
+{
+  return const_cast<APFloat *>(this)->significandParts();
+}
+
+integerPart *
+APFloat::significandParts()
+{
+  assert(category == fcNormal);
+
+  if(partCount() > 1)
+    return significand.parts;
+  else
+    return &significand.part;
+}
+
+/* Combine the effect of two lost fractions.  */
+lostFraction
+APFloat::combineLostFractions(lostFraction moreSignificant,
+			      lostFraction lessSignificant)
+{
+  if(lessSignificant != lfExactlyZero) {
+    if(moreSignificant == lfExactlyZero)
+      moreSignificant = lfLessThanHalf;
+    else if(moreSignificant == lfExactlyHalf)
+      moreSignificant = lfMoreThanHalf;
+  }
+
+  return moreSignificant;
+}
+
+void
+APFloat::zeroSignificand()
+{
+  category = fcNormal;
+  APInt::tcSet(significandParts(), 0, partCount());
+}
+
+/* Increment an fcNormal floating point number's significand.  */
+void
+APFloat::incrementSignificand()
+{
+  integerPart carry;
+
+  carry = APInt::tcIncrement(significandParts(), partCount());
+
+  /* Our callers should never cause us to overflow.  */
+  assert(carry == 0);
+}
+
+/* Add the significand of the RHS.  Returns the carry flag.  */
+integerPart
+APFloat::addSignificand(const APFloat &rhs)
+{
+  integerPart *parts;
+
+  parts = significandParts();
+
+  assert(semantics == rhs.semantics);
+  assert(exponent == rhs.exponent);
+
+  return APInt::tcAdd(parts, rhs.significandParts(), 0, partCount());
+}
+
+/* Subtract the significand of the RHS with a borrow flag.  Returns
+   the borrow flag.  */
+integerPart
+APFloat::subtractSignificand(const APFloat &rhs, integerPart borrow)
+{
+  integerPart *parts;
+
+  parts = significandParts();
+
+  assert(semantics == rhs.semantics);
+  assert(exponent == rhs.exponent);
+
+  return APInt::tcSubtract(parts, rhs.significandParts(), borrow,
+			   partCount());
+}
+
+/* Multiply the significand of the RHS.  If ADDEND is non-NULL, add it
+   on to the full-precision result of the multiplication.  Returns the
+   lost fraction.  */
+lostFraction
+APFloat::multiplySignificand(const APFloat &rhs, const APFloat *addend)
+{
+  unsigned int omsb;	// One, not zero, based MSB.
+  unsigned int partsCount, newPartsCount, precision;
+  integerPart *lhsSignificand;
+  integerPart scratch[4];
+  integerPart *fullSignificand;
+  lostFraction lost_fraction;
+
+  assert(semantics == rhs.semantics);
+
+  precision = semantics->precision;
+  newPartsCount = partCountForBits(precision * 2);
+
+  if(newPartsCount > 4)
+    fullSignificand = new integerPart[newPartsCount];
+  else
+    fullSignificand = scratch;
+
+  lhsSignificand = significandParts();
+  partsCount = partCount();
+
+  APInt::tcFullMultiply(fullSignificand, lhsSignificand,
+			rhs.significandParts(), partsCount);
+
+  lost_fraction = lfExactlyZero;
+  omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
+  exponent += rhs.exponent;
+
+  if(addend) {
+    Significand savedSignificand = significand;
+    const fltSemantics *savedSemantics = semantics;
+    fltSemantics extendedSemantics;
+    opStatus status;
+    unsigned int extendedPrecision;
+
+    /* Normalize our MSB.  */
+    extendedPrecision = precision + precision - 1;
+    if(omsb != extendedPrecision)
+      {
+	APInt::tcShiftLeft(fullSignificand, newPartsCount,
+			   extendedPrecision - omsb);
+	exponent -= extendedPrecision - omsb;
+      }
+
+    /* Create new semantics.  */
+    extendedSemantics = *semantics;
+    extendedSemantics.precision = extendedPrecision;
+
+    if(newPartsCount == 1)
+      significand.part = fullSignificand[0];
+    else
+      significand.parts = fullSignificand;
+    semantics = &extendedSemantics;
+
+    APFloat extendedAddend(*addend);
+    status = extendedAddend.convert(extendedSemantics, rmTowardZero);
+    assert(status == opOK);
+    lost_fraction = addOrSubtractSignificand(extendedAddend, false);
+
+    /* Restore our state.  */
+    if(newPartsCount == 1)
+      fullSignificand[0] = significand.part;
+    significand = savedSignificand;
+    semantics = savedSemantics;
+
+    omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
+  }
+
+  exponent -= (precision - 1);
+
+  if(omsb > precision) {
+    unsigned int bits, significantParts;
+    lostFraction lf;
+
+    bits = omsb - precision;
+    significantParts = partCountForBits(omsb);
+    lf = shiftRight(fullSignificand, significantParts, bits);
+    lost_fraction = combineLostFractions(lf, lost_fraction);
+    exponent += bits;
+  }
+
+  APInt::tcAssign(lhsSignificand, fullSignificand, partsCount);
+
+  if(newPartsCount > 4)
+    delete [] fullSignificand;
+
+  return lost_fraction;
+}
+
+/* Multiply the significands of LHS and RHS to DST.  */
+lostFraction
+APFloat::divideSignificand(const APFloat &rhs)
+{
+  unsigned int bit, i, partsCount;
+  const integerPart *rhsSignificand;
+  integerPart *lhsSignificand, *dividend, *divisor;
+  integerPart scratch[4];
+  lostFraction lost_fraction;
+
+  assert(semantics == rhs.semantics);
+
+  lhsSignificand = significandParts();
+  rhsSignificand = rhs.significandParts();
+  partsCount = partCount();
+
+  if(partsCount > 2)
+    dividend = new integerPart[partsCount * 2];
+  else
+    dividend = scratch;
+
+  divisor = dividend + partsCount;
+
+  /* Copy the dividend and divisor as they will be modified in-place.  */
+  for(i = 0; i < partsCount; i++) {
+    dividend[i] = lhsSignificand[i];
+    divisor[i] = rhsSignificand[i];
+    lhsSignificand[i] = 0;
+  }
+
+  exponent -= rhs.exponent;
+
+  unsigned int precision = semantics->precision;
+
+  /* Normalize the divisor.  */
+  bit = precision - APInt::tcMSB(divisor, partsCount) - 1;
+  if(bit) {
+    exponent += bit;
+    APInt::tcShiftLeft(divisor, partsCount, bit);
+  }
+
+  /* Normalize the dividend.  */
+  bit = precision - APInt::tcMSB(dividend, partsCount) - 1;
+  if(bit) {
+    exponent -= bit;
+    APInt::tcShiftLeft(dividend, partsCount, bit);
+  }
+
+  if(APInt::tcCompare(dividend, divisor, partsCount) < 0) {
+    exponent--;
+    APInt::tcShiftLeft(dividend, partsCount, 1);
+    assert(APInt::tcCompare(dividend, divisor, partsCount) >= 0);
+  }
+
+  /* Long division.  */
+  for(bit = precision; bit; bit -= 1) {
+    if(APInt::tcCompare(dividend, divisor, partsCount) >= 0) {
+      APInt::tcSubtract(dividend, divisor, 0, partsCount);
+      APInt::tcSetBit(lhsSignificand, bit - 1);
+    }
+
+    APInt::tcShiftLeft(dividend, partsCount, 1);
+  }
+
+  /* Figure out the lost fraction.  */
+  int cmp = APInt::tcCompare(dividend, divisor, partsCount);
+
+  if(cmp > 0)
+    lost_fraction = lfMoreThanHalf;
+  else if(cmp == 0)
+    lost_fraction = lfExactlyHalf;
+  else if(APInt::tcIsZero(dividend, partsCount))
+    lost_fraction = lfExactlyZero;
+  else
+    lost_fraction = lfLessThanHalf;
+
+  if(partsCount > 2)
+    delete [] dividend;
+
+  return lost_fraction;
+}
+
+unsigned int
+APFloat::significandMSB() const
+{
+  return APInt::tcMSB(significandParts(), partCount());
+}
+
+unsigned int
+APFloat::significandLSB() const
+{
+  return APInt::tcLSB(significandParts(), partCount());
+}
+
+/* Note that a zero result is NOT normalized to fcZero.  */
+lostFraction
+APFloat::shiftSignificandRight(unsigned int bits)
+{
+  /* Our exponent should not overflow.  */
+  assert((exponent_t) (exponent + bits) >= exponent);
+
+  exponent += bits;
+
+  return shiftRight(significandParts(), partCount(), bits);
+}
+
+/* Shift the significand left BITS bits, subtract BITS from its exponent.  */
+void
+APFloat::shiftSignificandLeft(unsigned int bits)
+{
+  assert(bits < semantics->precision);
+
+  if(bits) {
+    unsigned int partsCount = partCount();
+
+    APInt::tcShiftLeft(significandParts(), partsCount, bits);
+    exponent -= bits;
+
+    assert(!APInt::tcIsZero(significandParts(), partsCount));
+  }
+}
+
+APFloat::cmpResult
+APFloat::compareAbsoluteValue(const APFloat &rhs) const
+{
+  int compare;
+
+  assert(semantics == rhs.semantics);
+  assert(category == fcNormal);
+  assert(rhs.category == fcNormal);
+
+  compare = exponent - rhs.exponent;
+
+  /* If exponents are equal, do an unsigned bignum comparison of the
+     significands.  */
+  if(compare == 0)
+    compare = APInt::tcCompare(significandParts(), rhs.significandParts(),
+			       partCount());
+
+  if(compare > 0)
+    return cmpGreaterThan;
+  else if(compare < 0)
+    return cmpLessThan;
+  else
+    return cmpEqual;
+}
+
+/* Handle overflow.  Sign is preserved.  We either become infinity or
+   the largest finite number.  */
+APFloat::opStatus
+APFloat::handleOverflow(roundingMode rounding_mode)
+{
+  /* Infinity?  */
+  if(rounding_mode == rmNearestTiesToEven
+     || rounding_mode == rmNearestTiesToAway
+     || (rounding_mode == rmTowardPositive && !sign)
+     || (rounding_mode == rmTowardNegative && sign))
+    {
+      category = fcInfinity;
+      return (opStatus) (opOverflow | opInexact);
+    }
+
+  /* Otherwise we become the largest finite number.  */
+  category = fcNormal;
+  exponent = semantics->maxExponent;
+  APInt::tcSetLeastSignificantBits(significandParts(), partCount(),
+				   semantics->precision);
+
+  return opInexact;
+}
+
+/* This routine must work for fcZero of both signs, and fcNormal
+   numbers.  */
+bool
+APFloat::roundAwayFromZero(roundingMode rounding_mode,
+			   lostFraction lost_fraction)
+{
+  /* QNaNs and infinities should not have lost fractions.  */
+  assert(category == fcNormal || category == fcZero);
+
+  /* Our caller has already handled this case.  */
+  assert(lost_fraction != lfExactlyZero);
+
+  switch(rounding_mode) {
+  default:
+    assert(0);
+
+  case rmNearestTiesToAway:
+    return lost_fraction == lfExactlyHalf || lost_fraction == lfMoreThanHalf;
+
+  case rmNearestTiesToEven:
+    if(lost_fraction == lfMoreThanHalf)
+      return true;
+
+    /* Our zeroes don't have a significand to test.  */
+    if(lost_fraction == lfExactlyHalf && category != fcZero)
+      return significandParts()[0] & 1;
+
+    return false;
+
+  case rmTowardZero:
+    return false;
+
+  case rmTowardPositive:
+    return sign == false;
+
+  case rmTowardNegative:
+    return sign == true;
+  }
+}
+
+APFloat::opStatus
+APFloat::normalize(roundingMode rounding_mode,
+		   lostFraction lost_fraction)
+{
+  unsigned int omsb;		/* One, not zero, based MSB.  */
+  int exponentChange;
+
+  if(category != fcNormal)
+    return opOK;
+
+  /* Before rounding normalize the exponent of fcNormal numbers.  */
+  omsb = significandMSB() + 1;
+
+  if(omsb) {
+    /* OMSB is numbered from 1.  We want to place it in the integer
+       bit numbered PRECISON if possible, with a compensating change in
+       the exponent.  */
+    exponentChange = omsb - semantics->precision;
+
+    /* If the resulting exponent is too high, overflow according to
+       the rounding mode.  */
+    if(exponent + exponentChange > semantics->maxExponent)
+      return handleOverflow(rounding_mode);
+
+    /* Subnormal numbers have exponent minExponent, and their MSB
+       is forced based on that.  */
+    if(exponent + exponentChange < semantics->minExponent)
+      exponentChange = semantics->minExponent - exponent;
+
+    /* Shifting left is easy as we don't lose precision.  */
+    if(exponentChange < 0) {
+      assert(lost_fraction == lfExactlyZero);
+
+      shiftSignificandLeft(-exponentChange);
+
+      return opOK;
+    }
+
+    if(exponentChange > 0) {
+      lostFraction lf;
+
+      /* Shift right and capture any new lost fraction.  */
+      lf = shiftSignificandRight(exponentChange);
+
+      lost_fraction = combineLostFractions(lf, lost_fraction);
+
+      /* Keep OMSB up-to-date.  */
+      if(omsb > (unsigned) exponentChange)
+	omsb -= (unsigned) exponentChange;
+      else
+	omsb = 0;
+    }
+  }
+
+  /* Now round the number according to rounding_mode given the lost
+     fraction.  */
+
+  /* As specified in IEEE 754, since we do not trap we do not report
+     underflow for exact results.  */
+  if(lost_fraction == lfExactlyZero) {
+    /* Canonicalize zeroes.  */
+    if(omsb == 0)
+      category = fcZero;
+
+    return opOK;
+  }
+
+  /* Increment the significand if we're rounding away from zero.  */
+  if(roundAwayFromZero(rounding_mode, lost_fraction)) {
+    if(omsb == 0)
+      exponent = semantics->minExponent;
+
+    incrementSignificand();
+    omsb = significandMSB() + 1;
+
+    /* Did the significand increment overflow?  */
+    if(omsb == (unsigned) semantics->precision + 1) {
+      /* Renormalize by incrementing the exponent and shifting our
+	 significand right one.  However if we already have the
+	 maximum exponent we overflow to infinity.  */
+      if(exponent == semantics->maxExponent) {
+	category = fcInfinity;
+
+	return (opStatus) (opOverflow | opInexact);
+      }
+
+      shiftSignificandRight(1);
+
+      return opInexact;
+    }
+  }
+
+  /* The normal case - we were and are not denormal, and any
+     significand increment above didn't overflow.  */
+  if(omsb == semantics->precision)
+    return opInexact;
+
+  /* We have a non-zero denormal.  */
+  assert(omsb < semantics->pr