MCKL
Monte Carlo Kernel Library
threefry_sse2_2x32.hpp
Go to the documentation of this file.
1 //============================================================================
2 // MCKL/include/mckl/random/internal/threefry_sse2_2x32.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_2X32_HPP
33 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_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_2x32.hpp>
39 
40 MCKL_PUSH_GCC_WARNING("-Wignored-attributes")
41 
42 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_KBOX(N) \
43  xmmt0 = _mm_set1_epi32( \
44  static_cast<int>(ThreefryKBox<T, K, N>::template key<0>(par))); \
45  xmmt1 = _mm_set1_epi32( \
46  static_cast<int>(ThreefryKBox<T, K, N>::template key<1>(par))); \
47  xmms0 = _mm_add_epi32(xmms0, xmmt0); \
48  xmms1 = _mm_add_epi32(xmms1, xmmt1); \
49  xmms2 = _mm_add_epi32(xmms2, xmmt0); \
50  xmms3 = _mm_add_epi32(xmms3, xmmt1); \
51  xmms4 = _mm_add_epi32(xmms4, xmmt0); \
52  xmms5 = _mm_add_epi32(xmms5, xmmt1); \
53  xmms6 = _mm_add_epi32(xmms6, xmmt0); \
54  xmms7 = _mm_add_epi32(xmms7, xmmt1);
55 
56 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_RBOX(N) \
57  { \
58  constexpr int L0 = Constants::rotate::value[0][(N - 1) % 8]; \
59  constexpr int R0 = 32 - L0; \
60  xmms0 = _mm_add_epi32(xmms0, xmms1); \
61  xmms2 = _mm_add_epi32(xmms2, xmms3); \
62  xmms4 = _mm_add_epi32(xmms4, xmms5); \
63  xmms6 = _mm_add_epi32(xmms6, xmms7); \
64  xmmt1 = _mm_slli_epi32(xmms1, L0); \
65  xmmt3 = _mm_slli_epi32(xmms3, L0); \
66  xmmt5 = _mm_slli_epi32(xmms5, L0); \
67  xmmt7 = _mm_slli_epi32(xmms7, L0); \
68  xmms1 = _mm_srli_epi32(xmms1, R0); \
69  xmms3 = _mm_srli_epi32(xmms3, R0); \
70  xmms5 = _mm_srli_epi32(xmms5, R0); \
71  xmms7 = _mm_srli_epi32(xmms7, R0); \
72  xmmt1 = _mm_or_si128(xmms1, xmmt1); \
73  xmmt3 = _mm_or_si128(xmms3, xmmt3); \
74  xmmt5 = _mm_or_si128(xmms5, xmmt5); \
75  xmmt7 = _mm_or_si128(xmms7, xmmt7); \
76  xmms1 = _mm_xor_si128(xmms0, xmmt1); \
77  xmms3 = _mm_xor_si128(xmms2, xmmt3); \
78  xmms5 = _mm_xor_si128(xmms4, xmmt5); \
79  xmms7 = _mm_xor_si128(xmms6, xmmt7); \
80  }
81 
82 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_CYCLE_4(N) \
83  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_RBOX(N * 8 + 1); \
84  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_RBOX(N * 8 + 2); \
85  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_RBOX(N * 8 + 3); \
86  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_RBOX(N * 8 + 4); \
87  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_KBOX(N * 8 + 4);
88 
89 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_CYCLE_8(N) \
90  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_RBOX(N * 8 + 1); \
91  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_RBOX(N * 8 + 2); \
92  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_RBOX(N * 8 + 3); \
93  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_RBOX(N * 8 + 4); \
94  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_KBOX(N * 8 + 4); \
95  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_RBOX(N * 8 + 5); \
96  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_RBOX(N * 8 + 6); \
97  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_RBOX(N * 8 + 7); \
98  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_RBOX(N * 8 + 8); \
99  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_KBOX(N * 8 + 8);
100 
101 namespace mckl {
102 
103 namespace internal {
104 
105 template <typename T>
106 class Threefry2x32GeneratorSSE2Impl
107 {
108  static_assert(std::numeric_limits<T>::digits == 32,
109  "**Threefry2x32GeneratorSSE2Impl** used with T other than a "
110  "32-bit unsigned integers");
111 
112  static constexpr std::size_t K = 2;
113  static constexpr std::size_t Rounds = 20;
114  using Constants = ThreefryConstants<T, 2>;
115 
116  public:
117  static void eval(
118  const void *plain, void *cipher, const std::array<T, K + 4> &par)
119  {
120  Threefry2x32GeneratorGenericImpl<T>::eval(plain, cipher, par);
121  }
122 
123  template <typename ResultType>
124  static void eval(std::array<std::uint64_t, 1> &ctr, ResultType *r,
125  const std::array<T, K + 4> &par)
126  {
127  Threefry2x32GeneratorGenericImpl<T>::eval(ctr, r, par);
128  }
129 
130  template <typename ResultType>
131  static void eval(std::array<std::uint64_t, 1> &ctr, std::size_t n,
132  ResultType *r, const std::array<T, K + 4> &par)
133  {
134  constexpr std::size_t R = sizeof(T) * K / sizeof(ResultType);
135 
136  const std::size_t n0 =
137  static_cast<std::size_t>(std::min(static_cast<std::uint64_t>(n),
138  std::numeric_limits<std::uint64_t>::max() - ctr.front()));
139 
140  eval_kernel(ctr, n0, r, par);
141  n -= n0;
142  r += n0 * R;
143 
144  if (n != 0) {
145  eval(ctr, r, par);
146  n -= 1;
147  r += R;
148  }
149 
150  eval_kernel(ctr, n, r, par);
151  }
152 
153  private:
154  template <typename ResultType>
155  static void eval_kernel(std::array<std::uint64_t, 1> &ctr, std::size_t n,
156  ResultType *r, const std::array<T, K + 4> &par)
157  {
158  constexpr std::size_t S = 8;
159  constexpr std::size_t N = sizeof(__m128i) * S / (sizeof(T) * K);
160 
161  __m128i xmmc =
162  _mm_set1_epi64x(static_cast<MCKL_INT64>(std::get<0>(ctr)));
163  ctr.front() += n;
164 
165  __m128i *rptr = reinterpret_cast<__m128i *>(r);
166  while (n != 0) {
167  __m128i xmms0 = _mm_add_epi64(xmmc, _mm_set_epi64x(0x02, 0x01));
168  __m128i xmms1 = _mm_add_epi64(xmmc, _mm_set_epi64x(0x04, 0x03));
169  __m128i xmms2 = _mm_add_epi64(xmmc, _mm_set_epi64x(0x06, 0x05));
170  __m128i xmms3 = _mm_add_epi64(xmmc, _mm_set_epi64x(0x08, 0x07));
171  __m128i xmms4 = _mm_add_epi64(xmmc, _mm_set_epi64x(0x0A, 0x09));
172  __m128i xmms5 = _mm_add_epi64(xmmc, _mm_set_epi64x(0x0C, 0x0B));
173  __m128i xmms6 = _mm_add_epi64(xmmc, _mm_set_epi64x(0x0E, 0x0D));
174  __m128i xmms7 = _mm_add_epi64(xmmc, _mm_set_epi64x(0x10, 0x0F));
175  xmmc = _mm_add_epi64(xmmc, _mm_set1_epi64x(0x10));
176 
177  __m128i xmmt0;
178  __m128i xmmt1;
179  __m128i xmmt2;
180  __m128i xmmt3;
181  __m128i xmmt4;
182  __m128i xmmt5;
183  __m128i xmmt6;
184  __m128i xmmt7;
185 
186  xmmt1 = xmms1;
187  xmms1 = xmms2;
188  xmms2 = xmms4;
189  xmms4 = xmmt1;
190 
191  xmmt3 = xmms3;
192  xmms3 = xmms6;
193  xmms6 = xmms5;
194  xmms5 = xmmt3;
195 
196  xmmt0 = _mm_unpacklo_epi32(xmms0, xmms1);
197  xmmt1 = _mm_unpacklo_epi32(xmms2, xmms3);
198  xmmt2 = _mm_unpackhi_epi32(xmms0, xmms1);
199  xmmt3 = _mm_unpackhi_epi32(xmms2, xmms3);
200  xmmt4 = _mm_unpacklo_epi32(xmms4, xmms5);
201  xmmt5 = _mm_unpacklo_epi32(xmms6, xmms7);
202  xmmt6 = _mm_unpackhi_epi32(xmms4, xmms5);
203  xmmt7 = _mm_unpackhi_epi32(xmms6, xmms7);
204  xmms0 = _mm_unpacklo_epi64(xmmt0, xmmt1);
205  xmms1 = _mm_unpackhi_epi64(xmmt0, xmmt1);
206  xmms2 = _mm_unpacklo_epi64(xmmt2, xmmt3);
207  xmms3 = _mm_unpackhi_epi64(xmmt2, xmmt3);
208  xmms4 = _mm_unpacklo_epi64(xmmt4, xmmt5);
209  xmms5 = _mm_unpackhi_epi64(xmmt4, xmmt5);
210  xmms6 = _mm_unpacklo_epi64(xmmt6, xmmt7);
211  xmms7 = _mm_unpackhi_epi64(xmmt6, xmmt7);
212 
213  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_KBOX(0)
214  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_CYCLE_8(0)
215  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_CYCLE_8(1)
216  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_CYCLE_4(2)
217 
218  xmmt0 = _mm_unpacklo_epi32(xmms0, xmms1);
219  xmmt1 = _mm_unpacklo_epi32(xmms2, xmms3);
220  xmmt2 = _mm_unpackhi_epi32(xmms0, xmms1);
221  xmmt3 = _mm_unpackhi_epi32(xmms2, xmms3);
222  xmmt4 = _mm_unpacklo_epi32(xmms4, xmms5);
223  xmmt5 = _mm_unpacklo_epi32(xmms6, xmms7);
224  xmmt6 = _mm_unpackhi_epi32(xmms4, xmms5);
225  xmmt7 = _mm_unpackhi_epi32(xmms6, xmms7);
226  xmms0 = _mm_unpacklo_epi64(xmmt0, xmmt1);
227  xmms1 = _mm_unpackhi_epi64(xmmt0, xmmt1);
228  xmms2 = _mm_unpacklo_epi64(xmmt2, xmmt3);
229  xmms3 = _mm_unpackhi_epi64(xmmt2, xmmt3);
230  xmms4 = _mm_unpacklo_epi64(xmmt4, xmmt5);
231  xmms5 = _mm_unpackhi_epi64(xmmt4, xmmt5);
232  xmms6 = _mm_unpacklo_epi64(xmmt6, xmmt7);
233  xmms7 = _mm_unpackhi_epi64(xmmt6, xmmt7);
234 
235  xmmt2 = xmms2;
236  xmms2 = xmms1;
237  xmms1 = xmms4;
238  xmms4 = xmmt2;
239 
240  xmmt6 = xmms6;
241  xmms6 = xmms3;
242  xmms3 = xmms5;
243  xmms5 = xmmt6;
244 
245  if (n >= N) {
246  n -= N;
247  _mm_storeu_si128(rptr++, xmms0);
248  _mm_storeu_si128(rptr++, xmms1);
249  _mm_storeu_si128(rptr++, xmms2);
250  _mm_storeu_si128(rptr++, xmms3);
251  _mm_storeu_si128(rptr++, xmms4);
252  _mm_storeu_si128(rptr++, xmms5);
253  _mm_storeu_si128(rptr++, xmms6);
254  _mm_storeu_si128(rptr++, xmms7);
255  } else {
256  std::array<__m128i, S> s;
257  std::get<0>(s) = xmms0;
258  std::get<1>(s) = xmms1;
259  std::get<2>(s) = xmms2;
260  std::get<3>(s) = xmms3;
261  std::get<4>(s) = xmms4;
262  std::get<5>(s) = xmms5;
263  std::get<6>(s) = xmms6;
264  std::get<7>(s) = xmms7;
265  std::memcpy(rptr, s.data(), n * sizeof(T) * K);
266  break;
267  }
268  }
269  }
270 }; // class Threefry2x32GeneratorSSE2Impl
271 
272 } // namespace internal
273 
274 } // namespace mckl
275 
276 MCKL_POP_GCC_WARNING
277 
278 #endif // MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_HPP
MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_CYCLE_8
#define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_CYCLE_8(N)
Definition: threefry_sse2_2x32.hpp:89
threefry_generic_2x32.hpp
MCKL_PUSH_GCC_WARNING
#define MCKL_PUSH_GCC_WARNING(warning)
Definition: compiler.h:78
mckl::internal::Threefry2x32GeneratorSSE2Impl::eval
static void eval(std::array< std::uint64_t, 1 > &ctr, std::size_t n, ResultType *r, const std::array< T, K+4 > &par)
Definition: threefry_sse2_2x32.hpp:131
mckl::internal::Threefry2x32GeneratorGenericImpl::eval
static void eval(const void *plain, void *cipher, const std::array< T, K+4 > &par)
Definition: threefry_generic_2x32.hpp:98
MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_KBOX
#define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_KBOX(N)
Definition: threefry_sse2_2x32.hpp:42
mckl::internal::Threefry2x32GeneratorSSE2Impl::eval
static void eval(std::array< std::uint64_t, 1 > &ctr, ResultType *r, const std::array< T, K+4 > &par)
Definition: threefry_sse2_2x32.hpp:124
mckl::ThreefryConstants
Default Threefry constants.
Definition: threefry_constants.hpp:147
mckl::internal::Threefry2x32GeneratorSSE2Impl::eval
static void eval(const void *plain, void *cipher, const std::array< T, K+4 > &par)
Definition: threefry_sse2_2x32.hpp:117
MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_CYCLE_4
#define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X32_CYCLE_4(N)
Definition: threefry_sse2_2x32.hpp:82
threefry_constants.hpp
mckl
Definition: mcmc.hpp:40
threefry_common.hpp
MCKL_POP_GCC_WARNING
#define MCKL_POP_GCC_WARNING
Definition: compiler.h:79
Generated by   doxygen 1.8.13