MCKL
Monte Carlo Kernel Library
threefry_sse2_16x64.hpp
Go to the documentation of this file.
1 //============================================================================
2 // MCKL/include/mckl/random/internal/threefry_sse2_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_SSE2_16X64_HPP
33 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_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_SSE2_16X64_KBOX(N) \
44  xmmt0 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
45  ThreefryKBox<T, K, N>::template key<0x0>(par))); \
46  xmmt1 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
47  ThreefryKBox<T, K, N>::template key<0x1>(par))); \
48  xmmt2 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
49  ThreefryKBox<T, K, N>::template key<0x2>(par))); \
50  xmmt3 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
51  ThreefryKBox<T, K, N>::template key<0x3>(par))); \
52  xmmt4 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
53  ThreefryKBox<T, K, N>::template key<0x4>(par))); \
54  xmmt5 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
55  ThreefryKBox<T, K, N>::template key<0x5>(par))); \
56  xmmt6 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
57  ThreefryKBox<T, K, N>::template key<0x6>(par))); \
58  xmmt7 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
59  ThreefryKBox<T, K, N>::template key<0x7>(par))); \
60  xmmt8 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
61  ThreefryKBox<T, K, N>::template key<0x8>(par))); \
62  xmmt9 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
63  ThreefryKBox<T, K, N>::template key<0x9>(par))); \
64  xmmtA = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
65  ThreefryKBox<T, K, N>::template key<0xA>(par))); \
66  xmmtB = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
67  ThreefryKBox<T, K, N>::template key<0xB>(par))); \
68  xmmtC = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
69  ThreefryKBox<T, K, N>::template key<0xC>(par))); \
70  xmmtD = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
71  ThreefryKBox<T, K, N>::template key<0xD>(par))); \
72  xmmtE = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
73  ThreefryKBox<T, K, N>::template key<0xE>(par))); \
74  xmmtF = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
75  ThreefryKBox<T, K, N>::template key<0xF>(par))); \
76  xmms0 = _mm_add_epi64(xmms0, xmmt0); \
77  xmms1 = _mm_add_epi64(xmms1, xmmt1); \
78  xmms2 = _mm_add_epi64(xmms2, xmmt2); \
79  xmms3 = _mm_add_epi64(xmms3, xmmt3); \
80  xmms4 = _mm_add_epi64(xmms4, xmmt4); \
81  xmms5 = _mm_add_epi64(xmms5, xmmt5); \
82  xmms6 = _mm_add_epi64(xmms6, xmmt6); \
83  xmms7 = _mm_add_epi64(xmms7, xmmt7); \
84  xmms8 = _mm_add_epi64(xmms8, xmmt8); \
85  xmms9 = _mm_add_epi64(xmms9, xmmt9); \
86  xmmsA = _mm_add_epi64(xmmsA, xmmtA); \
87  xmmsB = _mm_add_epi64(xmmsB, xmmtB); \
88  xmmsC = _mm_add_epi64(xmmsC, xmmtC); \
89  xmmsD = _mm_add_epi64(xmmsD, xmmtD); \
90  xmmsE = _mm_add_epi64(xmmsE, xmmtE); \
91  xmmsF = _mm_add_epi64(xmmsF, xmmtF);
92 
93 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_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  xmmt0 = _mm_add_epi64(xmms0, xmms1); \
112  xmmt2 = _mm_add_epi64(xmms2, xmms3); \
113  xmmt4 = _mm_add_epi64(xmms4, xmms5); \
114  xmmt6 = _mm_add_epi64(xmms6, xmms7); \
115  xmmt8 = _mm_add_epi64(xmms8, xmms9); \
116  xmmtA = _mm_add_epi64(xmmsA, xmmsB); \
117  xmmtC = _mm_add_epi64(xmmsC, xmmsD); \
118  xmmtE = _mm_add_epi64(xmmsE, xmmsF); \
119  xmmt1 = _mm_slli_epi64(xmms1, L0); \
120  xmmt3 = _mm_slli_epi64(xmms3, L1); \
121  xmmt5 = _mm_slli_epi64(xmms5, L2); \
122  xmmt7 = _mm_slli_epi64(xmms7, L3); \
123  xmmt9 = _mm_slli_epi64(xmms9, L4); \
124  xmmtB = _mm_slli_epi64(xmmsB, L5); \
125  xmmtD = _mm_slli_epi64(xmmsD, L6); \
126  xmmtF = _mm_slli_epi64(xmmsF, L7); \
127  xmms1 = _mm_srli_epi64(xmms1, R0); \
128  xmms3 = _mm_srli_epi64(xmms3, R1); \
129  xmms5 = _mm_srli_epi64(xmms5, R2); \
130  xmms7 = _mm_srli_epi64(xmms7, R3); \
131  xmms9 = _mm_srli_epi64(xmms9, R4); \
132  xmmsB = _mm_srli_epi64(xmmsB, R5); \
133  xmmsD = _mm_srli_epi64(xmmsD, R6); \
134  xmmsF = _mm_srli_epi64(xmmsF, R7); \
135  xmmt1 = _mm_or_si128(xmms1, xmmt1); \
136  xmmt3 = _mm_or_si128(xmms3, xmmt3); \
137  xmmt5 = _mm_or_si128(xmms5, xmmt5); \
138  xmmt7 = _mm_or_si128(xmms7, xmmt7); \
139  xmmt9 = _mm_or_si128(xmms9, xmmt9); \
140  xmmtB = _mm_or_si128(xmmsB, xmmtB); \
141  xmmtD = _mm_or_si128(xmmsD, xmmtD); \
142  xmmtF = _mm_or_si128(xmmsF, xmmtF); \
143  xmms0 = xmmt0; \
144  xmms2 = xmmt2; \
145  xmms4 = xmmt6; \
146  xmms6 = xmmt4; \
147  xmms8 = xmmtA; \
148  xmmsA = xmmtC; \
149  xmmsC = xmmtE; \
150  xmmsE = xmmt8; \
151  xmms1 = _mm_xor_si128(xmmt8, xmmt9); \
152  xmms3 = _mm_xor_si128(xmmtC, xmmtD); \
153  xmms5 = _mm_xor_si128(xmmtA, xmmtB); \
154  xmms7 = _mm_xor_si128(xmmtE, xmmtF); \
155  xmms9 = _mm_xor_si128(xmmt6, xmmt7); \
156  xmmsB = _mm_xor_si128(xmmt2, xmmt3); \
157  xmmsD = _mm_xor_si128(xmmt4, xmmt5); \
158  xmmsF = _mm_xor_si128(xmmt0, xmmt1); \
159  }
160 
161 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_4(N) \
162  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_RBOX(N * 8 + 1); \
163  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_RBOX(N * 8 + 2); \
164  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_RBOX(N * 8 + 3); \
165  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_RBOX(N * 8 + 4); \
166  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_KBOX(N * 8 + 4);
167 
168 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_8(N) \
169  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_RBOX(N * 8 + 1); \
170  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_RBOX(N * 8 + 2); \
171  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_RBOX(N * 8 + 3); \
172  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_RBOX(N * 8 + 4); \
173  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_KBOX(N * 8 + 4); \
174  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_RBOX(N * 8 + 5); \
175  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_RBOX(N * 8 + 6); \
176  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_RBOX(N * 8 + 7); \
177  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_RBOX(N * 8 + 8); \
178  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_KBOX(N * 8 + 8);
179 
180 namespace mckl {
181 
182 namespace internal {
183 
184 template <typename T>
185 class Threefry16x64GeneratorSSE2Impl
186 {
187  static_assert(std::numeric_limits<T>::digits == 64,
188  "**Threefry16x64GeneratorSSE2Impl** 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(__m128i) * S / (sizeof(T) * K);
239 
240  __m128i xmmc0 =
241  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<0x1>(ctr)),
242  static_cast<MCKL_INT64>(std::get<0x0>(ctr)));
243  __m128i xmmc1 =
244  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<0x3>(ctr)),
245  static_cast<MCKL_INT64>(std::get<0x2>(ctr)));
246  __m128i xmmc2 =
247  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<0x5>(ctr)),
248  static_cast<MCKL_INT64>(std::get<0x4>(ctr)));
249  __m128i xmmc3 =
250  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<0x7>(ctr)),
251  static_cast<MCKL_INT64>(std::get<0x6>(ctr)));
252  __m128i xmmc4 =
253  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<0x9>(ctr)),
254  static_cast<MCKL_INT64>(std::get<0x8>(ctr)));
255  __m128i xmmc5 =
256  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<0xB>(ctr)),
257  static_cast<MCKL_INT64>(std::get<0xA>(ctr)));
258  __m128i xmmc6 =
259  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<0xD>(ctr)),
260  static_cast<MCKL_INT64>(std::get<0xC>(ctr)));
261  __m128i xmmc7 =
262  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<0xF>(ctr)),
263  static_cast<MCKL_INT64>(std::get<0xE>(ctr)));
264  ctr.front() += n;
265 
266  __m128i *rptr = reinterpret_cast<__m128i *>(r);
267  while (n != 0) {
268  __m128i xmmt0 = _mm_add_epi64(xmmc0, _mm_set_epi64x(0, 1));
269  __m128i xmmt8 = _mm_add_epi64(xmmc0, _mm_set_epi64x(0, 2));
270  xmmc0 = xmmt8;
271 
272  __m128i xmmt1 = xmmc1;
273  __m128i xmmt2 = xmmc2;
274  __m128i xmmt3 = xmmc3;
275  __m128i xmmt4 = xmmc4;
276  __m128i xmmt5 = xmmc5;
277  __m128i xmmt6 = xmmc6;
278  __m128i xmmt7 = xmmc7;
279  __m128i xmmt9 = xmmc1;
280  __m128i xmmtA = xmmc2;
281  __m128i xmmtB = xmmc3;
282  __m128i xmmtC = xmmc4;
283  __m128i xmmtD = xmmc5;
284  __m128i xmmtE = xmmc6;
285  __m128i xmmtF = xmmc7;
286 
287  __m128i xmms0;
288  __m128i xmms1;
289  __m128i xmms2;
290  __m128i xmms3;
291  __m128i xmms4;
292  __m128i xmms5;
293  __m128i xmms6;
294  __m128i xmms7;
295  __m128i xmms8;
296  __m128i xmms9;
297  __m128i xmmsA;
298  __m128i xmmsB;
299  __m128i xmmsC;
300  __m128i xmmsD;
301  __m128i xmmsE;
302  __m128i xmmsF;
303 
304  xmms2 = xmmt2;
305  xmmt2 = xmmt1;
306  xmmt1 = xmmt8;
307  xmmt8 = xmmt4;
308  xmmt4 = xmms2;
309 
310  xmms6 = xmmt6;
311  xmmt6 = xmmt3;
312  xmmt3 = xmmt9;
313  xmmt9 = xmmtC;
314  xmmtC = xmms6;
315 
316  xmmsE = xmmtE;
317  xmmtE = xmmt7;
318  xmmt7 = xmmtB;
319  xmmtB = xmmtD;
320  xmmtD = xmmsE;
321 
322  xmms5 = xmmt5;
323  xmmt5 = xmmtA;
324  xmmtA = xmms5;
325 
326  xmms0 = _mm_unpacklo_epi64(xmmt0, xmmt1);
327  xmms1 = _mm_unpackhi_epi64(xmmt0, xmmt1);
328  xmms2 = _mm_unpacklo_epi64(xmmt2, xmmt3);
329  xmms3 = _mm_unpackhi_epi64(xmmt2, xmmt3);
330  xmms4 = _mm_unpacklo_epi64(xmmt4, xmmt5);
331  xmms5 = _mm_unpackhi_epi64(xmmt4, xmmt5);
332  xmms6 = _mm_unpacklo_epi64(xmmt6, xmmt7);
333  xmms7 = _mm_unpackhi_epi64(xmmt6, xmmt7);
334  xmms8 = _mm_unpacklo_epi64(xmmt8, xmmt9);
335  xmms9 = _mm_unpackhi_epi64(xmmt8, xmmt9);
336  xmmsA = _mm_unpacklo_epi64(xmmtA, xmmtB);
337  xmmsB = _mm_unpackhi_epi64(xmmtA, xmmtB);
338  xmmsC = _mm_unpacklo_epi64(xmmtC, xmmtD);
339  xmmsD = _mm_unpackhi_epi64(xmmtC, xmmtD);
340  xmmsE = _mm_unpacklo_epi64(xmmtE, xmmtF);
341  xmmsF = _mm_unpackhi_epi64(xmmtE, xmmtF);
342 
343  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_KBOX(0);
344  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_8(0);
345  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_8(1);
346  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_4(2);
347 
348  xmmt0 = _mm_unpacklo_epi64(xmms0, xmms1);
349  xmmt1 = _mm_unpackhi_epi64(xmms0, xmms1);
350  xmmt2 = _mm_unpacklo_epi64(xmms2, xmms3);
351  xmmt3 = _mm_unpackhi_epi64(xmms2, xmms3);
352  xmmt4 = _mm_unpacklo_epi64(xmms4, xmms5);
353  xmmt5 = _mm_unpackhi_epi64(xmms4, xmms5);
354  xmmt6 = _mm_unpacklo_epi64(xmms6, xmms7);
355  xmmt7 = _mm_unpackhi_epi64(xmms6, xmms7);
356  xmmt8 = _mm_unpacklo_epi64(xmms8, xmms9);
357  xmmt9 = _mm_unpackhi_epi64(xmms8, xmms9);
358  xmmtA = _mm_unpacklo_epi64(xmmsA, xmmsB);
359  xmmtB = _mm_unpackhi_epi64(xmmsA, xmmsB);
360  xmmtC = _mm_unpacklo_epi64(xmmsC, xmmsD);
361  xmmtD = _mm_unpackhi_epi64(xmmsC, xmmsD);
362  xmmtE = _mm_unpacklo_epi64(xmmsE, xmmsF);
363  xmmtF = _mm_unpackhi_epi64(xmmsE, xmmsF);
364 
365  xmms1 = xmmt1;
366  xmmt1 = xmmt2;
367  xmmt2 = xmmt4;
368  xmmt4 = xmmt8;
369  xmmt8 = xmms1;
370 
371  xmms3 = xmmt3;
372  xmmt3 = xmmt6;
373  xmmt6 = xmmtC;
374  xmmtC = xmmt9;
375  xmmt9 = xmms3;
376 
377  xmms7 = xmmt7;
378  xmmt7 = xmmtE;
379  xmmtE = xmmtD;
380  xmmtD = xmmtB;
381  xmmtB = xmms7;
382 
383  xmms5 = xmmt5;
384  xmmt5 = xmmtA;
385  xmmtA = xmms5;
386 
387  if (n >= N) {
388  n -= N;
389  _mm_storeu_si128(rptr++, xmmt0);
390  _mm_storeu_si128(rptr++, xmmt1);
391  _mm_storeu_si128(rptr++, xmmt2);
392  _mm_storeu_si128(rptr++, xmmt3);
393  _mm_storeu_si128(rptr++, xmmt4);
394  _mm_storeu_si128(rptr++, xmmt5);
395  _mm_storeu_si128(rptr++, xmmt6);
396  _mm_storeu_si128(rptr++, xmmt7);
397  _mm_storeu_si128(rptr++, xmmt8);
398  _mm_storeu_si128(rptr++, xmmt9);
399  _mm_storeu_si128(rptr++, xmmtA);
400  _mm_storeu_si128(rptr++, xmmtB);
401  _mm_storeu_si128(rptr++, xmmtC);
402  _mm_storeu_si128(rptr++, xmmtD);
403  _mm_storeu_si128(rptr++, xmmtE);
404  _mm_storeu_si128(rptr++, xmmtF);
405  } else {
406  std::array<__m128i, S> s;
407  std::get<0x0>(s) = xmmt0;
408  std::get<0x1>(s) = xmmt1;
409  std::get<0x2>(s) = xmmt2;
410  std::get<0x3>(s) = xmmt3;
411  std::get<0x4>(s) = xmmt4;
412  std::get<0x5>(s) = xmmt5;
413  std::get<0x6>(s) = xmmt6;
414  std::get<0x7>(s) = xmmt7;
415  std::get<0x8>(s) = xmmt8;
416  std::get<0x9>(s) = xmmt9;
417  std::get<0xA>(s) = xmmtA;
418  std::get<0xB>(s) = xmmtB;
419  std::get<0xC>(s) = xmmtC;
420  std::get<0xD>(s) = xmmtD;
421  std::get<0xE>(s) = xmmtE;
422  std::get<0xF>(s) = xmmtF;
423  std::memcpy(rptr, s.data(), n * sizeof(T) * K);
424  break;
425  }
426  }
427  }
428 }; // class Threefry16x64GeneratorSSE2Impl
429 
430 template <typename T>
431 class Threefish1024GeneratorSSE2Impl
432 {
433  static_assert(std::numeric_limits<T>::digits == 64,
434  "**Threefish1024GeneratorSSE2Impl** used with T other than a 64-bit "
435  "unsigned integers");
436 
437  static constexpr std::size_t K = 16;
438  static constexpr std::size_t Rounds = 80;
439  using Constants = ThreefryConstants<T, 16>;
440 
441  public:
442  static void eval(
443  const void *plain, void *cipher, const std::array<T, K + 4> &par)
444  {
445  Threefish1024GeneratorGenericImpl<T>::eval(plain, cipher, par);
446  }
447 
448  template <typename ResultType>
449  static void eval(std::array<std::uint64_t, 16> &ctr, ResultType *r,
450  const std::array<T, K + 4> &par)
451  {
452  Threefish1024GeneratorGenericImpl<T>::eval(ctr, r, par);
453  }
454 
455  template <typename ResultType>
456  static void eval(std::array<std::uint64_t, 16> &ctr, std::size_t n,
457  ResultType *r, const std::array<T, K + 4> &par)
458  {
459  constexpr std::size_t R = sizeof(T) * K / sizeof(ResultType);
460 
461  const std::size_t n0 =
462  static_cast<std::size_t>(std::min(static_cast<std::uint64_t>(n),
463  std::numeric_limits<std::uint64_t>::max() - ctr.front()));
464 
465  eval_kernel(ctr, n0, r, par);
466  n -= n0;
467  r += n0 * R;
468 
469  if (n != 0) {
470  eval(ctr, r, par);
471  n -= 1;
472  r += R;
473  }
474 
475  eval_kernel(ctr, n, r, par);
476  }
477 
478  private:
479  template <typename ResultType>
480  static void eval_kernel(std::array<std::uint64_t, 16> &ctr, std::size_t n,
481  ResultType *r, const std::array<T, K + 4> &par)
482  {
483  constexpr std::size_t S = 16;
484  constexpr std::size_t N = sizeof(__m128i) * S / (sizeof(T) * K);
485 
486  __m128i xmmc0 =
487  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<0x1>(ctr)),
488  static_cast<MCKL_INT64>(std::get<0x0>(ctr)));
489  __m128i xmmc1 =
490  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<0x3>(ctr)),
491  static_cast<MCKL_INT64>(std::get<0x2>(ctr)));
492  __m128i xmmc2 =
493  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<0x5>(ctr)),
494  static_cast<MCKL_INT64>(std::get<0x4>(ctr)));
495  __m128i xmmc3 =
496  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<0x7>(ctr)),
497  static_cast<MCKL_INT64>(std::get<0x6>(ctr)));
498  __m128i xmmc4 =
499  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<0x9>(ctr)),
500  static_cast<MCKL_INT64>(std::get<0x8>(ctr)));
501  __m128i xmmc5 =
502  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<0xB>(ctr)),
503  static_cast<MCKL_INT64>(std::get<0xA>(ctr)));
504  __m128i xmmc6 =
505  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<0xD>(ctr)),
506  static_cast<MCKL_INT64>(std::get<0xC>(ctr)));
507  __m128i xmmc7 =
508  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<0xF>(ctr)),
509  static_cast<MCKL_INT64>(std::get<0xE>(ctr)));
510  ctr.front() += n;
511 
512  __m128i *rptr = reinterpret_cast<__m128i *>(r);
513  while (n != 0) {
514  __m128i xmmt0 = _mm_add_epi64(xmmc0, _mm_set_epi64x(0, 1));
515  __m128i xmmt8 = _mm_add_epi64(xmmc0, _mm_set_epi64x(0, 2));
516  xmmc0 = xmmt8;
517 
518  __m128i xmmt1 = xmmc1;
519  __m128i xmmt2 = xmmc2;
520  __m128i xmmt3 = xmmc3;
521  __m128i xmmt4 = xmmc4;
522  __m128i xmmt5 = xmmc5;
523  __m128i xmmt6 = xmmc6;
524  __m128i xmmt7 = xmmc7;
525  __m128i xmmt9 = xmmc1;
526  __m128i xmmtA = xmmc2;
527  __m128i xmmtB = xmmc3;
528  __m128i xmmtC = xmmc4;
529  __m128i xmmtD = xmmc5;
530  __m128i xmmtE = xmmc6;
531  __m128i xmmtF = xmmc7;
532 
533  __m128i xmms0;
534  __m128i xmms1;
535  __m128i xmms2;
536  __m128i xmms3;
537  __m128i xmms4;
538  __m128i xmms5;
539  __m128i xmms6;
540  __m128i xmms7;
541  __m128i xmms8;
542  __m128i xmms9;
543  __m128i xmmsA;
544  __m128i xmmsB;
545  __m128i xmmsC;
546  __m128i xmmsD;
547  __m128i xmmsE;
548  __m128i xmmsF;
549 
550  xmms2 = xmmt2;
551  xmmt2 = xmmt1;
552  xmmt1 = xmmt8;
553  xmmt8 = xmmt4;
554  xmmt4 = xmms2;
555 
556  xmms6 = xmmt6;
557  xmmt6 = xmmt3;
558  xmmt3 = xmmt9;
559  xmmt9 = xmmtC;
560  xmmtC = xmms6;
561 
562  xmmsE = xmmtE;
563  xmmtE = xmmt7;
564  xmmt7 = xmmtB;
565  xmmtB = xmmtD;
566  xmmtD = xmmsE;
567 
568  xmms5 = xmmt5;
569  xmmt5 = xmmtA;
570  xmmtA = xmms5;
571 
572  xmms0 = _mm_unpacklo_epi64(xmmt0, xmmt1);
573  xmms1 = _mm_unpackhi_epi64(xmmt0, xmmt1);
574  xmms2 = _mm_unpacklo_epi64(xmmt2, xmmt3);
575  xmms3 = _mm_unpackhi_epi64(xmmt2, xmmt3);
576  xmms4 = _mm_unpacklo_epi64(xmmt4, xmmt5);
577  xmms5 = _mm_unpackhi_epi64(xmmt4, xmmt5);
578  xmms6 = _mm_unpacklo_epi64(xmmt6, xmmt7);
579  xmms7 = _mm_unpackhi_epi64(xmmt6, xmmt7);
580  xmms8 = _mm_unpacklo_epi64(xmmt8, xmmt9);
581  xmms9 = _mm_unpackhi_epi64(xmmt8, xmmt9);
582  xmmsA = _mm_unpacklo_epi64(xmmtA, xmmtB);
583  xmmsB = _mm_unpackhi_epi64(xmmtA, xmmtB);
584  xmmsC = _mm_unpacklo_epi64(xmmtC, xmmtD);
585  xmmsD = _mm_unpackhi_epi64(xmmtC, xmmtD);
586  xmmsE = _mm_unpacklo_epi64(xmmtE, xmmtF);
587  xmmsF = _mm_unpackhi_epi64(xmmtE, xmmtF);
588 
589  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_KBOX(0);
590  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_8(0);
591  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_8(1);
592  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_8(2);
593  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_8(3);
594  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_8(4);
595  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_8(5);
596  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_8(6);
597  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_8(7);
598  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_8(8);
599  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_8(9);
600 
601  xmmt0 = _mm_unpacklo_epi64(xmms0, xmms1);
602  xmmt1 = _mm_unpackhi_epi64(xmms0, xmms1);
603  xmmt2 = _mm_unpacklo_epi64(xmms2, xmms3);
604  xmmt3 = _mm_unpackhi_epi64(xmms2, xmms3);
605  xmmt4 = _mm_unpacklo_epi64(xmms4, xmms5);
606  xmmt5 = _mm_unpackhi_epi64(xmms4, xmms5);
607  xmmt6 = _mm_unpacklo_epi64(xmms6, xmms7);
608  xmmt7 = _mm_unpackhi_epi64(xmms6, xmms7);
609  xmmt8 = _mm_unpacklo_epi64(xmms8, xmms9);
610  xmmt9 = _mm_unpackhi_epi64(xmms8, xmms9);
611  xmmtA = _mm_unpacklo_epi64(xmmsA, xmmsB);
612  xmmtB = _mm_unpackhi_epi64(xmmsA, xmmsB);
613  xmmtC = _mm_unpacklo_epi64(xmmsC, xmmsD);
614  xmmtD = _mm_unpackhi_epi64(xmmsC, xmmsD);
615  xmmtE = _mm_unpacklo_epi64(xmmsE, xmmsF);
616  xmmtF = _mm_unpackhi_epi64(xmmsE, xmmsF);
617 
618  xmms1 = xmmt1;
619  xmmt1 = xmmt2;
620  xmmt2 = xmmt4;
621  xmmt4 = xmmt8;
622  xmmt8 = xmms1;
623 
624  xmms3 = xmmt3;
625  xmmt3 = xmmt6;
626  xmmt6 = xmmtC;
627  xmmtC = xmmt9;
628  xmmt9 = xmms3;
629 
630  xmms7 = xmmt7;
631  xmmt7 = xmmtE;
632  xmmtE = xmmtD;
633  xmmtD = xmmtB;
634  xmmtB = xmms7;
635 
636  xmms5 = xmmt5;
637  xmmt5 = xmmtA;
638  xmmtA = xmms5;
639 
640  if (n >= N) {
641  n -= N;
642  _mm_storeu_si128(rptr++, xmmt0);
643  _mm_storeu_si128(rptr++, xmmt1);
644  _mm_storeu_si128(rptr++, xmmt2);
645  _mm_storeu_si128(rptr++, xmmt3);
646  _mm_storeu_si128(rptr++, xmmt4);
647  _mm_storeu_si128(rptr++, xmmt5);
648  _mm_storeu_si128(rptr++, xmmt6);
649  _mm_storeu_si128(rptr++, xmmt7);
650  _mm_storeu_si128(rptr++, xmmt8);
651  _mm_storeu_si128(rptr++, xmmt9);
652  _mm_storeu_si128(rptr++, xmmtA);
653  _mm_storeu_si128(rptr++, xmmtB);
654  _mm_storeu_si128(rptr++, xmmtC);
655  _mm_storeu_si128(rptr++, xmmtD);
656  _mm_storeu_si128(rptr++, xmmtE);
657  _mm_storeu_si128(rptr++, xmmtF);
658  } else {
659  std::array<__m128i, S> s;
660  std::get<0x0>(s) = xmmt0;
661  std::get<0x1>(s) = xmmt1;
662  std::get<0x2>(s) = xmmt2;
663  std::get<0x3>(s) = xmmt3;
664  std::get<0x4>(s) = xmmt4;
665  std::get<0x5>(s) = xmmt5;
666  std::get<0x6>(s) = xmmt6;
667  std::get<0x7>(s) = xmmt7;
668  std::get<0x8>(s) = xmmt8;
669  std::get<0x9>(s) = xmmt9;
670  std::get<0xA>(s) = xmmtA;
671  std::get<0xB>(s) = xmmtB;
672  std::get<0xC>(s) = xmmtC;
673  std::get<0xD>(s) = xmmtD;
674  std::get<0xE>(s) = xmmtE;
675  std::get<0xF>(s) = xmmtF;
676  std::memcpy(rptr, s.data(), n * sizeof(T) * K);
677  break;
678  }
679  }
680  }
681 }; // class Threefish1024GeneratorSSE2Impl
682 
683 } // namespace internal
684 
685 } // namespace mckl
686 
687 MCKL_POP_GCC_WARNING
688 
689 #endif // MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_HPP
MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_8
#define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_8(N)
Definition: threefry_sse2_16x64.hpp:168
MCKL_PUSH_GCC_WARNING
#define MCKL_PUSH_GCC_WARNING(warning)
Definition: compiler.h:78
mckl::internal::Threefish1024GeneratorSSE2Impl::eval
static void eval(const void *plain, void *cipher, const std::array< T, K+4 > &par)
Definition: threefry_sse2_16x64.hpp:442
mckl::internal::Threefry16x64GeneratorSSE2Impl::eval
static void eval(const void *plain, void *cipher, const std::array< T, K+4 > &par)
Definition: threefry_sse2_16x64.hpp:196
mckl::internal::Threefry16x64GeneratorSSE2Impl::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_sse2_16x64.hpp:210
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::Threefish1024GeneratorSSE2Impl::eval
static void eval(std::array< std::uint64_t, 16 > &ctr, ResultType *r, const std::array< T, K+4 > &par)
Definition: threefry_sse2_16x64.hpp:449
mckl::internal::Threefry16x64GeneratorSSE2Impl::eval
static void eval(std::array< std::uint64_t, 16 > &ctr, ResultType *r, const std::array< T, K+4 > &par)
Definition: threefry_sse2_16x64.hpp:203
mckl::ThreefryConstants
Default Threefry constants.
Definition: threefry_constants.hpp:147
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_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_4
#define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_CYCLE_4(N)
Definition: threefry_sse2_16x64.hpp:161
mckl
Definition: mcmc.hpp:40
threefry_common.hpp
mckl::internal::Threefish1024GeneratorSSE2Impl::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_sse2_16x64.hpp:456
MCKL_POP_GCC_WARNING
#define MCKL_POP_GCC_WARNING
Definition: compiler.h:79
MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_KBOX
#define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_16X64_KBOX(N)
Definition: threefry_sse2_16x64.hpp:43
Generated by   doxygen 1.8.13