libstdc++: Use 128-bit arithmetic for std::linear_congruential_engine [PR87744]

For 32-bit targets without __int128 we need to implement the LCG
transition function by hand using 64-bit types.

We can also slightly simplify the __mod function by using if-constexpr
unconditionally, disabling -Wc++17-extensions warnings with diagnostic
pragmas.

libstdc++-v3/ChangeLog:

	PR libstdc++/87744
	* include/bits/random.h [!__SIZEOF_INT128__] (_Select_uint_least_t):
	Define specialization for 64-bit generators with
	non-power-of-two modulus and large constants.
	(__mod): Use if constexpr unconditionally.
	* testsuite/26_numerics/random/pr60037-neg.cc: Adjust dg-error
	line number.
	* testsuite/26_numerics/random/linear_congruential_engine/87744.cc:
	New test.

diff --git a/libstdc++-v3/include/bits/random.h b/libstdc++-v3/include/bits/random.h
index 0fbd092e7efb8db08bd3127fff44aa9f017be499..5fda21af882c9c9993168b40842bd87ebaf7ea17 100644 (file)
--- a/libstdc++-v3/include/bits/random.h
+++ b/libstdc++-v3/include/bits/random.h
@@ -64,6 +64,9 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
   // Implementation-space details.
   namespace __detail
   {
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wc++17-extensions"
+
     template<typename _UIntType, size_t __w,
             bool = __w < static_cast<size_t>
                          (std::numeric_limits<_UIntType>::digits)>
@@ -102,6 +105,108 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
     template<int __s>
       struct _Select_uint_least_t<__s, 1>
       { __extension__ using type = unsigned __int128; };
+#elif __has_builtin(__builtin_add_overflow) \
+    && __has_builtin(__builtin_sub_overflow) \
+    && defined __UINT64_TYPE__
+    template<int __s>
+      struct _Select_uint_least_t<__s, 1>
+      {
+       // This is NOT a general-purpose 128-bit integer type.
+       // It only supports (type(a) * x + c) % m as needed by __mod.
+       struct type
+       {
+         explicit
+         type(uint64_t __a) noexcept : _M_lo(__a), _M_hi(0) { }
+
+         // pre: __l._M_hi == 0
+         friend type
+         operator*(type __l, uint64_t __x) noexcept
+         {
+           // Split 64-bit values __l._M_lo and __x into high and low 32-bit
+           // limbs and multiply those individually.
+           // l * x = (l0 + l1) * (x0 + x1) = l0x0 + l0x1 + l1x0 + l1x1
+
+           constexpr uint64_t __mask = 0xffffffff;
+           uint64_t __ll[2] = { __l._M_lo >> 32, __l._M_lo & __mask };
+           uint64_t __xx[2] = { __x >> 32, __x & __mask };
+           uint64_t __l0x0 = __ll[0] * __xx[0];
+           uint64_t __l0x1 = __ll[0] * __xx[1];
+           uint64_t __l1x0 = __ll[1] * __xx[0];
+           uint64_t __l1x1 = __ll[1] * __xx[1];
+           // These bits are the low half of __l._M_hi
+           // and the high half of __l._M_lo.
+           uint64_t __mid
+             = (__l0x1 & __mask) + (__l1x0 & __mask) + (__l1x1 >> 32);
+           __l._M_hi = __l0x0 + (__l0x1 >> 32) + (__l1x0 >> 32) + (__mid >> 32);
+           __l._M_lo = (__mid << 32) + (__l1x1 & __mask);
+           return __l;
+         }
+
+         friend type
+         operator+(type __l, uint64_t __c) noexcept
+         {
+           __l._M_hi += __builtin_add_overflow(__l._M_lo, __c, &__l._M_lo);
+           return __l;
+         }
+
+         friend type
+         operator%(type __l, uint64_t __m) noexcept
+         {
+           if (__builtin_expect(__l._M_hi == 0, 0))
+             {
+               __l._M_lo %= __m;
+               return __l;
+             }
+
+           int __shift = __builtin_clzll(__m) + 64
+                           - __builtin_clzll(__l._M_hi);
+           type __x(0);
+           if (__shift >= 64)
+             {
+               __x._M_hi = __m << (__shift - 64);
+               __x._M_lo = 0;
+             }
+           else
+             {
+               __x._M_hi = __m >> (64 - __shift);
+               __x._M_lo = __m << __shift;
+             }
+
+           while (__l._M_hi != 0 || __l._M_lo >= __m)
+             {
+               if (__x <= __l)
+                 {
+                   __l._M_hi -= __x._M_hi;
+                   __l._M_hi -= __builtin_sub_overflow(__l._M_lo, __x._M_lo,
+                                                       &__l._M_lo);
+                 }
+               __x._M_lo = (__x._M_lo >> 1) | (__x._M_hi << 63);
+               __x._M_hi >>= 1;
+             }
+           return __l;
+         }
+
+         // pre: __l._M_hi == 0
+         explicit operator uint64_t() const noexcept
+         { return _M_lo; }
+
+         friend bool operator<(const type& __l, const type& __r) noexcept
+         {
+           if (__l._M_hi < __r._M_hi)
+             return true;
+           else if (__l._M_hi == __r._M_hi)
+             return __l._M_lo < __r._M_lo;
+           else
+             return false;
+         }
+
+         friend bool operator<=(const type& __l, const type& __r) noexcept
+         { return !(__r < __l); }
+
+         uint64_t _M_lo;
+         uint64_t _M_hi;
+       };
+      };
 #endif
 
     // Assume a != 0, a < m, c < m, x < m.
@@ -149,14 +254,10 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
       inline _Tp
       __mod(_Tp __x)
       {
-       if _GLIBCXX17_CONSTEXPR (__a == 0)
+       if constexpr (__a == 0)
          return __c;
-       else
-         {
-           // _Mod must not be instantiated with a == 0
-           constexpr _Tp __a1 = __a ? __a : 1;
-           return _Mod<_Tp, __m, __a1, __c>::__calc(__x);
-         }
+       else // N.B. _Mod must not be instantiated with a == 0
+         return _Mod<_Tp, __m, __a, __c>::__calc(__x);
       }
 
     /*
@@ -216,6 +317,7 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
        __not_<is_convertible<_Sseq, _Res>>
       >;
 
+#pragma GCC diagnostic pop
   } // namespace __detail
   /// @endcond
 

diff --git a/libstdc++-v3/testsuite/26_numerics/random/linear_congruential_engine/87744.cc b/libstdc++-v3/testsuite/26_numerics/random/linear_congruential_engine/87744.cc
new file mode 100644 (file)
index 0000000..f8a8f9f
--- /dev/null
+++ b/libstdc++-v3/testsuite/26_numerics/random/linear_congruential_engine/87744.cc
@@ -0,0 +1,22 @@
+// { dg-do run { target c++11 } }
+// PR libstdc++/87744 Some valid instantiations of linear_congruential_engine
+// produce compiler errors when __int128 isn't available
+
+#include <random>
+#include <testsuite_hooks.h>
+
+int main()
+{
+  std::linear_congruential_engine<unsigned long long int,
+                                 864691128455135232ULL, 12347ULL,
+                                 4052555153018976267ULL> gen;
+  gen();
+  VERIFY( gen() == 20120904253298372 );
+  VERIFY( gen() == 499698276788149646 );
+
+  std::linear_congruential_engine<unsigned long long, 6364136223846793005ULL,
+                                 1ULL, (-1ULL >> 1)> gen2;
+  for (int i = 0; i < 99; ++i)
+    gen2();
+  VERIFY( gen2() == 5913590678204212798 );
+}

diff --git a/libstdc++-v3/testsuite/26_numerics/random/pr60037-neg.cc b/libstdc++-v3/testsuite/26_numerics/random/pr60037-neg.cc
index 59cf84adb48d32e5f08ac920a792e9efb9edd3cb..4c24e56cea275735dfff072d9ab73bd74fe47173 100644 (file)
--- a/libstdc++-v3/testsuite/26_numerics/random/pr60037-neg.cc
+++ b/libstdc++-v3/testsuite/26_numerics/random/pr60037-neg.cc
@@ -10,6 +10,6 @@ std::__detail::_Adaptor<std::mt19937, unsigned long> aurng(urng);
 auto x = std::generate_canonical<std::size_t,
                        std::numeric_limits<std::size_t>::digits>(urng);
 
-// { dg-error "static assertion failed: template argument must be a floating point type" "" { target *-*-* } 169 }
+// { dg-error "static assertion failed: template argument must be a floating point type" "" { target *-*-* } 271 }
 
 // { dg-error "static assertion failed: template argument must be a floating point type" "" { target *-*-* } 3351 }