Basic shuffle implementation
Resolve linux compilation, add benchmark
Add shuffle_copy, tidy up and comment
doxygen
Address review comments
Silence warnings
Guard c++11
MODIFY
dependencies/cub
dependencies/cub
-Subproject commit 8bdfee73970eac30c2c776956b5aab3fc87b604e
+Subproject commit f7ad39d345ffe970ae586245eeca6a38581a95a9

MODIFY
internal/benchmark/bench.cu
internal/benchmark/bench.cu
#include <thrust/scan.h>
#include <thrust/version.h>
+#if THRUST_CPP_DIALECT >= 2011
+#include <thrust/shuffle.h>
+#include <random>
+#endif
+
#include <algorithm>
#include <numeric>
@@ -710,6 +715,22 @@ struct copy_trial_base : trial_base<TrialKind>
}
};
+#if THRUST_CPP_DIALECT >= 2011
+template <typename Container, typename TrialKind = regular_trial>
+struct shuffle_trial_base : trial_base<TrialKind>
+{
+ Container input;
+ std::default_random_engine g;
+
+ void setup(uint64_t elements)
+ {
+ input.resize(elements);
+
+ randomize(input);
+ }
+};
+#endif
+
///////////////////////////////////////////////////////////////////////////////
template <typename T>
@@ -905,6 +926,35 @@ struct copy_tester
#endif
};
+#if THRUST_CPP_DIALECT >= 2011
+template <typename T>
+struct shuffle_tester
+{
+ static char const* test_name() { return "shuffle"; }
+
+ struct std_trial : shuffle_trial_base<std::vector<T>, baseline_trial>
+ {
+ void operator()()
+ {
+ std::shuffle(this->input.begin(), this->input.end(), this->g);
+ }
+ };
+
+ struct thrust_trial : shuffle_trial_base<thrust::device_vector<T> >
+ {
+ void operator()()
+ {
+ thrust::shuffle(this->input.begin(), this->input.end(), this->g);
+ #if THRUST_DEVICE_SYSTEM == THRUST_DEVICE_SYSTEM_CUDA
+ cudaError_t err = cudaDeviceSynchronize();
+ if (err != cudaSuccess)
+ throw thrust::error_code(err, thrust::cuda_category());
+ #endif
+ }
+ };
+};
+#endif
+
///////////////////////////////////////////////////////////////////////////////
template <
@@ -956,6 +1006,16 @@ void run_core_primitives_experiments_for_type()
, BaselineTrials
, RegularTrials
>::run_experiment();
+
+#if THRUST_CPP_DIALECT >= 2011
+ experiment_driver<
+ shuffle_tester
+ , ElementMetaType
+ , Elements / sizeof(typename ElementMetaType::type)
+ , BaselineTrials
+ , RegularTrials
+ >::run_experiment();
+#endif
}
///////////////////////////////////////////////////////////////////////////////

ADD
testing/shuffle.cu
+#include <thrust/detail/config.h>
+
+#if THRUST_CPP_DIALECT >= 2011
+#include <thrust/random.h>
+#include <thrust/shuffle.h>
+#include <thrust/sort.h>
+#include <unittest/unittest.h>
+#include <map>
+
+template <typename Vector>
+void TestShuffleSimple() {
+ Vector data(5);
+ data[0] = 0;
+ data[1] = 1;
+ data[2] = 2;
+ data[3] = 3;
+ data[4] = 4;
+ Vector shuffled(data.begin(), data.end());
+ thrust::default_random_engine g(2);
+ thrust::shuffle(shuffled.begin(), shuffled.end(), g);
+ thrust::sort(shuffled.begin(), shuffled.end());
+ // Check all of our data is present
+ // This only tests for strange conditions like duplicated elements
+ ASSERT_EQUAL(shuffled, data);
+}
+DECLARE_VECTOR_UNITTEST(TestShuffleSimple);
+
+template <typename Vector>
+void TestShuffleCopySimple() {
+ Vector data(5);
+ data[0] = 0;
+ data[1] = 1;
+ data[2] = 2;
+ data[3] = 3;
+ data[4] = 4;
+ Vector shuffled(5);
+ thrust::default_random_engine g(2);
+ thrust::shuffle_copy(data.begin(), data.end(), shuffled.begin(), g);
+ g.seed(2);
+ thrust::shuffle(data.begin(), data.end(), g);
+ ASSERT_EQUAL(shuffled, data);
+}
+DECLARE_VECTOR_UNITTEST(TestShuffleCopySimple);
+
+template <typename T>
+void TestHostDeviceIdentical(size_t m) {
+ thrust::host_vector<T> host_result(m);
+ thrust::host_vector<T> device_result(m);
+ thrust::sequence(host_result.begin(), host_result.end(), 0llu);
+ thrust::sequence(device_result.begin(), device_result.end(), 0llu);
+
+ thrust::default_random_engine host_g(183);
+ thrust::default_random_engine device_g(183);
+
+ thrust::shuffle(host_result.begin(), host_result.end(), host_g);
+ thrust::shuffle(device_result.begin(), device_result.end(), device_g);
+
+ ASSERT_EQUAL(device_result, host_result);
+}
+DECLARE_VARIABLE_UNITTEST(TestHostDeviceIdentical);
+
+// Individual input keys should be permuted to output locations with uniform
+// probability. Perform chi-squared test with confidence 99.9%.
+template <typename Vector>
+void TestShuffleKeyPosition() {
+ typedef typename Vector::value_type T;
+ size_t m = 20;
+ size_t num_samples = 100;
+ thrust::host_vector<size_t> index_sum(m, 0);
+ thrust::host_vector<T> sequence(m);
+ thrust::sequence(sequence.begin(), sequence.end(), T(0));
+
+ for (size_t i = 0; i < num_samples; i++) {
+ Vector shuffled(sequence.begin(), sequence.end());
+ thrust::default_random_engine g(i);
+ thrust::shuffle(shuffled.begin(), shuffled.end(), g);
+ thrust::host_vector<T> tmp(shuffled.begin(), shuffled.end());
+
+ for (auto j = 0ull; j < m; j++) {
+ index_sum[tmp[j]] += j;
+ }
+ }
+ double expected_average_position = static_cast<double>(m - 1) / 2;
+ double chi_squared = 0.0;
+ for (auto j = 0ull; j < m; j++) {
+ double average_position = static_cast<double>(index_sum[j]) / num_samples;
+ chi_squared += std::pow(expected_average_position - average_position, 2) /
+ expected_average_position;
+ }
+ // Tabulated chi-squared critical value for m-1=19 degrees of freedom
+ // and 99.9% confidence
+ double confidence_threshold = 43.82;
+ ASSERT_LESS(chi_squared, confidence_threshold);
+}
+DECLARE_INTEGRAL_VECTOR_UNITTEST(TestShuffleKeyPosition);
+
+struct vector_compare {
+ template <typename VectorT>
+ bool operator()(const VectorT& a, const VectorT& b) const {
+ for (auto i = 0ull; i < a.size(); i++) {
+ if (a[i] < b[i]) return true;
+ if (a[i] > b[i]) return false;
+ }
+ return false;
+ }
+};
+
+// Brute force check permutations are uniformly distributed on small input
+// Uses a chi-squared test indicating 99% confidence the output is uniformly
+// random
+template <typename Vector>
+void TestShuffleUniformPermutation() {
+ typedef typename Vector::value_type T;
+
+ size_t m = 5;
+ size_t num_samples = 1000;
+ size_t total_permutations = 1 * 2 * 3 * 4 * 5;
+ std::map<thrust::host_vector<T>, size_t, vector_compare> permutation_counts;
+ Vector sequence(m);
+ thrust::sequence(sequence.begin(), sequence.end(), T(0));
+ thrust::default_random_engine g(17);
+ for (auto i = 0ull; i < num_samples; i++) {
+ thrust::shuffle(sequence.begin(), sequence.end(), g);
+ thrust::host_vector<T> tmp(sequence.begin(), sequence.end());
+ permutation_counts[tmp]++;
+ }
+
+ ASSERT_EQUAL(permutation_counts.size(), total_permutations);
+
+ double chi_squared = 0.0;
+ double expected_count = static_cast<double>(num_samples) / total_permutations;
+ for (auto kv : permutation_counts) {
+ chi_squared += std::pow(expected_count - kv.second, 2) / expected_count;
+ }
+ // Tabulated chi-squared critical value for 119 degrees of freedom (5! - 1)
+ // and 99% confidence
+ double confidence_threshold = 157.8;
+ ASSERT_LESS(chi_squared, confidence_threshold);
+}
+DECLARE_VECTOR_UNITTEST(TestShuffleUniformPermutation);
+#endif

ADD
thrust/detail/shuffle.inl
+/*
+ * Copyright 2008-2020 NVIDIA Corporation
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/*! \file shuffle.inl
+ * \brief Inline file for shuffle.h.
+ */
+
+#include <thrust/detail/config.h>
+#include <thrust/detail/cpp11_required.h>
+
+#if THRUST_CPP_DIALECT >= 2011
+
+#include <thrust/iterator/iterator_traits.h>
+#include <thrust/shuffle.h>
+#include <thrust/system/detail/generic/select_system.h>
+#include <thrust/system/detail/generic/shuffle.h>
+
+namespace thrust {
+
+__thrust_exec_check_disable__
+template <typename DerivedPolicy, typename RandomIterator, typename URBG>
+__host__ __device__ void shuffle(
+ const thrust::detail::execution_policy_base<DerivedPolicy>& exec,
+ RandomIterator first, RandomIterator last, URBG&& g) {
+ using thrust::system::detail::generic::shuffle;
+ return shuffle(
+ thrust::detail::derived_cast(thrust::detail::strip_const(exec)),
+ first, last, g);
+}
+
+template <typename RandomIterator, typename URBG>
+__host__ __device__ void shuffle(RandomIterator first, RandomIterator last,
+ URBG&& g) {
+ using thrust::system::detail::generic::select_system;
+
+ typedef typename thrust::iterator_system<RandomIterator>::type System;
+ System system;
+
+ return thrust::shuffle(select_system(system), first, last, g);
+}
+
+__thrust_exec_check_disable__
+template <typename DerivedPolicy, typename RandomIterator,
+ typename OutputIterator, typename URBG>
+__host__ __device__ void shuffle_copy(
+ const thrust::detail::execution_policy_base<DerivedPolicy>& exec,
+ RandomIterator first, RandomIterator last, OutputIterator result,
+ URBG&& g) {
+ using thrust::system::detail::generic::shuffle_copy;
+ return shuffle_copy(
+ thrust::detail::derived_cast(thrust::detail::strip_const(exec)),
+ first, last, result, g);
+}
+
+template <typename RandomIterator, typename OutputIterator, typename URBG>
+__host__ __device__ void shuffle_copy(RandomIterator first, RandomIterator last,
+ OutputIterator result, URBG&& g) {
+ using thrust::system::detail::generic::select_system;
+
+ typedef typename thrust::iterator_system<RandomIterator>::type System1;
+ typedef typename thrust::iterator_system<OutputIterator>::type System2;
+
+ System1 system1;
+ System2 system2;
+
+ return thrust::shuffle_copy(select_system(system1, system2), first, last,
+ result, g);
+}
+
+} // namespace thrust
+
+#endif

ADD
thrust/shuffle.h
+/*
+ * Copyright 2008-2020 NVIDIA Corporation
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/*! \file shuffle.h
+ * \brief Reorders range by a uniform random permutation
+ */
+
+#pragma once
+
+#include <thrust/detail/config.h>
+#include <thrust/detail/cpp11_required.h>
+
+#if THRUST_CPP_DIALECT >= 2011
+
+#include <thrust/detail/config.h>
+#include <thrust/detail/execution_policy.h>
+
+namespace thrust {
+
+/*! \addtogroup reordering
+* \ingroup algorithms
+*
+* \addtogroup shuffling
+* \ingroup reordering
+* \{
+*/
+
+
+/*! \p shuffle reorders the elements <tt>[first, last)</tt> by a uniform pseudorandom permutation, defined by
+ * random engine \p g.
+ *
+ * The algorithm's execution is parallelized as determined by \p exec.
+ *
+ * \param exec The execution policy to use for parallelization.
+ * \param first The beginning of the sequence to shuffle.
+ * \param last The end of the sequence to shuffle.
+ * \param g A UniformRandomBitGenerator
+ *
+ * \tparam DerivedPolicy The name of the derived execution policy.
+ * \tparam RandomIterator is a random access iterator
+ * \tparam URBG is a uniform random bit generator
+ *
+ * The following code snippet demonstrates how to use \p shuffle to create a random permutation
+ * using the \p thrust::host execution policy for parallelization:
+ *
+ * \code
+ * #include <thrust/shuffle.h>
+ * #include <thrust/random.h>
+ * #include <thrust/execution_policy.h>
+ * int A[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
+ * const int N = sizeof(A)/sizeof(int);
+ * thrust::default_random_engine g;
+ * thrust::shuffle(thrust::host, A, A + N, g);
+ * // A is now {6, 5, 8, 7, 2, 1, 4, 3, 10, 9}
+ * \endcode
+ *
+ * \see \p shuffle_copy
+ */
+template <typename DerivedPolicy, typename RandomIterator, typename URBG>
+__host__ __device__ void shuffle(
+ const thrust::detail::execution_policy_base<DerivedPolicy>& exec,
+ RandomIterator first, RandomIterator last, URBG&& g);
+
+/*! \p shuffle reorders the elements <tt>[first, last)</tt> by a uniform pseudorandom permutation, defined by
+ * random engine \p g.
+ *
+ * \param first The beginning of the sequence to shuffle.
+ * \param last The end of the sequence to shuffle.
+ * \param g A UniformRandomBitGenerator
+ *
+ * \tparam RandomIterator is a random access iterator
+ * \tparam URBG is a uniform random bit generator
+ *
+ * The following code snippet demonstrates how to use \p shuffle to create a random permutation.
+ *
+ * \code
+ * #include <thrust/shuffle.h>
+ * #include <thrust/random.h>
+ * int A[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
+ * const int N = sizeof(A)/sizeof(int);
+ * thrust::default_random_engine g;
+ * thrust::shuffle(A, A + N, g);
+ * // A is now {6, 5, 8, 7, 2, 1, 4, 3, 10, 9}
+ * \endcode
+ *
+ * \see \p shuffle_copy
+ */
+template <typename RandomIterator, typename URBG>
+__host__ __device__ void shuffle(RandomIterator first, RandomIterator last,
+ URBG&& g);
+
+/*! shuffle_copy differs from shuffle only in that the reordered sequence is written to different output sequences, rather than in place.
+ * \p shuffle_copy reorders the elements <tt>[first, last)</tt> by a uniform pseudorandom permutation, defined by
+ * random engine \p g.
+ *
+ * The algorithm's execution is parallelized as determined by \p exec.
+
+ * \param exec The execution policy to use for parallelization.
+ * \param first The beginning of the sequence to shuffle.
+ * \param last The end of the sequence to shuffle.
+ * \param result Destination of shuffled sequence
+ * \param g A UniformRandomBitGenerator
+ *
+ * \tparam DerivedPolicy The name of the derived execution policy.
+ * \tparam RandomIterator is a random access iterator
+ * \tparam OutputIterator is a model of <a href="http://www.sgi.com/tech/stl/OutputIterator.html">Output Iterator</a>.
+ * \tparam URBG is a uniform random bit generator
+ *
+ * The following code snippet demonstrates how to use \p shuffle_copy to create a random permutation.
+ *
+ * \code
+ * #include <thrust/shuffle.h>
+ * #include <thrust/random.h>
+ * #include <thrust/execution_policy.h>
+ * int A[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
+ * int result[10];
+ * const int N = sizeof(A)/sizeof(int);
+ * thrust::default_random_engine g;
+ * thrust::shuffle_copy(thrust::host, A, A + N, result, g);
+ * // result is now {6, 5, 8, 7, 2, 1, 4, 3, 10, 9}
+ * \endcode
+ *
+ * \see \p shuffle
+ */
+template <typename DerivedPolicy, typename RandomIterator,
+ typename OutputIterator, typename URBG>
+__host__ __device__ void shuffle_copy(
+ const thrust::detail::execution_policy_base<DerivedPolicy>& exec,
+ RandomIterator first, RandomIterator last, OutputIterator result, URBG&& g);
+
+/*! shuffle_copy differs from shuffle only in that the reordered sequence is written to different output sequences, rather than in place.
+ *\p shuffle_copy reorders the elements <tt>[first, last)</tt> by a uniform pseudorandom permutation, defined by
+ * random engine \p g.
+ *
+ * \param first The beginning of the sequence to shuffle.
+ * \param last The end of the sequence to shuffle.
+ * \param result Destination of shuffled sequence
+ * \param g A UniformRandomBitGenerator
+ *
+ * \tparam RandomIterator is a random access iterator
+ * \tparam OutputIterator is a model of <a href="http://www.sgi.com/tech/stl/OutputIterator.html">Output Iterator</a>.
+ * \tparam URBG is a uniform random bit generator
+ *
+ * The following code snippet demonstrates how to use \p shuffle_copy to create a random permutation.
+ *
+ * \code
+ * #include <thrust/shuffle.h>
+ * #include <thrust/random.h>
+ * int A[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
+ * int result[10];
+ * const int N = sizeof(A)/sizeof(int);
+ * thrust::default_random_engine g;
+ * thrust::shuffle_copy(A, A + N, result, g);
+ * // result is now {6, 5, 8, 7, 2, 1, 4, 3, 10, 9}
+ * \endcode
+ *
+ * \see \p shuffle
+ */
+template <typename RandomIterator, typename OutputIterator, typename URBG>
+__host__ __device__ void shuffle_copy(RandomIterator first, RandomIterator last,
+ OutputIterator result, URBG&& g);
+
+} // namespace thrust
+
+#include <thrust/detail/shuffle.inl>
+#endif

ADD
thrust/system/detail/generic/shuffle.h
+/*
+ * Copyright 2008-2020 NVIDIA Corporation
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/*! \file shuffle.h
+ * \brief Generic implementations of shuffle functions.
+ */
+
+#pragma once
+
+
+#include <thrust/detail/config.h>
+#include <thrust/detail/cpp11_required.h>
+
+#if THRUST_CPP_DIALECT >= 2011
+
+#include <thrust/system/detail/generic/tag.h>
+
+namespace thrust {
+namespace system {
+namespace detail {
+namespace generic {
+
+template <typename ExecutionPolicy, typename RandomIterator, typename URBG>
+__host__ __device__ void shuffle(
+ thrust::execution_policy<ExecutionPolicy>& exec, RandomIterator first,
+ RandomIterator last, URBG&& g);
+
+template <typename ExecutionPolicy, typename RandomIterator,
+ typename OutputIterator, typename URBG>
+__host__ __device__ void shuffle_copy(
+ thrust::execution_policy<ExecutionPolicy>& exec, RandomIterator first,
+ RandomIterator last, OutputIterator result, URBG&& g);
+
+} // end namespace generic
+} // end namespace detail
+} // end namespace system
+} // end namespace thrust
+
+#include <thrust/system/detail/generic/shuffle.inl>
+
+#endif

ADD
thrust/system/detail/generic/shuffle.inl
+/*
+ * Copyright 2008-20120 NVIDIA Corporation
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include <thrust/detail/config.h>
+#include <thrust/detail/cpp11_required.h>
+
+#if THRUST_CPP_DIALECT >= 2011
+
+#include <thrust/detail/temporary_array.h>
+#include <thrust/iterator/discard_iterator.h>
+#include <thrust/iterator/transform_iterator.h>
+#include <thrust/iterator/transform_output_iterator.h>
+#include <thrust/random.h>
+#include <thrust/scan.h>
+#include <thrust/system/detail/generic/shuffle.h>
+
+namespace thrust {
+namespace system {
+namespace detail {
+namespace generic {
+
+// An implementation of a Feistel cipher for operating on 64 bit keys
+class feistel_bijection {
+ private:
+ struct round_state {
+ uint32_t left;
+ uint32_t right;
+ };
+
+ public:
+ template <class URBG>
+ __host__ __device__ feistel_bijection(uint64_t m, URBG&& g) {
+ uint64_t total_bits = get_cipher_bits(m);
+ // Half bits rounded down
+ left_side_bits = total_bits / 2;
+ left_side_mask = (1ull << left_side_bits) - 1;
+ // Half the bits rounded up
+ right_side_bits = total_bits - left_side_bits;
+ right_side_mask = (1ull << right_side_bits) - 1;
+
+ for (uint64_t i = 0; i < num_rounds; i++) {
+ key[i] = g();
+ }
+ }
+
+ __host__ __device__ uint64_t nearest_power_of_two() const {
+ return 1ull << (left_side_bits + right_side_bits);
+ }
+ __host__ __device__ uint64_t operator()(const uint64_t val) const {
+ // Extract the right and left sides of the input
+ uint32_t left = (uint32_t)(val >> right_side_bits);
+ uint32_t right = (uint32_t)(val & right_side_mask);
+ round_state state = {left, right};
+
+ for (uint64_t i = 0; i < num_rounds; i++) {
+ state = do_round(state, i);
+ }
+
+ // Check we have the correct number of bits on each side
+ assert((state.left >> left_side_bits) == 0);
+ assert((state.right >> right_side_bits) == 0);
+
+ // Combine the left and right sides together to get result
+ return state.left << right_side_bits | state.right;
+ }
+
+ private:
+ // Find the nearest power of two
+ __host__ __device__ uint64_t get_cipher_bits(uint64_t m) {
+ uint64_t i = 0;
+ while (m != 0) {
+ i++;
+ m >>= 1;
+ }
+ return i;
+ }
+
+ // Round function, a 'pseudorandom function' whos output is indistinguishable
+ // from random for each key value input. This is not cryptographically secure
+ // but sufficient for generating permutations. We hash the value with the
+ // tau88 engine and combine it with the random bits of the key (provided by
+ // the user-defined engine).
+ __host__ __device__ uint32_t round_function(uint64_t value,
+ const uint64_t key) const {
+ uint64_t value_hash = thrust::random::taus88(value)();
+ return (value_hash ^ key) & left_side_mask;
+ }
+
+ __host__ __device__ round_state do_round(const round_state state,
+ const uint64_t round) const {
+ const uint32_t new_left = state.right & left_side_mask;
+ const uint32_t round_function_res =
+ state.left ^ round_function(state.right, key[round]);
+ if (right_side_bits != left_side_bits) {
+ // Upper bit of the old right becomes lower bit of new right if we have
+ // odd length feistel
+ const uint32_t new_right =
+ (round_function_res << 1ull) | state.right >> left_side_bits;
+ return {new_left, new_right};
+ }
+ return {new_left, round_function_res};
+ }
+
+ static const uint64_t num_rounds = 8;
+ uint64_t right_side_bits;
+ uint64_t left_side_bits;
+ uint64_t right_side_mask;
+ uint64_t left_side_mask;
+ uint64_t key[num_rounds];
+};
+
+struct key_flag_tuple {
+ uint64_t key;
+ uint64_t flag;
+};
+
+// scan only flags
+struct key_flag_scan_op {
+ __host__ __device__ key_flag_tuple operator()(const key_flag_tuple& a,
+ const key_flag_tuple& b) {
+ return {b.key, a.flag + b.flag};
+ }
+};
+
+struct construct_key_flag_op {
+ uint64_t m;
+ feistel_bijection bijection;
+ __host__ __device__ construct_key_flag_op(uint64_t m,
+ feistel_bijection bijection)
+ : m(m), bijection(bijection) {}
+ __host__ __device__ key_flag_tuple operator()(uint64_t idx) {
+ auto gather_key = bijection(idx);
+ return key_flag_tuple{gather_key, (gather_key < m) ? 1ull : 0ull};
+ }
+};
+
+template <typename InputIterT, typename OutputIterT>
+struct write_output_op {
+ uint64_t m;
+ InputIterT in;
+ OutputIterT out;
+ // flag contains inclusive scan of valid keys
+ // perform gather using valid keys
+ __thrust_exec_check_disable__
+ __host__ __device__ size_t operator()(key_flag_tuple x) {
+ if (x.key < m) {
+ // -1 because inclusive scan
+ out[x.flag - 1] = in[x.key];
+ }
+ return 0; // Discarded
+ }
+};
+
+template <typename ExecutionPolicy, typename RandomIterator, typename URBG>
+__host__ __device__ void shuffle(
+ thrust::execution_policy<ExecutionPolicy>& exec, RandomIterator first,
+ RandomIterator last, URBG&& g) {
+ typedef
+ typename thrust::iterator_traits<RandomIterator>::value_type InputType;
+
+ // copy input to temp buffer
+ thrust::detail::temporary_array<InputType, ExecutionPolicy> temp(exec, first,
+ last);
+ thrust::shuffle_copy(exec, temp.begin(), temp.end(), first, g);
+}
+
+template <typename ExecutionPolicy, typename RandomIterator,
+ typename OutputIterator, typename URBG>
+__host__ __device__ void shuffle_copy(
+ thrust::execution_policy<ExecutionPolicy>& exec, RandomIterator first,
+ RandomIterator last, OutputIterator result, URBG&& g) {
+ // m is the length of the input
+ // we have an available bijection of length n via a feistel cipher
+ size_t m = last - first;
+ feistel_bijection bijection(m, g);
+ uint64_t n = bijection.nearest_power_of_two();
+
+ // perform stream compaction over length n bijection to get length m
+ // pseudorandom bijection over the original input
+ thrust::counting_iterator<uint64_t> indices(0);
+ thrust::transform_iterator<construct_key_flag_op, decltype(indices),
+ key_flag_tuple>
+ key_flag_it(indices, construct_key_flag_op(m, bijection));
+ write_output_op<RandomIterator, decltype(result)> write_functor{m, first,
+ result};
+ auto gather_output_it = thrust::make_transform_output_iterator(
+ thrust::discard_iterator<size_t>(), write_functor);
+ // the feistel_bijection outputs a stream of permuted indices in range [0,n)
+ // flag each value < m and compact it, so we have a set of permuted indices in
+ // range [0,m) each thread gathers an input element according to its
+ // pseudorandom permuted index
+ thrust::inclusive_scan(exec, key_flag_it, key_flag_it + n, gather_output_it,
+ key_flag_scan_op());
+}
+
+} // end namespace generic
+} // end namespace detail
+} // end namespace system
+} // end namespace thrust
+#endif

Basic shuffle implementation

Resolve linux compilation, add benchmark

Add shuffle_copy, tidy up and comment

doxygen

Address review comments

Silence warnings

Guard c++11