[PATCH 1/4] libstdc++: Factor out uses of __int128 into a type alias

[PATCH 1/4] libstdc++: Factor out uses of __int128 into a type alias

[Patrick Palka]

Since Ryu has an alias uint128_t for this same purpose, it seems best
for us to use this name as well, so as to minimize the amount of local
modifications we'd need to make to our copy of Ryu.  (In a subsequent
patch, we're going to remove Ryu's aliases so that it uses the one
defined in floating_to_chars.cc.)

libstdc++-v3/ChangeLog:

	* src/c++17/floating_to_chars.cc (uint128_t): New conditionally
	defined alias of unsigned __int128.
	(floating_type_traits_binary128::mantissa_t): Use uint128_t
	instead of unsigned __int128.
	[LONG_DOUBLE_KIND == LDK_IBM128]
	(floating_type_traits<long double>::mantissa_t): Likewise.
	(get_ieee_repr): Likewise.  Make casts from uint_t to mantissa_t
	and uint32_t explicit.  Simplify the extraction of mantissa,
	exponent and sign bit.
---
 libstdc++-v3/src/c++17/floating_to_chars.cc | 25 +++++++++++++--------
 1 file changed, 16 insertions(+), 9 deletions(-)

diff --git a/libstdc++-v3/src/c++17/floating_to_chars.cc b/libstdc++-v3/src/c++17/floating_to_chars.cc
index 611747bb99e..da3fbaa1ed1 100644
--- a/libstdc++-v3/src/c++17/floating_to_chars.cc
+++ b/libstdc++-v3/src/c++17/floating_to_chars.cc
@@ -98,6 +98,10 @@ using F128_type = void;
 
 namespace
 {
+#if defined __SIZEOF_INT128__
+  using uint128_t = unsigned __int128;
+#endif
+
   namespace ryu
   {
 #include "ryu/common.h"
@@ -171,7 +175,7 @@ namespace
     static constexpr int mantissa_bits = 112;
     static constexpr int exponent_bits = 15;
     static constexpr bool has_implicit_leading_bit = true;
-    using mantissa_t = unsigned __int128;
+    using mantissa_t = uint128_t;
     using shortest_scientific_t = ryu::floating_decimal_128;
 
     static constexpr uint64_t pow10_adjustment_tab[]
@@ -367,7 +371,7 @@ namespace
       static constexpr int mantissa_bits = 105;
       static constexpr int exponent_bits = 11;
       static constexpr bool has_implicit_leading_bit = true;
-      using mantissa_t = unsigned __int128;
+      using mantissa_t = uint128_t;
       using shortest_scientific_t = ryu::floating_decimal_128;
 
       static constexpr uint64_t pow10_adjustment_tab[]
@@ -393,6 +397,7 @@ namespace
     ieee_t<T>
     get_ieee_repr(const T value)
     {
+      using mantissa_t = typename floating_type_traits<T>::mantissa_t;
       constexpr int mantissa_bits = floating_type_traits<T>::mantissa_bits;
       constexpr int exponent_bits = floating_type_traits<T>::exponent_bits;
       constexpr int total_bits = mantissa_bits + exponent_bits + 1;
@@ -404,7 +409,7 @@ namespace
 	  return uint64_t{};
 #ifdef __SIZEOF_INT128__
 	else if constexpr (total_bits <= 128)
-	  return (unsigned __int128){};
+	  return uint128_t{};
 #endif
       };
       using uint_t = decltype(get_uint_t());
@@ -412,10 +417,13 @@ namespace
       memcpy(&value_bits, &value, sizeof(value));
 
       ieee_t<T> ieee_repr;
-      ieee_repr.mantissa = value_bits & ((uint_t{1} << mantissa_bits) - 1u);
+      ieee_repr.mantissa
+	= static_cast<mantissa_t>(value_bits & ((uint_t{1} << mantissa_bits) - 1u));
+      value_bits >>= mantissa_bits;
       ieee_repr.biased_exponent
-	= (value_bits >> mantissa_bits) & ((uint_t{1} << exponent_bits) - 1u);
-      ieee_repr.sign = (value_bits >> (mantissa_bits + exponent_bits)) & 1;
+	= static_cast<uint32_t>(value_bits & ((uint_t{1} << exponent_bits) - 1u));
+      value_bits >>= exponent_bits;
+      ieee_repr.sign = (value_bits & 1) != 0;
       return ieee_repr;
     }
 
@@ -430,7 +438,6 @@ namespace
       // mantissa (plus an implicit leading bit).  We use the exponent and sign
       // of the high part, and we merge the mantissa of the high part with the
       // mantissa (and the implicit leading bit) of the low part.
-      using uint_t = unsigned __int128;
       uint64_t value_bits[2] = {};
       memcpy(value_bits, &value, sizeof(value_bits));
 
@@ -478,8 +485,8 @@ namespace
 	}
 
       ieee_t<long double> ieee_repr;
-      ieee_repr.mantissa = ((uint_t{mantissa_hi} << 64)
-			    | (uint_t{mantissa_lo} << 4)) >> 11;
+      ieee_repr.mantissa = ((uint128_t{mantissa_hi} << 64)
+			    | (uint128_t{mantissa_lo} << 4)) >> 11;
       ieee_repr.biased_exponent = exponent_hi;
       ieee_repr.sign = sign_hi;
       return ieee_repr;
-- 
2.31.0.rc2

[PATCH 2/4] libstdc++: Add LOCAL_PATCHES file to Ryu sources

[Patrick Palka]

This file keeps track of the local modifications we've made to our copy
of Ryu.

libstdc++-v3/ChangeLog:

	* src/c++17/ryu/LOCAL_PATCHES: New file.
---
 libstdc++-v3/src/c++17/ryu/LOCAL_PATCHES | 1 +
 1 file changed, 1 insertion(+)
 create mode 100644 libstdc++-v3/src/c++17/ryu/LOCAL_PATCHES

diff --git a/libstdc++-v3/src/c++17/ryu/LOCAL_PATCHES b/libstdc++-v3/src/c++17/ryu/LOCAL_PATCHES
new file mode 100644
index 00000000000..c26633ee6c2
--- /dev/null
+++ b/libstdc++-v3/src/c++17/ryu/LOCAL_PATCHES
@@ -0,0 +1 @@
+r11-6248
-- 
2.31.0.rc2

[PATCH 3/4] libstdc++: Remove Ryu's uint128_t aliases

[Patrick Palka]

This makes Ryu consistently use the uint128_t alias defined in
floating_to_chars.cc.

libstdc++-v3/ChangeLog:

	* src/c++17/ryu/LOCAL_PATCHES: Update.
	* src/c++17/ryu/d2s_intrinsics.h: Don't define uint128_t.
	* src/c++17/ryu/generic_128.h: Likewise.
	* src/c++17/ryu/ryu_generic_128.h (struct floating_decimal_128):
	Use uint128_t instead of __uint128_t.
	(generic_binary_to_decimal): Likewise.
---
 libstdc++-v3/src/c++17/ryu/LOCAL_PATCHES     | 1 +
 libstdc++-v3/src/c++17/ryu/d2s_intrinsics.h  | 4 ----
 libstdc++-v3/src/c++17/ryu/generic_128.h     | 3 ---
 libstdc++-v3/src/c++17/ryu/ryu_generic_128.h | 4 ++--
 4 files changed, 3 insertions(+), 9 deletions(-)

diff --git a/libstdc++-v3/src/c++17/ryu/LOCAL_PATCHES b/libstdc++-v3/src/c++17/ryu/LOCAL_PATCHES
index c26633ee6c2..a72168c1f71 100644
--- a/libstdc++-v3/src/c++17/ryu/LOCAL_PATCHES
+++ b/libstdc++-v3/src/c++17/ryu/LOCAL_PATCHES
@@ -1 +1,2 @@
 r11-6248
+r11-????
diff --git a/libstdc++-v3/src/c++17/ryu/d2s_intrinsics.h b/libstdc++-v3/src/c++17/ryu/d2s_intrinsics.h
index fa993e6fad6..bbac4dfd48f 100644
--- a/libstdc++-v3/src/c++17/ryu/d2s_intrinsics.h
+++ b/libstdc++-v3/src/c++17/ryu/d2s_intrinsics.h
@@ -28,10 +28,6 @@
 #define HAS_64_BIT_INTRINSICS
 #endif
 
-#if defined(HAS_UINT128)
-typedef __uint128_t uint128_t;
-#endif
-
 #if defined(HAS_64_BIT_INTRINSICS)
 
 
diff --git a/libstdc++-v3/src/c++17/ryu/generic_128.h b/libstdc++-v3/src/c++17/ryu/generic_128.h
index 88e96776664..d3ca398fbdb 100644
--- a/libstdc++-v3/src/c++17/ryu/generic_128.h
+++ b/libstdc++-v3/src/c++17/ryu/generic_128.h
@@ -17,9 +17,6 @@
 #ifndef RYU_GENERIC128_H
 #define RYU_GENERIC128_H
 
-
-typedef __uint128_t uint128_t;
-
 #define FLOAT_128_POW5_INV_BITCOUNT 249
 #define FLOAT_128_POW5_BITCOUNT 249
 #define POW5_TABLE_SIZE 56
diff --git a/libstdc++-v3/src/c++17/ryu/ryu_generic_128.h b/libstdc++-v3/src/c++17/ryu/ryu_generic_128.h
index f5d476343b6..2afbf274e11 100644
--- a/libstdc++-v3/src/c++17/ryu/ryu_generic_128.h
+++ b/libstdc++-v3/src/c++17/ryu/ryu_generic_128.h
@@ -37,7 +37,7 @@ extern "C" {
 
 // A floating decimal representing (-1)^s * m * 10^e.
 struct floating_decimal_128 {
-  __uint128_t mantissa;
+  uint128_t mantissa;
   int32_t exponent;
   bool sign;
 };
@@ -52,7 +52,7 @@ struct floating_decimal_128 long_double_to_fd128(long double d);
 // Converts the given binary floating point number to the shortest decimal floating point number
 // that still accurately represents it.
 struct floating_decimal_128 generic_binary_to_decimal(
-    const __uint128_t bits, const uint32_t mantissaBits, const uint32_t exponentBits, const bool explicitLeadingBit);
+    const uint128_t bits, const uint32_t mantissaBits, const uint32_t exponentBits, const bool explicitLeadingBit);
 
 // Converts the given decimal floating point number to a string, writing to result, and returning
 // the number characters written. Does not terminate the buffer with a 0. In the worst case, this
-- 
2.31.0.rc2

[PATCH 4/4] libstdc++: Add fallback 128-bit integer class type and use it

[Patrick Palka]

This implements a minimal integer class type that emulates 128-bit
unsigned arithmetic using a pair of 64-bit integers, which the
floating-point std::to_chars implementation then uses as a drop-in
replacement for unsigned __int128 on targets that lack the latter.
This allows us to fully support formatting of large long double types
on targets that lack __int128.

Since Ryu performs 128-bit division/modulus only by 2, 5 and 10, the
integer class type supports only these divisors rather than supporting
general division/modulus.

Tested on x86, x86_64, ppc64le, ppc64be and aarch64, with and without
performing the equivalent of -U__SIZEOF_INT128__ in floating_to_chars.cc
(so that we also test using the class type on targets when __int128 is
available).

libstdc++-v3/ChangeLog:

	* src/c++17/floating_to_chars.cc: Simplify the file as if
	__SIZEOF_INT128__ is always defined.
	[!defined __SIZEOF_INT128__]: Include "uint128_t.h".  Define
	a to_chars overload for the uint128_t class type.
	* src/c++17/uint128_t.h: New file.
	* testsuite/20_util/to_chars/long_double.cc: No longer expect an
	execution FAIL on targets that have a large long double type
	but lack __int128.
---
 libstdc++-v3/src/c++17/floating_to_chars.cc   |  58 ++--
 libstdc++-v3/src/c++17/uint128_t.h            | 297 ++++++++++++++++++
 .../testsuite/20_util/to_chars/long_double.cc |   1 -
 3 files changed, 332 insertions(+), 24 deletions(-)
 create mode 100644 libstdc++-v3/src/c++17/uint128_t.h

diff --git a/libstdc++-v3/src/c++17/floating_to_chars.cc b/libstdc++-v3/src/c++17/floating_to_chars.cc
index da3fbaa1ed1..86f4401e134 100644
--- a/libstdc++-v3/src/c++17/floating_to_chars.cc
+++ b/libstdc++-v3/src/c++17/floating_to_chars.cc
@@ -64,25 +64,19 @@ extern "C" int __sprintfieee128(char*, const char*, ...);
 
 #if __LDBL_MANT_DIG__ == __DBL_MANT_DIG__
 # define LONG_DOUBLE_KIND LDK_BINARY64
-#elif defined(__SIZEOF_INT128__)
-// The Ryu routines need a 128-bit integer type in order to do shortest
-// formatting of types larger than 64-bit double, so without __int128 we can't
-// support any large long double format.  This is the case for e.g. i386.
-# if __LDBL_MANT_DIG__ == 64
+#elif __LDBL_MANT_DIG__ == 64
 #  define LONG_DOUBLE_KIND LDK_FLOAT80
-# elif __LDBL_MANT_DIG__ == 113
-#  define LONG_DOUBLE_KIND LDK_BINARY128
-# elif __LDBL_MANT_DIG__ == 106
-#  define LONG_DOUBLE_KIND LDK_IBM128
-# endif
-# if defined _GLIBCXX_USE_FLOAT128 && __FLT128_MANT_DIG__ == 113
-// Define overloads of std::to_chars for __float128.
-#  define FLOAT128_TO_CHARS 1
-# endif
+#elif __LDBL_MANT_DIG__ == 113
+# define LONG_DOUBLE_KIND LDK_BINARY128
+#elif __LDBL_MANT_DIG__ == 106
+# define LONG_DOUBLE_KIND LDK_IBM128
+#else
+# define LONG_DOUBLE_KIND LDK_UNSUPPORTED
 #endif
 
-#if !defined(LONG_DOUBLE_KIND)
-# define LONG_DOUBLE_KIND LDK_UNSUPPORTED
+#if defined _GLIBCXX_USE_FLOAT128 && __FLT128_MANT_DIG__ == 113
+// Define overloads of std::to_chars for __float128.
+# define FLOAT128_TO_CHARS 1
 #endif
 
 // For now we only support __float128 when it's the powerpc64 __ieee128 type.
@@ -100,6 +94,8 @@ namespace
 {
 #if defined __SIZEOF_INT128__
   using uint128_t = unsigned __int128;
+#else
+# include "uint128_t.h"
 #endif
 
   namespace ryu
@@ -114,7 +110,6 @@ namespace
 #include "ryu/d2fixed.c"
 #include "ryu/f2s.c"
 
-#ifdef __SIZEOF_INT128__
     namespace generic128
     {
       // Put the generic Ryu bits in their own namespace to avoid name conflicts.
@@ -129,7 +124,6 @@ namespace
     int
     to_chars(const floating_decimal_128 v, char* const result)
     { return generic128::generic_to_chars(v, result); }
-#endif
   } // namespace ryu
 
   // A traits class that contains pertinent information about the binary
@@ -407,10 +401,8 @@ namespace
 	  return uint32_t{};
 	else if constexpr (total_bits <= 64)
 	  return uint64_t{};
-#ifdef __SIZEOF_INT128__
 	else if constexpr (total_bits <= 128)
 	  return uint128_t{};
-#endif
       };
       using uint_t = decltype(get_uint_t());
       uint_t value_bits = 0;
@@ -503,7 +495,6 @@ namespace
 	return ryu::floating_to_fd32(value);
       else if constexpr (std::is_same_v<T, double>)
 	return ryu::floating_to_fd64(value);
-#ifdef __SIZEOF_INT128__
       else if constexpr (std::is_same_v<T, long double>
 			 || std::is_same_v<T, F128_type>)
 	{
@@ -519,7 +510,6 @@ namespace
 						mantissa_bits, exponent_bits,
 						!has_implicit_leading_bit);
 	}
-#endif
     }
 
   // This subroutine returns true if the shortest scientific form fd is a
@@ -558,10 +548,32 @@ namespace
   get_mantissa_length(const ryu::floating_decimal_64 fd)
   { return ryu::decimalLength17(fd.mantissa); }
 
-#ifdef __SIZEOF_INT128__
   int
   get_mantissa_length(const ryu::floating_decimal_128 fd)
   { return ryu::generic128::decimalLength(fd.mantissa); }
+
+#if !defined __SIZEOF_INT128__
+  // An implementation of base-10 std::to_chars for uint128_t on targets that
+  // lack __int128.
+  std::to_chars_result
+  to_chars(char* first, char* last, uint128_t x)
+  {
+    const int len = ryu::generic128::decimalLength(x);
+    if (last - first < len)
+      return {last, std::errc::value_too_large};
+    if (x == 0)
+      {
+	*first++ = '0';
+	return {first, std::errc{}};
+      }
+    for (int i = 0; i < len; ++i)
+      {
+	first[len - 1 - i] = '0' + static_cast<char>(x % 10);
+	x /= 10;
+      }
+    __glibcxx_assert(x == 0);
+    return {first + len, std::errc{}};
+  }
 #endif
 } // anon namespace
 
diff --git a/libstdc++-v3/src/c++17/uint128_t.h b/libstdc++-v3/src/c++17/uint128_t.h
new file mode 100644
index 00000000000..90ebae2ffd2
--- /dev/null
+++ b/libstdc++-v3/src/c++17/uint128_t.h
@@ -0,0 +1,297 @@
+// A relatiely minimal unsigned 128-bit integer class type, used by the
+// floating-point std::to_chars implementation on targets that lack __int128.
+
+// Copyright (C) 2021 Free Software Foundation, Inc.
+//
+// This file is part of the GNU ISO C++ Library.  This library is free
+// software; you can redistribute it and/or modify it under the
+// terms of the GNU General Public License as published by the
+// Free Software Foundation; either version 3, or (at your option)
+// any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+
+// Under Section 7 of GPL version 3, you are granted additional
+// permissions described in the GCC Runtime Library Exception, version
+// 3.1, as published by the Free Software Foundation.
+
+// You should have received a copy of the GNU General Public License and
+// a copy of the GCC Runtime Library Exception along with this program;
+// see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
+// <http://www.gnu.org/licenses/>.
+
+struct uint128_t
+{
+#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
+  uint64_t lo, hi;
+#else
+  uint64_t hi, lo;
+#endif
+
+  uint128_t() = default;
+
+  constexpr
+  uint128_t(uint64_t lo, uint64_t hi = 0)
+#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
+    : lo(lo), hi(hi)
+#else
+    : hi(hi), lo(lo)
+#endif
+  { }
+
+  constexpr explicit
+  operator bool() const
+  { return *this != 0; }
+
+  template<typename T, typename = std::enable_if_t<std::is_integral_v<T>>>
+    constexpr explicit
+    operator T() const
+    {
+      static_assert(sizeof(T) <= sizeof(uint64_t));
+      return static_cast<T>(lo);
+    }
+
+  friend constexpr uint128_t
+  operator&(uint128_t x, const uint128_t y)
+  {
+    x.lo &= y.lo;
+    x.hi &= y.hi;
+    return x;
+  }
+
+  friend constexpr uint128_t
+  operator|(uint128_t x, const uint128_t y)
+  {
+    x.lo |= y.lo;
+    x.hi |= y.hi;
+    return x;
+  }
+
+  friend constexpr uint128_t
+  operator<<(uint128_t x, const uint128_t y)
+  {
+    __glibcxx_assert(y < 128);
+    // TODO: Convince GCC to use shldq on x86 here.
+    if (y.lo >= 64)
+      {
+	x.hi = x.lo << (y.lo - 64);
+	x.lo = 0;
+      }
+    else if (y.lo != 0)
+      {
+	x.hi <<= y.lo;
+	x.hi |= x.lo >> (64 - y.lo);
+	x.lo <<= y.lo;
+      }
+    return x;
+  }
+
+  friend constexpr uint128_t
+  operator>>(uint128_t x, const uint128_t y)
+  {
+    __glibcxx_assert(y < 128);
+    // TODO: Convince GCC to use shrdq on x86 here.
+    if (y.lo >= 64)
+      {
+	x.lo = x.hi >> (y.lo - 64);
+	x.hi = 0;
+      }
+    else if (y.lo != 0)
+      {
+	x.lo >>= y.lo;
+	x.lo |= x.hi << (64 - y.lo);
+	x.hi >>= y.lo;
+      }
+    return x;
+  }
+
+  constexpr uint128_t
+  operator~() const
+  { return {~lo, ~hi}; }
+
+  constexpr uint128_t
+  operator-() const
+  { return operator~() + 1; }
+
+  friend constexpr uint128_t
+  operator+(uint128_t x, const uint128_t y)
+  {
+    x.hi += __builtin_add_overflow(x.lo, y.lo, &x.lo);
+    x.hi += y.hi;
+    return x;
+  }
+
+  friend constexpr uint128_t
+  operator-(uint128_t x, const uint128_t y)
+  {
+    x.hi -= __builtin_sub_overflow(x.lo, y.lo, &x.lo);
+    x.hi -= y.hi;
+    return x;
+  }
+
+  static constexpr uint128_t
+  umul64_64_128(const uint64_t x, const uint64_t y)
+  {
+    const uint64_t xl = x & 0xffffffff;
+    const uint64_t xh = x >> 32;
+    const uint64_t yl = y & 0xffffffff;
+    const uint64_t yh = y >> 32;
+    const uint64_t ll = xl * yl;
+    const uint64_t lh = xl * yh;
+    const uint64_t hl = xh * yl;
+    const uint64_t hh = xh * yh;
+    const uint64_t m = (ll >> 32) + lh + (hl & 0xffffffff);
+    const uint64_t l = (ll & 0xffffffff ) | (m << 32);
+    const uint64_t h = (m >> 32) + (hl >> 32) + hh;
+    return {l, h};
+  }
+
+  friend constexpr uint128_t
+  operator*(uint128_t x, const uint128_t y)
+  {
+    uint128_t z = umul64_64_128(x.lo, y.lo);
+    z.hi += x.hi*y.lo + y.hi*x.lo;
+    return z;
+  }
+
+  friend constexpr uint128_t
+  operator/(const uint128_t x, const uint128_t y)
+  {
+    // Ryu performs 128-bit division only by 5 and 10, so that's what we
+    // implement.  The strategy here is to relate division of x with that of
+    // x.hi and x.lo separately.
+    __glibcxx_assert(y == 5 || y == 10);
+    // The following implements division by 5 and 10.  In either case, we
+    // first compute division by 5:
+    //   x/5 = (x.hi*2^64 + x.lo)/5
+    //       = (x.hi*(2^64-1) + x.hi + x.lo)/5
+    //       = x.hi*((2^64-1)/5) + (x.hi + x.lo)/5 since CST=(2^64-1)/5 is exact
+    //       = x.hi*CST + x.hi/5 + x.lo/5 + ((x.lo%5) + (x.hi%5) >= 5)
+    // We go a step further and replace the last adjustment term with a
+    // lookup table, which we encode as a binary literal.  This seems to
+    // yield smaller code on x86 at least.
+    constexpr auto cst = ~uint64_t(0) / 5;
+    uint128_t q = uint128_t{x.hi}*cst + uint128_t{x.hi/5 + x.lo/5};
+    constexpr auto lookup = 0b111100000u;
+    q += (lookup >> ((x.hi % 5) + (x.lo % 5))) & 1;
+    if (y == 10)
+      q >>= 1;
+    return q;
+  }
+
+  friend constexpr uint128_t
+  operator%(uint128_t x, const uint128_t y)
+  {
+    // Ryu performs 128-bit modulus only by 2, 5 and 10, so that's what we
+    // implement.  The strategy here is to relate modulus of x with that of
+    // x.hi and x.lo separately.
+    if (y == 2)
+      return x & 1;
+    __glibcxx_assert(y == 5 || y == 10);
+    // The following implements modulus by 5 and 10.  In either case,
+    // we first compute modulus by 5:
+    //   x (mod 5) = x.hi*2^64 + x.lo (mod 5)
+    //             = x.hi + x.lo (mod 5) since 2^64 ≡ 1 (mod 5)
+    // So the straightforward implementation would be
+    //   ((x.hi % 5) + (x.lo % 5)) % 5
+    // But we go a step further and replace the outermost % with a
+    // lookup table:
+    //           = {0,1,2,3,4,0,1,2,3}[(x.hi % 5) + (x.lo % 5)] (mod 5)
+    // which we encode as an octal literal.
+    constexpr auto lookup = 0321043210u;
+    auto r = (lookup >> 3*((x.hi % 5) + (x.lo % 5))) & 7;
+    if (y == 10)
+      // x % 10 = (x % 5)      if x / 5 is even
+      //          (x % 5) + 5  if x / 5 is odd
+      // The compiler should be able to CSE the below computation of x/5 and
+      // the above modulus operations with a nearby inlined computation of x/10.
+      r += ((x / 5).lo & 1) * 5;
+    return r;
+  }
+
+  friend constexpr bool
+  operator==(const uint128_t x, const uint128_t y)
+  { return x.hi == y.hi && x.lo == y.lo; }
+
+  friend constexpr bool
+  operator<(const uint128_t x, const uint128_t y)
+  { return x.hi < y.hi || (x.hi == y.hi && x.lo < y.lo); }
+
+  friend constexpr auto
+  __bit_width(const uint128_t x)
+  {
+    if (auto w = std::__bit_width(x.hi))
+      return w + 64;
+    else
+      return std::__bit_width(x.lo);
+  }
+
+  friend constexpr auto
+  __countr_zero(const uint128_t x)
+  {
+    auto c = std::__countr_zero(x.lo);
+    if (c == 64)
+      return 64 + std::__countr_zero(x.hi);
+    else
+      return c;
+  }
+
+  constexpr uint128_t&
+  operator--()
+  { return *this -= 1; }
+
+  constexpr uint128_t&
+  operator++()
+  { return *this += 1; }
+
+  constexpr uint128_t&
+  operator+=(const uint128_t y)
+  { return *this = *this + y; }
+
+  constexpr uint128_t&
+  operator-=(const uint128_t y)
+  { return *this = *this - y; }
+
+  constexpr uint128_t&
+  operator*=(const uint128_t y)
+  { return *this = *this * y; }
+
+  constexpr uint128_t&
+  operator<<=(const uint128_t y)
+  { return *this = *this << y; }
+
+  constexpr uint128_t&
+  operator>>=(const uint128_t y)
+  { return *this = *this >> y; }
+
+  constexpr uint128_t&
+  operator|=(const uint128_t y)
+  { return *this = *this | y; }
+
+  constexpr uint128_t&
+  operator&=(const uint128_t y)
+  { return *this = *this & y; }
+
+  constexpr uint128_t&
+  operator%=(const uint128_t y)
+  { return *this = *this % y; }
+
+  constexpr uint128_t&
+  operator/=(const uint128_t y)
+  { return *this = *this / y; }
+
+  friend constexpr bool
+  operator!=(const uint128_t x, const uint128_t y)
+  { return !(x == y); }
+
+  friend constexpr bool
+  operator>(const uint128_t x, const uint128_t y)
+  { return y < x; }
+
+  friend constexpr bool
+  operator>=(const uint128_t x, const uint128_t y)
+  { return !(x < y); }
+};
diff --git a/libstdc++-v3/testsuite/20_util/to_chars/long_double.cc b/libstdc++-v3/testsuite/20_util/to_chars/long_double.cc
index da847ae5401..5c1f7136f21 100644
--- a/libstdc++-v3/testsuite/20_util/to_chars/long_double.cc
+++ b/libstdc++-v3/testsuite/20_util/to_chars/long_double.cc
@@ -18,7 +18,6 @@
 // <charconv> is supported in C++14 as a GNU extension, but this test uses C++17
 // hexadecimal floating-point literals.
 // { dg-do run { target c++17 } }
-// { dg-xfail-run-if "Ryu needs __int128" { large_long_double && { ! int128 } } }
 // { dg-require-effective-target ieee-floats }
 
 #include <charconv>
-- 
2.31.0.rc2

Re: [PATCH 4/4] libstdc++: Add fallback 128-bit integer class type and use it

[Daniel Krügler]

Am Do., 11. März 2021 um 18:17 Uhr schrieb Patrick Palka via Libstdc++
<libstdc++@gcc.gnu.org>:
>
> This implements a minimal integer class type that emulates 128-bit
> unsigned arithmetic using a pair of 64-bit integers, which the
> floating-point std::to_chars implementation then uses as a drop-in
> replacement for unsigned __int128 on targets that lack the latter.
> This allows us to fully support formatting of large long double types
> on targets that lack __int128.
>
> Since Ryu performs 128-bit division/modulus only by 2, 5 and 10, the
> integer class type supports only these divisors rather than supporting
> general division/modulus.
>
> Tested on x86, x86_64, ppc64le, ppc64be and aarch64, with and without
> performing the equivalent of -U__SIZEOF_INT128__ in floating_to_chars.cc
> (so that we also test using the class type on targets when __int128 is
> available).
>
> libstdc++-v3/ChangeLog:
>
>         * src/c++17/floating_to_chars.cc: Simplify the file as if
>         __SIZEOF_INT128__ is always defined.
>         [!defined __SIZEOF_INT128__]: Include "uint128_t.h".  Define
>         a to_chars overload for the uint128_t class type.
>         * src/c++17/uint128_t.h: New file.
>         * testsuite/20_util/to_chars/long_double.cc: No longer expect an
>         execution FAIL on targets that have a large long double type
>         but lack __int128.
> ---
>  libstdc++-v3/src/c++17/floating_to_chars.cc   |  58 ++--
>  libstdc++-v3/src/c++17/uint128_t.h            | 297 ++++++++++++++++++
>  .../testsuite/20_util/to_chars/long_double.cc |   1 -
>  3 files changed, 332 insertions(+), 24 deletions(-)
>  create mode 100644 libstdc++-v3/src/c++17/uint128_t.h
>
> diff --git a/libstdc++-v3/src/c++17/floating_to_chars.cc b/libstdc++-v3/src/c++17/floating_to_chars.cc
> index da3fbaa1ed1..86f4401e134 100644
> --- a/libstdc++-v3/src/c++17/floating_to_chars.cc
> +++ b/libstdc++-v3/src/c++17/floating_to_chars.cc
> @@ -64,25 +64,19 @@ extern "C" int __sprintfieee128(char*, const char*, ...);
>
>  #if __LDBL_MANT_DIG__ == __DBL_MANT_DIG__
>  # define LONG_DOUBLE_KIND LDK_BINARY64
> -#elif defined(__SIZEOF_INT128__)
> -// The Ryu routines need a 128-bit integer type in order to do shortest
> -// formatting of types larger than 64-bit double, so without __int128 we can't
> -// support any large long double format.  This is the case for e.g. i386.
> -# if __LDBL_MANT_DIG__ == 64
> +#elif __LDBL_MANT_DIG__ == 64
>  #  define LONG_DOUBLE_KIND LDK_FLOAT80
> -# elif __LDBL_MANT_DIG__ == 113
> -#  define LONG_DOUBLE_KIND LDK_BINARY128
> -# elif __LDBL_MANT_DIG__ == 106
> -#  define LONG_DOUBLE_KIND LDK_IBM128
> -# endif
> -# if defined _GLIBCXX_USE_FLOAT128 && __FLT128_MANT_DIG__ == 113
> -// Define overloads of std::to_chars for __float128.
> -#  define FLOAT128_TO_CHARS 1
> -# endif
> +#elif __LDBL_MANT_DIG__ == 113
> +# define LONG_DOUBLE_KIND LDK_BINARY128
> +#elif __LDBL_MANT_DIG__ == 106
> +# define LONG_DOUBLE_KIND LDK_IBM128
> +#else
> +# define LONG_DOUBLE_KIND LDK_UNSUPPORTED
>  #endif
>
> -#if !defined(LONG_DOUBLE_KIND)
> -# define LONG_DOUBLE_KIND LDK_UNSUPPORTED
> +#if defined _GLIBCXX_USE_FLOAT128 && __FLT128_MANT_DIG__ == 113
> +// Define overloads of std::to_chars for __float128.
> +# define FLOAT128_TO_CHARS 1
>  #endif
>
>  // For now we only support __float128 when it's the powerpc64 __ieee128 type.
> @@ -100,6 +94,8 @@ namespace
>  {
>  #if defined __SIZEOF_INT128__
>    using uint128_t = unsigned __int128;
> +#else
> +# include "uint128_t.h"
>  #endif
>
>    namespace ryu
> @@ -114,7 +110,6 @@ namespace
>  #include "ryu/d2fixed.c"
>  #include "ryu/f2s.c"
>
> -#ifdef __SIZEOF_INT128__
>      namespace generic128
>      {
>        // Put the generic Ryu bits in their own namespace to avoid name conflicts.
> @@ -129,7 +124,6 @@ namespace
>      int
>      to_chars(const floating_decimal_128 v, char* const result)
>      { return generic128::generic_to_chars(v, result); }
> -#endif
>    } // namespace ryu
>
>    // A traits class that contains pertinent information about the binary
> @@ -407,10 +401,8 @@ namespace
>           return uint32_t{};
>         else if constexpr (total_bits <= 64)
>           return uint64_t{};
> -#ifdef __SIZEOF_INT128__
>         else if constexpr (total_bits <= 128)
>           return uint128_t{};
> -#endif
>        };
>        using uint_t = decltype(get_uint_t());
>        uint_t value_bits = 0;
> @@ -503,7 +495,6 @@ namespace
>         return ryu::floating_to_fd32(value);
>        else if constexpr (std::is_same_v<T, double>)
>         return ryu::floating_to_fd64(value);
> -#ifdef __SIZEOF_INT128__
>        else if constexpr (std::is_same_v<T, long double>
>                          || std::is_same_v<T, F128_type>)
>         {
> @@ -519,7 +510,6 @@ namespace
>                                                 mantissa_bits, exponent_bits,
>                                                 !has_implicit_leading_bit);
>         }
> -#endif
>      }
>
>    // This subroutine returns true if the shortest scientific form fd is a
> @@ -558,10 +548,32 @@ namespace
>    get_mantissa_length(const ryu::floating_decimal_64 fd)
>    { return ryu::decimalLength17(fd.mantissa); }
>
> -#ifdef __SIZEOF_INT128__
>    int
>    get_mantissa_length(const ryu::floating_decimal_128 fd)
>    { return ryu::generic128::decimalLength(fd.mantissa); }
> +
> +#if !defined __SIZEOF_INT128__
> +  // An implementation of base-10 std::to_chars for uint128_t on targets that
> +  // lack __int128.
> +  std::to_chars_result
> +  to_chars(char* first, char* last, uint128_t x)
> +  {
> +    const int len = ryu::generic128::decimalLength(x);
> +    if (last - first < len)
> +      return {last, std::errc::value_too_large};
> +    if (x == 0)
> +      {
> +       *first++ = '0';
> +       return {first, std::errc{}};
> +      }
> +    for (int i = 0; i < len; ++i)
> +      {
> +       first[len - 1 - i] = '0' + static_cast<char>(x % 10);
> +       x /= 10;
> +      }
> +    __glibcxx_assert(x == 0);
> +    return {first + len, std::errc{}};
> +  }
>  #endif
>  } // anon namespace
>
> diff --git a/libstdc++-v3/src/c++17/uint128_t.h b/libstdc++-v3/src/c++17/uint128_t.h
> new file mode 100644
> index 00000000000..90ebae2ffd2
> --- /dev/null
> +++ b/libstdc++-v3/src/c++17/uint128_t.h
> @@ -0,0 +1,297 @@
> +// A relatiely minimal unsigned 128-bit integer class type, used by the
> +// floating-point std::to_chars implementation on targets that lack __int128.
> +
> +// Copyright (C) 2021 Free Software Foundation, Inc.
> +//
> +// This file is part of the GNU ISO C++ Library.  This library is free
> +// software; you can redistribute it and/or modify it under the
> +// terms of the GNU General Public License as published by the
> +// Free Software Foundation; either version 3, or (at your option)
> +// any later version.
> +
> +// This library is distributed in the hope that it will be useful,
> +// but WITHOUT ANY WARRANTY; without even the implied warranty of
> +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
> +// GNU General Public License for more details.
> +
> +// Under Section 7 of GPL version 3, you are granted additional
> +// permissions described in the GCC Runtime Library Exception, version
> +// 3.1, as published by the Free Software Foundation.
> +
> +// You should have received a copy of the GNU General Public License and
> +// a copy of the GCC Runtime Library Exception along with this program;
> +// see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
> +// <http://www.gnu.org/licenses/>.
> +
> +struct uint128_t
> +{
> +#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
> +  uint64_t lo, hi;
> +#else
> +  uint64_t hi, lo;
> +#endif
> +
> +  uint128_t() = default;
> +
> +  constexpr
> +  uint128_t(uint64_t lo, uint64_t hi = 0)
> +#if __BYTE_ORDER__ == __ORDER_LITTLE_ENDIAN__
> +    : lo(lo), hi(hi)
> +#else
> +    : hi(hi), lo(lo)
> +#endif
> +  { }
> +
> +  constexpr explicit
> +  operator bool() const
> +  { return *this != 0; }
> +
> +  template<typename T, typename = std::enable_if_t<std::is_integral_v<T>>>
> +    constexpr explicit
> +    operator T() const
> +    {
> +      static_assert(sizeof(T) <= sizeof(uint64_t));
> +      return static_cast<T>(lo);
> +    }
> +
> +  friend constexpr uint128_t
> +  operator&(uint128_t x, const uint128_t y)
> +  {
> +    x.lo &= y.lo;
> +    x.hi &= y.hi;
> +    return x;
> +  }
> +
> +  friend constexpr uint128_t
> +  operator|(uint128_t x, const uint128_t y)
> +  {
> +    x.lo |= y.lo;
> +    x.hi |= y.hi;
> +    return x;
> +  }
> +
> +  friend constexpr uint128_t
> +  operator<<(uint128_t x, const uint128_t y)
> +  {
> +    __glibcxx_assert(y < 128);
> +    // TODO: Convince GCC to use shldq on x86 here.
> +    if (y.lo >= 64)
> +      {
> +       x.hi = x.lo << (y.lo - 64);
> +       x.lo = 0;
> +      }
> +    else if (y.lo != 0)
> +      {
> +       x.hi <<= y.lo;
> +       x.hi |= x.lo >> (64 - y.lo);
> +       x.lo <<= y.lo;
> +      }
> +    return x;
> +  }
> +
> +  friend constexpr uint128_t
> +  operator>>(uint128_t x, const uint128_t y)
> +  {
> +    __glibcxx_assert(y < 128);
> +    // TODO: Convince GCC to use shrdq on x86 here.
> +    if (y.lo >= 64)
> +      {
> +       x.lo = x.hi >> (y.lo - 64);
> +       x.hi = 0;
> +      }
> +    else if (y.lo != 0)
> +      {
> +       x.lo >>= y.lo;
> +       x.lo |= x.hi << (64 - y.lo);
> +       x.hi >>= y.lo;
> +      }
> +    return x;
> +  }
> +
> +  constexpr uint128_t
> +  operator~() const
> +  { return {~lo, ~hi}; }
> +
> +  constexpr uint128_t
> +  operator-() const
> +  { return operator~() + 1; }
> +
> +  friend constexpr uint128_t
> +  operator+(uint128_t x, const uint128_t y)
> +  {
> +    x.hi += __builtin_add_overflow(x.lo, y.lo, &x.lo);
> +    x.hi += y.hi;
> +    return x;
> +  }
> +
> +  friend constexpr uint128_t
> +  operator-(uint128_t x, const uint128_t y)
> +  {
> +    x.hi -= __builtin_sub_overflow(x.lo, y.lo, &x.lo);
> +    x.hi -= y.hi;
> +    return x;
> +  }
> +
> +  static constexpr uint128_t
> +  umul64_64_128(const uint64_t x, const uint64_t y)
> +  {
> +    const uint64_t xl = x & 0xffffffff;
> +    const uint64_t xh = x >> 32;
> +    const uint64_t yl = y & 0xffffffff;
> +    const uint64_t yh = y >> 32;
> +    const uint64_t ll = xl * yl;
> +    const uint64_t lh = xl * yh;
> +    const uint64_t hl = xh * yl;
> +    const uint64_t hh = xh * yh;
> +    const uint64_t m = (ll >> 32) + lh + (hl & 0xffffffff);
> +    const uint64_t l = (ll & 0xffffffff ) | (m << 32);
> +    const uint64_t h = (m >> 32) + (hl >> 32) + hh;
> +    return {l, h};
> +  }
> +
> +  friend constexpr uint128_t
> +  operator*(uint128_t x, const uint128_t y)
> +  {
> +    uint128_t z = umul64_64_128(x.lo, y.lo);
> +    z.hi += x.hi*y.lo + y.hi*x.lo;
> +    return z;
> +  }
> +
> +  friend constexpr uint128_t
> +  operator/(const uint128_t x, const uint128_t y)
> +  {
> +    // Ryu performs 128-bit division only by 5 and 10, so that's what we
> +    // implement.  The strategy here is to relate division of x with that of
> +    // x.hi and x.lo separately.
> +    __glibcxx_assert(y == 5 || y == 10);
> +    // The following implements division by 5 and 10.  In either case, we
> +    // first compute division by 5:
> +    //   x/5 = (x.hi*2^64 + x.lo)/5
> +    //       = (x.hi*(2^64-1) + x.hi + x.lo)/5
> +    //       = x.hi*((2^64-1)/5) + (x.hi + x.lo)/5 since CST=(2^64-1)/5 is exact
> +    //       = x.hi*CST + x.hi/5 + x.lo/5 + ((x.lo%5) + (x.hi%5) >= 5)
> +    // We go a step further and replace the last adjustment term with a
> +    // lookup table, which we encode as a binary literal.  This seems to
> +    // yield smaller code on x86 at least.
> +    constexpr auto cst = ~uint64_t(0) / 5;
> +    uint128_t q = uint128_t{x.hi}*cst + uint128_t{x.hi/5 + x.lo/5};
> +    constexpr auto lookup = 0b111100000u;
> +    q += (lookup >> ((x.hi % 5) + (x.lo % 5))) & 1;
> +    if (y == 10)
> +      q >>= 1;
> +    return q;
> +  }
> +
> +  friend constexpr uint128_t
> +  operator%(uint128_t x, const uint128_t y)
> +  {
> +    // Ryu performs 128-bit modulus only by 2, 5 and 10, so that's what we
> +    // implement.  The strategy here is to relate modulus of x with that of
> +    // x.hi and x.lo separately.
> +    if (y == 2)
> +      return x & 1;
> +    __glibcxx_assert(y == 5 || y == 10);
> +    // The following implements modulus by 5 and 10.  In either case,
> +    // we first compute modulus by 5:
> +    //   x (mod 5) = x.hi*2^64 + x.lo (mod 5)
> +    //             = x.hi + x.lo (mod 5) since 2^64 ≡ 1 (mod 5)
> +    // So the straightforward implementation would be
> +    //   ((x.hi % 5) + (x.lo % 5)) % 5
> +    // But we go a step further and replace the outermost % with a
> +    // lookup table:
> +    //           = {0,1,2,3,4,0,1,2,3}[(x.hi % 5) + (x.lo % 5)] (mod 5)
> +    // which we encode as an octal literal.
> +    constexpr auto lookup = 0321043210u;
> +    auto r = (lookup >> 3*((x.hi % 5) + (x.lo % 5))) & 7;
> +    if (y == 10)
> +      // x % 10 = (x % 5)      if x / 5 is even
> +      //          (x % 5) + 5  if x / 5 is odd
> +      // The compiler should be able to CSE the below computation of x/5 and
> +      // the above modulus operations with a nearby inlined computation of x/10.
> +      r += ((x / 5).lo & 1) * 5;
> +    return r;
> +  }
> +
> +  friend constexpr bool
> +  operator==(const uint128_t x, const uint128_t y)
> +  { return x.hi == y.hi && x.lo == y.lo; }
> +
> +  friend constexpr bool
> +  operator<(const uint128_t x, const uint128_t y)
> +  { return x.hi < y.hi || (x.hi == y.hi && x.lo < y.lo); }
> +
> +  friend constexpr auto
> +  __bit_width(const uint128_t x)
> +  {
> +    if (auto w = std::__bit_width(x.hi))
> +      return w + 64;
> +    else
> +      return std::__bit_width(x.lo);
> +  }
> +
> +  friend constexpr auto
> +  __countr_zero(const uint128_t x)
> +  {
> +    auto c = std::__countr_zero(x.lo);
> +    if (c == 64)
> +      return 64 + std::__countr_zero(x.hi);
> +    else
> +      return c;
> +  }
> +
> +  constexpr uint128_t&
> +  operator--()
> +  { return *this -= 1; }
> +
> +  constexpr uint128_t&
> +  operator++()
> +  { return *this += 1; }
> +
> +  constexpr uint128_t&
> +  operator+=(const uint128_t y)
> +  { return *this = *this + y; }
> +
> +  constexpr uint128_t&
> +  operator-=(const uint128_t y)
> +  { return *this = *this - y; }
> +
> +  constexpr uint128_t&
> +  operator*=(const uint128_t y)
> +  { return *this = *this * y; }
> +
> +  constexpr uint128_t&
> +  operator<<=(const uint128_t y)
> +  { return *this = *this << y; }
> +
> +  constexpr uint128_t&
> +  operator>>=(const uint128_t y)
> +  { return *this = *this >> y; }
> +
> +  constexpr uint128_t&
> +  operator|=(const uint128_t y)
> +  { return *this = *this | y; }
> +
> +  constexpr uint128_t&
> +  operator&=(const uint128_t y)
> +  { return *this = *this & y; }
> +
> +  constexpr uint128_t&
> +  operator%=(const uint128_t y)
> +  { return *this = *this % y; }
> +
> +  constexpr uint128_t&
> +  operator/=(const uint128_t y)
> +  { return *this = *this / y; }
> +
> +  friend constexpr bool
> +  operator!=(const uint128_t x, const uint128_t y)
> +  { return !(x == y); }
> +
> +  friend constexpr bool
> +  operator>(const uint128_t x, const uint128_t y)
> +  { return y < x; }
> +
> +  friend constexpr bool
> +  operator>=(const uint128_t x, const uint128_t y)
> +  { return !(x < y); }
> +};
> diff --git a/libstdc++-v3/testsuite/20_util/to_chars/long_double.cc b/libstdc++-v3/testsuite/20_util/to_chars/long_double.cc
> index da847ae5401..5c1f7136f21 100644
> --- a/libstdc++-v3/testsuite/20_util/to_chars/long_double.cc
> +++ b/libstdc++-v3/testsuite/20_util/to_chars/long_double.cc
> @@ -18,7 +18,6 @@
>  // <charconv> is supported in C++14 as a GNU extension, but this test uses C++17
>  // hexadecimal floating-point literals.
>  // { dg-do run { target c++17 } }
> -// { dg-xfail-run-if "Ryu needs __int128" { large_long_double && { ! int128 } } }
>  // { dg-require-effective-target ieee-floats }

It seems to me that basically all uint128_t operations should be
unconditionally noexcept. Currently none of them is so (The constexpr
keyword is a red herring in this regard).

Is it worth considering to add (conditional) support for operator<=>
and corresponding simplifications of comparison operators?

Thanks,

- Daniel

Re: [PATCH 4/4] libstdc++: Add fallback 128-bit integer class type and use it

[Jonathan Wakely]

On 11/03/21 19:09 +0100, Daniel Krügler via Libstdc++ wrote:
>It seems to me that basically all uint128_t operations should be
>unconditionally noexcept. Currently none of them is so (The constexpr
>keyword is a red herring in this regard).

None of this code is exposed to users (so nobody is using the noexcept
operator to heck if it throws), and is all inline, so the compiler can
see that it doesn't throw. We could make it noexcept, but I don't
think it makes any difference to anything. That could be done as a
follow-up patch.

>Is it worth considering to add (conditional) support for operator<=>
>and corresponding simplifications of comparison operators?

This is (currently) compiled with -std=gnu++17 so that's not
available. We could of course change that, but enabling experimental
C++20 features in the implementation of our less-experimental C++17
features has non-zero risk.

I definitely don't want conditional support for operator<=>
complicating the code. Either we compile it with -std=gnu++17 or we
compile it with -std=gnu++20, but we don't need to check which it is
with #if.

Re: [PATCH 4/4] libstdc++: Add fallback 128-bit integer class type and use it

[Jonathan Wakely]

On 11/03/21 12:16 -0500, Patrick Palka via Libstdc++ wrote:
>--- /dev/null
>+++ b/libstdc++-v3/src/c++17/uint128_t.h
>@@ -0,0 +1,297 @@
>+// A relatiely minimal unsigned 128-bit integer class type, used by the

"relatively"

Apart from that, all four patches in the series are OK for trunk.