From 8f25d96f558f979cea18b5d6ca6aa8db71edab05 Mon Sep 17 00:00:00 2001 From: wantonsushi Date: Sun, 21 Jun 2026 21:27:33 -0400 Subject: [PATCH] Speed up scalar Sobol sampling via Ahmed 2024 The scalar Sobol path evaluates each dimension with a per-bit matrix-vector product over 16 columns. Ahmed (EGSR 2024) computes the same result more efficiently using a sequence of shift-mask-xor steps, avoiding the loop over the columns. Replace the scalar sobolReversedIndex loop with closed-form steps for dimensions 0 to 3, steps emitted by a new matrices CLI tool, and cite the technique in the comments and docs. Leave the SIMD paths unchanged, as experiments showed Ahmed to be slower since the SIMD path also avoids the loop by packing the 16 columns into a wide register. Add test to validate the new path against old classic implementation. Additionally, as suggested by Chris Kulla, on CPU we reverse the shared input step masks and right shift on CPU, which removes the need for reversal in every draw. GPU keeps left shifts as experiments showed they are faster than right shifts. Signed-off-by: wantonsushi --- CHANGELOG.md | 3 + README.md | 7 ++- include/oqmc/owen.h | 112 +++++++++++++++++++++++++++++++++---- include/oqmc/sobol.h | 4 +- src/tests/owen.cpp | 69 +++++++++++++++++++++++ src/tools/cli/matrices.cpp | 97 ++++++++++++++++++++++++++++++++ 6 files changed, 279 insertions(+), 13 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f0923c6..e08f23b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added ### Changed + +- Sped up scalar paths of `oqmc::sobolReversedIndex` using the construction method from Ahmed 2024. + ### Deprecated ### Removed ### Fixed diff --git a/README.md b/README.md index 25e0d16..29f8de7 100644 --- a/README.md +++ b/README.md @@ -975,7 +975,8 @@ The implementation uses an elegant construction by Burley [^1] for an Owen scrambled Sobol sequence. This also includes performance improvements such as limiting the index to 16 bits, pre-inverting the input and output matrices, and making use of CPU vector intrinsics. You need to select a `OPENQMC_ARCH_TYPE` to -make use of the performance from vector intrinsics for a given architecture. +make use of the performance from vector intrinsics for a given architecture. The +scalar path uses the closed-form 2D Sobol construction by Ahmed [^5]. @@ -1470,6 +1471,10 @@ https://dx.doi.org/10.1007/s11075-011-9482-5 Lessons Learned and Improvements when Building Screen-Space Samplers with Blue-Noise Error Distribution. ACM SIGGRAPH 2021 Talks. https://dx.doi.org/10.1145/3450623.3464645 +[^5]: Abdalla G. M. Ahmed. 2024. +An Implementation Algorithm of 2D Sobol Sequence: Fast, Elegant, and Compact. EGSR. +https://dx.doi.org/10.2312/sr.20241147 + [^6]: Melissa E. O'Neill. 2014. PCG: A Family of Simple Fast Space-Efficient Statistically Good Algorithms for Random Number Generation. https://www.pcg-random.org diff --git a/include/oqmc/owen.h b/include/oqmc/owen.h index 507f0ed..64b58b8 100644 --- a/include/oqmc/owen.h +++ b/include/oqmc/owen.h @@ -48,6 +48,107 @@ OQMC_HOST_DEVICE inline std::uint16_t sobolReversedIndex(std::uint16_t index, assert(dimension >= 0); assert(dimension <= 3); +#if defined(OQMC_ARCH_SCALAR) + + // Each matrix factors into shift-mask-xor steps (Ahmed 2024, eq. 18), + // dimension 1 being the Pascal matrix (Listing 19). Steps emitted by the + // matrices cli tool in src/tools/cli/matrices.cpp. + +#if defined(__CUDA_ARCH__) + + // Prefer left shifts on GPU since they are more efficient than right + // shifts. Reversal must be applied to the result. + switch(dimension) + { + case 1: + index ^= static_cast((index & 0x00ff) << 8); + index ^= static_cast((index & 0x0f0f) << 4); + index ^= static_cast((index & 0x3333) << 2); + index ^= static_cast((index & 0x5555) << 1); + break; + case 2: + index ^= static_cast((index & 0x0003) << 14); + index ^= static_cast((index & 0x0004) << 13); + index ^= static_cast((index & 0x000f) << 12); + index ^= static_cast((index & 0x0033) << 10); + index ^= static_cast((index & 0x0055) << 9); + index ^= static_cast((index & 0x0030) << 8); + index ^= static_cast((index & 0x0303) << 6); + index ^= static_cast((index & 0x0505) << 5); + index ^= static_cast((index & 0x0cf3) << 4); + index ^= static_cast((index & 0x1111) << 3); + index ^= static_cast((index & 0x0f0f) << 2); + index ^= static_cast((index & 0x6666) << 1); + break; + case 3: + index ^= static_cast((index & 0x000f) << 12); + index ^= static_cast((index & 0x0030) << 10); + index ^= static_cast((index & 0x0050) << 9); + index ^= static_cast((index & 0x007f) << 8); + index ^= static_cast((index & 0x03f0) << 6); + index ^= static_cast((index & 0x0410) << 5); + index ^= static_cast((index & 0x07e0) << 4); + index ^= static_cast((index & 0x1c71) << 3); + index ^= static_cast((index & 0x1c71) << 2); + index ^= static_cast((index & 0x4924) << 1); + break; + default: + break; + } + + return reverseBits16(index); + +#else + + // Reversed masks and right shifts put the reversal on the shared input. + // This optimization hoists the reversal out of a draw + // (Chris Kulla, PR #97). + index = reverseBits16(index); + + switch(dimension) + { + case 1: + index ^= static_cast((index & 0xff00) >> 8); + index ^= static_cast((index & 0xf0f0) >> 4); + index ^= static_cast((index & 0xcccc) >> 2); + index ^= static_cast((index & 0xaaaa) >> 1); + break; + case 2: + index ^= static_cast((index & 0xc000) >> 14); + index ^= static_cast((index & 0x2000) >> 13); + index ^= static_cast((index & 0xf000) >> 12); + index ^= static_cast((index & 0xcc00) >> 10); + index ^= static_cast((index & 0xaa00) >> 9); + index ^= static_cast((index & 0x0c00) >> 8); + index ^= static_cast((index & 0xc0c0) >> 6); + index ^= static_cast((index & 0xa0a0) >> 5); + index ^= static_cast((index & 0xcf30) >> 4); + index ^= static_cast((index & 0x8888) >> 3); + index ^= static_cast((index & 0xf0f0) >> 2); + index ^= static_cast((index & 0x6666) >> 1); + break; + case 3: + index ^= static_cast((index & 0xf000) >> 12); + index ^= static_cast((index & 0x0c00) >> 10); + index ^= static_cast((index & 0x0a00) >> 9); + index ^= static_cast((index & 0xfe00) >> 8); + index ^= static_cast((index & 0x0fc0) >> 6); + index ^= static_cast((index & 0x0820) >> 5); + index ^= static_cast((index & 0x07e0) >> 4); + index ^= static_cast((index & 0x8e38) >> 3); + index ^= static_cast((index & 0x8e38) >> 2); + index ^= static_cast((index & 0x2492) >> 1); + break; + default: + break; + } + + return index; + +#endif + +#else + if(dimension == 0) { return reverseBits16(index); @@ -242,17 +343,6 @@ OQMC_HOST_DEVICE inline std::uint16_t sobolReversedIndex(std::uint16_t index, return vgetq_lane_u16(bits, 0); #endif -#if defined(OQMC_ARCH_SCALAR) - std::uint16_t sample = 0; - for(int i = 0; i < 16; ++i) - { - if((index & masks[i]) != 0) - { - sample ^= matrix[i]; - } - } - - return sample; #endif } diff --git a/include/oqmc/sobol.h b/include/oqmc/sobol.h index 08dce09..034ec82 100644 --- a/include/oqmc/sobol.h +++ b/include/oqmc/sobol.h @@ -99,7 +99,9 @@ void SobolImpl::drawRnd(std::uint32_t rnd[Size]) const /// includes performance improvements such as limiting the index to 16 bits, /// pre-inverting the input and output matrices, and making use of CPU vector /// intrinsics. You need to select an `OPENQMC_ARCH_TYPE` to make use of the -/// performance from vector intrinsics for a given architecture. +/// performance from vector intrinsics for a given architecture. The scalar path +/// uses the construction method by Abdalla G. M. Ahmed in 'An Implementation +/// Algorithm of 2D Sobol Sequence Fast, Elegant, and Compact' (EGSR 2024). /// /// This sampler has no cache initialisation cost, it generates all samples on /// the fly without touching memory. However the cost per draw sample call is diff --git a/src/tests/owen.cpp b/src/tests/owen.cpp index 4a5b377..a788fe9 100644 --- a/src/tests/owen.cpp +++ b/src/tests/owen.cpp @@ -4,6 +4,7 @@ #include "hypothesis.h" #include #include +#include #include @@ -142,4 +143,72 @@ TEST(OwenTest, ShirleyRemapping) } } +TEST(OwenTest, SobolReversedIndex) +{ + // clang-format off + constexpr std::uint16_t masks[16] = { + 0b0000000000000001, + 0b0000000000000010, + 0b0000000000000100, + 0b0000000000001000, + 0b0000000000010000, + 0b0000000000100000, + 0b0000000001000000, + 0b0000000010000000, + 0b0000000100000000, + 0b0000001000000000, + 0b0000010000000000, + 0b0000100000000000, + 0b0001000000000000, + 0b0010000000000000, + 0b0100000000000000, + 0b1000000000000000, + }; + + constexpr std::uint16_t directions[4][16] = { + {0x8000, 0x4000, 0x2000, 0x1000, 0x0800, 0x0400, 0x0200, 0x0100, + 0x0080, 0x0040, 0x0020, 0x0010, 0x0008, 0x0004, 0x0002, 0x0001}, + {0xffff, 0x5555, 0x3333, 0x1111, 0x0f0f, 0x0505, 0x0303, 0x0101, + 0x00ff, 0x0055, 0x0033, 0x0011, 0x000f, 0x0005, 0x0003, 0x0001}, + {0xaa09, 0x7706, 0x3903, 0x1601, 0x09aa, 0x0677, 0x0339, 0x0116, + 0x00a3, 0x0071, 0x003a, 0x0017, 0x0009, 0x0006, 0x0003, 0x0001}, + {0xa0c3, 0x4041, 0x302d, 0x101e, 0x0b67, 0x079a, 0x02a4, 0x011b, + 0x00c9, 0x0045, 0x002e, 0x001f, 0x000a, 0x0004, 0x0003, 0x0001}, + }; + // clang-format on + + // Reference code from the classic scalar implementation. Verifies that + // Ahmed 2024 closed form output matches (PR #97). + const auto reference = [&](std::uint16_t index, int dimension) { + if(dimension == 0) + { + return oqmc::reverseBits16(index); + } + + const auto matrix = directions[dimension]; + + std::uint16_t sample = 0; + for(int i = 0; i < 16; ++i) + { + if((index & masks[i]) != 0) + { + sample ^= matrix[i]; + } + } + + return sample; + }; + + for(int dimension = 0; dimension < 4; ++dimension) + { + for(std::uint32_t index = 0; index < (1u << 16); ++index) + { + const auto value = static_cast(index); + + EXPECT_EQ(oqmc::sobolReversedIndex(value, dimension), + reference(value, dimension)); + } + } +} + } // namespace diff --git a/src/tools/cli/matrices.cpp b/src/tools/cli/matrices.cpp index 38c9f07..30d94ff 100644 --- a/src/tools/cli/matrices.cpp +++ b/src/tools/cli/matrices.cpp @@ -32,10 +32,13 @@ #include +#include #include #include #include #include +#include +#include constexpr auto numDimensions = 1024; constexpr auto size = 52; @@ -13379,6 +13382,98 @@ void printMatrices(int dimensionSize, int indexSize) } } +// In the bit reversed basis, the generator matrix is unit lower-triangular, so +// GF(2) elimination factors it into shift-mask-xor steps (Ahmed 2024, eq. 18), +// one per sub-diagonal. Reversed, the steps form the program for that +// dimension. Programs print twice: left shifts for the device, and the +// conjugate form with reversed masks and right shifts for the host. + +constexpr int wordBits = 16; + +std::uint16_t reversedRow(int dimension, int row) +{ + assert(dimension >= 0 && dimension < numDimensions); + assert(row >= 0 && row < wordBits); + + std::uint16_t bits = 0; + for(int column = 0; column < wordBits; ++column) + { + const auto reversedColumn = wordBits - 1 - column; + const auto reversedRowIndex = wordBits - 1 - row; + + const auto value = + oqmc::reverseBits32(matrices[dimension][reversedColumn]); + const auto reversed = std::bitset(value); + if(reversed[reversedRowIndex]) + { + bits |= static_cast(1u << column); + } + } + return bits; +} + +void printPrograms(int dimensionSize, bool conjugate) +{ + assert(dimensionSize >= 0); + + for(int dimension = 0; dimension < dimensionSize; ++dimension) + { + std::array rows; + for(int row = 0; row < wordBits; ++row) + { + rows[row] = reversedRow(dimension, row); + } + + std::vector> program; + for(int shift = 1; shift < wordBits; ++shift) + { + std::uint16_t mask = 0; + for(int row = 0; row + shift < wordBits; ++row) + { + if(((rows[row + shift] >> row) & 1) != 0) + { + mask |= static_cast(1u << row); + } + } + if(mask == 0) + { + continue; + } + for(int row = 0; row + shift < wordBits; ++row) + { + if(((mask >> row) & 1) != 0) + { + rows[row + shift] ^= rows[row]; + } + } + program.emplace_back(mask, shift); + } + + for(int row = 0; row < wordBits; ++row) + { + assert(rows[row] == static_cast(1u << row)); + } + + std::printf("// dimension %d\n", dimension); + for(auto step = program.rbegin(); step != program.rend(); ++step) + { + if(conjugate) + { + std::printf("index ^= static_cast((index & " + "0x%04x) >> %d);\n", + oqmc::reverseBits16(step->first), step->second); + } + else + { + std::printf("index ^= static_cast((index & " + "0x%04x) << %d);\n", + step->first, step->second); + } + } + std::printf("\n"); + } +} + int main() { constexpr auto samplePrecision = 16; @@ -13386,6 +13481,8 @@ int main() constexpr auto indexSize = 16; printMatrices(dimensionSize, indexSize); + printPrograms(dimensionSize, false); + printPrograms(dimensionSize, true); return 0; }