MCKL
Monte Carlo Kernel Library
threefry_avx2_16x64.hpp
Go to the documentation of this file.
1 //============================================================================
2 // MCKL/include/mckl/random/internal/threefry_avx2_16x64.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_16X64_HPP
33 #define MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_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_16x64.hpp>
39 #include <mckl/random/increment.hpp>
40 
41 MCKL_PUSH_GCC_WARNING("-Wignored-attributes")
42 
43 #define MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_KBOX(N) \
44  ymmt0 = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
45  ThreefryKBox<T, K, N>::template key<0x0>(par))); \
46  ymmt1 = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
47  ThreefryKBox<T, K, N>::template key<0x1>(par))); \
48  ymmt2 = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
49  ThreefryKBox<T, K, N>::template key<0x2>(par))); \
50  ymmt3 = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
51  ThreefryKBox<T, K, N>::template key<0x3>(par))); \
52  ymmt4 = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
53  ThreefryKBox<T, K, N>::template key<0x4>(par))); \
54  ymmt5 = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
55  ThreefryKBox<T, K, N>::template key<0x5>(par))); \
56  ymmt6 = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
57  ThreefryKBox<T, K, N>::template key<0x6>(par))); \
58  ymmt7 = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
59  ThreefryKBox<T, K, N>::template key<0x7>(par))); \
60  ymmt8 = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
61  ThreefryKBox<T, K, N>::template key<0x8>(par))); \
62  ymmt9 = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
63  ThreefryKBox<T, K, N>::template key<0x9>(par))); \
64  ymmtA = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
65  ThreefryKBox<T, K, N>::template key<0xA>(par))); \
66  ymmtB = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
67  ThreefryKBox<T, K, N>::template key<0xB>(par))); \
68  ymmtC = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
69  ThreefryKBox<T, K, N>::template key<0xC>(par))); \
70  ymmtD = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
71  ThreefryKBox<T, K, N>::template key<0xD>(par))); \
72  ymmtE = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
73  ThreefryKBox<T, K, N>::template key<0xE>(par))); \
74  ymmtF = _mm256_set1_epi64x(static_cast<MCKL_INT64>( \
75  ThreefryKBox<T, K, N>::template key<0xF>(par))); \
76  ymms0 = _mm256_add_epi64(ymms0, ymmt0); \
77  ymms1 = _mm256_add_epi64(ymms1, ymmt1); \
78  ymms2 = _mm256_add_epi64(ymms2, ymmt2); \
79  ymms3 = _mm256_add_epi64(ymms3, ymmt3); \
80  ymms4 = _mm256_add_epi64(ymms4, ymmt4); \
81  ymms5 = _mm256_add_epi64(ymms5, ymmt5); \
82  ymms6 = _mm256_add_epi64(ymms6, ymmt6); \
83  ymms7 = _mm256_add_epi64(ymms7, ymmt7); \
84  ymms8 = _mm256_add_epi64(ymms8, ymmt8); \
85  ymms9 = _mm256_add_epi64(ymms9, ymmt9); \
86  ymmsA = _mm256_add_epi64(ymmsA, ymmtA); \
87  ymmsB = _mm256_add_epi64(ymmsB, ymmtB); \
88  ymmsC = _mm256_add_epi64(ymmsC, ymmtC); \
89  ymmsD = _mm256_add_epi64(ymmsD, ymmtD); \
90  ymmsE = _mm256_add_epi64(ymmsE, ymmtE); \
91  ymmsF = _mm256_add_epi64(ymmsF, ymmtF);
92 
93 #define MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_RBOX(N) \
94  { \
95  constexpr int L0 = Constants::rotate::value[0][(N - 1) % 8]; \
96  constexpr int L1 = Constants::rotate::value[1][(N - 1) % 8]; \
97  constexpr int L2 = Constants::rotate::value[2][(N - 1) % 8]; \
98  constexpr int L3 = Constants::rotate::value[3][(N - 1) % 8]; \
99  constexpr int L4 = Constants::rotate::value[4][(N - 1) % 8]; \
100  constexpr int L5 = Constants::rotate::value[5][(N - 1) % 8]; \
101  constexpr int L6 = Constants::rotate::value[6][(N - 1) % 8]; \
102  constexpr int L7 = Constants::rotate::value[7][(N - 1) % 8]; \
103  constexpr int R0 = 64 - L0; \
104  constexpr int R1 = 64 - L1; \
105  constexpr int R2 = 64 - L2; \
106  constexpr int R3 = 64 - L3; \
107  constexpr int R4 = 64 - L4; \
108  constexpr int R5 = 64 - L5; \
109  constexpr int R6 = 64 - L6; \
110  constexpr int R7 = 64 - L7; \
111  ymmt0 = _mm256_add_epi64(ymms0, ymms1); \
112  ymmt2 = _mm256_add_epi64(ymms2, ymms3); \
113  ymmt4 = _mm256_add_epi64(ymms4, ymms5); \
114  ymmt6 = _mm256_add_epi64(ymms6, ymms7); \
115  ymmt8 = _mm256_add_epi64(ymms8, ymms9); \
116  ymmtA = _mm256_add_epi64(ymmsA, ymmsB); \
117  ymmtC = _mm256_add_epi64(ymmsC, ymmsD); \
118  ymmtE = _mm256_add_epi64(ymmsE, ymmsF); \
119  ymmt1 = _mm256_slli_epi64(ymms1, L0); \
120  ymmt3 = _mm256_slli_epi64(ymms3, L1); \
121  ymmt5 = _mm256_slli_epi64(ymms5, L2); \
122  ymmt7 = _mm256_slli_epi64(ymms7, L3); \
123  ymmt9 = _mm256_slli_epi64(ymms9, L4); \
124  ymmtB = _mm256_slli_epi64(ymmsB, L5); \
125  ymmtD = _mm256_slli_epi64(ymmsD, L6); \
126  ymmtF = _mm256_slli_epi64(ymmsF, L7); \
127  ymms1 = _mm256_srli_epi64(ymms1, R0); \
128  ymms3 = _mm256_srli_epi64(ymms3, R1); \
129  ymms5 = _mm256_srli_epi64(ymms5, R2); \
130  ymms7 = _mm256_srli_epi64(ymms7, R3); \
131  ymms9 = _mm256_srli_epi64(ymms9, R4); \
132  ymmsB = _mm256_srli_epi64(ymmsB, R5); \
133  ymmsD = _mm256_srli_epi64(ymmsD, R6); \
134  ymmsF = _mm256_srli_epi64(ymmsF, R7); \
135  ymmt1 = _mm256_or_si256(ymms1, ymmt1); \
136  ymmt3 = _mm256_or_si256(ymms3, ymmt3); \
137  ymmt5 = _mm256_or_si256(ymms5, ymmt5); \
138  ymmt7 = _mm256_or_si256(ymms7, ymmt7); \
139  ymmt9 = _mm256_or_si256(ymms9, ymmt9); \
140  ymmtB = _mm256_or_si256(ymmsB, ymmtB); \
141  ymmtD = _mm256_or_si256(ymmsD, ymmtD); \
142  ymmtF = _mm256_or_si256(ymmsF, ymmtF); \
143  ymms0 = ymmt0; \
144  ymms2 = ymmt2; \
145  ymms4 = ymmt6; \
146  ymms6 = ymmt4; \
147  ymms8 = ymmtA; \
148  ymmsA = ymmtC; \
149  ymmsC = ymmtE; \
150  ymmsE = ymmt8; \
151  ymms1 = _mm256_xor_si256(ymmt8, ymmt9); \
152  ymms3 = _mm256_xor_si256(ymmtC, ymmtD); \
153  ymms5 = _mm256_xor_si256(ymmtA, ymmtB); \
154  ymms7 = _mm256_xor_si256(ymmtE, ymmtF); \
155  ymms9 = _mm256_xor_si256(ymmt6, ymmt7); \
156  ymmsB = _mm256_xor_si256(ymmt2, ymmt3); \
157  ymmsD = _mm256_xor_si256(ymmt4, ymmt5); \
158  ymmsF = _mm256_xor_si256(ymmt0, ymmt1); \
159  }
160 
161 #define MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_4(N) \
162  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_RBOX(N * 8 + 1); \
163  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_RBOX(N * 8 + 2); \
164  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_RBOX(N * 8 + 3); \
165  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_RBOX(N * 8 + 4); \
166  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_KBOX(N * 8 + 4);
167 
168 #define MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_8(N) \
169  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_RBOX(N * 8 + 1); \
170  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_RBOX(N * 8 + 2); \
171  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_RBOX(N * 8 + 3); \
172  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_RBOX(N * 8 + 4); \
173  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_KBOX(N * 8 + 4); \
174  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_RBOX(N * 8 + 5); \
175  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_RBOX(N * 8 + 6); \
176  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_RBOX(N * 8 + 7); \
177  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_RBOX(N * 8 + 8); \
178  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_KBOX(N * 8 + 8);
179 
180 namespace mckl {
181 
182 namespace internal {
183 
184 template <typename T>
185 class Threefry16x64GeneratorAVX2Impl
186 {
187  static_assert(std::numeric_limits<T>::digits == 64,
188  "**Threefry16x64GeneratorAVX2Impl** used with T other than a 64-bit "
189  "unsigned integers");
190 
191  static constexpr std::size_t K = 16;
192  static constexpr std::size_t Rounds = 20;
193  using Constants = ThreefryConstants<T, 16>;
194 
195  public:
196  static void eval(
197  const void *plain, void *cipher, const std::array<T, K + 4> &par)
198  {
199  Threefry16x64GeneratorGenericImpl<T>::eval(plain, cipher, par);
200  }
201 
202  template <typename ResultType>
203  static void eval(std::array<std::uint64_t, 16> &ctr, ResultType *r,
204  const std::array<T, K + 4> &par)
205  {
206  Threefry16x64GeneratorGenericImpl<T>::eval(ctr, r, par);
207  }
208 
209  template <typename ResultType>
210  static void eval(std::array<std::uint64_t, 16> &ctr, std::size_t n,
211  ResultType *r, const std::array<T, K + 4> &par)
212  {
213  constexpr std::size_t R = sizeof(T) * K / sizeof(ResultType);
214 
215  const std::size_t n0 =
216  static_cast<std::size_t>(std::min(static_cast<std::uint64_t>(n),
217  std::numeric_limits<std::uint64_t>::max() - ctr.front()));
218 
219  eval_kernel(ctr, n0, r, par);
220  n -= n0;
221  r += n0 * R;
222 
223  if (n != 0) {
224  eval(ctr, r, par);
225  n -= 1;
226  r += R;
227  }
228 
229  eval_kernel(ctr, n, r, par);
230  }
231 
232  private:
233  template <typename ResultType>
234  static void eval_kernel(std::array<std::uint64_t, 16> &ctr, std::size_t n,
235  ResultType *r, const std::array<T, K + 4> &par)
236  {
237  constexpr std::size_t S = 16;
238  constexpr std::size_t N = sizeof(__m256i) * S / (sizeof(T) * K);
239 
240  __m256i ymmc0 =
241  _mm256_set_epi64x(static_cast<MCKL_INT64>(std::get<0x3>(ctr)),
242  static_cast<MCKL_INT64>(std::get<0x2>(ctr)),
243  static_cast<MCKL_INT64>(std::get<0x1>(ctr)),
244  static_cast<MCKL_INT64>(std::get<0x0>(ctr)));
245  __m256i ymmc1 =
246  _mm256_set_epi64x(static_cast<MCKL_INT64>(std::get<0x7>(ctr)),
247  static_cast<MCKL_INT64>(std::get<0x6>(ctr)),
248  static_cast<MCKL_INT64>(std::get<0x5>(ctr)),
249  static_cast<MCKL_INT64>(std::get<0x4>(ctr)));
250  __m256i ymmc2 =
251  _mm256_set_epi64x(static_cast<MCKL_INT64>(std::get<0xB>(ctr)),
252  static_cast<MCKL_INT64>(std::get<0xA>(ctr)),
253  static_cast<MCKL_INT64>(std::get<0x9>(ctr)),
254  static_cast<MCKL_INT64>(std::get<0x8>(ctr)));
255  __m256i ymmc3 =
256  _mm256_set_epi64x(static_cast<MCKL_INT64>(std::get<0xF>(ctr)),
257  static_cast<MCKL_INT64>(std::get<0xE>(ctr)),
258  static_cast<MCKL_INT64>(std::get<0xD>(ctr)),
259  static_cast<MCKL_INT64>(std::get<0xC>(ctr)));
260  ctr.front() += n;
261 
262  __m256i *rptr = reinterpret_cast<__m256i *>(r);
263  while (n != 0) {
264  __m256i ymms0 =
265  _mm256_add_epi64(ymmc0, _mm256_set_epi64x(0, 0, 0, 1));
266  __m256i ymms4 =
267  _mm256_add_epi64(ymmc0, _mm256_set_epi64x(0, 0, 0, 2));
268  __m256i ymms8 =
269  _mm256_add_epi64(ymmc0, _mm256_set_epi64x(0, 0, 0, 3));
270  __m256i ymmsC =
271  _mm256_add_epi64(ymmc0, _mm256_set_epi64x(0, 0, 0, 4));
272  ymmc0 = ymmsC;
273 
274  __m256i ymms1 = ymmc1;
275  __m256i ymms2 = ymmc2;
276  __m256i ymms3 = ymmc3;
277  __m256i ymms5 = ymmc1;
278  __m256i ymms6 = ymmc2;
279  __m256i ymms7 = ymmc3;
280  __m256i ymms9 = ymmc1;
281  __m256i ymmsA = ymmc2;
282  __m256i ymmsB = ymmc3;
283  __m256i ymmsD = ymmc1;
284  __m256i ymmsE = ymmc2;
285  __m256i ymmsF = ymmc3;
286 
287  __m256i ymmt0;
288  __m256i ymmt1;
289  __m256i ymmt2;
290  __m256i ymmt3;
291  __m256i ymmt4;
292  __m256i ymmt5;
293  __m256i ymmt6;
294  __m256i ymmt7;
295  __m256i ymmt8;
296  __m256i ymmt9;
297  __m256i ymmtA;
298  __m256i ymmtB;
299  __m256i ymmtC;
300  __m256i ymmtD;
301  __m256i ymmtE;
302  __m256i ymmtF;
303 
304  ymmt1 = ymms1;
305  ymmt2 = ymms2;
306  ymmt3 = ymms3;
307  ymmt6 = ymms6;
308  ymmt7 = ymms7;
309  ymmtB = ymmsB;
310  ymms1 = ymms4;
311  ymms2 = ymms8;
312  ymms3 = ymmsC;
313  ymms6 = ymms9;
314  ymms7 = ymmsD;
315  ymmsB = ymmsE;
316  ymms4 = ymmt1;
317  ymms8 = ymmt2;
318  ymmsC = ymmt3;
319  ymms9 = ymmt6;
320  ymmsD = ymmt7;
321  ymmsE = ymmtB;
322 
323  ymmt0 = _mm256_unpacklo_epi64(ymms0, ymms1);
324  ymmt1 = _mm256_unpacklo_epi64(ymms2, ymms3);
325  ymmt2 = _mm256_unpackhi_epi64(ymms0, ymms1);
326  ymmt3 = _mm256_unpackhi_epi64(ymms2, ymms3);
327  ymmt4 = _mm256_unpacklo_epi64(ymms4, ymms5);
328  ymmt5 = _mm256_unpacklo_epi64(ymms6, ymms7);
329  ymmt6 = _mm256_unpackhi_epi64(ymms4, ymms5);
330  ymmt7 = _mm256_unpackhi_epi64(ymms6, ymms7);
331  ymmt8 = _mm256_unpacklo_epi64(ymms8, ymms9);
332  ymmt9 = _mm256_unpacklo_epi64(ymmsA, ymmsB);
333  ymmtA = _mm256_unpackhi_epi64(ymms8, ymms9);
334  ymmtB = _mm256_unpackhi_epi64(ymmsA, ymmsB);
335  ymmtC = _mm256_unpacklo_epi64(ymmsC, ymmsD);
336  ymmtD = _mm256_unpacklo_epi64(ymmsE, ymmsF);
337  ymmtE = _mm256_unpackhi_epi64(ymmsC, ymmsD);
338  ymmtF = _mm256_unpackhi_epi64(ymmsE, ymmsF);
339  ymms0 = _mm256_permute2x128_si256(ymmt0, ymmt1, 0x20);
340  ymms1 = _mm256_permute2x128_si256(ymmt2, ymmt3, 0x20);
341  ymms2 = _mm256_permute2x128_si256(ymmt0, ymmt1, 0x31);
342  ymms3 = _mm256_permute2x128_si256(ymmt2, ymmt3, 0x31);
343  ymms4 = _mm256_permute2x128_si256(ymmt4, ymmt5, 0x20);
344  ymms5 = _mm256_permute2x128_si256(ymmt6, ymmt7, 0x20);
345  ymms6 = _mm256_permute2x128_si256(ymmt4, ymmt5, 0x31);
346  ymms7 = _mm256_permute2x128_si256(ymmt6, ymmt7, 0x31);
347  ymms8 = _mm256_permute2x128_si256(ymmt8, ymmt9, 0x20);
348  ymms9 = _mm256_permute2x128_si256(ymmtA, ymmtB, 0x20);
349  ymmsA = _mm256_permute2x128_si256(ymmt8, ymmt9, 0x31);
350  ymmsB = _mm256_permute2x128_si256(ymmtA, ymmtB, 0x31);
351  ymmsC = _mm256_permute2x128_si256(ymmtC, ymmtD, 0x20);
352  ymmsD = _mm256_permute2x128_si256(ymmtE, ymmtF, 0x20);
353  ymmsE = _mm256_permute2x128_si256(ymmtC, ymmtD, 0x31);
354  ymmsF = _mm256_permute2x128_si256(ymmtE, ymmtF, 0x31);
355 
356  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_KBOX(0);
357  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_8(0);
358  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_8(1);
359  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_4(2);
360 
361  ymmt0 = _mm256_unpacklo_epi64(ymms0, ymms1);
362  ymmt1 = _mm256_unpacklo_epi64(ymms2, ymms3);
363  ymmt2 = _mm256_unpackhi_epi64(ymms0, ymms1);
364  ymmt3 = _mm256_unpackhi_epi64(ymms2, ymms3);
365  ymmt4 = _mm256_unpacklo_epi64(ymms4, ymms5);
366  ymmt5 = _mm256_unpacklo_epi64(ymms6, ymms7);
367  ymmt6 = _mm256_unpackhi_epi64(ymms4, ymms5);
368  ymmt7 = _mm256_unpackhi_epi64(ymms6, ymms7);
369  ymmt8 = _mm256_unpacklo_epi64(ymms8, ymms9);
370  ymmt9 = _mm256_unpacklo_epi64(ymmsA, ymmsB);
371  ymmtA = _mm256_unpackhi_epi64(ymms8, ymms9);
372  ymmtB = _mm256_unpackhi_epi64(ymmsA, ymmsB);
373  ymmtC = _mm256_unpacklo_epi64(ymmsC, ymmsD);
374  ymmtD = _mm256_unpacklo_epi64(ymmsE, ymmsF);
375  ymmtE = _mm256_unpackhi_epi64(ymmsC, ymmsD);
376  ymmtF = _mm256_unpackhi_epi64(ymmsE, ymmsF);
377  ymms0 = _mm256_permute2x128_si256(ymmt0, ymmt1, 0x20);
378  ymms1 = _mm256_permute2x128_si256(ymmt2, ymmt3, 0x20);
379  ymms2 = _mm256_permute2x128_si256(ymmt0, ymmt1, 0x31);
380  ymms3 = _mm256_permute2x128_si256(ymmt2, ymmt3, 0x31);
381  ymms4 = _mm256_permute2x128_si256(ymmt4, ymmt5, 0x20);
382  ymms5 = _mm256_permute2x128_si256(ymmt6, ymmt7, 0x20);
383  ymms6 = _mm256_permute2x128_si256(ymmt4, ymmt5, 0x31);
384  ymms7 = _mm256_permute2x128_si256(ymmt6, ymmt7, 0x31);
385  ymms8 = _mm256_permute2x128_si256(ymmt8, ymmt9, 0x20);
386  ymms9 = _mm256_permute2x128_si256(ymmtA, ymmtB, 0x20);
387  ymmsA = _mm256_permute2x128_si256(ymmt8, ymmt9, 0x31);
388  ymmsB = _mm256_permute2x128_si256(ymmtA, ymmtB, 0x31);
389  ymmsC = _mm256_permute2x128_si256(ymmtC, ymmtD, 0x20);
390  ymmsD = _mm256_permute2x128_si256(ymmtE, ymmtF, 0x20);
391  ymmsE = _mm256_permute2x128_si256(ymmtC, ymmtD, 0x31);
392  ymmsF = _mm256_permute2x128_si256(ymmtE, ymmtF, 0x31);
393 
394  ymmt1 = ymms1;
395  ymmt2 = ymms2;
396  ymmt3 = ymms3;
397  ymmt6 = ymms6;
398  ymmt7 = ymms7;
399  ymmtB = ymmsB;
400  ymms1 = ymms4;
401  ymms2 = ymms8;
402  ymms3 = ymmsC;
403  ymms6 = ymms9;
404  ymms7 = ymmsD;
405  ymmsB = ymmsE;
406  ymms4 = ymmt1;
407  ymms8 = ymmt2;
408  ymmsC = ymmt3;
409  ymms9 = ymmt6;
410  ymmsD = ymmt7;
411  ymmsE = ymmtB;
412 
413  if (n >= N) {
414  n -= N;
415  _mm256_storeu_si256(rptr++, ymms0);
416  _mm256_storeu_si256(rptr++, ymms1);
417  _mm256_storeu_si256(rptr++, ymms2);
418  _mm256_storeu_si256(rptr++, ymms3);
419  _mm256_storeu_si256(rptr++, ymms4);
420  _mm256_storeu_si256(rptr++, ymms5);
421  _mm256_storeu_si256(rptr++, ymms6);
422  _mm256_storeu_si256(rptr++, ymms7);
423  _mm256_storeu_si256(rptr++, ymms8);
424  _mm256_storeu_si256(rptr++, ymms9);
425  _mm256_storeu_si256(rptr++, ymmsA);
426  _mm256_storeu_si256(rptr++, ymmsB);
427  _mm256_storeu_si256(rptr++, ymmsC);
428  _mm256_storeu_si256(rptr++, ymmsD);
429  _mm256_storeu_si256(rptr++, ymmsE);
430  _mm256_storeu_si256(rptr++, ymmsF);
431  } else {
432  std::array<__m256i, S> s;
433  std::get<0x0>(s) = ymms0;
434  std::get<0x1>(s) = ymms1;
435  std::get<0x2>(s) = ymms2;
436  std::get<0x3>(s) = ymms3;
437  std::get<0x4>(s) = ymms4;
438  std::get<0x5>(s) = ymms5;
439  std::get<0x6>(s) = ymms6;
440  std::get<0x7>(s) = ymms7;
441  std::get<0x8>(s) = ymms8;
442  std::get<0x9>(s) = ymms9;
443  std::get<0xA>(s) = ymmsA;
444  std::get<0xB>(s) = ymmsB;
445  std::get<0xC>(s) = ymmsC;
446  std::get<0xD>(s) = ymmsD;
447  std::get<0xE>(s) = ymmsE;
448  std::get<0xF>(s) = ymmsF;
449  std::memcpy(rptr, s.data(), n * sizeof(T) * K);
450  break;
451  }
452  }
453  }
454 }; // class Threefry16x64GeneratorAVX2Impl
455 
456 template <typename T>
457 class Threefish1024GeneratorAVX2Impl
458 {
459  static_assert(std::numeric_limits<T>::digits == 64,
460  "**Threefish1024GeneratorAVX2Impl** used with T other than a 64-bit "
461  "unsigned integers");
462 
463  static constexpr std::size_t K = 16;
464  static constexpr std::size_t Rounds = 80;
465  using Constants = ThreefryConstants<T, 16>;
466 
467  public:
468  static void eval(
469  const void *plain, void *cipher, const std::array<T, K + 4> &par)
470  {
471  Threefish1024GeneratorGenericImpl<T>::eval(plain, cipher, par);
472  }
473 
474  template <typename ResultType>
475  static void eval(std::array<std::uint64_t, 16> &ctr, ResultType *r,
476  const std::array<T, K + 4> &par)
477  {
478  Threefish1024GeneratorGenericImpl<T>::eval(ctr, r, par);
479  }
480 
481  template <typename ResultType>
482  static void eval(std::array<std::uint64_t, 16> &ctr, std::size_t n,
483  ResultType *r, const std::array<T, K + 4> &par)
484  {
485  constexpr std::size_t R = sizeof(T) * K / sizeof(ResultType);
486 
487  const std::size_t n0 =
488  static_cast<std::size_t>(std::min(static_cast<std::uint64_t>(n),
489  std::numeric_limits<std::uint64_t>::max() - ctr.front()));
490 
491  eval_kernel(ctr, n0, r, par);
492  n -= n0;
493  r += n0 * R;
494 
495  if (n != 0) {
496  eval(ctr, r, par);
497  n -= 1;
498  r += R;
499  }
500 
501  eval_kernel(ctr, n, r, par);
502  }
503 
504  private:
505  template <typename ResultType>
506  static void eval_kernel(std::array<std::uint64_t, 16> &ctr, std::size_t n,
507  ResultType *r, const std::array<T, K + 4> &par)
508  {
509  constexpr std::size_t S = 16;
510  constexpr std::size_t N = sizeof(__m256i) * S / (sizeof(T) * K);
511 
512  __m256i ymmc0 =
513  _mm256_set_epi64x(static_cast<MCKL_INT64>(std::get<0x3>(ctr)),
514  static_cast<MCKL_INT64>(std::get<0x2>(ctr)),
515  static_cast<MCKL_INT64>(std::get<0x1>(ctr)),
516  static_cast<MCKL_INT64>(std::get<0x0>(ctr)));
517  __m256i ymmc1 =
518  _mm256_set_epi64x(static_cast<MCKL_INT64>(std::get<0x7>(ctr)),
519  static_cast<MCKL_INT64>(std::get<0x6>(ctr)),
520  static_cast<MCKL_INT64>(std::get<0x5>(ctr)),
521  static_cast<MCKL_INT64>(std::get<0x4>(ctr)));
522  __m256i ymmc2 =
523  _mm256_set_epi64x(static_cast<MCKL_INT64>(std::get<0xB>(ctr)),
524  static_cast<MCKL_INT64>(std::get<0xA>(ctr)),
525  static_cast<MCKL_INT64>(std::get<0x9>(ctr)),
526  static_cast<MCKL_INT64>(std::get<0x8>(ctr)));
527  __m256i ymmc3 =
528  _mm256_set_epi64x(static_cast<MCKL_INT64>(std::get<0xF>(ctr)),
529  static_cast<MCKL_INT64>(std::get<0xE>(ctr)),
530  static_cast<MCKL_INT64>(std::get<0xD>(ctr)),
531  static_cast<MCKL_INT64>(std::get<0xC>(ctr)));
532  ctr.front() += n;
533 
534  __m256i *rptr = reinterpret_cast<__m256i *>(r);
535  while (n != 0) {
536  __m256i ymms0 =
537  _mm256_add_epi64(ymmc0, _mm256_set_epi64x(0, 0, 0, 1));
538  __m256i ymms4 =
539  _mm256_add_epi64(ymmc0, _mm256_set_epi64x(0, 0, 0, 2));
540  __m256i ymms8 =
541  _mm256_add_epi64(ymmc0, _mm256_set_epi64x(0, 0, 0, 3));
542  __m256i ymmsC =
543  _mm256_add_epi64(ymmc0, _mm256_set_epi64x(0, 0, 0, 4));
544  ymmc0 = ymmsC;
545 
546  __m256i ymms1 = ymmc1;
547  __m256i ymms2 = ymmc2;
548  __m256i ymms3 = ymmc3;
549  __m256i ymms5 = ymmc1;
550  __m256i ymms6 = ymmc2;
551  __m256i ymms7 = ymmc3;
552  __m256i ymms9 = ymmc1;
553  __m256i ymmsA = ymmc2;
554  __m256i ymmsB = ymmc3;
555  __m256i ymmsD = ymmc1;
556  __m256i ymmsE = ymmc2;
557  __m256i ymmsF = ymmc3;
558 
559  __m256i ymmt0;
560  __m256i ymmt1;
561  __m256i ymmt2;
562  __m256i ymmt3;
563  __m256i ymmt4;
564  __m256i ymmt5;
565  __m256i ymmt6;
566  __m256i ymmt7;
567  __m256i ymmt8;
568  __m256i ymmt9;
569  __m256i ymmtA;
570  __m256i ymmtB;
571  __m256i ymmtC;
572  __m256i ymmtD;
573  __m256i ymmtE;
574  __m256i ymmtF;
575 
576  ymmt1 = ymms1;
577  ymmt2 = ymms2;
578  ymmt3 = ymms3;
579  ymmt6 = ymms6;
580  ymmt7 = ymms7;
581  ymmtB = ymmsB;
582  ymms1 = ymms4;
583  ymms2 = ymms8;
584  ymms3 = ymmsC;
585  ymms6 = ymms9;
586  ymms7 = ymmsD;
587  ymmsB = ymmsE;
588  ymms4 = ymmt1;
589  ymms8 = ymmt2;
590  ymmsC = ymmt3;
591  ymms9 = ymmt6;
592  ymmsD = ymmt7;
593  ymmsE = ymmtB;
594 
595  ymmt0 = _mm256_unpacklo_epi64(ymms0, ymms1);
596  ymmt1 = _mm256_unpacklo_epi64(ymms2, ymms3);
597  ymmt2 = _mm256_unpackhi_epi64(ymms0, ymms1);
598  ymmt3 = _mm256_unpackhi_epi64(ymms2, ymms3);
599  ymmt4 = _mm256_unpacklo_epi64(ymms4, ymms5);
600  ymmt5 = _mm256_unpacklo_epi64(ymms6, ymms7);
601  ymmt6 = _mm256_unpackhi_epi64(ymms4, ymms5);
602  ymmt7 = _mm256_unpackhi_epi64(ymms6, ymms7);
603  ymmt8 = _mm256_unpacklo_epi64(ymms8, ymms9);
604  ymmt9 = _mm256_unpacklo_epi64(ymmsA, ymmsB);
605  ymmtA = _mm256_unpackhi_epi64(ymms8, ymms9);
606  ymmtB = _mm256_unpackhi_epi64(ymmsA, ymmsB);
607  ymmtC = _mm256_unpacklo_epi64(ymmsC, ymmsD);
608  ymmtD = _mm256_unpacklo_epi64(ymmsE, ymmsF);
609  ymmtE = _mm256_unpackhi_epi64(ymmsC, ymmsD);
610  ymmtF = _mm256_unpackhi_epi64(ymmsE, ymmsF);
611  ymms0 = _mm256_permute2x128_si256(ymmt0, ymmt1, 0x20);
612  ymms1 = _mm256_permute2x128_si256(ymmt2, ymmt3, 0x20);
613  ymms2 = _mm256_permute2x128_si256(ymmt0, ymmt1, 0x31);
614  ymms3 = _mm256_permute2x128_si256(ymmt2, ymmt3, 0x31);
615  ymms4 = _mm256_permute2x128_si256(ymmt4, ymmt5, 0x20);
616  ymms5 = _mm256_permute2x128_si256(ymmt6, ymmt7, 0x20);
617  ymms6 = _mm256_permute2x128_si256(ymmt4, ymmt5, 0x31);
618  ymms7 = _mm256_permute2x128_si256(ymmt6, ymmt7, 0x31);
619  ymms8 = _mm256_permute2x128_si256(ymmt8, ymmt9, 0x20);
620  ymms9 = _mm256_permute2x128_si256(ymmtA, ymmtB, 0x20);
621  ymmsA = _mm256_permute2x128_si256(ymmt8, ymmt9, 0x31);
622  ymmsB = _mm256_permute2x128_si256(ymmtA, ymmtB, 0x31);
623  ymmsC = _mm256_permute2x128_si256(ymmtC, ymmtD, 0x20);
624  ymmsD = _mm256_permute2x128_si256(ymmtE, ymmtF, 0x20);
625  ymmsE = _mm256_permute2x128_si256(ymmtC, ymmtD, 0x31);
626  ymmsF = _mm256_permute2x128_si256(ymmtE, ymmtF, 0x31);
627 
628  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_KBOX(0);
629  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_8(0);
630  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_8(1);
631  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_8(2);
632  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_8(3);
633  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_8(4);
634  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_8(5);
635  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_8(6);
636  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_8(7);
637  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_8(8);
638  MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_8(9);
639 
640  ymmt0 = _mm256_unpacklo_epi64(ymms0, ymms1);
641  ymmt1 = _mm256_unpacklo_epi64(ymms2, ymms3);
642  ymmt2 = _mm256_unpackhi_epi64(ymms0, ymms1);
643  ymmt3 = _mm256_unpackhi_epi64(ymms2, ymms3);
644  ymmt4 = _mm256_unpacklo_epi64(ymms4, ymms5);
645  ymmt5 = _mm256_unpacklo_epi64(ymms6, ymms7);
646  ymmt6 = _mm256_unpackhi_epi64(ymms4, ymms5);
647  ymmt7 = _mm256_unpackhi_epi64(ymms6, ymms7);
648  ymmt8 = _mm256_unpacklo_epi64(ymms8, ymms9);
649  ymmt9 = _mm256_unpacklo_epi64(ymmsA, ymmsB);
650  ymmtA = _mm256_unpackhi_epi64(ymms8, ymms9);
651  ymmtB = _mm256_unpackhi_epi64(ymmsA, ymmsB);
652  ymmtC = _mm256_unpacklo_epi64(ymmsC, ymmsD);
653  ymmtD = _mm256_unpacklo_epi64(ymmsE, ymmsF);
654  ymmtE = _mm256_unpackhi_epi64(ymmsC, ymmsD);
655  ymmtF = _mm256_unpackhi_epi64(ymmsE, ymmsF);
656  ymms0 = _mm256_permute2x128_si256(ymmt0, ymmt1, 0x20);
657  ymms1 = _mm256_permute2x128_si256(ymmt2, ymmt3, 0x20);
658  ymms2 = _mm256_permute2x128_si256(ymmt0, ymmt1, 0x31);
659  ymms3 = _mm256_permute2x128_si256(ymmt2, ymmt3, 0x31);
660  ymms4 = _mm256_permute2x128_si256(ymmt4, ymmt5, 0x20);
661  ymms5 = _mm256_permute2x128_si256(ymmt6, ymmt7, 0x20);
662  ymms6 = _mm256_permute2x128_si256(ymmt4, ymmt5, 0x31);
663  ymms7 = _mm256_permute2x128_si256(ymmt6, ymmt7, 0x31);
664  ymms8 = _mm256_permute2x128_si256(ymmt8, ymmt9, 0x20);
665  ymms9 = _mm256_permute2x128_si256(ymmtA, ymmtB, 0x20);
666  ymmsA = _mm256_permute2x128_si256(ymmt8, ymmt9, 0x31);
667  ymmsB = _mm256_permute2x128_si256(ymmtA, ymmtB, 0x31);
668  ymmsC = _mm256_permute2x128_si256(ymmtC, ymmtD, 0x20);
669  ymmsD = _mm256_permute2x128_si256(ymmtE, ymmtF, 0x20);
670  ymmsE = _mm256_permute2x128_si256(ymmtC, ymmtD, 0x31);
671  ymmsF = _mm256_permute2x128_si256(ymmtE, ymmtF, 0x31);
672 
673  ymmt1 = ymms1;
674  ymmt2 = ymms2;
675  ymmt3 = ymms3;
676  ymmt6 = ymms6;
677  ymmt7 = ymms7;
678  ymmtB = ymmsB;
679  ymms1 = ymms4;
680  ymms2 = ymms8;
681  ymms3 = ymmsC;
682  ymms6 = ymms9;
683  ymms7 = ymmsD;
684  ymmsB = ymmsE;
685  ymms4 = ymmt1;
686  ymms8 = ymmt2;
687  ymmsC = ymmt3;
688  ymms9 = ymmt6;
689  ymmsD = ymmt7;
690  ymmsE = ymmtB;
691 
692  if (n >= N) {
693  n -= N;
694  _mm256_storeu_si256(rptr++, ymms0);
695  _mm256_storeu_si256(rptr++, ymms1);
696  _mm256_storeu_si256(rptr++, ymms2);
697  _mm256_storeu_si256(rptr++, ymms3);
698  _mm256_storeu_si256(rptr++, ymms4);
699  _mm256_storeu_si256(rptr++, ymms5);
700  _mm256_storeu_si256(rptr++, ymms6);
701  _mm256_storeu_si256(rptr++, ymms7);
702  _mm256_storeu_si256(rptr++, ymms8);
703  _mm256_storeu_si256(rptr++, ymms9);
704  _mm256_storeu_si256(rptr++, ymmsA);
705  _mm256_storeu_si256(rptr++, ymmsB);
706  _mm256_storeu_si256(rptr++, ymmsC);
707  _mm256_storeu_si256(rptr++, ymmsD);
708  _mm256_storeu_si256(rptr++, ymmsE);
709  _mm256_storeu_si256(rptr++, ymmsF);
710  } else {
711  std::array<__m256i, S> s;
712  std::get<0x0>(s) = ymms0;
713  std::get<0x1>(s) = ymms1;
714  std::get<0x2>(s) = ymms2;
715  std::get<0x3>(s) = ymms3;
716  std::get<0x4>(s) = ymms4;
717  std::get<0x5>(s) = ymms5;
718  std::get<0x6>(s) = ymms6;
719  std::get<0x7>(s) = ymms7;
720  std::get<0x8>(s) = ymms8;
721  std::get<0x9>(s) = ymms9;
722  std::get<0xA>(s) = ymmsA;
723  std::get<0xB>(s) = ymmsB;
724  std::get<0xC>(s) = ymmsC;
725  std::get<0xD>(s) = ymmsD;
726  std::get<0xE>(s) = ymmsE;
727  std::get<0xF>(s) = ymmsF;
728  std::memcpy(rptr, s.data(), n * sizeof(T) * K);
729  break;
730  }
731  }
732  }
733 }; // class Threefish1024GeneratorAVX2Impl
734 
735 } // namespace internal
736 
737 } // namespace mckl
738 
739 MCKL_POP_GCC_WARNING
740 
741 #endif // MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_HPP
MCKL_PUSH_GCC_WARNING
#define MCKL_PUSH_GCC_WARNING(warning)
Definition: compiler.h:78
MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_4
#define MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_4(N)
Definition: threefry_avx2_16x64.hpp:161
mckl::internal::Threefish1024GeneratorAVX2Impl::eval
static void eval(std::array< std::uint64_t, 16 > &ctr, std::size_t n, ResultType *r, const std::array< T, K+4 > &par)
Definition: threefry_avx2_16x64.hpp:482
MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_8
#define MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_CYCLE_8(N)
Definition: threefry_avx2_16x64.hpp:168
threefry_generic_16x64.hpp
mckl::internal::Threefish1024GeneratorGenericImpl::eval
static void eval(const void *plain, void *cipher, const std::array< T, K+4 > &par)
Definition: threefry_generic_16x64.hpp:276
mckl::internal::Threefish1024GeneratorAVX2Impl::eval
static void eval(std::array< std::uint64_t, 16 > &ctr, ResultType *r, const std::array< T, K+4 > &par)
Definition: threefry_avx2_16x64.hpp:475
mckl::internal::Threefry16x64GeneratorAVX2Impl::eval
static void eval(std::array< std::uint64_t, 16 > &ctr, std::size_t n, ResultType *r, const std::array< T, K+4 > &par)
Definition: threefry_avx2_16x64.hpp:210
mckl::internal::Threefry16x64GeneratorAVX2Impl::eval
static void eval(const void *plain, void *cipher, const std::array< T, K+4 > &par)
Definition: threefry_avx2_16x64.hpp:196
mckl::ThreefryConstants
Default Threefry constants.
Definition: threefry_constants.hpp:147
MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_KBOX
#define MCKL_RANDOM_INTERNAL_THREEFRY_AVX2_16X64_KBOX(N)
Definition: threefry_avx2_16x64.hpp:43
mckl::internal::Threefry16x64GeneratorGenericImpl::eval
static void eval(const void *plain, void *cipher, const std::array< T, K+4 > &par)
Definition: threefry_generic_16x64.hpp:215
threefry_constants.hpp
mckl
Definition: mcmc.hpp:40
threefry_common.hpp
mckl::internal::Threefry16x64GeneratorAVX2Impl::eval
static void eval(std::array< std::uint64_t, 16 > &ctr, ResultType *r, const std::array< T, K+4 > &par)
Definition: threefry_avx2_16x64.hpp:203
mckl::internal::Threefish1024GeneratorAVX2Impl::eval
static void eval(const void *plain, void *cipher, const std::array< T, K+4 > &par)
Definition: threefry_avx2_16x64.hpp:468
MCKL_POP_GCC_WARNING
#define MCKL_POP_GCC_WARNING
Definition: compiler.h:79
Generated by   doxygen 1.8.13