MCKL
Monte Carlo Kernel Library
threefry_avx2_2x64.hpp
Go to the documentation of this file.
1 //============================================================================
2 // MCKL/include/mckl/random/internal/threefry_avx2_2x64.hpp
3 //----------------------------------------------------------------------------
4 // MCKL: Monte Carlo Kernel Library
5 //----------------------------------------------------------------------------
6 // Copyright (c) 2013-2018, Yan Zhou
7 // All rights reserved.
8 //
9 // Redistribution and use in source and binary forms, with or without
10 // modification, are permitted provided that the following conditions are met:
11 //
12 // Redistributions of source code must retain the above copyright notice,
13 // this list of conditions and the following disclaimer.
14 //
15 // Redistributions in binary form must reproduce the above copyright notice,
16 // this list of conditions and the following disclaimer in the documentation
17 // and/or other materials provided with the distribution.
18 //
19 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS AS IS
20 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 // POSSIBILITY OF SUCH DAMAGE.
30 //============================================================================
31 
32 #ifndef MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_HPP
33 #define MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_HPP
34 
35 #include <mckl/random/internal/common.hpp>
36 #include <mckl/random/internal/threefry_common.hpp>
37 #include <mckl/random/internal/threefry_constants.hpp>
38 #include <mckl/random/internal/threefry_generic_2x64.hpp>
39 #include <mckl/random/increment.hpp>
40 
41 MCKL_PUSH_GCC_WARNING("-Wignored-attributes")
42 
43 #define MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_KBOX(N) \
44  ymmt0 = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
45  ThreefryKBox<T, K, N>::template key<0>(par))); \
46  ymmt1 = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
47  ThreefryKBox<T, K, N>::template key<1>(par))); \
48  ymms0 = _mm256_add_epi64(ymms0, ymmt0); \
49  ymms1 = _mm256_add_epi64(ymms1, ymmt1); \
50  ymms2 = _mm256_add_epi64(ymms2, ymmt0); \
51  ymms3 = _mm256_add_epi64(ymms3, ymmt1); \
52  ymms4 = _mm256_add_epi64(ymms4, ymmt0); \
53  ymms5 = _mm256_add_epi64(ymms5, ymmt1); \
54  ymms6 = _mm256_add_epi64(ymms6, ymmt0); \
55  ymms7 = _mm256_add_epi64(ymms7, ymmt1);
56 
57 #define MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_RBOX(N) \
58  { \
59  constexpr int L0 = Constants::rotate::value[0][(N - 1) % 8]; \
60  constexpr int R0 = 64 - L0; \
61  ymms0 = _mm256_add_epi64(ymms0, ymms1); \
62  ymms2 = _mm256_add_epi64(ymms2, ymms3); \
63  ymms4 = _mm256_add_epi64(ymms4, ymms5); \
64  ymms6 = _mm256_add_epi64(ymms6, ymms7); \
65  ymmt1 = _mm256_slli_epi64(ymms1, L0); \
66  ymmt3 = _mm256_slli_epi64(ymms3, L0); \
67  ymmt5 = _mm256_slli_epi64(ymms5, L0); \
68  ymmt7 = _mm256_slli_epi64(ymms7, L0); \
69  ymms1 = _mm256_srli_epi64(ymms1, R0); \
70  ymms3 = _mm256_srli_epi64(ymms3, R0); \
71  ymms5 = _mm256_srli_epi64(ymms5, R0); \
72  ymms7 = _mm256_srli_epi64(ymms7, R0); \
73  ymmt1 = _mm256_or_si256(ymms1, ymmt1); \
74  ymmt3 = _mm256_or_si256(ymms3, ymmt3); \
75  ymmt5 = _mm256_or_si256(ymms5, ymmt5); \
76  ymmt7 = _mm256_or_si256(ymms7, ymmt7); \
77  ymms1 = _mm256_xor_si256(ymms0, ymmt1); \
78  ymms3 = _mm256_xor_si256(ymms2, ymmt3); \
79  ymms5 = _mm256_xor_si256(ymms4, ymmt5); \
80  ymms7 = _mm256_xor_si256(ymms6, ymmt7); \
81  }
82 
83 #define MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_CYCLE_4(N) \
84  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_RBOX(N * 8 + 1); \
85  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_RBOX(N * 8 + 2); \
86  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_RBOX(N * 8 + 3); \
87  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_RBOX(N * 8 + 4); \
88  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_KBOX(N * 8 + 4);
89 
90 #define MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_CYCLE_8(N) \
91  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_RBOX(N * 8 + 1); \
92  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_RBOX(N * 8 + 2); \
93  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_RBOX(N * 8 + 3); \
94  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_RBOX(N * 8 + 4); \
95  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_KBOX(N * 8 + 4); \
96  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_RBOX(N * 8 + 5); \
97  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_RBOX(N * 8 + 6); \
98  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_RBOX(N * 8 + 7); \
99  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_RBOX(N * 8 + 8); \
100  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_KBOX(N * 8 + 8);
101 
102 namespace mckl {
103 
104 namespace internal {
105 
106 template <typename T>
107 class Threefry2x64GeneratorAVX2Impl
108 {
109  static_assert(std::numeric_limits<T>::digits == 64,
110  "**Threefry2x64GeneratorAVX2Impl** used with T other than a 64-bit "
111  "unsigned integers");
112 
113  static constexpr std::size_t K = 2;
114  static constexpr std::size_t Rounds = 20;
115  using Constants = ThreefryConstants<T, 2>;
116 
117  public:
118  static void eval(
119  const void *plain, void *cipher, const std::array<T, K + 4> &par)
120  {
121  Threefry2x64GeneratorGenericImpl<T>::eval(plain, cipher, par);
122  }
123 
124  template <typename ResultType>
125  static void eval(std::array<std::uint64_t, 2> &ctr, ResultType *r,
126  const std::array<T, K + 4> &par)
127  {
128  Threefry2x64GeneratorGenericImpl<T>::eval(ctr, r, par);
129  }
130 
131  template <typename ResultType>
132  static void eval(std::array<std::uint64_t, 2> &ctr, std::size_t n,
133  ResultType *r, const std::array<T, K + 4> &par)
134  {
135  constexpr std::size_t R = sizeof(T) * K / sizeof(ResultType);
136 
137  const std::size_t n0 =
138  static_cast<std::size_t>(std::min(static_cast<std::uint64_t>(n),
139  std::numeric_limits<std::uint64_t>::max() - ctr.front()));
140 
141  eval_kernel(ctr, n0, r, par);
142  n -= n0;
143  r += n0 * R;
144 
145  if (n != 0) {
146  eval(ctr, r, par);
147  n -= 1;
148  r += R;
149  }
150 
151  eval_kernel(ctr, n, r, par);
152  }
153 
154  private:
155  template <typename ResultType>
156  static void eval_kernel(std::array<std::uint64_t, 2> &ctr, std::size_t n,
157  ResultType *r, const std::array<T, K + 4> &par)
158  {
159  constexpr std::size_t S = 8;
160  constexpr std::size_t N = sizeof(__m256i) * S / (sizeof(T) * K);
161 
162  __m256i ymmc =
163  _mm256_set_epi64x(static_cast<MCKL_INT64>(std::get<1>(ctr)),
164  static_cast<MCKL_INT64>(std::get<0>(ctr)),
165  static_cast<MCKL_INT64>(std::get<1>(ctr)),
166  static_cast<MCKL_INT64>(std::get<0>(ctr)));
167  ctr.front() += n;
168 
169  __m256i *rptr = reinterpret_cast<__m256i *>(r);
170  while (n != 0) {
171  __m256i ymms0 =
172  _mm256_add_epi64(ymmc, _mm256_set_epi64x(0, 0x02, 0, 0x01));
173  __m256i ymms1 =
174  _mm256_add_epi64(ymmc, _mm256_set_epi64x(0, 0x04, 0, 0x03));
175  __m256i ymms2 =
176  _mm256_add_epi64(ymmc, _mm256_set_epi64x(0, 0x06, 0, 0x05));
177  __m256i ymms3 =
178  _mm256_add_epi64(ymmc, _mm256_set_epi64x(0, 0x08, 0, 0x07));
179  __m256i ymms4 =
180  _mm256_add_epi64(ymmc, _mm256_set_epi64x(0, 0x0A, 0, 0x09));
181  __m256i ymms5 =
182  _mm256_add_epi64(ymmc, _mm256_set_epi64x(0, 0x0C, 0, 0x0B));
183  __m256i ymms6 =
184  _mm256_add_epi64(ymmc, _mm256_set_epi64x(0, 0x0E, 0, 0x0D));
185  __m256i ymms7 =
186  _mm256_add_epi64(ymmc, _mm256_set_epi64x(0, 0x10, 0, 0x0F));
187  ymmc = _mm256_add_epi64(ymmc, _mm256_set_epi64x(0, 0x10, 0, 0x10));
188 
189  __m256i ymmt0;
190  __m256i ymmt1;
191  __m256i ymmt2;
192  __m256i ymmt3;
193  __m256i ymmt4;
194  __m256i ymmt5;
195  __m256i ymmt6;
196  __m256i ymmt7;
197 
198  ymmt1 = ymms1;
199  ymms1 = ymms2;
200  ymms2 = ymms4;
201  ymms4 = ymmt1;
202 
203  ymmt3 = ymms3;
204  ymms3 = ymms6;
205  ymms6 = ymms5;
206  ymms5 = ymmt3;
207 
208  ymmt0 = _mm256_unpacklo_epi64(ymms0, ymms1);
209  ymmt1 = _mm256_unpacklo_epi64(ymms2, ymms3);
210  ymmt2 = _mm256_unpackhi_epi64(ymms0, ymms1);
211  ymmt3 = _mm256_unpackhi_epi64(ymms2, ymms3);
212  ymmt4 = _mm256_unpacklo_epi64(ymms4, ymms5);
213  ymmt5 = _mm256_unpacklo_epi64(ymms6, ymms7);
214  ymmt6 = _mm256_unpackhi_epi64(ymms4, ymms5);
215  ymmt7 = _mm256_unpackhi_epi64(ymms6, ymms7);
216  ymms0 = _mm256_permute2x128_si256(ymmt0, ymmt1, 0x20);
217  ymms1 = _mm256_permute2x128_si256(ymmt2, ymmt3, 0x20);
218  ymms2 = _mm256_permute2x128_si256(ymmt0, ymmt1, 0x31);
219  ymms3 = _mm256_permute2x128_si256(ymmt2, ymmt3, 0x31);
220  ymms4 = _mm256_permute2x128_si256(ymmt4, ymmt5, 0x20);
221  ymms5 = _mm256_permute2x128_si256(ymmt6, ymmt7, 0x20);
222  ymms6 = _mm256_permute2x128_si256(ymmt4, ymmt5, 0x31);
223  ymms7 = _mm256_permute2x128_si256(ymmt6, ymmt7, 0x31);
224 
225  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_KBOX(0)
226  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_CYCLE_8(0)
227  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_CYCLE_8(1)
228  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_CYCLE_4(2)
229 
230  ymmt0 = _mm256_unpacklo_epi64(ymms0, ymms1);
231  ymmt1 = _mm256_unpacklo_epi64(ymms2, ymms3);
232  ymmt2 = _mm256_unpackhi_epi64(ymms0, ymms1);
233  ymmt3 = _mm256_unpackhi_epi64(ymms2, ymms3);
234  ymmt4 = _mm256_unpacklo_epi64(ymms4, ymms5);
235  ymmt5 = _mm256_unpacklo_epi64(ymms6, ymms7);
236  ymmt6 = _mm256_unpackhi_epi64(ymms4, ymms5);
237  ymmt7 = _mm256_unpackhi_epi64(ymms6, ymms7);
238  ymms0 = _mm256_permute2x128_si256(ymmt0, ymmt1, 0x20);
239  ymms1 = _mm256_permute2x128_si256(ymmt2, ymmt3, 0x20);
240  ymms2 = _mm256_permute2x128_si256(ymmt0, ymmt1, 0x31);
241  ymms3 = _mm256_permute2x128_si256(ymmt2, ymmt3, 0x31);
242  ymms4 = _mm256_permute2x128_si256(ymmt4, ymmt5, 0x20);
243  ymms5 = _mm256_permute2x128_si256(ymmt6, ymmt7, 0x20);
244  ymms6 = _mm256_permute2x128_si256(ymmt4, ymmt5, 0x31);
245  ymms7 = _mm256_permute2x128_si256(ymmt6, ymmt7, 0x31);
246 
247  ymmt2 = ymms2;
248  ymms2 = ymms1;
249  ymms1 = ymms4;
250  ymms4 = ymmt2;
251 
252  ymmt6 = ymms6;
253  ymms6 = ymms3;
254  ymms3 = ymms5;
255  ymms5 = ymmt6;
256 
257  if (n >= N) {
258  n -= N;
259  _mm256_storeu_si256(rptr++, ymms0);
260  _mm256_storeu_si256(rptr++, ymms1);
261  _mm256_storeu_si256(rptr++, ymms2);
262  _mm256_storeu_si256(rptr++, ymms3);
263  _mm256_storeu_si256(rptr++, ymms4);
264  _mm256_storeu_si256(rptr++, ymms5);
265  _mm256_storeu_si256(rptr++, ymms6);
266  _mm256_storeu_si256(rptr++, ymms7);
267  } else {
268  std::array<__m256i, S> s;
269  std::get<0>(s) = ymms0;
270  std::get<1>(s) = ymms1;
271  std::get<2>(s) = ymms2;
272  std::get<3>(s) = ymms3;
273  std::get<4>(s) = ymms4;
274  std::get<5>(s) = ymms5;
275  std::get<6>(s) = ymms6;
276  std::get<7>(s) = ymms7;
277  std::memcpy(rptr, s.data(), n * sizeof(T) * K);
278  break;
279  }
280  }
281  }
282 }; // class Threefry2x64GeneratorAVX2Impl
283 
284 } // namespace internal
285 
286 } // namespace mckl
287 
288 MCKL_POP_GCC_WARNING
289 
290 #endif // MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_HPP
MCKL_PUSH_GCC_WARNING
#define MCKL_PUSH_GCC_WARNING(warning)
Definition: compiler.h:78
threefry_generic_2x64.hpp
mckl::internal::Threefry2x64GeneratorAVX2Impl::eval
static void eval(std::array< std::uint64_t, 2 > &ctr, std::size_t n, ResultType *r, const std::array< T, K+4 > &par)
Definition: threefry_avx2_2x64.hpp:132
MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_KBOX
#define MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_KBOX(N)
Definition: threefry_avx2_2x64.hpp:43
MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_CYCLE_8
#define MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_CYCLE_8(N)
Definition: threefry_avx2_2x64.hpp:90
mckl::ThreefryConstants
Default Threefry constants.
Definition: threefry_constants.hpp:147
mckl::internal::Threefry2x64GeneratorAVX2Impl::eval
static void eval(const void *plain, void *cipher, const std::array< T, K+4 > &par)
Definition: threefry_avx2_2x64.hpp:118
threefry_constants.hpp
mckl::internal::Threefry2x64GeneratorAVX2Impl::eval
static void eval(std::array< std::uint64_t, 2 > &ctr, ResultType *r, const std::array< T, K+4 > &par)
Definition: threefry_avx2_2x64.hpp:125
mckl
Definition: mcmc.hpp:40
threefry_common.hpp
MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_CYCLE_4
#define MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_2X64_CYCLE_4(N)
Definition: threefry_avx2_2x64.hpp:83
mckl::internal::Threefry2x64GeneratorGenericImpl::eval
static void eval(const void *plain, void *cipher, const std::array< T, K+4 > &par)
Definition: threefry_generic_2x64.hpp:98
MCKL_POP_GCC_WARNING
#define MCKL_POP_GCC_WARNING
Definition: compiler.h:79
Generated by   doxygen 1.8.13