[libcxx-commits] [libcxx] [libc++] Extend the scope of radix sorting inside std::stable_sort to floating-point types (PR #129452)

[Дмитрий Изволов via libcxx-commits]

================
@@ -298,6 +301,84 @@ _LIBCPP_HIDE_FROM_ABI constexpr auto __shift_to_unsigned(_Ip __n) {
   return static_cast<make_unsigned_t<_Ip> >(__n ^ __min_value);
 }
 
+template <size_t _Size>
+struct __unsigned_integer_of_size {};
+
+template <>
+struct __unsigned_integer_of_size<1> {
+  using type = uint8_t;
+};
+
+template <>
+struct __unsigned_integer_of_size<2> {
+  using type = uint16_t;
+};
+
+template <>
+struct __unsigned_integer_of_size<4> {
+  using type = uint32_t;
+};
+
+template <>
+struct __unsigned_integer_of_size<8> {
+  using type = uint64_t;
+};
+
+template <>
+struct __unsigned_integer_of_size<16> {
+#  if _LIBCPP_HAS_INT128
+  using type = __int128;
+#  endif
+};
+
+template <size_t _Size>
+using __unsigned_integer_of_size_t _LIBCPP_NODEBUG = typename __unsigned_integer_of_size<_Size>::type;
+
+template <class _Sc>
+using __unsigned_representation_for_t _LIBCPP_NODEBUG = __unsigned_integer_of_size_t<sizeof(_Sc)>;
+
+// The function `__to_ordered_integral` is defined for integers and IEEE 754 floating-point numbers.
+// Returns an integer representation such that for any `x` and `y` such that `x < y`, the expression
+// `__to_ordered_integral(x) < __to_ordered_integral(y)` is true, where `x`, `y` are integers or IEEE 754 floats.
+template <class _Integral, enable_if_t< is_integral<_Integral>::value, int> = 0>
+_LIBCPP_HIDE_FROM_ABI constexpr auto __to_ordered_integral(_Integral __n) {
+  return __n;
+}
+
+// An overload for IEEE 754 floating-point numbers
+
+// For the floats conforming to IEEE 754 (IEC 559) standard, we know that:
+// 1. The bit representation of positive floats directly reflects their order:
+//    When comparing floats by magnitude, the number with the larger exponent is greater, and if the exponents are
+//    equal, the one with the larger mantissa is greater.
+// 2. The bit representation of negative floats reflects their reverse order (for the same reasons).
+// 3. The most significant bit (sign bit) is zero for positive floats and one for negative floats. Therefore, in the raw
+//    bit representation, any negative number will be greater than any positive number.
+
+// The only exception from this rule is `NaN`, which is unordered by definition.
+
+// Based on the above, to obtain correctly ordered integral representation of floating-point numbers, we need to:
+// 1. Invert the bit representation (including the sign bit) of negative floats to switch from reverse order to direct
+//    order;
+// 2. Invert the sign bit for positive floats.
+
+// Thus, in final integral representation, we have reversed the order for negative floats and made all negative floats
+// smaller than all positive numbers (by inverting the sign bit).
+template <class _Floating, enable_if_t< numeric_limits<_Floating>::is_iec559, int> = 0>
----------------
izvolov wrote:

Had to tinker with `long double`.  
First, it turned out that `is_iec559` can be true not only for 128-bit `long double` but also for 80-bit ones.  
At the same time, the actual size of an 80-bit `long double` (what `sizeof` returns) can be 10, 12, or 16.  
So, in general, you can’t determine its bit width based on the type’s size.  

Tried two approaches to determine the bit width:  

1. **Using `numeric_limits`:**  
   Tried calculating it as the sum of the mantissa, exponent, and sign bit sizes.  
   The problem here is that for all types except `long double`, the mantissa size includes the hidden bit.  
   For example, `std::numeric_limits<float>::digits == 24`, even though the mantissa is actually stored in 23 bits.  
   But for 80-bit `long double`, this isn’t the case: its mantissa is stored in 64 bits, and `std::numeric_limits<long double>::digits == 64`.  

   As a result, the sum `std::numeric_limits<Float>::digits + bit_log2(static_cast<std::size_t>(std::numeric_limits<Float>::max_exponent) + 1)`  
   gives the correct bit count for all floating-point types except `long double`.  
   And it seems impossible to adjust this formula for the special case of `long double` because, on some architectures, discrepancies might arise again.  

2. **Determining padding size:**  
   Online sources claimed that if an 80-bit `long double` is represented as a 12-byte or 16-byte value,  
   the unused bits are padded with zeros.  
   Decided to use this assumption and compute the bit width based on the padding size.  
   The idea was to take a negative number (its most significant bit—the sign bit—would be `1`) and count the number of leading zero bits after the sign bit.  
   This didn’t work either because I encountered cases where the padding wasn’t necessarily zeros but could contain arbitrary data.  

Also I've encountered a problem that `bit_cast` from `long double` to `uint128` is not working:
```
  # | /__w/llvm-project/llvm-project/build/generic-cxx26/libcxx/test-suite-install/include/c++/v1/__bit/bit_cast.h:26:10: note: indeterminate value can only initialize an object of type 'unsigned char' or 'std::byte'; 'unsigned __int128' is invalid
  # |    26 |   return __builtin_bit_cast(_ToType, __from);
  # |       |          ^
  # | /__w/llvm-project/llvm-project/build/generic-cxx26/libcxx/test-suite-install/include/c++/v1/__algorithm/radix_sort.h:353:29: note: in call to '__bit_cast<unsigned __int128, long double>(-numeric_limits<long double>::infinity())'
  # |   353 |          std::__countl_zero(std::__bit_cast<__integral_type>(-numeric_limits<_Floating>::infinity()));
  # |       |                             ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
```

Came to the conclusion that `long double` should be excluded for now, leaving only `float` and `double`.  
Meanwhile, the new floating-point types from `<stdfloat>` should work "out of the box" since they are defined in compliance with IEEE 754 from the start.  


https://github.com/llvm/llvm-project/pull/129452