MCKL
Monte Carlo Kernel Library
threefry_sse2_2x64.hpp
Go to the documentation of this file.
1 //============================================================================
2 // MCKL/include/mckl/random/internal/threefry_sse2_2x64.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_2X64_HPP
33 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_HPP
34 
40 
41 MCKL_PUSH_GCC_WARNING("-Wignored-attributes")
42 
43 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_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  xmms0 = _mm_add_epi64(xmms0, xmmt0); \
49  xmms1 = _mm_add_epi64(xmms1, xmmt1); \
50  xmms2 = _mm_add_epi64(xmms2, xmmt0); \
51  xmms3 = _mm_add_epi64(xmms3, xmmt1); \
52  xmms4 = _mm_add_epi64(xmms4, xmmt0); \
53  xmms5 = _mm_add_epi64(xmms5, xmmt1); \
54  xmms6 = _mm_add_epi64(xmms6, xmmt0); \
55  xmms7 = _mm_add_epi64(xmms7, xmmt1);
56 
57 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_RBOX(N) \
58  { \
59  constexpr int L0 = Constants::rotate::value[0][(N - 1) % 8]; \
60  constexpr int R0 = 64 - L0; \
61  xmms0 = _mm_add_epi64(xmms0, xmms1); \
62  xmms2 = _mm_add_epi64(xmms2, xmms3); \
63  xmms4 = _mm_add_epi64(xmms4, xmms5); \
64  xmms6 = _mm_add_epi64(xmms6, xmms7); \
65  xmmt1 = _mm_slli_epi64(xmms1, L0); \
66  xmmt3 = _mm_slli_epi64(xmms3, L0); \
67  xmmt5 = _mm_slli_epi64(xmms5, L0); \
68  xmmt7 = _mm_slli_epi64(xmms7, L0); \
69  xmms1 = _mm_srli_epi64(xmms1, R0); \
70  xmms3 = _mm_srli_epi64(xmms3, R0); \
71  xmms5 = _mm_srli_epi64(xmms5, R0); \
72  xmms7 = _mm_srli_epi64(xmms7, R0); \
73  xmmt1 = _mm_or_si128(xmms1, xmmt1); \
74  xmmt3 = _mm_or_si128(xmms3, xmmt3); \
75  xmmt5 = _mm_or_si128(xmms5, xmmt5); \
76  xmmt7 = _mm_or_si128(xmms7, xmmt7); \
77  xmms1 = _mm_xor_si128(xmms0, xmmt1); \
78  xmms3 = _mm_xor_si128(xmms2, xmmt3); \
79  xmms5 = _mm_xor_si128(xmms4, xmmt5); \
80  xmms7 = _mm_xor_si128(xmms6, xmmt7); \
81  }
82 
83 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_CYCLE_4(N) \
84  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_RBOX(N * 8 + 1); \
85  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_RBOX(N * 8 + 2); \
86  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_RBOX(N * 8 + 3); \
87  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_RBOX(N * 8 + 4); \
88  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_KBOX(N * 8 + 4);
89 
90 #define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_CYCLE_8(N) \
91  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_RBOX(N * 8 + 1); \
92  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_RBOX(N * 8 + 2); \
93  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_RBOX(N * 8 + 3); \
94  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_RBOX(N * 8 + 4); \
95  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_KBOX(N * 8 + 4); \
96  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_RBOX(N * 8 + 5); \
97  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_RBOX(N * 8 + 6); \
98  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_RBOX(N * 8 + 7); \
99  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_RBOX(N * 8 + 8); \
100  MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_KBOX(N * 8 + 8);
101 
102 namespace mckl {
103 
104 namespace internal {
105 
106 template <typename T>
108 {
109  static_assert(std::numeric_limits<T>::digits == 64,
110  "**Threefry2x64GeneratorSSE2Impl** used with T other than a 64-bit "
111  "unsigned integers");
112 
113  static constexpr std::size_t K = 2;
114  static constexpr std::size_t Rounds = 20;
116 
117  public:
118  static void eval(
119  const void *plain, void *cipher, const std::array<T, K + 4> &par)
120  {
121  Threefry2x64GeneratorGenericImpl<T>::eval(plain, cipher, par);
122  }
123 
124  template <typename ResultType>
125  static void eval(std::array<std::uint64_t, 2> &ctr, ResultType *r,
126  const std::array<T, K + 4> &par)
127  {
129  }
130 
131  template <typename ResultType>
132  static void eval(std::array<std::uint64_t, 2> &ctr, std::size_t n,
133  ResultType *r, const std::array<T, K + 4> &par)
134  {
135  constexpr std::size_t R = sizeof(T) * K / sizeof(ResultType);
136 
137  const std::size_t n0 =
138  static_cast<std::size_t>(std::min(static_cast<std::uint64_t>(n),
139  std::numeric_limits<std::uint64_t>::max() - ctr.front()));
140 
141  eval_kernel(ctr, n0, r, par);
142  n -= n0;
143  r += n0 * R;
144 
145  if (n != 0) {
146  eval(ctr, r, par);
147  n -= 1;
148  r += R;
149  }
150 
151  eval_kernel(ctr, n, r, par);
152  }
153 
154  private:
155  template <typename ResultType>
156  static void eval_kernel(std::array<std::uint64_t, 2> &ctr, std::size_t n,
157  ResultType *r, const std::array<T, K + 4> &par)
158  {
159  constexpr std::size_t S = 8;
160  constexpr std::size_t N = sizeof(__m128i) * S / (sizeof(T) * K);
161 
162  __m128i xmmc =
163  _mm_set_epi64x(static_cast<MCKL_INT64>(std::get<1>(ctr)),
164  static_cast<MCKL_INT64>(std::get<0>(ctr)));
165  ctr.front() += n;
166 
167  __m128i *rptr = reinterpret_cast<__m128i *>(r);
168  while (n != 0) {
169  __m128i xmmt0 = _mm_add_epi64(xmmc, _mm_set_epi64x(0, 1));
170  __m128i xmmt1 = _mm_add_epi64(xmmc, _mm_set_epi64x(0, 2));
171  __m128i xmmt2 = _mm_add_epi64(xmmc, _mm_set_epi64x(0, 3));
172  __m128i xmmt3 = _mm_add_epi64(xmmc, _mm_set_epi64x(0, 4));
173  __m128i xmmt4 = _mm_add_epi64(xmmc, _mm_set_epi64x(0, 5));
174  __m128i xmmt5 = _mm_add_epi64(xmmc, _mm_set_epi64x(0, 6));
175  __m128i xmmt6 = _mm_add_epi64(xmmc, _mm_set_epi64x(0, 7));
176  __m128i xmmt7 = _mm_add_epi64(xmmc, _mm_set_epi64x(0, 8));
177  xmmc = xmmt7;
178 
179  __m128i xmms0 = _mm_unpacklo_epi64(xmmt0, xmmt1);
180  __m128i xmms1 = _mm_unpackhi_epi64(xmmt0, xmmt1);
181  __m128i xmms2 = _mm_unpacklo_epi64(xmmt2, xmmt3);
182  __m128i xmms3 = _mm_unpackhi_epi64(xmmt2, xmmt3);
183  __m128i xmms4 = _mm_unpacklo_epi64(xmmt4, xmmt5);
184  __m128i xmms5 = _mm_unpackhi_epi64(xmmt4, xmmt5);
185  __m128i xmms6 = _mm_unpacklo_epi64(xmmt6, xmmt7);
186  __m128i xmms7 = _mm_unpackhi_epi64(xmmt6, xmmt7);
187 
192 
193  xmmt0 = _mm_unpacklo_epi64(xmms0, xmms1);
194  xmmt1 = _mm_unpackhi_epi64(xmms0, xmms1);
195  xmmt2 = _mm_unpacklo_epi64(xmms2, xmms3);
196  xmmt3 = _mm_unpackhi_epi64(xmms2, xmms3);
197  xmmt4 = _mm_unpacklo_epi64(xmms4, xmms5);
198  xmmt5 = _mm_unpackhi_epi64(xmms4, xmms5);
199  xmmt6 = _mm_unpacklo_epi64(xmms6, xmms7);
200  xmmt7 = _mm_unpackhi_epi64(xmms6, xmms7);
201 
202  if (n >= N) {
203  n -= N;
204  _mm_storeu_si128(rptr++, xmmt0);
205  _mm_storeu_si128(rptr++, xmmt1);
206  _mm_storeu_si128(rptr++, xmmt2);
207  _mm_storeu_si128(rptr++, xmmt3);
208  _mm_storeu_si128(rptr++, xmmt4);
209  _mm_storeu_si128(rptr++, xmmt5);
210  _mm_storeu_si128(rptr++, xmmt6);
211  _mm_storeu_si128(rptr++, xmmt7);
212  } else {
213  std::array<__m128i, S> s;
214  std::get<0>(s) = xmmt0;
215  std::get<1>(s) = xmmt1;
216  std::get<2>(s) = xmmt2;
217  std::get<3>(s) = xmmt3;
218  std::get<4>(s) = xmmt4;
219  std::get<5>(s) = xmmt5;
220  std::get<6>(s) = xmmt6;
221  std::get<7>(s) = xmmt7;
222  std::memcpy(rptr, s.data(), n * sizeof(T) * K);
223  break;
224  }
225  }
226  }
227 }; // class Threefry2x64GeneratorSSE2Impl
228 
229 } // namespace internal
230 
231 } // namespace mckl
232 
234 
235 #endif // MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_HPP
#define MCKL_PUSH_GCC_WARNING(warning)
Definition: compiler.h:78
static void eval(std::array< std::uint64_t, 2 > &ctr, std::size_t n, ResultType *r, const std::array< T, K+4 > &par)
#define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_CYCLE_8(N)
static void eval(const void *plain, void *cipher, const std::array< T, K+4 > &par)
Default Threefry constants.
#define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_CYCLE_4(N)
Definition: mcmc.hpp:40
static void eval(std::array< std::uint64_t, 2 > &ctr, ResultType *r, const std::array< T, K+4 > &par)
#define MCKL_RANDOM_INTERNAL_THREEFRY_SSE2_2X64_KBOX(N)
static void eval(const void *plain, void *cipher, const std::array< T, K+4 > &par)
#define MCKL_POP_GCC_WARNING
Definition: compiler.h:79