[committed] libstdc++: Prefer double to long double in std::shuffle_order_engine

[Jonathan Wakely]

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.

Tested powerpc64le-linux. Committed to trunk.

-------------- next part --------------
commit 60d9f254876a00260992b2f37639ef4d82d9db8f
Author: Jonathan Wakely <jwakely@redhat.com>
Date:   Sat Oct 31 07:16:47 2020

    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.

diff --git a/libstdc++-v3/include/bits/random.h b/libstdc++-v3/include/bits/random.h
index 0be1191e07de..32537f80df72 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 bf39a51559bb..3205442f2f69 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 0b5f597040b7..ba252ef34fe4 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 }