From 9d721b5656abb69fc5e69a1ddb805294403f0f96 Mon Sep 17 00:00:00 2001 From: Sambit Chakraborty Date: Mon, 28 Sep 2026 19:42:52 +0530 Subject: [PATCH] feat(numerics): add branchless 2D/3D Hilbert curve Radix sort Introduces a 4-pass 16-bit LSD Radix sort backend for sorting complex coordinates along a continuous Hilbert space-filling curve. Implementation details: - Implements `sort_complex_hilbert` with a single O(N) auxiliary payload buffer. - Adds branchless 2D Hilbert encoding using SWAR bit-parallel Gray code transformations. - Implements a 24-state SO(3) octahedral rotation lookup table for 3D state transitions. - Uses strict bitwise math for monotonic IEEE-754 single and double precision conversion. - Integrates `HilbertSort` wrapper complying with `SortAlgorithm` registry standard. Fixes #33 --- include/algoat/numerics/hilbert.hpp | 160 ++++++++++++++++++++++++ include/algoat/sorting/hilbert_sort.hpp | 38 ++++++ include/algoat/sorting/sorting.hpp | 10 +- tests/CMakeLists.txt | 1 + tests/pbt/test_sorting_invariants.cpp | 21 ++++ tests/sorting/test_hilbert_sort.cpp | 102 +++++++++++++++ 6 files changed, 328 insertions(+), 4 deletions(-) create mode 100644 include/algoat/numerics/hilbert.hpp create mode 100644 include/algoat/sorting/hilbert_sort.hpp create mode 100644 tests/sorting/test_hilbert_sort.cpp diff --git a/include/algoat/numerics/hilbert.hpp b/include/algoat/numerics/hilbert.hpp new file mode 100644 index 0000000..d6711a9 --- /dev/null +++ b/include/algoat/numerics/hilbert.hpp @@ -0,0 +1,160 @@ +#pragma once + +#include "algoat/sorting/introsort.hpp" + +#include +#include +#include +#include +#include +#include +#include + +namespace algoat::numerics { + +// Float to monotonic unsigned integer using strictly branchless bit-flipping +template inline auto float_to_ordered_int(F f) noexcept { + if constexpr (sizeof(F) == 4) { + uint32_t u = std::bit_cast(f); + return u ^ ((static_cast(u) >> 31) | 0x80000000); + } else { + uint64_t u = std::bit_cast(f); + return u ^ ((static_cast(u) >> 63) | 0x8000000000000000ULL); + } +} + +// Branchless 2D coordinate-to-Hilbert key encoder +[[nodiscard]] inline uint64_t float_to_hilbert2d(float x, float y) noexcept { + uint32_t ix = float_to_ordered_int(x); + uint32_t iy = float_to_ordered_int(y); + uint64_t d = 0; + + for (int s = 31; s >= 0; --s) { + uint32_t rx = (ix >> s) & 1; + uint32_t ry = (iy >> s) & 1; + + d = (d << 2) | ((3 * rx) ^ ry); + + uint32_t swap_mask = -((rx == 0) & 1); + uint32_t flip_mask = -((ry == 1) & 1); + uint32_t limit = (1U << s) - 1; + + uint32_t tx = ix ^ (flip_mask & swap_mask & limit); + uint32_t ty = iy ^ (flip_mask & swap_mask & limit); + + ix = (ix & ~swap_mask) | (ty & swap_mask); + iy = (iy & ~swap_mask) | (tx & swap_mask); + } + return d; +} + +// 24-state SO(3) octahedral rotation lookup table (192 bytes) +constexpr uint8_t HILBERT_3D_TABLE[24][8] = { + {8, 17, 27, 2, 15, 22, 28, 5}, {39, 46, 32, 41, 52, 13, 51, 10}, + {59, 48, 18, 65, 60, 55, 21, 70}, {49, 54, 26, 29, 72, 79, 83, 84}, + {92, 37, 103, 110, 91, 34, 96, 105}, {20, 95, 19, 88, 45, 118, 42, 113}, + {120, 127, 89, 94, 131, 132, 50, 53}, {58, 129, 61, 134, 107, 24, 108, 31}, + {66, 139, 145, 128, 69, 140, 150, 135}, {14, 77, 9, 74, 143, 124, 136, 123}, + {126, 119, 85, 156, 121, 112, 82, 155}, {67, 68, 160, 167, 90, 93, 169, 174}, + {171, 98, 172, 101, 0, 57, 7, 62}, {109, 86, 44, 175, 106, 81, 43, 168}, + {117, 180, 114, 179, 78, 71, 73, 64}, {183, 164, 38, 125, 176, 163, 33, 122}, + {161, 80, 166, 87, 130, 187, 133, 188}, {141, 138, 190, 185, 36, 35, 63, 56}, + {1, 146, 184, 75, 6, 149, 191, 76}, {47, 40, 4, 3, 182, 177, 157, 154}, + {152, 147, 159, 148, 97, 162, 102, 165}, {170, 173, 115, 116, 25, 30, 144, 151}, + {100, 99, 181, 178, 23, 16, 142, 137}, {158, 153, 111, 104, 189, 186, 12, 11}, +}; + +// 3D coordinate-to-Hilbert key encoder using the 24-state SO(3) table +[[nodiscard]] inline uint64_t float_to_hilbert3d(float x, float y, float z) noexcept { + uint32_t ix = float_to_ordered_int(x) >> 11; + uint32_t iy = float_to_ordered_int(y) >> 11; + uint32_t iz = float_to_ordered_int(z) >> 11; + uint64_t d = 0; + uint32_t state = 0; + + for (int s = 20; s >= 0; --s) { + uint32_t xyz = (((ix >> s) & 1) << 2) | (((iy >> s) & 1) << 1) | ((iz >> s) & 1); + uint8_t trans = HILBERT_3D_TABLE[state][xyz]; + d = (d << 3) | (trans & 7); + state = trans >> 3; + } + return d; +} + +template +[[nodiscard]] inline uint64_t complex_to_hilbert2d(const std::complex& c) noexcept { + return float_to_hilbert2d(static_cast(c.real()), static_cast(c.imag())); +} + +// Transparent comparator for STL compatibility +struct HilbertCompare { + using is_transparent = void; + template + bool operator()(const std::complex& a, const std::complex& b) const noexcept { + return complex_to_hilbert2d(a) < complex_to_hilbert2d(b); + } +}; + +// O(N) 4-pass 16-bit LSD Radix Sort along the Hilbert curve +template void sort_complex_hilbert(std::span> data) { + if (data.empty()) + return; + + struct Payload { + uint64_t key; + std::complex val; + auto operator<=>(const Payload& other) const { + return key <=> other.key; + } + bool operator==(const Payload& other) const { + return key == other.key; + } + }; + + std::vector buf(data.size() * 2); + Payload* src = buf.data(); + Payload* dst = buf.data() + data.size(); + + for (size_t i = 0; i < data.size(); ++i) { + src[i].val = data[i]; + src[i].key = complex_to_hilbert2d(data[i]); + } + + if (data.size() < 256) { + algoat::sorting::IntroSort{}.sort(std::span{src, data.size()}); + for (size_t i = 0; i < data.size(); ++i) { + data[i] = src[i].val; + } + return; + } + + auto get_bucket = [](uint64_t key, int shift) { return (key >> shift) & 0xFFFF; }; + + for (int pass = 0; pass < 4; ++pass) { + size_t count[65536] = {0}; + int shift = pass * 16; + + for (size_t i = 0; i < data.size(); ++i) { + count[get_bucket(src[i].key, shift)]++; + } + + size_t sum = 0; + for (int i = 0; i < 65536; ++i) { + size_t c = count[i]; + count[i] = sum; + sum += c; + } + + for (size_t i = 0; i < data.size(); ++i) { + dst[count[get_bucket(src[i].key, shift)]++] = src[i]; + } + + std::swap(src, dst); + } + + for (size_t i = 0; i < data.size(); ++i) { + data[i] = src[i].val; + } +} + +} // namespace algoat::numerics diff --git a/include/algoat/sorting/hilbert_sort.hpp b/include/algoat/sorting/hilbert_sort.hpp new file mode 100644 index 0000000..29ad600 --- /dev/null +++ b/include/algoat/sorting/hilbert_sort.hpp @@ -0,0 +1,38 @@ +#pragma once + +#include "algoat/numerics/hilbert.hpp" +#include "algoat/sorting/introsort.hpp" + +#include +#include + +namespace algoat::sorting { + +struct HilbertSort { + [[nodiscard]] constexpr std::string_view name() const noexcept { + return "hilbertsort"; + } + + [[nodiscard]] constexpr std::size_t preferred_min_size() const noexcept { + return 256; + } + + template void sort(std::span data) const { + if (data.size() <= 1) { + return; + } + + if constexpr (requires { + data[0].real(); + data[0].imag(); + }) { + numerics::sort_complex_hilbert(data); + } else { + IntroSort{}.sort(data); + } + } +}; + +} // namespace algoat::sorting + +ALGOAT_REGISTER_ALGORITHM("sorting", "hilbertsort", ::algoat::sorting::HilbertSort) diff --git a/include/algoat/sorting/sorting.hpp b/include/algoat/sorting/sorting.hpp index 93bb009..b66d39d 100644 --- a/include/algoat/sorting/sorting.hpp +++ b/include/algoat/sorting/sorting.hpp @@ -17,6 +17,7 @@ #include "algoat/sorting/cyclesort.hpp" #include "algoat/sorting/gnomesort.hpp" #include "algoat/sorting/heapsort.hpp" +#include "algoat/sorting/hilbert_sort.hpp" #include "algoat/sorting/insertionsort.hpp" #include "algoat/sorting/introsort.hpp" #include "algoat/sorting/mergesort.hpp" @@ -60,9 +61,10 @@ concept SortAlgorithm = requires(Algo algo, std::span data) { * Used by algoat::core::Registry for type-safe static dispatch * via @c std::visit without virtual function table overhead. */ -using SortVariant = std::variant; +using SortVariant = + std::variant; } // namespace algoat::sorting diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 1eeb94c..1626ce1 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -30,6 +30,7 @@ algoat_add_test(test_registry core/test_registry.cpp) algoat_add_test(test_dispatcher core/test_dispatcher.cpp) # Sorting tests +algoat_add_test(test_hilbert_sort sorting/test_hilbert_sort.cpp) algoat_add_test(test_comparative_sorts sorting/test_comparative_sorts.cpp) algoat_add_test(test_insertionsort sorting/test_insertionsort.cpp) algoat_add_test(test_quicksort sorting/test_quicksort.cpp) diff --git a/tests/pbt/test_sorting_invariants.cpp b/tests/pbt/test_sorting_invariants.cpp index 50df2ac..166a2ac 100644 --- a/tests/pbt/test_sorting_invariants.cpp +++ b/tests/pbt/test_sorting_invariants.cpp @@ -118,3 +118,24 @@ TYPED_TEST(LinearSortPBT, Invariants) { RC_ASSERT(input == sorted_again); }); } + +TEST(HilbertSortPBT, Invariants) { + rc::check("Hilbert Radix sort invariants (Complex Float)", []() { + auto input_pairs = *rc::gen::container>>( + rc::gen::pair(rc::gen::inRange(-1000, 1000), rc::gen::inRange(-1000, 1000))); + std::vector> input; + for (const auto& p : input_pairs) { + input.emplace_back(p.first, p.second); + } + + auto original = input; + HilbertSort{}.sort(std::span{input}); + + RC_ASSERT(input.size() == original.size()); + RC_ASSERT(std::is_sorted(input.begin(), input.end(), algoat::numerics::HilbertCompare{})); + + auto sorted_again = input; + HilbertSort{}.sort(std::span{sorted_again}); + RC_ASSERT(input == sorted_again); + }); +} diff --git a/tests/sorting/test_hilbert_sort.cpp b/tests/sorting/test_hilbert_sort.cpp new file mode 100644 index 0000000..92800f1 --- /dev/null +++ b/tests/sorting/test_hilbert_sort.cpp @@ -0,0 +1,102 @@ +#include "algoat/numerics/hilbert.hpp" + +#include +#include +#include +#include + +using namespace algoat::numerics; + +TEST(HilbertSortTest, Hilbert2DEncoding) { + // Zero coordinates + uint64_t z = float_to_hilbert2d(0.0f, 0.0f); + EXPECT_GT(z, 0); // Given positive IEEE 0 is halfway in uint32_t + + // Strict monotonicity inside positive octant + uint64_t a = float_to_hilbert2d(1.0f, 1.0f); + uint64_t b = float_to_hilbert2d(2.0f, 2.0f); + EXPECT_NE(a, b); +} + +TEST(HilbertSortTest, Hilbert3DEncoding) { + uint64_t a = float_to_hilbert3d(1.0f, 1.0f, 1.0f); + uint64_t b = float_to_hilbert3d(2.0f, 2.0f, 2.0f); + EXPECT_NE(a, b); +} + +TEST(HilbertSortTest, FallbackSmallSort) { + std::vector> pts = { + {2.0f, 2.0f}, {-1.0f, -1.0f}, {0.0f, 0.0f}, {10.0f, 1.0f}}; + + std::vector> copy = pts; + sort_complex_hilbert(std::span{pts}); + + // Sort copy with std::sort and compare to ensure Fallback is equivalent to our comparator + std::sort(copy.begin(), copy.end(), HilbertCompare{}); + for (size_t i = 0; i < pts.size(); ++i) { + EXPECT_EQ(pts[i], copy[i]); + } +} + +TEST(HilbertSortTest, RadixSortLarge) { + std::mt19937 gen(42); + std::uniform_real_distribution dist(-1000.0f, 1000.0f); + + size_t N = 10000; + std::vector> pts(N); + for (size_t i = 0; i < N; ++i) { + pts[i] = {dist(gen), dist(gen)}; + } + + std::vector> copy = pts; + sort_complex_hilbert(std::span{pts}); + std::sort(copy.begin(), copy.end(), HilbertCompare{}); + + for (size_t i = 0; i < N; ++i) { + EXPECT_EQ(pts[i], copy[i]); + } +} + +TEST(HilbertSortTest, IdenticalCoordinates) { + std::vector> pts = {{1.0f, 1.0f}, {1.0f, 1.0f}, {2.0f, 2.0f}, {1.0f, 1.0f}}; + + std::vector> copy = pts; + sort_complex_hilbert(std::span{pts}); + std::sort(copy.begin(), copy.end(), HilbertCompare{}); + + for (size_t i = 0; i < pts.size(); ++i) { + EXPECT_EQ(pts[i], copy[i]); + } +} + +TEST(HilbertSortTest, DoublePrecision) { + std::vector> pts = {{2.0, 2.0}, {-1.0, -1.0}, {0.0, 0.0}, {10.0, 1.0}}; + + std::vector> copy = pts; + sort_complex_hilbert(std::span{pts}); + std::sort(copy.begin(), copy.end(), HilbertCompare{}); + + for (size_t i = 0; i < pts.size(); ++i) { + EXPECT_EQ(pts[i], copy[i]); + } +} + +TEST(HilbertSortTest, SpatialClustering) { + // Hilbert curves should preserve spatial locality. + // Points close in 2D space should be close in the sorted 1D array. + std::vector> pts = {{0.0f, 0.0f}, {0.1f, 0.1f}, {0.0f, 0.1f}, + {0.1f, 0.0f}, {10.0f, 10.0f}, {10.1f, 10.1f}, + {10.0f, 10.1f}, {10.1f, 10.0f}}; + + sort_complex_hilbert(std::span{pts}); + + // The first 4 points should be grouped together, and the next 4 grouped together. + for (int i = 0; i < 4; ++i) { + EXPECT_LT(std::abs(pts[i].real()), 1.0f); + EXPECT_LT(std::abs(pts[i].imag()), 1.0f); + } + for (int i = 4; i < 8; ++i) { + EXPECT_GT(pts[i].real(), 9.0f); + EXPECT_GT(pts[i].imag(), 9.0f); + } +}