libstdc++: Prefer double to long double in std::shuffle_order_engine

The transition algorithm for std::shuffle_order_engine uses long double
to ensure that the value (max() - min() + 1) can be accurately
represented, to avoid bias in the shuffling. However, when the base
engine's range is small enough we can avoid slower long double
arithmetic by using double. For example, long double is unnecessary for
any base engine returning 32-bit values.

This makes std::knuth_b::operator() about 15% faster on x86_64, and
probably even more on targets where long double uses soft-float.

libstdc++-v3/ChangeLog:

	* include/bits/random.h (independent_bit_engine): Fix typo
	in comment.
	(shuffle_order_engine): Fix incorrect description in comment.
	* include/bits/random.tcc (__representable_as_double
	(__p1_representable_as_double): New helper functions.
	(shuffle_order_engine::operator()): Use double for calculation
	if (max() - min() + 1) is representable as double.
	* testsuite/26_numerics/random/pr60037-neg.cc: Adjust dg-error
	line number.

3 files changed, 43 insertions, 7 deletions

diff --git a/libstdc++-v3/include/bits/random.h b/libstdc++-v3/include/bits/random.h
index 0be1191..32537f8 100644
--- a/libstdc++-v3/include/bits/random.h
+++ b/libstdc++-v3/include/bits/random.h
@@ -1099,7 +1099,7 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
 
   /**
    * Produces random numbers by combining random numbers from some base
-   * engine to produce random numbers with a specifies number of bits @p __w.
+   * engine to produce random numbers with a specified number of bits @p __w.
    */
   template<typename _RandomNumberEngine, size_t __w, typename _UIntType>
     class independent_bits_engine
@@ -1316,9 +1316,11 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
 
 
   /**
-   * @brief Produces random numbers by combining random numbers from some
-   * base engine to produce random numbers with a specifies number of bits
-   * @p __k.
+   * @brief Produces random numbers by reordering random numbers from some
+   * base engine.
+   *
+   * The values from the base engine are stored in a sequence of size @p __k
+   * and shuffled by an algorithm that depends on those values.
    */
   template<typename _RandomNumberEngine, size_t __k>
     class shuffle_order_engine
diff --git a/libstdc++-v3/include/bits/random.tcc b/libstdc++-v3/include/bits/random.tcc
index bf39a51..3205442 100644
--- a/libstdc++-v3/include/bits/random.tcc
+++ b/libstdc++-v3/include/bits/random.tcc
@@ -804,13 +804,47 @@ _GLIBCXX_BEGIN_NAMESPACE_VERSION
     constexpr size_t
     shuffle_order_engine<_RandomNumberEngine, __k>::table_size;
 
+  namespace __detail
+  {
+    // Determine whether an integer is representable as double.
+    template<typename _Tp>
+      constexpr bool
+      __representable_as_double(_Tp __x) noexcept
+      {
+	static_assert(numeric_limits<_Tp>::is_integer);
+	static_assert(!numeric_limits<_Tp>::is_signed);
+	// All integers <= 2^53 are representable.
+	return (__x <= (1ull << __DBL_MANT_DIG__))
+	  // Between 2^53 and 2^54 only even numbers are representable.
+	  || (!(__x & 1) && __detail::__representable_as_double(__x >> 1));
+      }
+
+    // Determine whether x+1 is representable as double.
+    template<typename _Tp>
+      constexpr bool
+      __p1_representable_as_double(_Tp __x) noexcept
+      {
+	static_assert(numeric_limits<_Tp>::is_integer);
+	static_assert(!numeric_limits<_Tp>::is_signed);
+	return numeric_limits<_Tp>::digits < __DBL_MANT_DIG__
+	  || (bool(__x + 1u) // return false if x+1 wraps around to zero
+	      && __detail::__representable_as_double(__x + 1u));
+      }
+  }
+
   template<typename _RandomNumberEngine, size_t __k>
     typename shuffle_order_engine<_RandomNumberEngine, __k>::result_type
     shuffle_order_engine<_RandomNumberEngine, __k>::
     operator()()
     {
-      size_t __j = __k * ((_M_y - _M_b.min())
-			  / (_M_b.max() - _M_b.min() + 1.0L));
+      constexpr result_type __range = max() - min();
+      size_t __j = __k;
+      const result_type __y = _M_y - min();
+      // Avoid using slower long double arithmetic if possible.
+      if _GLIBCXX17_CONSTEXPR (__detail::__p1_representable_as_double(__range))
+	__j *= __y / (__range + 1.0);
+      else
+	__j *= __y / (__range + 1.0L);
       _M_y = _M_v[__j];
       _M_v[__j] = _M_b();
 
diff --git a/libstdc++-v3/testsuite/26_numerics/random/pr60037-neg.cc b/libstdc++-v3/testsuite/26_numerics/random/pr60037-neg.cc
index 0b5f597..ba252ef 100644
--- a/libstdc++-v3/testsuite/26_numerics/random/pr60037-neg.cc
+++ b/libstdc++-v3/testsuite/26_numerics/random/pr60037-neg.cc
@@ -12,4 +12,4 @@ auto x = std::generate_canonical<std::size_t,
 
 // { dg-error "static assertion failed: template argument must be a floating point type" "" { target *-*-* } 167 }
 
-// { dg-error "static assertion failed: template argument must be a floating point type" "" { target *-*-* } 3312 }
+// { dg-error "static assertion failed: template argument must be a floating point type" "" { target *-*-* } 3346 }