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