MCKL
Monte Carlo Kernel Library
threefry_sse2_8x64.hpp
Go to the documentation of this file.
1 //============================================================================
2 // MCKL/include/mckl/random/internal/threefry_sse2_8x64.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_8X64_HPP
33 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_HPP
34 
40 
41 MCKL_PUSH_GCC_WARNING("-Wignored-attributes")
42 
43 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_KBOX(N) \
44  xmmt0 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
45  ThreefryKBox<T, K, N>::template key<0>(par))); \
46  xmmt1 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
47  ThreefryKBox<T, K, N>::template key<1>(par))); \
48  xmmt2 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
49  ThreefryKBox<T, K, N>::template key<2>(par))); \
50  xmmt3 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
51  ThreefryKBox<T, K, N>::template key<3>(par))); \
52  xmmt4 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
53  ThreefryKBox<T, K, N>::template key<4>(par))); \
54  xmmt5 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
55  ThreefryKBox<T, K, N>::template key<5>(par))); \
56  xmmt6 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
57  ThreefryKBox<T, K, N>::template key<6>(par))); \
58  xmmt7 = _mm_set1_epi64x(static_cast<MCKL_INT64>( \
59  ThreefryKBox<T, K, N>::template key<7>(par))); \
60  xmms0 = _mm_add_epi64(xmms0, xmmt0); \
61  xmms1 = _mm_add_epi64(xmms1, xmmt1); \
62  xmms2 = _mm_add_epi64(xmms2, xmmt2); \
63  xmms3 = _mm_add_epi64(xmms3, xmmt3); \
64  xmms4 = _mm_add_epi64(xmms4, xmmt4); \
65  xmms5 = _mm_add_epi64(xmms5, xmmt5); \
66  xmms6 = _mm_add_epi64(xmms6, xmmt6); \
67  xmms7 = _mm_add_epi64(xmms7, xmmt7);
68 
69 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_RBOX(N) \
70  { \
71  constexpr int L0 = Constants::rotate::value[0][(N - 1) % 8]; \
72  constexpr int L1 = Constants::rotate::value[1][(N - 1) % 8]; \
73  constexpr int L2 = Constants::rotate::value[2][(N - 1) % 8]; \
74  constexpr int L3 = Constants::rotate::value[3][(N - 1) % 8]; \
75  constexpr int R0 = 64 - L0; \
76  constexpr int R1 = 64 - L1; \
77  constexpr int R2 = 64 - L2; \
78  constexpr int R3 = 64 - L3; \
79  xmmt0 = _mm_add_epi64(xmms0, xmms1); \
80  xmmt2 = _mm_add_epi64(xmms2, xmms3); \
81  xmmt4 = _mm_add_epi64(xmms4, xmms5); \
82  xmmt6 = _mm_add_epi64(xmms6, xmms7); \
83  xmmt1 = _mm_slli_epi64(xmms1, L0); \
84  xmmt3 = _mm_slli_epi64(xmms3, L1); \
85  xmmt5 = _mm_slli_epi64(xmms5, L2); \
86  xmmt7 = _mm_slli_epi64(xmms7, L3); \
87  xmms1 = _mm_srli_epi64(xmms1, R0); \
88  xmms3 = _mm_srli_epi64(xmms3, R1); \
89  xmms5 = _mm_srli_epi64(xmms5, R2); \
90  xmms7 = _mm_srli_epi64(xmms7, R3); \
91  xmmt1 = _mm_or_si128(xmms1, xmmt1); \
92  xmmt3 = _mm_or_si128(xmms3, xmmt3); \
93  xmmt5 = _mm_or_si128(xmms5, xmmt5); \
94  xmmt7 = _mm_or_si128(xmms7, xmmt7); \
95  xmms0 = xmmt2; \
96  xmms2 = xmmt4; \
97  xmms4 = xmmt6; \
98  xmms6 = xmmt0; \
99  xmms1 = _mm_xor_si128(xmmt0, xmmt1); \
100  xmms3 = _mm_xor_si128(xmmt6, xmmt7); \
101  xmms5 = _mm_xor_si128(xmmt4, xmmt5); \
102  xmms7 = _mm_xor_si128(xmmt2, xmmt3); \
103  }
104 
105 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_CYCLE_4(N) \
106  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_RBOX(N * 8 + 1); \
107  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_RBOX(N * 8 + 2); \
108  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_RBOX(N * 8 + 3); \
109  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_RBOX(N * 8 + 4); \
110  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_KBOX(N * 8 + 4);
111 
112 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_CYCLE_8(N) \
113  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_RBOX(N * 8 + 1); \
114  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_RBOX(N * 8 + 2); \
115  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_RBOX(N * 8 + 3); \
116  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_RBOX(N * 8 + 4); \
117  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_KBOX(N * 8 + 4); \
118  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_RBOX(N * 8 + 5); \
119  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_RBOX(N * 8 + 6); \
120  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_RBOX(N * 8 + 7); \
121  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_RBOX(N * 8 + 8); \
122  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_KBOX(N * 8 + 8);
123 
124 namespace mckl {
125 
126 namespace internal {
127 
128 template <typename T>
130 {
131  static_assert(std::numeric_limits<T>::digits == 64,
132  "**Threefry8x64GeneratorSSE2Impl** used with T other than a 64-bit "
133  "unsigned integers");
134 
135  static constexpr std::size_t K = 8;
136  static constexpr std::size_t Rounds = 20;
138 
139  public:
140  static void eval(
141  const void *plain, void *cipher, const std::array<T, K + 4> &par)
142  {
143  Threefry8x64GeneratorGenericImpl<T>::eval(plain, cipher, par);
144  }
145 
146  template <typename ResultType>
147  static void eval(std::array<std::uint64_t, 8> &ctr, ResultType *r,
148  const std::array<T, K + 4> &par)
149  {
151  }
152 
153  template <typename ResultType>
154  static void eval(std::array<std::uint64_t, 8> &ctr, std::size_t n,
155  ResultType *r, const std::array<T, K + 4> &par)
156  {
157  constexpr std::size_t R = sizeof(T) * K / sizeof(ResultType);
158 
159  const std::size_t n0 =
160  static_cast<std::size_t>(std::min(static_cast<std::uint64_t>(n),
161  std::numeric_limits<std::uint64_t>::max() - ctr.front()));
162 
163  eval_kernel(ctr, n0, r, par);
164  n -= n0;
165  r += n0 * R;
166 
167  if (n != 0) {
168  eval(ctr, r, par);
169  n -= 1;
170  r += R;
171  }
172 
173  eval_kernel(ctr, n, r, par);
174  }
175 
176  private:
177  template <typename ResultType>
178  static void eval_kernel(std::array<std::uint64_t, 8> &ctr, std::size_t n,
179  ResultType *r, const std::array<T, K + 4> &par)
180  {
181  constexpr std::size_t S = 8;
182  constexpr std::size_t N = sizeof(__m128i) * S / (sizeof(T) * K);
183 
184  __m128i xmmc0 =
185  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<1>(ctr)),
186  static_cast<MCKL_INT64>(std::get<0>(ctr)));
187  __m128i xmmc1 =
188  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<3>(ctr)),
189  static_cast<MCKL_INT64>(std::get<2>(ctr)));
190  __m128i xmmc2 =
191  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<5>(ctr)),
192  static_cast<MCKL_INT64>(std::get<4>(ctr)));
193  __m128i xmmc3 =
194  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<7>(ctr)),
195  static_cast<MCKL_INT64>(std::get<6>(ctr)));
196  ctr.front() += n;
197 
198  __m128i *rptr = reinterpret_cast<__m128i *>(r);
199  while (n != 0) {
200  __m128i xmmt0 = _mm_add_epi64(xmmc0, _mm_set_epi64x(0, 1));
201  __m128i xmmt4 = _mm_add_epi64(xmmc0, _mm_set_epi64x(0, 2));
202  xmmc0 = xmmt4;
203 
204  __m128i xmmt1 = xmmc1;
205  __m128i xmmt2 = xmmc2;
206  __m128i xmmt3 = xmmc3;
207  __m128i xmmt5 = xmmc1;
208  __m128i xmmt6 = xmmc2;
209  __m128i xmmt7 = xmmc3;
210 
211  __m128i xmms0;
212  __m128i xmms1;
213  __m128i xmms2;
214  __m128i xmms3;
215  __m128i xmms4;
216  __m128i xmms5;
217  __m128i xmms6;
218  __m128i xmms7;
219 
220  xmms2 = xmmt2;
221  xmmt2 = xmmt1;
222  xmmt1 = xmmt4;
223  xmmt4 = xmms2;
224 
225  xmms6 = xmmt6;
226  xmmt6 = xmmt3;
227  xmmt3 = xmmt5;
228  xmmt5 = xmms6;
229 
230  xmms0 = _mm_unpacklo_epi64(xmmt0, xmmt1);
231  xmms1 = _mm_unpackhi_epi64(xmmt0, xmmt1);
232  xmms2 = _mm_unpacklo_epi64(xmmt2, xmmt3);
233  xmms3 = _mm_unpackhi_epi64(xmmt2, xmmt3);
234  xmms4 = _mm_unpacklo_epi64(xmmt4, xmmt5);
235  xmms5 = _mm_unpackhi_epi64(xmmt4, xmmt5);
236  xmms6 = _mm_unpacklo_epi64(xmmt6, xmmt7);
237  xmms7 = _mm_unpackhi_epi64(xmmt6, xmmt7);
238 
243 
244  xmmt0 = _mm_unpacklo_epi64(xmms0, xmms1);
245  xmmt1 = _mm_unpackhi_epi64(xmms0, xmms1);
246  xmmt2 = _mm_unpacklo_epi64(xmms2, xmms3);
247  xmmt3 = _mm_unpackhi_epi64(xmms2, xmms3);
248  xmmt4 = _mm_unpacklo_epi64(xmms4, xmms5);
249  xmmt5 = _mm_unpackhi_epi64(xmms4, xmms5);
250  xmmt6 = _mm_unpacklo_epi64(xmms6, xmms7);
251  xmmt7 = _mm_unpackhi_epi64(xmms6, xmms7);
252 
253  xmms1 = xmmt1;
254  xmmt1 = xmmt2;
255  xmmt2 = xmmt4;
256  xmmt4 = xmms1;
257 
258  xmms3 = xmmt3;
259  xmmt3 = xmmt6;
260  xmmt6 = xmmt5;
261  xmmt5 = xmms3;
262 
263  if (n >= N) {
264  n -= N;
265  _mm_storeu_si128(rptr++, xmmt0);
266  _mm_storeu_si128(rptr++, xmmt1);
267  _mm_storeu_si128(rptr++, xmmt2);
268  _mm_storeu_si128(rptr++, xmmt3);
269  _mm_storeu_si128(rptr++, xmmt4);
270  _mm_storeu_si128(rptr++, xmmt5);
271  _mm_storeu_si128(rptr++, xmmt6);
272  _mm_storeu_si128(rptr++, xmmt7);
273  } else {
274  std::array<__m128i, S> s;
275  std::get<0>(s) = xmmt0;
276  std::get<1>(s) = xmmt1;
277  std::get<2>(s) = xmmt2;
278  std::get<3>(s) = xmmt3;
279  std::get<4>(s) = xmmt4;
280  std::get<5>(s) = xmmt5;
281  std::get<6>(s) = xmmt6;
282  std::get<7>(s) = xmmt7;
283  std::memcpy(rptr, s.data(), n * sizeof(T) * K);
284  break;
285  }
286  }
287  }
288 }; // class Threefry8x64GeneratorSSE2Impl
289 
290 template <typename T>
292 {
293  static_assert(std::numeric_limits<T>::digits == 64,
294  "**Threefish512GeneratorSSE2Impl** used with T other than a 64-bit "
295  "unsigned integers");
296 
297  static constexpr std::size_t K = 8;
298  static constexpr std::size_t Rounds = 20;
300 
301  public:
302  static void eval(
303  const void *plain, void *cipher, const std::array<T, K + 4> &par)
304  {
305  Threefish512GeneratorGenericImpl<T>::eval(plain, cipher, par);
306  }
307 
308  template <typename ResultType>
309  static void eval(std::array<std::uint64_t, 8> &ctr, ResultType *r,
310  const std::array<T, K + 4> &par)
311  {
313  }
314 
315  template <typename ResultType>
316  static void eval(std::array<std::uint64_t, 8> &ctr, std::size_t n,
317  ResultType *r, const std::array<T, K + 4> &par)
318  {
319  constexpr std::size_t R = sizeof(T) * K / sizeof(ResultType);
320 
321  const std::size_t n0 =
322  static_cast<std::size_t>(std::min(static_cast<std::uint64_t>(n),
323  std::numeric_limits<std::uint64_t>::max() - ctr.front()));
324 
325  eval_kernel(ctr, n0, r, par);
326  n -= n0;
327  r += n0 * R;
328 
329  if (n != 0) {
330  eval(ctr, r, par);
331  n -= 1;
332  r += R;
333  }
334 
335  eval_kernel(ctr, n, r, par);
336  }
337 
338  template <typename ResultType>
339  static void eval_kernel(std::array<std::uint64_t, 8> &ctr, std::size_t n,
340  ResultType *r, const std::array<T, K + 4> &par)
341  {
342  constexpr std::size_t S = 8;
343  constexpr std::size_t N = sizeof(__m128i) * S / (sizeof(T) * K);
344 
345  __m128i xmmc0 =
346  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<1>(ctr)),
347  static_cast<MCKL_INT64>(std::get<0>(ctr)));
348  __m128i xmmc1 =
349  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<3>(ctr)),
350  static_cast<MCKL_INT64>(std::get<2>(ctr)));
351  __m128i xmmc2 =
352  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<5>(ctr)),
353  static_cast<MCKL_INT64>(std::get<4>(ctr)));
354  __m128i xmmc3 =
355  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<7>(ctr)),
356  static_cast<MCKL_INT64>(std::get<6>(ctr)));
357  ctr.front() += n;
358 
359  __m128i *rptr = reinterpret_cast<__m128i *>(r);
360  while (n != 0) {
361  __m128i xmmt0 = _mm_add_epi64(xmmc0, _mm_set_epi64x(0, 1));
362  __m128i xmmt4 = _mm_add_epi64(xmmc0, _mm_set_epi64x(0, 2));
363  xmmc0 = xmmt4;
364 
365  __m128i xmmt1 = xmmc1;
366  __m128i xmmt2 = xmmc2;
367  __m128i xmmt3 = xmmc3;
368  __m128i xmmt5 = xmmc1;
369  __m128i xmmt6 = xmmc2;
370  __m128i xmmt7 = xmmc3;
371 
372  __m128i xmms0;
373  __m128i xmms1;
374  __m128i xmms2;
375  __m128i xmms3;
376  __m128i xmms4;
377  __m128i xmms5;
378  __m128i xmms6;
379  __m128i xmms7;
380 
381  xmms2 = xmmt2;
382  xmmt2 = xmmt1;
383  xmmt1 = xmmt4;
384  xmmt4 = xmms2;
385 
386  xmms6 = xmmt6;
387  xmmt6 = xmmt3;
388  xmmt3 = xmmt5;
389  xmmt5 = xmms6;
390 
391  xmms0 = _mm_unpacklo_epi64(xmmt0, xmmt1);
392  xmms1 = _mm_unpackhi_epi64(xmmt0, xmmt1);
393  xmms2 = _mm_unpacklo_epi64(xmmt2, xmmt3);
394  xmms3 = _mm_unpackhi_epi64(xmmt2, xmmt3);
395  xmms4 = _mm_unpacklo_epi64(xmmt4, xmmt5);
396  xmms5 = _mm_unpackhi_epi64(xmmt4, xmmt5);
397  xmms6 = _mm_unpacklo_epi64(xmmt6, xmmt7);
398  xmms7 = _mm_unpackhi_epi64(xmmt6, xmmt7);
399 
410 
411  xmmt0 = _mm_unpacklo_epi64(xmms0, xmms1);
412  xmmt1 = _mm_unpackhi_epi64(xmms0, xmms1);
413  xmmt2 = _mm_unpacklo_epi64(xmms2, xmms3);
414  xmmt3 = _mm_unpackhi_epi64(xmms2, xmms3);
415  xmmt4 = _mm_unpacklo_epi64(xmms4, xmms5);
416  xmmt5 = _mm_unpackhi_epi64(xmms4, xmms5);
417  xmmt6 = _mm_unpacklo_epi64(xmms6, xmms7);
418  xmmt7 = _mm_unpackhi_epi64(xmms6, xmms7);
419 
420  xmms1 = xmmt1;
421  xmmt1 = xmmt2;
422  xmmt2 = xmmt4;
423  xmmt4 = xmms1;
424 
425  xmms3 = xmmt3;
426  xmmt3 = xmmt6;
427  xmmt6 = xmmt5;
428  xmmt5 = xmms3;
429 
430  if (n >= N) {
431  n -= N;
432  _mm_storeu_si128(rptr++, xmmt0);
433  _mm_storeu_si128(rptr++, xmmt1);
434  _mm_storeu_si128(rptr++, xmmt2);
435  _mm_storeu_si128(rptr++, xmmt3);
436  _mm_storeu_si128(rptr++, xmmt4);
437  _mm_storeu_si128(rptr++, xmmt5);
438  _mm_storeu_si128(rptr++, xmmt6);
439  _mm_storeu_si128(rptr++, xmmt7);
440  } else {
441  std::array<__m128i, S> s;
442  std::get<0>(s) = xmmt0;
443  std::get<1>(s) = xmmt1;
444  std::get<2>(s) = xmmt2;
445  std::get<3>(s) = xmmt3;
446  std::get<4>(s) = xmmt4;
447  std::get<5>(s) = xmmt5;
448  std::get<6>(s) = xmmt6;
449  std::get<7>(s) = xmmt7;
450  std::memcpy(rptr, s.data(), n * sizeof(T) * K);
451  break;
452  }
453  }
454  }
455 }; // class Threefish512GeneratorSSE2Impl
456 
457 } // namespace internal
458 
459 } // namespace mckl
460 
462 
463 #endif // MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_HPP
#define MCKL_PUSH_GCC_WARNING(warning)
Definition: compiler.h:78
static void eval(const void *plain, void *cipher, const std::array< T, K+4 > &par)
static void eval(std::array< std::uint64_t, 8 > &ctr, ResultType *r, const std::array< T, K+4 > &par)
static void eval(const void *plain, void *cipher, const std::array< T, K+4 > &par)
static void eval(std::array< std::uint64_t, 8 > &ctr, std::size_t n, ResultType *r, const std::array< T, K+4 > &par)
#define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_CYCLE_8(N)
#define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_CYCLE_4(N)
Default Threefry constants.
static void eval(std::array< std::uint64_t, 8 > &ctr, std::size_t n, ResultType *r, const std::array< T, K+4 > &par)
Definition: mcmc.hpp:40
static void eval(std::array< std::uint64_t, 8 > &ctr, ResultType *r, const std::array< T, K+4 > &par)
static void eval(const void *plain, void *cipher, const std::array< T, K+4 > &par)
static void eval_kernel(std::array< std::uint64_t, 8 > &ctr, std::size_t n, ResultType *r, const std::array< T, K+4 > &par)
static void eval(const void *plain, void *cipher, const std::array< T, K+4 > &par)
#define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_8X64_KBOX(N)
#define MCKL_POP_GCC_WARNING
Definition: compiler.h:79