[APFloat] Enlarge ExponentType to 32bit integer

[ekatz]

Summary

Enlarge the size of ExponentType from 16bit integer to 32bit. This is required to prevent exponent overflow/underflow.

Note that IEEEFloat size doesn't change in 64bit or 32bit compilation targets.

Fix PR34851.

[arsenm]

llvm/unittests/ADT/APFloatTest.cpp

2612

EXPECT_EQ

2913

EXPECT_EQ

[ekatz]

Changed to EXPECT_EQ as requested.

[efriedma]

llvm/lib/Support/APFloat.cpp

1033

These static_casts don't do anything.

1237

EXPONENT_MIN is -32768, which is more than twice as small as the smallest exponent in any format we support, right? If the exponent is smaller than that, can we "cheat", and pretend the result came out to -32768? The result is essentially zero anyway.

Do we need to handle EXPONENT_MAX?

llvm/unittests/ADT/APFloatTest.cpp

2910

Maybe check what happens with other rounding modes, where the result isn't zero?

[ekatz]

This is a simpler take for the solution.

Also, added tests for other rounding modes, as requested.

[efriedma]

llvm/include/llvm/ADT/APFloat.h

146 ↗

(On Diff #229752)

What's the effect of this on the size of APFloat on a 64-bit target? On a 32-bit target?

llvm/lib/Support/APFloat.cpp

1015

Does this have some functional effect? If it does, is there a test?

[ekatz]

llvm/include/llvm/ADT/APFloat.h

146 ↗

(On Diff #229752)

APFloat remains the same size. I added a note about that in the description. Maybe it is not clear enough.

detail::IEEEFloat is the only record that has a field of ExponentType type. It remains the size of 24 bytes and alignment of 8 bytes (in both 32 bit and 64 bit targets).

llvm/lib/Support/APFloat.cpp

1015

Now that I think of it, I'll just separate it into another patch. It is a bug, revealed from the tests I've made.

[ekatz]

Removed changes unrelated to topic.

[efriedma]

LGTM

llvm/include/llvm/ADT/APFloat.h

146 ↗

(On Diff #229752)

Okay. I missed the comment in the updated commit message.

llvm/lib/Support/APFloat.cpp

 #include "llvm/Config/llvm-config.h"
 #include "llvm/Support/Debug.h"
 #include "llvm/Support/Error.h"
 #include "llvm/Support/MathExtras.h"
 #include "llvm/Support/raw_ostream.h"
 #include <cstring>
 #include <limits.h>
 
+#define EXPONENT_MIN ExponentType(1 << (sizeof(ExponentType) * CHAR_BIT - 1))
+
 #define APFLOAT_DISPATCH_ON_SEMANTICS(METHOD_CALL)                             \
   do {                                                                         \
     if (usesLayout<IEEEFloat>(getSemantics()))                                 \
       return U.IEEE.METHOD_CALL;                                               \
     if (usesLayout<DoubleAPFloat>(getSemantics()))                             \
       return U.Double.METHOD_CALL;                                             \
     llvm_unreachable("Unexpected semantics");                                  \
   } while (false)

 }
 
 /* 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 IEEEFloat::multiplySignificand(const IEEEFloat &rhs,
                                             const IEEEFloat *addend) {
   unsigned int omsb;        // One, not zero, based MSB.
-  unsigned int partsCount, newPartsCount, precision;
+  unsigned int partsCount, newPartsCount, precision, bits;
+  int newExponent;
   integerPart *lhsSignificand;
   integerPart scratch[4];
   integerPart *fullSignificand;
   lostFraction lost_fraction;
   bool ignored;
 
   assert(semantics == rhs.semantics);
 
   precision = semantics->precision;
 
   // Allocate space for twice as many bits as the original significand, plus one
   // extra bit for the addition to overflow into.
   newPartsCount = partCountForBits(precision * 2 + 1);
 
   if (newPartsCount > 4)
     fullSignificand = new integerPart[newPartsCount];
   else
     fullSignificand = scratch;
 
   lhsSignificand = significandParts();
   partsCount = partCount();
 
   APInt::tcFullMultiply(fullSignificand, lhsSignificand,
                         rhs.significandParts(), partsCount, partsCount);
 
   lost_fraction = lfExactlyZero;
   omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
-  exponent += rhs.exponent;
+  newExponent = static_cast<int>(exponent) + static_cast<int>(rhs.exponent);
 
   // Assume the operands involved in the multiplication are single-precision
   // FP, and the two multiplicants are:
   //   *this = a23 . a22 ... a0 * 2^e1
   //     rhs = b23 . b22 ... b0 * 2^e2
   // the result of multiplication is:
   //   *this = c48 c47 c46 . c45 ... c0 * 2^(e1+e2)
   // Note that there are three significant bits at the left-hand side of the
   // radix point: two for the multiplication, and an overflow bit for the
   // addition (that will always be zero at this point). Move the radix point
   // toward left by two bits, and adjust exponent accordingly.
-  exponent += 2;
+  newExponent += 2;
 
   if (addend && addend->isNonZero()) {
     // The intermediate result of the multiplication has "2 * precision"
     // signicant bit; adjust the addend to be consistent with mul result.
     //
     Significand savedSignificand = significand;
     const fltSemantics *savedSemantics = semantics;
     fltSemantics extendedSemantics;
     opStatus status;
     unsigned int extendedPrecision;
 
     // Normalize our MSB to one below the top bit to allow for overflow.
     extendedPrecision = 2 * precision + 1;
     if (omsb != extendedPrecision - 1) {
       assert(extendedPrecision > omsb);
       APInt::tcShiftLeft(fullSignificand, newPartsCount,
                          (extendedPrecision - 1) - omsb);
-      exponent -= (extendedPrecision - 1) - omsb;
+      newExponent -= (extendedPrecision - 1) - omsb;
+    }
+
+    // Make sure the new exponent fits into ExponentType.
+    if (newExponent < EXPONENT_MIN) {
+      unsigned int significantParts;
+      lostFraction lf;
+
+      significantParts = partCountForBits(omsb);
+      lf = shiftRight(fullSignificand, significantParts,
+                      int(EXPONENT_MIN) - newExponent);
+      lost_fraction = combineLostFractions(lf, lost_fraction);
+      newExponent = EXPONENT_MIN;
     }
+    exponent = static_cast<ExponentType>(newExponent);
 
     /* Create new semantics.  */
     extendedSemantics = *semantics;
     extendedSemantics.precision = extendedPrecision;
 
     if (newPartsCount == 1)
       significand.part = fullSignificand[0];
     else
if (addend && addend->isNonZero()) {
     lost_fraction = addOrSubtractSignificand(extendedAddend, false);
 
     /* Restore our state.  */
     if (newPartsCount == 1)
       fullSignificand[0] = significand.part;
     significand = savedSignificand;
     semantics = savedSemantics;
 
+    newExponent = exponent;
     omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
   }
 
   // Convert the result having "2 * precision" significant-bits back to the one
   // having "precision" significant-bits. First, move the radix point from
   // poision "2*precision - 1" to "precision - 1". The exponent need to be
   // adjusted by "2*precision - 1" - "precision - 1" = "precision".
-  exponent -= precision + 1;
+  newExponent -= precision + 1;
 
   // In case MSB resides at the left-hand side of radix point, shift the
   // mantissa right by some amount to make sure the MSB reside right before
   // the radix point (i.e. "MSB . rest-significant-bits").
   //
   // Note that the result is not normalized when "omsb < precision". So, the
   // caller needs to call IEEEFloat::normalize() if normalized value is
   // expected.
-  if (omsb > precision) {
-    unsigned int bits, significantParts;
+  bits = 0;
+  if (omsb > precision)
+    bits += omsb - precision;
+
+  // Make sure the new exponent fits into ExponentType.
+  if (newExponent < EXPONENT_MIN) {
+    bits += int(EXPONENT_MIN) - newExponent;
+    newExponent = EXPONENT_MIN;
+  }
+
+  if (bits) {
+    unsigned int significantParts;
     lostFraction lf;
 
-    bits = omsb - precision;
     significantParts = partCountForBits(omsb);
     lf = shiftRight(fullSignificand, significantParts, bits);
     lost_fraction = combineLostFractions(lf, lost_fraction);
-    exponent += bits;
+    newExponent += bits;
   }
 
+  exponent = static_cast<ExponentType>(newExponent);
   APInt::tcAssign(lhsSignificand, fullSignificand, partsCount);
 
   if (newPartsCount > 4)
     delete [] fullSignificand;
 
   return lost_fraction;
 }
 
 /* Multiply the significands of LHS and RHS to DST.  */
 lostFraction IEEEFloat::divideSignificand(const IEEEFloat &rhs) {
   unsigned int bit, i, partsCount;
+  int newExponent;
   const integerPart *rhsSignificand;
   integerPart *lhsSignificand, *dividend, *divisor;
   integerPart scratch[4];
   lostFraction lost_fraction;
 
   assert(semantics == rhs.semantics);
 
   lhsSignificand = significandParts();
lostFraction IEEEFloat::divideSignificand(const IEEEFloat &rhs) {
 
   /* 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;
+  newExponent = exponent - rhs.exponent;
 
   unsigned int precision = semantics->precision;
 
   /* Normalize the divisor.  */
   bit = precision - APInt::tcMSB(divisor, partsCount) - 1;
   if (bit) {
-    exponent += bit;
+    newExponent += bit;
     APInt::tcShiftLeft(divisor, partsCount, bit);
   }
 
   /* Normalize the dividend.  */
   bit = precision - APInt::tcMSB(dividend, partsCount) - 1;
   if (bit) {
-    exponent -= bit;
+    newExponent -= bit;
     APInt::tcShiftLeft(dividend, partsCount, bit);
   }
 
   /* Ensure the dividend >= divisor initially for the loop below.
      Incidentally, this means that the division loop below is
      guaranteed to set the integer bit to one.  */
   if (APInt::tcCompare(dividend, divisor, partsCount) < 0) {
-    exponent--;
+    newExponent--;
     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);
lostFraction IEEEFloat::divideSignificand(const IEEEFloat &rhs) {
   else if (APInt::tcIsZero(dividend, partsCount))
     lost_fraction = lfExactlyZero;
   else
     lost_fraction = lfLessThanHalf;
 
   if (partsCount > 2)
     delete [] dividend;
 
+  // Make sure the new exponent fits into ExponentType.
+  if (newExponent < EXPONENT_MIN) {
+    unsigned int omsb, significantParts;
+    lostFraction lf;
+
+    omsb = APInt::tcMSB(lhsSignificand, partsCount) + 1;
+    significantParts = partCountForBits(omsb);
+    lf = shiftRight(lhsSignificand, significantParts,
+                    int(EXPONENT_MIN) - newExponent);
+    lost_fraction = combineLostFractions(lf, lost_fraction);
+    newExponent = EXPONENT_MIN;
+  }
+  exponent = static_cast<ExponentType>(newExponent);
+
   return lost_fraction;
 }
 
 unsigned int IEEEFloat::significandMSB() const {
   return APInt::tcMSB(significandParts(), partCount());
 }
 
 unsigned int IEEEFloat::significandLSB() const {

llvm/unittests/ADT/APFloatTest.cpp

for (size_t i = 0; i < NumTests; ++i) {
     APFloat::opStatus status = x.multiply(y, APFloat::rmNearestTiesToEven);
 
     APFloat result(APFloat::IEEEsingle(), SpecialCaseTests[i].result);
 
     EXPECT_TRUE(result.bitwiseIsEqual(x));
     EXPECT_TRUE((int)status == SpecialCaseTests[i].status);
     EXPECT_TRUE((int)x.getCategory() == SpecialCaseTests[i].category);
   }
+
+  // Test exponent underflow.
+  {
+    APFloat x(APFloat::x87DoubleExtended(), "3.645199531882474602528e-4951");
+    APFloat::opStatus status = x.multiply(x, APFloat::rmNearestTiesToEven);
+
+    APFloat zero = APFloat::getZero(APFloat::x87DoubleExtended());
+
+    EXPECT_TRUE(x.bitwiseIsEqual(zero));
+    EXPECT_TRUE((int)status == UnderflowStatus);
+  }
 }
 
 TEST(APFloatTest, divide) {
   // Test Special Cases against each other and normal values.
 
   // TODOS/NOTES:
   // 1. Since we perform only default exception handling all operations with
   // signaling NaNs should have a result that is a quiet NaN. Currently they
for (size_t i = 0; i < NumTests; ++i) {
     APFloat::opStatus status = x.divide(y, APFloat::rmNearestTiesToEven);
 
     APFloat result(APFloat::IEEEsingle(), SpecialCaseTests[i].result);
 
     EXPECT_TRUE(result.bitwiseIsEqual(x));
     EXPECT_TRUE((int)status == SpecialCaseTests[i].status);
     EXPECT_TRUE((int)x.getCategory() == SpecialCaseTests[i].category);
   }
+
+  // Test exponent underflow.
+  {
+    APFloat x(APFloat::x87DoubleExtended(), "3.645199531882474602528e-4951");
+    APFloat y(APFloat::x87DoubleExtended(), "1.18973149535723176502e+4932");
+    APFloat::opStatus status = x.divide(y, APFloat::rmNearestTiesToEven);
+
+    APFloat zero = APFloat::getZero(APFloat::x87DoubleExtended());
+
+    EXPECT_TRUE(x.bitwiseIsEqual(zero));
+    EXPECT_TRUE((int)status == UnderflowStatus);
+  }
 }
 
 TEST(APFloatTest, operatorOverloads) {
   // This is mostly testing that these operator overloads compile.
   APFloat One = APFloat(APFloat::IEEEsingle(), "0x1p+0");
   APFloat Two = APFloat(APFloat::IEEEsingle(), "0x2p+0");
   EXPECT_TRUE(Two.bitwiseIsEqual(One + One));
   EXPECT_TRUE(One.bitwiseIsEqual(Two - One));