libstdc++-v3/include/parallel/random_shuffle.h - gcc - Git at Google

 // -*- C++ -*-

 // Copyright (C) 2007-2020 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/>.

 /** @file parallel/random_shuffle.h
  *  @brief Parallel implementation of std::random_shuffle().
  *  This file is a GNU parallel extension to the Standard C++ Library.
  */

 // Written by Johannes Singler.

 #ifndef _GLIBCXX_PARALLEL_RANDOM_SHUFFLE_H
 #define _GLIBCXX_PARALLEL_RANDOM_SHUFFLE_H 1

 #include <limits>
 #include <bits/stl_numeric.h>
 #include <parallel/parallel.h>
 #include <parallel/random_number.h>

 namespace __gnu_parallel
 {
   /** @brief Type to hold the index of a bin.
     *
     *  Since many variables of this type are allocated, it should be
     *  chosen as small as possible.
     */
   typedef unsigned short _BinIndex;

   /** @brief Data known to every thread participating in
       __gnu_parallel::__parallel_random_shuffle(). */
   template<typename _RAIter>
     struct _DRandomShufflingGlobalData
     {
       typedef std::iterator_traits<_RAIter> _TraitsType;
       typedef typename _TraitsType::value_type _ValueType;
       typedef typename _TraitsType::difference_type _DifferenceType;

       /** @brief Begin iterator of the __source. */
       _RAIter& _M_source;

       /** @brief Temporary arrays for each thread. */
       _ValueType** _M_temporaries;

       /** @brief Two-dimensional array to hold the thread-bin distribution.
        *
        *  Dimensions (_M_num_threads + 1) __x (_M_num_bins + 1). */
       _DifferenceType** _M_dist;

       /** @brief Start indexes of the threads' __chunks. */
       _DifferenceType* _M_starts;

       /** @brief Number of the thread that will further process the
           corresponding bin. */
       _ThreadIndex* _M_bin_proc;

       /** @brief Number of bins to distribute to. */
       int _M_num_bins;

       /** @brief Number of bits needed to address the bins. */
       int _M_num_bits;

       /** @brief Constructor. */
       _DRandomShufflingGlobalData(_RAIter& __source)
       : _M_source(__source) { }
     };

   /** @brief Local data for a thread participating in
       __gnu_parallel::__parallel_random_shuffle().
     */
   template<typename _RAIter, typename _RandomNumberGenerator>
     struct _DRSSorterPU
     {
       /** @brief Number of threads participating in total. */
       int _M_num_threads;

       /** @brief Begin index for bins taken care of by this thread. */
       _BinIndex _M_bins_begin;

       /** @brief End index for bins taken care of by this thread. */
       _BinIndex __bins_end;

       /** @brief Random _M_seed for this thread. */
       uint32_t _M_seed;

       /** @brief Pointer to global data. */
       _DRandomShufflingGlobalData<_RAIter>* _M_sd;
     };

   /** @brief Generate a random number in @c [0,2^__logp).
     *  @param __logp Logarithm (basis 2) of the upper range __bound.
     *  @param __rng Random number generator to use.
     */
   template<typename _RandomNumberGenerator>
     inline int
     __random_number_pow2(int __logp, _RandomNumberGenerator& __rng)
     { return __rng.__genrand_bits(__logp); }

   /** @brief Random shuffle code executed by each thread.
     *  @param __pus Array of thread-local data records. */
   template<typename _RAIter, typename _RandomNumberGenerator>
     void
     __parallel_random_shuffle_drs_pu(_DRSSorterPU<_RAIter,
 				     _RandomNumberGenerator>* __pus)
     {
       typedef std::iterator_traits<_RAIter> _TraitsType;
       typedef typename _TraitsType::value_type _ValueType;
       typedef typename _TraitsType::difference_type _DifferenceType;

       _ThreadIndex __iam = omp_get_thread_num();
       _DRSSorterPU<_RAIter, _RandomNumberGenerator>* __d = &__pus[__iam];
       _DRandomShufflingGlobalData<_RAIter>* __sd = __d->_M_sd;

       // Indexing: _M_dist[bin][processor]
       _DifferenceType __length = (__sd->_M_starts[__iam + 1]
 				  - __sd->_M_starts[__iam]);
       _BinIndex* __oracles = new _BinIndex[__length];
       _DifferenceType* __dist = new _DifferenceType[__sd->_M_num_bins + 1];
       _BinIndex* __bin_proc = new _BinIndex[__sd->_M_num_bins];
       _ValueType** __temporaries = new _ValueType*[__d->_M_num_threads];

       // Compute oracles and count appearances.
       for (_BinIndex __b = 0; __b < __sd->_M_num_bins + 1; ++__b)
 	__dist[__b] = 0;
       int __num_bits = __sd->_M_num_bits;

       _RandomNumber __rng(__d->_M_seed);

       // First main loop.
       for (_DifferenceType __i = 0; __i < __length; ++__i)
 	{
           _BinIndex __oracle = __random_number_pow2(__num_bits, __rng);
           __oracles[__i] = __oracle;

           // To allow prefix (partial) sum.
           ++(__dist[__oracle + 1]);
 	}

       for (_BinIndex __b = 0; __b < __sd->_M_num_bins + 1; ++__b)
 	__sd->_M_dist[__b][__iam + 1] = __dist[__b];

 #     pragma omp barrier

 #     pragma omp single
       {
 	// Sum up bins, __sd->_M_dist[__s + 1][__d->_M_num_threads] now
 	// contains the total number of items in bin __s
 	for (_BinIndex __s = 0; __s < __sd->_M_num_bins; ++__s)
           __gnu_sequential::partial_sum(__sd->_M_dist[__s + 1],
 					__sd->_M_dist[__s + 1]
 					+ __d->_M_num_threads + 1,
 					__sd->_M_dist[__s + 1]);
       }

 #     pragma omp barrier

       _SequenceIndex __offset = 0, __global_offset = 0;
       for (_BinIndex __s = 0; __s < __d->_M_bins_begin; ++__s)
 	__global_offset += __sd->_M_dist[__s + 1][__d->_M_num_threads];

 #     pragma omp barrier

       for (_BinIndex __s = __d->_M_bins_begin; __s < __d->__bins_end; ++__s)
 	{
           for (int __t = 0; __t < __d->_M_num_threads + 1; ++__t)
             __sd->_M_dist[__s + 1][__t] += __offset;
           __offset = __sd->_M_dist[__s + 1][__d->_M_num_threads];
 	}

       __sd->_M_temporaries[__iam] = static_cast<_ValueType*>
 	(::operator new(sizeof(_ValueType) * __offset));

 #     pragma omp barrier

       // Draw local copies to avoid false sharing.
       for (_BinIndex __b = 0; __b < __sd->_M_num_bins + 1; ++__b)
 	__dist[__b] = __sd->_M_dist[__b][__iam];
       for (_BinIndex __b = 0; __b < __sd->_M_num_bins; ++__b)
 	__bin_proc[__b] = __sd->_M_bin_proc[__b];
       for (_ThreadIndex __t = 0; __t < __d->_M_num_threads; ++__t)
 	__temporaries[__t] = __sd->_M_temporaries[__t];

       _RAIter __source = __sd->_M_source;
       _DifferenceType __start = __sd->_M_starts[__iam];

       // Distribute according to oracles, second main loop.
       for (_DifferenceType __i = 0; __i < __length; ++__i)
 	{
           _BinIndex __target_bin = __oracles[__i];
           _ThreadIndex __target_p = __bin_proc[__target_bin];

           // Last column [__d->_M_num_threads] stays unchanged.
 	  ::new(&(__temporaries[__target_p][__dist[__target_bin + 1]++]))
               _ValueType(*(__source + __i + __start));
 	}

       delete[] __oracles;
       delete[] __dist;
       delete[] __bin_proc;
       delete[] __temporaries;

 #     pragma omp barrier

       // Shuffle bins internally.
       for (_BinIndex __b = __d->_M_bins_begin; __b < __d->__bins_end; ++__b)
 	{
           _ValueType* __begin =
 	    (__sd->_M_temporaries[__iam]
 	     + (__b == __d->_M_bins_begin
 		? 0 : __sd->_M_dist[__b][__d->_M_num_threads])),
 	    *__end = (__sd->_M_temporaries[__iam]
 		      + __sd->_M_dist[__b + 1][__d->_M_num_threads]);

           __sequential_random_shuffle(__begin, __end, __rng);
           std::copy(__begin, __end, __sd->_M_source + __global_offset
 		    + (__b == __d->_M_bins_begin
 		       ? 0 : __sd->_M_dist[__b][__d->_M_num_threads]));
 	}

       for (_SequenceIndex __i = 0; __i < __offset; ++__i)
 	__sd->_M_temporaries[__iam][__i].~_ValueType();
       ::operator delete(__sd->_M_temporaries[__iam]);
     }

   /** @brief Round up to the next greater power of 2.
     *  @param __x _Integer to round up */
   template<typename _Tp>
     _Tp
     __round_up_to_pow2(_Tp __x)
     {
       if (__x <= 1)
 	return 1;
       else
 	return (_Tp)1 << (__rd_log2(__x - 1) + 1);
     }

   /** @brief Main parallel random shuffle step.
     *  @param __begin Begin iterator of sequence.
     *  @param __end End iterator of sequence.
     *  @param __n Length of sequence.
     *  @param __num_threads Number of threads to use.
     *  @param __rng Random number generator to use.
     */
   template<typename _RAIter, typename _RandomNumberGenerator>
     void
     __parallel_random_shuffle_drs(_RAIter __begin, _RAIter __end,
 				  typename std::iterator_traits
 				  <_RAIter>::difference_type __n,
 				  _ThreadIndex __num_threads,
 				  _RandomNumberGenerator& __rng)
     {
       typedef std::iterator_traits<_RAIter> _TraitsType;
       typedef typename _TraitsType::value_type _ValueType;
       typedef typename _TraitsType::difference_type _DifferenceType;

       _GLIBCXX_CALL(__n)

       const _Settings& __s = _Settings::get();

       if (__num_threads > __n)
 	__num_threads = static_cast<_ThreadIndex>(__n);

       _BinIndex __num_bins, __num_bins_cache;

 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
       // Try the L1 cache first.

       // Must fit into L1.
       __num_bins_cache =
 	std::max<_DifferenceType>(1, __n / (__s.L1_cache_size_lb
 					    / sizeof(_ValueType)));
       __num_bins_cache = __round_up_to_pow2(__num_bins_cache);

       // No more buckets than TLB entries, power of 2
       // Power of 2 and at least one element per bin, at most the TLB size.
       __num_bins = std::min<_DifferenceType>(__n, __num_bins_cache);

 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
       // 2 TLB entries needed per bin.
       __num_bins = std::min<_DifferenceType>(__s.TLB_size / 2, __num_bins);
 #endif
       __num_bins = __round_up_to_pow2(__num_bins);

       if (__num_bins < __num_bins_cache)
 	{
 #endif
           // Now try the L2 cache
           // Must fit into L2
           __num_bins_cache = static_cast<_BinIndex>
 	    (std::max<_DifferenceType>(1, __n / (__s.L2_cache_size
 						 / sizeof(_ValueType))));
           __num_bins_cache = __round_up_to_pow2(__num_bins_cache);

           // No more buckets than TLB entries, power of 2.
           __num_bins = static_cast<_BinIndex>
 	    (std::min(__n, static_cast<_DifferenceType>(__num_bins_cache)));
           // Power of 2 and at least one element per bin, at most the TLB size.
 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
           // 2 TLB entries needed per bin.
           __num_bins = std::min(static_cast<_DifferenceType>(__s.TLB_size / 2),
 				__num_bins);
 #endif
             __num_bins = __round_up_to_pow2(__num_bins);
 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
 	}
 #endif

       __num_bins = __round_up_to_pow2(
                         std::max<_BinIndex>(__num_threads, __num_bins));

       if (__num_threads <= 1)
       {
         _RandomNumber __derived_rng(
                             __rng(std::numeric_limits<uint32_t>::max()));
 	__sequential_random_shuffle(__begin, __end, __derived_rng);
         return;
       }

       _DRandomShufflingGlobalData<_RAIter> __sd(__begin);
       _DRSSorterPU<_RAIter, _RandomNumber >* __pus;
       _DifferenceType* __starts;

 #     pragma omp parallel num_threads(__num_threads)
       {
 	_ThreadIndex __num_threads = omp_get_num_threads();
 #       pragma omp single
 	{
 	  __pus = new _DRSSorterPU<_RAIter, _RandomNumber>[__num_threads];

 	  __sd._M_temporaries = new _ValueType*[__num_threads];
 	  __sd._M_dist = new _DifferenceType*[__num_bins + 1];
 	  __sd._M_bin_proc = new _ThreadIndex[__num_bins];
 	  for (_BinIndex __b = 0; __b < __num_bins + 1; ++__b)
 	    __sd._M_dist[__b] = new _DifferenceType[__num_threads + 1];
 	  for (_BinIndex __b = 0; __b < (__num_bins + 1); ++__b)
 	    {
 	      __sd._M_dist[0][0] = 0;
 	      __sd._M_dist[__b][0] = 0;
 	    }
 	  __starts = __sd._M_starts = new _DifferenceType[__num_threads + 1];
 	  int __bin_cursor = 0;
 	  __sd._M_num_bins = __num_bins;
 	  __sd._M_num_bits = __rd_log2(__num_bins);

 	  _DifferenceType __chunk_length = __n / __num_threads,
 	                         __split = __n % __num_threads,
 	                         __start = 0;
 	  _DifferenceType __bin_chunk_length = __num_bins / __num_threads,
 	                         __bin_split = __num_bins % __num_threads;
 	  for (_ThreadIndex __i = 0; __i < __num_threads; ++__i)
 	    {
 	      __starts[__i] = __start;
 	      __start += (__i < __split
 			  ? (__chunk_length + 1) : __chunk_length);
 	      int __j = __pus[__i]._M_bins_begin = __bin_cursor;

 	      // Range of bins for this processor.
 	      __bin_cursor += (__i < __bin_split
 			       ? (__bin_chunk_length + 1)
 			       : __bin_chunk_length);
 	      __pus[__i].__bins_end = __bin_cursor;
 	      for (; __j < __bin_cursor; ++__j)
 		__sd._M_bin_proc[__j] = __i;
 	      __pus[__i]._M_num_threads = __num_threads;
 	      __pus[__i]._M_seed = __rng(std::numeric_limits<uint32_t>::max());
 	      __pus[__i]._M_sd = &__sd;
 	    }
 	  __starts[__num_threads] = __start;
 	} //single
           // Now shuffle in parallel.
 	__parallel_random_shuffle_drs_pu(__pus);
       }  // parallel

       delete[] __starts;
       delete[] __sd._M_bin_proc;
       for (int __s = 0; __s < (__num_bins + 1); ++__s)
 	delete[] __sd._M_dist[__s];
       delete[] __sd._M_dist;
       delete[] __sd._M_temporaries;

       delete[] __pus;
     }

   /** @brief Sequential cache-efficient random shuffle.
    *  @param __begin Begin iterator of sequence.
    *  @param __end End iterator of sequence.
    *  @param __rng Random number generator to use.
    */
   template<typename _RAIter, typename _RandomNumberGenerator>
     void
     __sequential_random_shuffle(_RAIter __begin, _RAIter __end,
 				_RandomNumberGenerator& __rng)
     {
       typedef std::iterator_traits<_RAIter> _TraitsType;
       typedef typename _TraitsType::value_type _ValueType;
       typedef typename _TraitsType::difference_type _DifferenceType;

       _DifferenceType __n = __end - __begin;
       const _Settings& __s = _Settings::get();

       _BinIndex __num_bins, __num_bins_cache;

 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
       // Try the L1 cache first, must fit into L1.
       __num_bins_cache = std::max<_DifferenceType>
 	(1, __n / (__s.L1_cache_size_lb / sizeof(_ValueType)));
       __num_bins_cache = __round_up_to_pow2(__num_bins_cache);

       // No more buckets than TLB entries, power of 2
       // Power of 2 and at least one element per bin, at most the TLB size
       __num_bins = std::min(__n, (_DifferenceType)__num_bins_cache);
 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
       // 2 TLB entries needed per bin
       __num_bins = std::min((_DifferenceType)__s.TLB_size / 2, __num_bins);
 #endif
       __num_bins = __round_up_to_pow2(__num_bins);

       if (__num_bins < __num_bins_cache)
 	{
 #endif
           // Now try the L2 cache, must fit into L2.
           __num_bins_cache = static_cast<_BinIndex>
 	    (std::max<_DifferenceType>(1, __n / (__s.L2_cache_size
 						 / sizeof(_ValueType))));
           __num_bins_cache = __round_up_to_pow2(__num_bins_cache);

           // No more buckets than TLB entries, power of 2
           // Power of 2 and at least one element per bin, at most the TLB size.
           __num_bins = static_cast<_BinIndex>
 	    (std::min(__n, static_cast<_DifferenceType>(__num_bins_cache)));

 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
           // 2 TLB entries needed per bin
           __num_bins = std::min<_DifferenceType>(__s.TLB_size / 2, __num_bins);
 #endif
           __num_bins = __round_up_to_pow2(__num_bins);
 #if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
 	}
 #endif

       int __num_bits = __rd_log2(__num_bins);

       if (__num_bins > 1)
 	{
           _ValueType* __target =
 	    static_cast<_ValueType*>(::operator new(sizeof(_ValueType) * __n));
           _BinIndex* __oracles = new _BinIndex[__n];
           _DifferenceType* __dist0 = new _DifferenceType[__num_bins + 1],
                 	 * __dist1 = new _DifferenceType[__num_bins + 1];

           for (int __b = 0; __b < __num_bins + 1; ++__b)
             __dist0[__b] = 0;

           _RandomNumber __bitrng(__rng(0xFFFFFFFF));

           for (_DifferenceType __i = 0; __i < __n; ++__i)
             {
               _BinIndex __oracle = __random_number_pow2(__num_bits, __bitrng);
               __oracles[__i] = __oracle;

               // To allow prefix (partial) sum.
               ++(__dist0[__oracle + 1]);
             }

           // Sum up bins.
           __gnu_sequential::partial_sum(__dist0, __dist0 + __num_bins + 1,
 					__dist0);

           for (int __b = 0; __b < __num_bins + 1; ++__b)
             __dist1[__b] = __dist0[__b];

           // Distribute according to oracles.
           for (_DifferenceType __i = 0; __i < __n; ++__i)
             ::new(&(__target[(__dist0[__oracles[__i]])++]))
 		_ValueType(*(__begin + __i));

           for (int __b = 0; __b < __num_bins; ++__b)
 	    __sequential_random_shuffle(__target + __dist1[__b],
 					__target + __dist1[__b + 1], __rng);

           // Copy elements back.
           std::copy(__target, __target + __n, __begin);

           delete[] __dist0;
           delete[] __dist1;
           delete[] __oracles;

 	  for (_DifferenceType __i = 0; __i < __n; ++__i)
 	    __target[__i].~_ValueType();
           ::operator delete(__target);
 	}
       else
 	__gnu_sequential::random_shuffle(__begin, __end, __rng);
     }

   /** @brief Parallel random public call.
    *  @param __begin Begin iterator of sequence.
    *  @param __end End iterator of sequence.
    *  @param __rng Random number generator to use.
    */
   template<typename _RAIter, typename _RandomNumberGenerator>
     inline void
     __parallel_random_shuffle(_RAIter __begin, _RAIter __end,
 			      _RandomNumberGenerator __rng = _RandomNumber())
     {
       typedef std::iterator_traits<_RAIter> _TraitsType;
       typedef typename _TraitsType::difference_type _DifferenceType;
       _DifferenceType __n = __end - __begin;
       __parallel_random_shuffle_drs(__begin, __end, __n,
 				    __get_max_threads(), __rng);
     }
 }

 #endif /* _GLIBCXX_PARALLEL_RANDOM_SHUFFLE_H */
	// -- C++ --

	// Copyright (C) 2007-2020 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/>.

	/** @file parallel/random_shuffle.h
	* @brief Parallel implementation of std::random_shuffle().
	* This file is a GNU parallel extension to the Standard C++ Library.
	*/

	// Written by Johannes Singler.

	#ifndef _GLIBCXX_PARALLEL_RANDOM_SHUFFLE_H
	#define _GLIBCXX_PARALLEL_RANDOM_SHUFFLE_H 1

	#include <limits>
	#include <bits/stl_numeric.h>
	#include <parallel/parallel.h>
	#include <parallel/random_number.h>

	namespace __gnu_parallel
	{
	/** @brief Type to hold the index of a bin.
	*
	* Since many variables of this type are allocated, it should be
	* chosen as small as possible.
	*/
	typedef unsigned short _BinIndex;

	/** @brief Data known to every thread participating in
	__gnu_parallel::__parallel_random_shuffle(). */
	template<typename _RAIter>
	struct _DRandomShufflingGlobalData
	{
	typedef std::iterator_traits<_RAIter> _TraitsType;
	typedef typename _TraitsType::value_type _ValueType;
	typedef typename _TraitsType::difference_type _DifferenceType;

	/** @brief Begin iterator of the __source. */
	_RAIter& _M_source;

	/** @brief Temporary arrays for each thread. */
	_ValueType** _M_temporaries;

	/** @brief Two-dimensional array to hold the thread-bin distribution.
	*
	* Dimensions (_M_num_threads + 1) __x (_M_num_bins + 1). */
	_DifferenceType** _M_dist;

	/** @brief Start indexes of the threads' __chunks. */
	_DifferenceType* _M_starts;

	/** @brief Number of the thread that will further process the
	corresponding bin. */
	_ThreadIndex* _M_bin_proc;

	/** @brief Number of bins to distribute to. */
	int _M_num_bins;

	/** @brief Number of bits needed to address the bins. */
	int _M_num_bits;

	/** @brief Constructor. */
	_DRandomShufflingGlobalData(_RAIter& __source)
	: _M_source(__source) { }
	};

	/** @brief Local data for a thread participating in
	__gnu_parallel::__parallel_random_shuffle().
	*/
	template<typename _RAIter, typename _RandomNumberGenerator>
	struct _DRSSorterPU
	{
	/** @brief Number of threads participating in total. */
	int _M_num_threads;

	/** @brief Begin index for bins taken care of by this thread. */
	_BinIndex _M_bins_begin;

	/** @brief End index for bins taken care of by this thread. */
	_BinIndex __bins_end;

	/** @brief Random _M_seed for this thread. */
	uint32_t _M_seed;

	/** @brief Pointer to global data. */
	_DRandomShufflingGlobalData<_RAIter>* _M_sd;
	};

	/** @brief Generate a random number in @c [0,2^__logp).
	* @param __logp Logarithm (basis 2) of the upper range __bound.
	* @param __rng Random number generator to use.
	*/
	template<typename _RandomNumberGenerator>
	inline int
	__random_number_pow2(int __logp, _RandomNumberGenerator& __rng)
	{ return __rng.__genrand_bits(__logp); }

	/** @brief Random shuffle code executed by each thread.
	* @param __pus Array of thread-local data records. */
	template<typename _RAIter, typename _RandomNumberGenerator>
	void
	__parallel_random_shuffle_drs_pu(_DRSSorterPU<_RAIter,
	_RandomNumberGenerator>* __pus)
	{
	typedef std::iterator_traits<_RAIter> _TraitsType;
	typedef typename _TraitsType::value_type _ValueType;
	typedef typename _TraitsType::difference_type _DifferenceType;

	_ThreadIndex __iam = omp_get_thread_num();
	_DRSSorterPU<_RAIter, _RandomNumberGenerator>* __d = &__pus[__iam];
	_DRandomShufflingGlobalData<_RAIter>* __sd = __d->_M_sd;

	// Indexing: _M_dist[bin][processor]
	_DifferenceType __length = (__sd->_M_starts[__iam + 1]
	- __sd->_M_starts[__iam]);
	_BinIndex* __oracles = new _BinIndex[__length];
	_DifferenceType* __dist = new _DifferenceType[__sd->_M_num_bins + 1];
	_BinIndex* __bin_proc = new _BinIndex[__sd->_M_num_bins];
	_ValueType** __temporaries = new _ValueType*[__d->_M_num_threads];

	// Compute oracles and count appearances.
	for (_BinIndex __b = 0; __b < __sd->_M_num_bins + 1; ++__b)
	__dist[__b] = 0;
	int __num_bits = __sd->_M_num_bits;

	_RandomNumber __rng(__d->_M_seed);

	// First main loop.
	for (_DifferenceType __i = 0; __i < __length; ++__i)
	{
	_BinIndex __oracle = __random_number_pow2(__num_bits, __rng);
	__oracles[__i] = __oracle;

	// To allow prefix (partial) sum.
	++(__dist[__oracle + 1]);
	}

	for (_BinIndex __b = 0; __b < __sd->_M_num_bins + 1; ++__b)
	__sd->_M_dist[__b][__iam + 1] = __dist[__b];

	# pragma omp barrier

	# pragma omp single
	{
	// Sum up bins, __sd->_M_dist[__s + 1][__d->_M_num_threads] now
	// contains the total number of items in bin __s
	for (_BinIndex __s = 0; __s < __sd->_M_num_bins; ++__s)
	__gnu_sequential::partial_sum(__sd->_M_dist[__s + 1],
	__sd->_M_dist[__s + 1]
	+ __d->_M_num_threads + 1,
	__sd->_M_dist[__s + 1]);
	}

	# pragma omp barrier

	_SequenceIndex __offset = 0, __global_offset = 0;
	for (_BinIndex __s = 0; __s < __d->_M_bins_begin; ++__s)
	__global_offset += __sd->_M_dist[__s + 1][__d->_M_num_threads];

	# pragma omp barrier

	for (_BinIndex __s = __d->_M_bins_begin; __s < __d->__bins_end; ++__s)
	{
	for (int __t = 0; __t < __d->_M_num_threads + 1; ++__t)
	__sd->_M_dist[__s + 1][__t] += __offset;
	__offset = __sd->_M_dist[__s + 1][__d->_M_num_threads];
	}

	__sd->_M_temporaries[__iam] = static_cast<_ValueType*>
	(::operator new(sizeof(_ValueType) * __offset));

	# pragma omp barrier

	// Draw local copies to avoid false sharing.
	for (_BinIndex __b = 0; __b < __sd->_M_num_bins + 1; ++__b)
	__dist[__b] = __sd->_M_dist[__b][__iam];
	for (_BinIndex __b = 0; __b < __sd->_M_num_bins; ++__b)
	__bin_proc[__b] = __sd->_M_bin_proc[__b];
	for (_ThreadIndex __t = 0; __t < __d->_M_num_threads; ++__t)
	__temporaries[__t] = __sd->_M_temporaries[__t];

	_RAIter __source = __sd->_M_source;
	_DifferenceType __start = __sd->_M_starts[__iam];

	// Distribute according to oracles, second main loop.
	for (_DifferenceType __i = 0; __i < __length; ++__i)
	{
	_BinIndex __target_bin = __oracles[__i];
	_ThreadIndex __target_p = __bin_proc[__target_bin];

	// Last column [__d->_M_num_threads] stays unchanged.
	::new(&(__temporaries[__target_p][__dist[__target_bin + 1]++]))
	_ValueType(*(__source + __i + __start));
	}

	delete[] __oracles;
	delete[] __dist;
	delete[] __bin_proc;
	delete[] __temporaries;

	# pragma omp barrier

	// Shuffle bins internally.
	for (_BinIndex __b = __d->_M_bins_begin; __b < __d->__bins_end; ++__b)
	{
	_ValueType* __begin =
	(__sd->_M_temporaries[__iam]
	+ (__b == __d->_M_bins_begin
	? 0 : __sd->_M_dist[__b][__d->_M_num_threads])),
	*__end = (__sd->_M_temporaries[__iam]
	+ __sd->_M_dist[__b + 1][__d->_M_num_threads]);

	__sequential_random_shuffle(__begin, __end, __rng);
	std::copy(__begin, __end, __sd->_M_source + __global_offset
	+ (__b == __d->_M_bins_begin
	? 0 : __sd->_M_dist[__b][__d->_M_num_threads]));
	}

	for (_SequenceIndex __i = 0; __i < __offset; ++__i)
	__sd->_M_temporaries[__iam][__i].~_ValueType();
	::operator delete(__sd->_M_temporaries[__iam]);
	}

	/** @brief Round up to the next greater power of 2.
	* @param __x _Integer to round up */
	template<typename _Tp>
	_Tp
	__round_up_to_pow2(_Tp __x)
	{
	if (__x <= 1)
	return 1;
	else
	return (_Tp)1 << (__rd_log2(__x - 1) + 1);
	}

	/** @brief Main parallel random shuffle step.
	* @param __begin Begin iterator of sequence.
	* @param __end End iterator of sequence.
	* @param __n Length of sequence.
	* @param __num_threads Number of threads to use.
	* @param __rng Random number generator to use.
	*/
	template<typename _RAIter, typename _RandomNumberGenerator>
	void
	__parallel_random_shuffle_drs(_RAIter __begin, _RAIter __end,
	typename std::iterator_traits
	<_RAIter>::difference_type __n,
	_ThreadIndex __num_threads,
	_RandomNumberGenerator& __rng)
	{
	typedef std::iterator_traits<_RAIter> _TraitsType;
	typedef typename _TraitsType::value_type _ValueType;
	typedef typename _TraitsType::difference_type _DifferenceType;

	_GLIBCXX_CALL(__n)

	const _Settings& __s = _Settings::get();

	if (__num_threads > __n)
	__num_threads = static_cast<_ThreadIndex>(__n);

	_BinIndex __num_bins, __num_bins_cache;

	#if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
	// Try the L1 cache first.

	// Must fit into L1.
	__num_bins_cache =
	std::max<_DifferenceType>(1, __n / (__s.L1_cache_size_lb
	/ sizeof(_ValueType)));
	__num_bins_cache = __round_up_to_pow2(__num_bins_cache);

	// No more buckets than TLB entries, power of 2
	// Power of 2 and at least one element per bin, at most the TLB size.
	__num_bins = std::min<_DifferenceType>(__n, __num_bins_cache);

	#if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
	// 2 TLB entries needed per bin.
	__num_bins = std::min<_DifferenceType>(__s.TLB_size / 2, __num_bins);
	#endif
	__num_bins = __round_up_to_pow2(__num_bins);

	if (__num_bins < __num_bins_cache)
	{
	#endif
	// Now try the L2 cache
	// Must fit into L2
	__num_bins_cache = static_cast<_BinIndex>
	(std::max<_DifferenceType>(1, __n / (__s.L2_cache_size
	/ sizeof(_ValueType))));
	__num_bins_cache = __round_up_to_pow2(__num_bins_cache);

	// No more buckets than TLB entries, power of 2.
	__num_bins = static_cast<_BinIndex>
	(std::min(__n, static_cast<_DifferenceType>(__num_bins_cache)));
	// Power of 2 and at least one element per bin, at most the TLB size.
	#if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
	// 2 TLB entries needed per bin.
	__num_bins = std::min(static_cast<_DifferenceType>(__s.TLB_size / 2),
	__num_bins);
	#endif
	__num_bins = __round_up_to_pow2(__num_bins);
	#if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
	}
	#endif

	__num_bins = __round_up_to_pow2(
	std::max<_BinIndex>(__num_threads, __num_bins));

	if (__num_threads <= 1)
	{
	_RandomNumber __derived_rng(
	__rng(std::numeric_limits<uint32_t>::max()));
	__sequential_random_shuffle(__begin, __end, __derived_rng);
	return;
	}

	_DRandomShufflingGlobalData<_RAIter> __sd(__begin);
	_DRSSorterPU<_RAIter, _RandomNumber >* __pus;
	_DifferenceType* __starts;

	# pragma omp parallel num_threads(__num_threads)
	{
	_ThreadIndex __num_threads = omp_get_num_threads();
	# pragma omp single
	{
	__pus = new _DRSSorterPU<_RAIter, _RandomNumber>[__num_threads];

	__sd._M_temporaries = new _ValueType*[__num_threads];
	__sd._M_dist = new _DifferenceType*[__num_bins + 1];
	__sd._M_bin_proc = new _ThreadIndex[__num_bins];
	for (_BinIndex __b = 0; __b < __num_bins + 1; ++__b)
	__sd._M_dist[__b] = new _DifferenceType[__num_threads + 1];
	for (_BinIndex __b = 0; __b < (__num_bins + 1); ++__b)
	{
	__sd._M_dist[0][0] = 0;
	__sd._M_dist[__b][0] = 0;
	}
	__starts = __sd._M_starts = new _DifferenceType[__num_threads + 1];
	int __bin_cursor = 0;
	__sd._M_num_bins = __num_bins;
	__sd._M_num_bits = __rd_log2(__num_bins);

	_DifferenceType __chunk_length = __n / __num_threads,
	__split = __n % __num_threads,
	__start = 0;
	_DifferenceType __bin_chunk_length = __num_bins / __num_threads,
	__bin_split = __num_bins % __num_threads;
	for (_ThreadIndex __i = 0; __i < __num_threads; ++__i)
	{
	__starts[__i] = __start;
	__start += (__i < __split
	? (__chunk_length + 1) : __chunk_length);
	int __j = __pus[__i]._M_bins_begin = __bin_cursor;

	// Range of bins for this processor.
	__bin_cursor += (__i < __bin_split
	? (__bin_chunk_length + 1)
	: __bin_chunk_length);
	__pus[__i].__bins_end = __bin_cursor;
	for (; __j < __bin_cursor; ++__j)
	__sd._M_bin_proc[__j] = __i;
	__pus[__i]._M_num_threads = __num_threads;
	__pus[__i]._M_seed = __rng(std::numeric_limits<uint32_t>::max());
	__pus[__i]._M_sd = &__sd;
	}
	__starts[__num_threads] = __start;
	} //single
	// Now shuffle in parallel.
	__parallel_random_shuffle_drs_pu(__pus);
	} // parallel

	delete[] __starts;
	delete[] __sd._M_bin_proc;
	for (int __s = 0; __s < (__num_bins + 1); ++__s)
	delete[] __sd._M_dist[__s];
	delete[] __sd._M_dist;
	delete[] __sd._M_temporaries;

	delete[] __pus;
	}

	/** @brief Sequential cache-efficient random shuffle.
	* @param __begin Begin iterator of sequence.
	* @param __end End iterator of sequence.
	* @param __rng Random number generator to use.
	*/
	template<typename _RAIter, typename _RandomNumberGenerator>
	void
	__sequential_random_shuffle(_RAIter __begin, _RAIter __end,
	_RandomNumberGenerator& __rng)
	{
	typedef std::iterator_traits<_RAIter> _TraitsType;
	typedef typename _TraitsType::value_type _ValueType;
	typedef typename _TraitsType::difference_type _DifferenceType;

	_DifferenceType __n = __end - __begin;
	const _Settings& __s = _Settings::get();

	_BinIndex __num_bins, __num_bins_cache;

	#if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
	// Try the L1 cache first, must fit into L1.
	__num_bins_cache = std::max<_DifferenceType>
	(1, __n / (__s.L1_cache_size_lb / sizeof(_ValueType)));
	__num_bins_cache = __round_up_to_pow2(__num_bins_cache);

	// No more buckets than TLB entries, power of 2
	// Power of 2 and at least one element per bin, at most the TLB size
	__num_bins = std::min(__n, (_DifferenceType)__num_bins_cache);
	#if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
	// 2 TLB entries needed per bin
	__num_bins = std::min((_DifferenceType)__s.TLB_size / 2, __num_bins);
	#endif
	__num_bins = __round_up_to_pow2(__num_bins);

	if (__num_bins < __num_bins_cache)
	{
	#endif
	// Now try the L2 cache, must fit into L2.
	__num_bins_cache = static_cast<_BinIndex>
	(std::max<_DifferenceType>(1, __n / (__s.L2_cache_size
	/ sizeof(_ValueType))));
	__num_bins_cache = __round_up_to_pow2(__num_bins_cache);

	// No more buckets than TLB entries, power of 2
	// Power of 2 and at least one element per bin, at most the TLB size.
	__num_bins = static_cast<_BinIndex>
	(std::min(__n, static_cast<_DifferenceType>(__num_bins_cache)));

	#if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_TLB
	// 2 TLB entries needed per bin
	__num_bins = std::min<_DifferenceType>(__s.TLB_size / 2, __num_bins);
	#endif
	__num_bins = __round_up_to_pow2(__num_bins);
	#if _GLIBCXX_RANDOM_SHUFFLE_CONSIDER_L1
	}
	#endif

	int __num_bits = __rd_log2(__num_bins);

	if (__num_bins > 1)
	{
	_ValueType* __target =
	static_cast<_ValueType>(::operator new(sizeof(_ValueType) __n));
	_BinIndex* __oracles = new _BinIndex[__n];
	_DifferenceType* __dist0 = new _DifferenceType[__num_bins + 1],
	* __dist1 = new _DifferenceType[__num_bins + 1];

	for (int __b = 0; __b < __num_bins + 1; ++__b)
	__dist0[__b] = 0;

	_RandomNumber __bitrng(__rng(0xFFFFFFFF));

	for (_DifferenceType __i = 0; __i < __n; ++__i)
	{
	_BinIndex __oracle = __random_number_pow2(__num_bits, __bitrng);
	__oracles[__i] = __oracle;

	// To allow prefix (partial) sum.
	++(__dist0[__oracle + 1]);
	}

	// Sum up bins.
	__gnu_sequential::partial_sum(__dist0, __dist0 + __num_bins + 1,
	__dist0);

	for (int __b = 0; __b < __num_bins + 1; ++__b)
	__dist1[__b] = __dist0[__b];

	// Distribute according to oracles.
	for (_DifferenceType __i = 0; __i < __n; ++__i)
	::new(&(__target[(__dist0[__oracles[__i]])++]))
	_ValueType(*(__begin + __i));

	for (int __b = 0; __b < __num_bins; ++__b)
	__sequential_random_shuffle(__target + __dist1[__b],
	__target + __dist1[__b + 1], __rng);

	// Copy elements back.
	std::copy(__target, __target + __n, __begin);

	delete[] __dist0;
	delete[] __dist1;
	delete[] __oracles;

	for (_DifferenceType __i = 0; __i < __n; ++__i)
	__target[__i].~_ValueType();
	::operator delete(__target);
	}
	else
	__gnu_sequential::random_shuffle(__begin, __end, __rng);
	}

	/** @brief Parallel random public call.
	* @param __begin Begin iterator of sequence.
	* @param __end End iterator of sequence.
	* @param __rng Random number generator to use.
	*/
	template<typename _RAIter, typename _RandomNumberGenerator>
	inline void
	__parallel_random_shuffle(_RAIter __begin, _RAIter __end,
	_RandomNumberGenerator __rng = _RandomNumber())
	{
	typedef std::iterator_traits<_RAIter> _TraitsType;
	typedef typename _TraitsType::difference_type _DifferenceType;
	_DifferenceType __n = __end - __begin;
	__parallel_random_shuffle_drs(__begin, __end, __n,
	__get_max_threads(), __rng);
	}
	}

	#endif /* _GLIBCXX_PARALLEL_RANDOM_SHUFFLE_H */