OpenQMC API
Loading...
Searching...
No Matches
owen.h
Go to the documentation of this file.
1// SPDX-License-Identifier: Apache-2.0
2// Copyright Contributors to the OpenQMC Project.
3
9
10#pragma once
11
12#include "arch.h"
13#include "gpu.h"
14#include "permute.h"
15#include "reverse.h"
16#include "rotate.h"
17
18#include <cassert>
19#include <cstdint>
20
21#if defined(OQMC_ARCH_AVX)
22#include <immintrin.h>
23#endif
24
25#if defined(OQMC_ARCH_SSE)
26#include <emmintrin.h>
27#endif
28
29#if defined(OQMC_ARCH_ARM)
30#include <arm_neon.h>
31#endif
32
33namespace oqmc
34{
35
45OQMC_HOST_DEVICE inline std::uint16_t sobolReversedIndex(std::uint16_t index,
46 int dimension)
47{
48 assert(dimension >= 0);
49 assert(dimension <= 3);
50
51#if defined(OQMC_ARCH_SCALAR)
52
53 // Each matrix factors into shift-mask-xor steps (Ahmed 2024, eq. 18),
54 // dimension 1 being the Pascal matrix (Listing 19). Steps emitted by the
55 // matrices cli tool in src/tools/cli/matrices.cpp.
56
57#if defined(__CUDA_ARCH__)
58
59 // Prefer left shifts on GPU since they are more efficient than right
60 // shifts. Reversal must be applied to the result.
61 switch(dimension)
62 {
63 case 1:
64 index ^= static_cast<std::uint16_t>((index & 0x00ff) << 8);
65 index ^= static_cast<std::uint16_t>((index & 0x0f0f) << 4);
66 index ^= static_cast<std::uint16_t>((index & 0x3333) << 2);
67 index ^= static_cast<std::uint16_t>((index & 0x5555) << 1);
68 break;
69 case 2:
70 index ^= static_cast<std::uint16_t>((index & 0x0003) << 14);
71 index ^= static_cast<std::uint16_t>((index & 0x0004) << 13);
72 index ^= static_cast<std::uint16_t>((index & 0x000f) << 12);
73 index ^= static_cast<std::uint16_t>((index & 0x0033) << 10);
74 index ^= static_cast<std::uint16_t>((index & 0x0055) << 9);
75 index ^= static_cast<std::uint16_t>((index & 0x0030) << 8);
76 index ^= static_cast<std::uint16_t>((index & 0x0303) << 6);
77 index ^= static_cast<std::uint16_t>((index & 0x0505) << 5);
78 index ^= static_cast<std::uint16_t>((index & 0x0cf3) << 4);
79 index ^= static_cast<std::uint16_t>((index & 0x1111) << 3);
80 index ^= static_cast<std::uint16_t>((index & 0x0f0f) << 2);
81 index ^= static_cast<std::uint16_t>((index & 0x6666) << 1);
82 break;
83 case 3:
84 index ^= static_cast<std::uint16_t>((index & 0x000f) << 12);
85 index ^= static_cast<std::uint16_t>((index & 0x0030) << 10);
86 index ^= static_cast<std::uint16_t>((index & 0x0050) << 9);
87 index ^= static_cast<std::uint16_t>((index & 0x007f) << 8);
88 index ^= static_cast<std::uint16_t>((index & 0x03f0) << 6);
89 index ^= static_cast<std::uint16_t>((index & 0x0410) << 5);
90 index ^= static_cast<std::uint16_t>((index & 0x07e0) << 4);
91 index ^= static_cast<std::uint16_t>((index & 0x1c71) << 3);
92 index ^= static_cast<std::uint16_t>((index & 0x1c71) << 2);
93 index ^= static_cast<std::uint16_t>((index & 0x4924) << 1);
94 break;
95 default:
96 break;
97 }
98
99 return reverseBits16(index);
100
101#else
102
103 // Reversed masks and right shifts put the reversal on the shared input.
104 // This optimization hoists the reversal out of a draw
105 // (Chris Kulla, PR #97).
106 index = reverseBits16(index);
107
108 switch(dimension)
109 {
110 case 1:
111 index ^= static_cast<std::uint16_t>((index & 0xff00) >> 8);
112 index ^= static_cast<std::uint16_t>((index & 0xf0f0) >> 4);
113 index ^= static_cast<std::uint16_t>((index & 0xcccc) >> 2);
114 index ^= static_cast<std::uint16_t>((index & 0xaaaa) >> 1);
115 break;
116 case 2:
117 index ^= static_cast<std::uint16_t>((index & 0xc000) >> 14);
118 index ^= static_cast<std::uint16_t>((index & 0x2000) >> 13);
119 index ^= static_cast<std::uint16_t>((index & 0xf000) >> 12);
120 index ^= static_cast<std::uint16_t>((index & 0xcc00) >> 10);
121 index ^= static_cast<std::uint16_t>((index & 0xaa00) >> 9);
122 index ^= static_cast<std::uint16_t>((index & 0x0c00) >> 8);
123 index ^= static_cast<std::uint16_t>((index & 0xc0c0) >> 6);
124 index ^= static_cast<std::uint16_t>((index & 0xa0a0) >> 5);
125 index ^= static_cast<std::uint16_t>((index & 0xcf30) >> 4);
126 index ^= static_cast<std::uint16_t>((index & 0x8888) >> 3);
127 index ^= static_cast<std::uint16_t>((index & 0xf0f0) >> 2);
128 index ^= static_cast<std::uint16_t>((index & 0x6666) >> 1);
129 break;
130 case 3:
131 index ^= static_cast<std::uint16_t>((index & 0xf000) >> 12);
132 index ^= static_cast<std::uint16_t>((index & 0x0c00) >> 10);
133 index ^= static_cast<std::uint16_t>((index & 0x0a00) >> 9);
134 index ^= static_cast<std::uint16_t>((index & 0xfe00) >> 8);
135 index ^= static_cast<std::uint16_t>((index & 0x0fc0) >> 6);
136 index ^= static_cast<std::uint16_t>((index & 0x0820) >> 5);
137 index ^= static_cast<std::uint16_t>((index & 0x07e0) >> 4);
138 index ^= static_cast<std::uint16_t>((index & 0x8e38) >> 3);
139 index ^= static_cast<std::uint16_t>((index & 0x8e38) >> 2);
140 index ^= static_cast<std::uint16_t>((index & 0x2492) >> 1);
141 break;
142 default:
143 break;
144 }
145
146 return index;
147
148#endif
149
150#else
151
152 if(dimension == 0)
153 {
154 return reverseBits16(index);
155 }
156
157 // Following matrices were produced using the matrices cli tool found in the
158 // source file src/tools/cli/matrices.cpp. This in turn uses matrices that
159 // were copied from MIT licensed code written by Leonhard Gruenschloss.
160
161 // clang-format off
162 constexpr std::uint16_t masks[16] = {
163 0b0000000000000001,
164 0b0000000000000010,
165 0b0000000000000100,
166 0b0000000000001000,
167 0b0000000000010000,
168 0b0000000000100000,
169 0b0000000001000000,
170 0b0000000010000000,
171 0b0000000100000000,
172 0b0000001000000000,
173 0b0000010000000000,
174 0b0000100000000000,
175 0b0001000000000000,
176 0b0010000000000000,
177 0b0100000000000000,
178 0b1000000000000000,
179 };
180
181 constexpr std::uint16_t directions[4][16] = {
182 {
183 0b1000000000000000,
184 0b0100000000000000,
185 0b0010000000000000,
186 0b0001000000000000,
187 0b0000100000000000,
188 0b0000010000000000,
189 0b0000001000000000,
190 0b0000000100000000,
191 0b0000000010000000,
192 0b0000000001000000,
193 0b0000000000100000,
194 0b0000000000010000,
195 0b0000000000001000,
196 0b0000000000000100,
197 0b0000000000000010,
198 0b0000000000000001,
199 },
200
201 {
202 0b1111111111111111,
203 0b0101010101010101,
204 0b0011001100110011,
205 0b0001000100010001,
206 0b0000111100001111,
207 0b0000010100000101,
208 0b0000001100000011,
209 0b0000000100000001,
210 0b0000000011111111,
211 0b0000000001010101,
212 0b0000000000110011,
213 0b0000000000010001,
214 0b0000000000001111,
215 0b0000000000000101,
216 0b0000000000000011,
217 0b0000000000000001,
218 },
219
220 {
221 0b1010101000001001,
222 0b0111011100000110,
223 0b0011100100000011,
224 0b0001011000000001,
225 0b0000100110101010,
226 0b0000011001110111,
227 0b0000001100111001,
228 0b0000000100010110,
229 0b0000000010100011,
230 0b0000000001110001,
231 0b0000000000111010,
232 0b0000000000010111,
233 0b0000000000001001,
234 0b0000000000000110,
235 0b0000000000000011,
236 0b0000000000000001,
237 },
238
239 {
240 0b1010000011000011,
241 0b0100000001000001,
242 0b0011000000101101,
243 0b0001000000011110,
244 0b0000101101100111,
245 0b0000011110011010,
246 0b0000001010100100,
247 0b0000000100011011,
248 0b0000000011001001,
249 0b0000000001000101,
250 0b0000000000101110,
251 0b0000000000011111,
252 0b0000000000001010,
253 0b0000000000000100,
254 0b0000000000000011,
255 0b0000000000000001,
256 },
257 };
258 // clang-format on
259
260 const auto matrix = directions[dimension];
261
262#if defined(OQMC_ARCH_AVX)
263 constexpr auto stride = 16;
264 const __m256i zero = _mm256_setzero_si256();
265
266 __m256i bits = zero;
267 for(int i = 0; i < 16; i += stride)
268 {
269 const auto maskPtr = reinterpret_cast<const __m256i*>(masks + i);
270 const __m256i mask = _mm256_loadu_si256(maskPtr);
271
272 const auto matrixPtr = reinterpret_cast<const __m256i*>(matrix + i);
273 const __m256i column = _mm256_loadu_si256(matrixPtr);
274
275 const __m256i masked = _mm256_and_si256(_mm256_set1_epi16(index), mask);
276 const __m256i cond = _mm256_cmpeq_epi16(masked, zero);
277
278 const __m256i xored = _mm256_xor_si256(bits, column);
279
280 bits = _mm256_or_si256(_mm256_and_si256(cond, bits),
281 _mm256_andnot_si256(cond, xored));
282 }
283
284 bits = _mm256_xor_si256(bits, _mm256_srli_si256(bits, 2));
285 bits = _mm256_xor_si256(bits, _mm256_srli_si256(bits, 4));
286 bits = _mm256_xor_si256(bits, _mm256_srli_si256(bits, 8));
287
288 return _mm256_extract_epi16(bits, 0) ^ _mm256_extract_epi16(bits, 8);
289#endif
290
291#if defined(OQMC_ARCH_SSE)
292 constexpr auto stride = 8;
293 const __m128i zero = _mm_setzero_si128();
294
295 __m128i bits = zero;
296 for(int i = 0; i < 16; i += stride)
297 {
298 const auto maskPtr = reinterpret_cast<const __m128i*>(masks + i);
299 const __m128i mask = _mm_loadu_si128(maskPtr);
300
301 const auto matrixPtr = reinterpret_cast<const __m128i*>(matrix + i);
302 const __m128i column = _mm_loadu_si128(matrixPtr);
303
304 const __m128i masked = _mm_and_si128(_mm_set1_epi16(index), mask);
305 const __m128i cond = _mm_cmpeq_epi16(masked, zero);
306
307 const __m128i xored = _mm_xor_si128(bits, column);
308
309 bits = _mm_or_si128(_mm_and_si128(cond, bits),
310 _mm_andnot_si128(cond, xored));
311 }
312
313 bits = _mm_xor_si128(bits, _mm_srli_si128(bits, 2));
314 bits = _mm_xor_si128(bits, _mm_srli_si128(bits, 4));
315 bits = _mm_xor_si128(bits, _mm_srli_si128(bits, 8));
316
317 return _mm_extract_epi16(bits, 0);
318#endif
319
320#if defined(OQMC_ARCH_ARM)
321 constexpr auto stride = 8;
322 const uint16x8_t zero = vdupq_n_u16(0);
323
324 uint16x8_t bits = zero;
325 for(int i = 0; i < 16; i += stride)
326 {
327 const uint16x8_t mask = vld1q_u16(masks + i);
328 const uint16x8_t column = vld1q_u16(matrix + i);
329
330 const uint16x8_t masked = vandq_u16(vdupq_n_u16(index), mask);
331 const uint16x8_t cond = vceqq_u16(masked, zero);
332
333 const uint16x8_t xored = veorq_u16(bits, column);
334
335 bits =
336 vorrq_u16(vandq_u16(cond, bits), vandq_u16(vmvnq_u16(cond), xored));
337 }
338
339 bits = veorq_u16(bits, vextq_u16(bits, zero, 1));
340 bits = veorq_u16(bits, vextq_u16(bits, zero, 2));
341 bits = veorq_u16(bits, vextq_u16(bits, zero, 4));
342
343 return vgetq_lane_u16(bits, 0);
344#endif
345
346#endif
347}
348
359OQMC_HOST_DEVICE constexpr std::uint32_t scrambleAndReverse(std::uint32_t value,
360 std::uint32_t seed)
361{
362 value = laineKarrasPermutation(value, seed);
363 value = reverseBits32(value);
364
365 return value;
366}
367
379template <int Depth>
380OQMC_HOST_DEVICE inline void shuffledScrambledSobol(std::uint32_t index,
381 std::uint32_t seed,
382 std::uint32_t sample[Depth])
383{
384 static_assert(Depth >= 1, "Pattern depth is greater or equal to one.");
385 static_assert(Depth <= 4, "Pattern depth is less or equal to four.");
386
387 index = reverseAndShuffle(index, seed);
388
389 for(int i = 0; i < Depth; ++i)
390 {
391 sample[i] = sobolReversedIndex(index >> 16, i);
392 sample[i] = scrambleAndReverse(sample[i], rotateBytes(seed, i));
393 }
394}
395
396} // namespace oqmc
OQMC_HOST_DEVICE
#define OQMC_HOST_DEVICE
Definition gpu.h:13
oqmc::reverseAndShuffle
constexpr std::uint32_t reverseAndShuffle(std::uint32_t value, std::uint32_t seed)
Reverse input bits and shuffle order.
Definition permute.h:54
oqmc::rotateBytes
constexpr std::uint32_t rotateBytes(std::uint32_t value, int distance)
Rotate bytes in an integer value.
Definition rotate.h:41
oqmc::decodeBits16
EncodeKey decodeBits16(std::uint16_t value)
Decode a value back into a key.
Definition encode.h:81
oqmc::reverseBits16
constexpr std::uint16_t reverseBits16(std::uint16_t value)
Reverse bits of an unsigned 16 bit integer.
Definition reverse.h:53
oqmc::laineKarrasPermutation
constexpr std::uint32_t laineKarrasPermutation(std::uint32_t value, std::uint32_t seed)
Laine and Karras style permutation.
Definition permute.h:34
oqmc::reverseBits32
constexpr std::uint32_t reverseBits32(std::uint32_t value)
Reverse bits of an unsigned 32 bit integer.
Definition reverse.h:25
oqmc
Definition bntables.h:21
oqmc::shuffledScrambledSobol
void shuffledScrambledSobol(std::uint32_t index, std::uint32_t seed, std::uint32_t sample[Depth])
Compute a randomised sobol sequence value.
Definition owen.h:380
oqmc::scrambleAndReverse
constexpr std::uint32_t scrambleAndReverse(std::uint32_t value, std::uint32_t seed)
Permute an input integer and reverse the bits.
Definition owen.h:359
oqmc::sobolReversedIndex
std::uint16_t sobolReversedIndex(std::uint16_t index, int dimension)
Compute sobol sequence value at an index with reversed bits.
Definition owen.h:45