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 "gpu.h"
13#include "permute.h"
14#include "reverse.h"
15#include "rotate.h"
16
17#include <cassert>
18#include <cstdint>
19
20namespace oqmc
21{
22
34OQMC_HOST_DEVICE inline std::uint16_t szReversedBasis(std::uint16_t index,
35 int dimension)
36{
37 assert(dimension >= 0);
38 assert(dimension <= 3);
39
40 // Each matrix factors into shift-mask-xor steps (Ahmed 2024, eq. 18).
41 // Steps emitted by the matrices cli tool in src/tools/cli/matrices.cpp.
42 switch(dimension)
43 {
44 case 1:
45 index ^= static_cast<std::uint16_t>((index & 0x00ff) << 8);
46 index ^= static_cast<std::uint16_t>((index & 0x0f0f) << 4);
47 index ^= static_cast<std::uint16_t>((index & 0x3333) << 2);
48 index ^= static_cast<std::uint16_t>((index & 0x5555) << 1);
49 break;
50 case 2:
51 index ^= static_cast<std::uint16_t>((index & 0x0055) << 9);
52 index ^= static_cast<std::uint16_t>((index & 0x0055) << 8);
53 index ^= static_cast<std::uint16_t>((index & 0x00aa) << 7);
54 index ^= static_cast<std::uint16_t>((index & 0x0202) << 6);
55 index ^= static_cast<std::uint16_t>((index & 0x0505) << 5);
56 index ^= static_cast<std::uint16_t>((index & 0x0a0a) << 4);
57 index ^= static_cast<std::uint16_t>((index & 0x1b1b) << 3);
58 index ^= static_cast<std::uint16_t>((index & 0x3333) << 2);
59 index ^= static_cast<std::uint16_t>((index & 0x7777) << 1);
60 break;
61 case 3:
62 index ^= static_cast<std::uint16_t>((index & 0x0055) << 9);
63 index ^= static_cast<std::uint16_t>((index & 0x00aa) << 8);
64 index ^= static_cast<std::uint16_t>((index & 0x00aa) << 7);
65 index ^= static_cast<std::uint16_t>((index & 0x0202) << 6);
66 index ^= static_cast<std::uint16_t>((index & 0x0505) << 5);
67 index ^= static_cast<std::uint16_t>((index & 0x0505) << 4);
68 index ^= static_cast<std::uint16_t>((index & 0x1b1b) << 3);
69 index ^= static_cast<std::uint16_t>((index & 0x2222) << 2);
70 index ^= static_cast<std::uint16_t>((index & 0x2222) << 1);
71 break;
72 default:
73 break;
74 }
75
76 return reverseBits16(index);
77}
78
88OQMC_HOST_DEVICE inline std::uint16_t szToReversed(std::uint16_t index,
89 int dimension)
90{
91 assert(dimension >= 0);
92 assert(dimension <= 3);
93
94 // The szReversedBasis programs, conjugated with reversed masks and right
95 // shifts to consume the index directly.
96 switch(dimension)
97 {
98 case 1:
99 index ^= static_cast<std::uint16_t>((index & 0xff00) >> 8);
100 index ^= static_cast<std::uint16_t>((index & 0xf0f0) >> 4);
101 index ^= static_cast<std::uint16_t>((index & 0xcccc) >> 2);
102 index ^= static_cast<std::uint16_t>((index & 0xaaaa) >> 1);
103 break;
104 case 2:
105 index ^= static_cast<std::uint16_t>((index & 0xaa00) >> 9);
106 index ^= static_cast<std::uint16_t>((index & 0xaa00) >> 8);
107 index ^= static_cast<std::uint16_t>((index & 0x5500) >> 7);
108 index ^= static_cast<std::uint16_t>((index & 0x4040) >> 6);
109 index ^= static_cast<std::uint16_t>((index & 0xa0a0) >> 5);
110 index ^= static_cast<std::uint16_t>((index & 0x5050) >> 4);
111 index ^= static_cast<std::uint16_t>((index & 0xd8d8) >> 3);
112 index ^= static_cast<std::uint16_t>((index & 0xcccc) >> 2);
113 index ^= static_cast<std::uint16_t>((index & 0xeeee) >> 1);
114 break;
115 case 3:
116 index ^= static_cast<std::uint16_t>((index & 0xaa00) >> 9);
117 index ^= static_cast<std::uint16_t>((index & 0x5500) >> 8);
118 index ^= static_cast<std::uint16_t>((index & 0x5500) >> 7);
119 index ^= static_cast<std::uint16_t>((index & 0x4040) >> 6);
120 index ^= static_cast<std::uint16_t>((index & 0xa0a0) >> 5);
121 index ^= static_cast<std::uint16_t>((index & 0xa0a0) >> 4);
122 index ^= static_cast<std::uint16_t>((index & 0xd8d8) >> 3);
123 index ^= static_cast<std::uint16_t>((index & 0x4444) >> 2);
124 index ^= static_cast<std::uint16_t>((index & 0x4444) >> 1);
125 break;
126 default:
127 break;
128 }
129
130 return index;
131}
132
143OQMC_HOST_DEVICE constexpr std::uint32_t scrambleAndReverse(std::uint32_t value,
144 std::uint32_t seed)
145{
146 value = laineKarrasPermutation(value, seed);
147 value = reverseBits32(value);
148
149 return value;
150}
151
163template <int Depth>
164OQMC_HOST_DEVICE inline void shuffledScrambledSobol(std::uint32_t index,
165 std::uint32_t seed,
166 std::uint32_t sample[Depth])
167{
168 static_assert(Depth >= 1, "Pattern depth is greater or equal to one.");
169 static_assert(Depth <= 4, "Pattern depth is less or equal to four.");
170
171 index = reverseAndShuffle(index, seed);
172
173#if defined(__CUDA_ARCH__)
174
175 for(int i = 0; i < Depth; ++i)
176 {
177 sample[i] = szReversedBasis(index >> 16, i);
178 sample[i] = scrambleAndReverse(sample[i], rotateBytes(seed, i));
179 }
180
181#else
182
183 // Reverse the shared index once for all dimensions (Kulla, PR #97), and
184 // evaluate dimension 3 from dimension 2 (Ahmed et al. 2025, eq. 13).
185 const auto natural = reverseBits16(index >> 16);
186
187 std::uint16_t value[4];
188 for(int i = 0; i < Depth; ++i)
189 {
190 value[i] =
191 i == 3 ? szToReversed(value[2], 1) : szToReversed(natural, i);
192 sample[i] = scrambleAndReverse(value[i], rotateBytes(seed, i));
193 }
194
195#endif
196}
197
198} // namespace oqmc
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::szToReversed
std::uint16_t szToReversed(std::uint16_t index, int dimension)
Compute SZ sequence value from a natural order index.
Definition owen.h:88
oqmc::szReversedBasis
std::uint16_t szReversedBasis(std::uint16_t index, int dimension)
Compute SZ sequence value in the reversed basis.
Definition owen.h:34
oqmc::shuffledScrambledSobol
void shuffledScrambledSobol(std::uint32_t index, std::uint32_t seed, std::uint32_t sample[Depth])
Compute a randomised SZ sequence value.
Definition owen.h:164
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:143