[libc] Correct rounding for hexadecimalStringToFloat with long inputs.

[lntue]

Summary

Update binaryExpTofloat so that it will round correctly for long inputs when converting hexadecimal strings to floating points.

[michaelrj]

libc/src/__support/str_to_float.h

522

This is unnecessary, since the FloatBits type will do it when setting the mantissa of the final value. It also doesn't really hurt anything to have it, but I generally prefer to get the full output bits for testing.

662–665

setting seenDigit to false isn't recommended, since it should just be a flag for if any digit has been seen. This particular if would cause numbers that end with alphanumeric characters that aren't in base to return as zero, for example 0x1z.

[lntue]

[libc] Correct rounding for hexadecimalStringToFloat with long inputs.

libc/src/__support/str_to_float.h

522

I think it is better for maintainability and consistency to have the outputs well-defined (as is expected results bits in the answer) instead of relying on the caller to do the cleanup without being documented. For instance, without this clean up, INF results might have outputMantissa == 0, 2^mantissaLength, or 2^(mantissaLength+1).

[michaelrj]

LGTM

libc/src/__support/str_to_float.h

522

That seems reasonable, in that case I should probably go fix the other functions to match, although that will require fixing a lot of tests.

[lntue]

[libc] Correct rounding for hexadecimalStringToFloat with long inputs.

libc/src/__support/str_to_float.h

 #include "src/__support/detailed_powers_of_ten.h"
 #include "src/__support/high_precision_decimal.h"
 #include "src/__support/str_to_integer.h"
 #include <errno.h>
 
 namespace __llvm_libc {
 namespace internal {
 
-// Shifts right and rounds according to the following rules:
-// 1) If the part being cut off is more than 2^(amountToShift - 1) then round
-// up
-// 2) If it is less than that number then round down
-// 3) If it is exactly that number, then round so that the final number will be
-// even
-template <class T>
-static inline T shiftRightAndRound(T numToShift, unsigned int amountToShift) {
-  T result = numToShift >> amountToShift;
-  T truncated = numToShift & ((1 << amountToShift) - 1);
-
-  if (truncated < (1 << (amountToShift - 1))) {
-    return result;
-  } else if (truncated > (1 << (amountToShift - 1))) {
-    return result + 1;
-  } else {
-    return result + (result & 1); // This rounds towards even.
-  }
-}
-
 template <class T> uint32_t inline leadingZeroes(T inputNumber) {
   // TODO(michaelrj): investigate the portability of using something like
   // __builtin_clz for specific types.
   constexpr uint32_t bitsInT = sizeof(T) * 8;
   if (inputNumber == 0) {
     return bitsInT;
   }
   uint32_t curGuess = bitsInT / 2;

 
 // Takes a mantissa and base 2 exponent and converts it into its closest
 // floating point type T equivalient. Since the exponent is already in the right
 // form, this is mostly just shifting and rounding. This is used for hexadecimal
 // numbers since a base 16 exponent multiplied by 4 is the base 2 exponent.
 template <class T>
 static inline void
 binaryExpToFloat(typename fputil::FPBits<T>::UIntType mantissa, int32_t exp2,
+                 bool truncated,
                  typename fputil::FPBits<T>::UIntType *outputMantissa,
                  uint32_t *outputExp2) {
   using BitsType = typename fputil::FPBits<T>::UIntType;
 
   // This is the number of leading zeroes a properly normalized float of type T
   // should have.
-  constexpr int32_t NORMALIZED_LEADING_ZEROES =
-      (sizeof(BitsType) * 8) - fputil::FloatProperties<T>::mantissaWidth - 1;
-  constexpr BitsType OVERFLOWED_MANTISSA =
-      BitsType(1) << (fputil::FloatProperties<T>::mantissaWidth + 1);
+  constexpr int32_t NUMBITS = sizeof(BitsType) * 8;
+  constexpr int32_t INF_EXP =
+      (1 << fputil::FloatProperties<T>::exponentWidth) - 1;
+
+  // Normalization step 1: Bring the leading bit to the highest bit of BitsType.
+  uint32_t amountToShiftLeft = leadingZeroes<BitsType>(mantissa);
+  mantissa <<= amountToShiftLeft;
 
-  // Normalization
-  int32_t amountToShift =
-      NORMALIZED_LEADING_ZEROES -
-      static_cast<int32_t>(leadingZeroes<BitsType>(mantissa));
-  if (amountToShift < 0) {
-    mantissa <<= -amountToShift;
-  } else {
-    mantissa = shiftRightAndRound(mantissa, amountToShift);
-    if (mantissa == OVERFLOWED_MANTISSA) {
-      mantissa >>= 1;
-      exp2 += 1;
-    }
-  }
-  exp2 += amountToShift;
-
-  // Account for the fact that the mantissa represented an integer
-  // previously, but now represents the fractional part of a normalized
-  // number.
-  exp2 += fputil::FloatProperties<T>::mantissaWidth;
+  // Keep exp2 representing the exponent of the lowest bit of BitsType.
+  exp2 -= amountToShiftLeft;
 
-  int32_t biasedExponent = exp2 + fputil::FPBits<T>::exponentBias;
-  // handle subnormals
-  if (biasedExponent <= 0) {
+  // biasedExponent represents the biased exponent of the most significant bit.
+  int32_t biasedExponent = exp2 + NUMBITS + fputil::FPBits<T>::exponentBias - 1;
 
-    // the most mantissa is currently normalized, meaning that the msb is
-    // one bit left of where the decimal point should go.
-    amountToShift = 1;
-    BitsType mantissaCopy = mantissa >> 1;
-    while (biasedExponent < 0 && mantissaCopy > 0) {
-      mantissaCopy = mantissaCopy >> 1;
-      ++amountToShift;
-      ++biasedExponent;
-    }
-    // If we cut off any bits to fit this number into a subnormal, then it's
-    // out of range for this size of float.
-    if ((mantissa & ((1 << amountToShift) - 1)) > 0) {
-      errno = ERANGE; // NOLINT
-    }
-    mantissa = shiftRightAndRound(mantissa, amountToShift);
-    if (mantissa == OVERFLOWED_MANTISSA) {
-      mantissa >>= 1;
-      exp2 += 1;
-    } else if (mantissa == 0) {
-      biasedExponent = 0;
+  // Handle numbers that're too large and get squashed to inf
+  if (biasedExponent >= INF_EXP) {
+    // This indicates an overflow, so we make the result INF and set errno.
+    *outputExp2 = (1 << fputil::FloatProperties<T>::exponentWidth) - 1;
+    *outputMantissa = 0;
+    return;
+
+  uint32_t amountToShiftRight =
+      NUMBITS - fputil::FloatProperties<T>::mantissaWidth - 1;
+
+  // Handle subnormals.
+  if (biasedExponent <= 0) {
+    amountToShiftRight += 1 - biasedExponent;
+
+    if (amountToShiftRight > NUMBITS) {
+      // Return 0 if the exponent is too small.
+      *outputMantissa = 0;
+      *outputExp2 = 0;
+      errno = ERANGE; // NOLINT
+      return;
     }
   }
-  // handle numbers that're too large and get squashed to inf
-  else if (biasedExponent >
-           (1 << fputil::FloatProperties<T>::exponentWidth) - 1) {
-    // This indicates an overflow, so we make the result INF and set errno.
-    biasedExponent = (1 << fputil::FloatProperties<T>::exponentWidth) - 1;
+
+  BitsType roundBitMask = BitsType(1) << (amountToShiftRight - 1);
+  BitsType stickyMask = roundBitMask - 1;
+  bool roundBit = mantissa & roundBitMask;
+  bool stickyBit = static_cast<bool>(mantissa & stickyMask) || truncated;
+
+  if (amountToShiftRight < NUMBITS) {
+    // Shift the mantissa and clear the implicit bit.
+    mantissa >>= amountToShiftRight;
+    mantissa &= fputil::FloatProperties<T>::mantissaMask;
+  } else {
     mantissa = 0;
+  }
+  bool leastSignificantBit = mantissa & BitsType(1);
+  // Perform rounding-to-nearest, tie-to-even.
+  if (roundBit && (leastSignificantBit || stickyBit)) {
+    ++mantissa;
+  }
+
+  if (mantissa > fputil::FloatProperties<T>::mantissaMask) {
+    // Rounding causes the exponent to increase.
+    ++biasedExponent;
+
+    if (biasedExponent == INF_EXP) {
       errno = ERANGE; // NOLINT
     }
-  *outputMantissa = mantissa;
+  }
+
+  if (biasedExponent == 0) {
+    errno = ERANGE; // NOLINT
+  }
+
+  *outputMantissa = mantissa & fputil::FloatProperties<T>::mantissaMask;
   *outputExp2 = biasedExponent;
 }
 
 // checks if the next 4 characters of the string pointer are the start of a
 // hexadecimal floating point number. Does not advance the string pointer.
 static inline bool is_float_hex_start(const char *__restrict src,
                                       const char decimalPoint) {
   if (!(*src == '0' && (*(src + 1) | 32) == 'x')) {
hexadecimalStringToFloat(const char *__restrict src, const char DECIMAL_POINT,
   // the format mantissa * (base ^ exponent)
 
   // The loop fills the mantissa with as many digits as it can hold
   const BitsType BITSTYPE_MAX_DIV_BY_BASE =
       __llvm_libc::cpp::NumericLimits<BitsType>::max() / BASE;
   while (true) {
     if (isalnum(*src)) {
       uint32_t digit = b36_char_to_int(*src);
-      if (digit >= BASE) {
-        seenDigit = false;
-        break;
-      }
+      if (digit < BASE)
         seenDigit = true;
+      else
+        break;
 
       if (mantissa < BITSTYPE_MAX_DIV_BY_BASE) {
         mantissa = (mantissa * BASE) + digit;
         if (afterDecimal)
           --exponent;
       } else {
         if (digit > 0)
           truncated = true;
if (*(src + 1) == '+' || *(src + 1) == '-' || isdigit(*(src + 1))) {
       exponent += add_to_exponent;
     }
   }
   *strEnd = const_cast<char *>(src);
   if (mantissa == 0) { // if we have a 0, then also 0 the exponent.
     *outputMantissa = 0;
     *outputExponent = 0;
   } else {
-    binaryExpToFloat<T>(mantissa, exponent, outputMantissa, outputExponent);
+    binaryExpToFloat<T>(mantissa, exponent, truncated, outputMantissa,
+                        outputExponent);
   }
   return true;
 }
 
 // Takes a pointer to a string and a pointer to a string pointer. This function
 // is used as the backend for all of the string to float functions.
 template <class T>
 static inline T strtofloatingpoint(const char *__restrict src,

libc/test/src/stdlib/strtof_test.cpp

 TEST_F(LlvmLibcStrToFTest, BasicHexadecimalTests) {
   runTest("0x1", 3, 0x3f800000);
   runTest("0x10", 4, 0x41800000);
   runTest("0x11", 4, 0x41880000);
   runTest("0x0.1234", 8, 0x3d91a000);
 }
 
 TEST_F(LlvmLibcStrToFTest, HexadecimalSubnormalTests) {
-  runTest("0x0.0000000000000000000000000000000002", 38, 0x4000);
+  runTest("0x0.0000000000000000000000000000000002", 38, 0x4000, ERANGE);
 
   // This is the largest subnormal number as represented in hex
-  runTest("0x0.00000000000000000000000000000003fffff8", 42, 0x7fffff);
+  runTest("0x0.00000000000000000000000000000003fffff8", 42, 0x7fffff, ERANGE);
 }
 
 TEST_F(LlvmLibcStrToFTest, HexadecimalSubnormalRoundingTests) {
   // This is the largest subnormal number that gets rounded down to 0 (as a
   // float)
   runTest("0x0.00000000000000000000000000000000000004", 42, 0x0, ERANGE);
 
   // This is slightly larger, and thus rounded up
   runTest("0x0.000000000000000000000000000000000000041", 43, 0x00000001,
           ERANGE);
 
   // These check that we're rounding to even properly
   runTest("0x0.0000000000000000000000000000000000000b", 42, 0x00000001, ERANGE);
   runTest("0x0.0000000000000000000000000000000000000c", 42, 0x00000002, ERANGE);
+
+  // These check that we're rounding to even properly even when the input bits
+  // are longer than the bit fields can contain.
+  runTest("0x1.000000000000000000000p-150", 30, 0x00000000, ERANGE);
+  runTest("0x1.000010000000000001000p-150", 30, 0x00000001, ERANGE);
+  runTest("0x1.000100000000000001000p-134", 30, 0x00008001, ERANGE);
+  runTest("0x1.FFFFFC000000000001000p-127", 30, 0x007FFFFF, ERANGE);
+  runTest("0x1.FFFFFE000000000000000p-127", 30, 0x00800000);
 }
 
 TEST_F(LlvmLibcStrToFTest, HexadecimalNormalRoundingTests) {
   // This also checks the round to even behavior by checking three adjacent
   // numbers.
   // This gets rounded down to even
   runTest("0x123456500", 11, 0x4f91a2b2);
   // This doesn't get rounded at all
   runTest("0x123456600", 11, 0x4f91a2b3);
   // This gets rounded up to even
   runTest("0x123456700", 11, 0x4f91a2b4);
+  // Correct rounding for long input
+  runTest("0x1.000001000000000000000", 25, 0x3f800000);
+  runTest("0x1.000001000000000000100", 25, 0x3f800001);
 }
 
 TEST_F(LlvmLibcStrToFTest, HexadecimalsWithRoundingProblems) {
   runTest("0xFFFFFFFF", 10, 0x4f800000);
 }
 
 TEST_F(LlvmLibcStrToFTest, HexadecimalOutOfRangeTests) {
   runTest("0x123456789123456789123456789123456789", 38, 0x7f800000, ERANGE);