From 4a95c670cc8616ca43b295b4444b186cebceeb84 Mon Sep 17 00:00:00 2001 From: mperikov Date: Sat, 25 Oct 2025 21:48:16 +0300 Subject: [PATCH] All SIMDs are moved to bits.h Alternative implementation without AVX512 Unittests for new functions --- include/bits.h | 167 +++++++- include/bitvector.h | 969 +++++++++++++++++++++---------------------- src/unittests.cpp | 988 +++++++++++++++++++++++--------------------- 3 files changed, 1156 insertions(+), 968 deletions(-) diff --git a/include/bits.h b/include/bits.h index b619bcb..bb2e8ce 100644 --- a/include/bits.h +++ b/include/bits.h @@ -4,10 +4,17 @@ #include #include #include +#include -/** - * TODO: compile/runtime checks fror SIMD - */ + +#if defined(__AVX512VPOPCNTDQ__) && defined(__AVX512F__) && defined(__AVX512BW__) +#define AVX512SUPP +#endif + + +inline uint64_t first_bits_mask(size_t num) { + return num >= 64 ? UINT64_MAX : ((1llu << num) - 1); +} /** * @brief Number of 1 bits in positions 0 .. count - 1 @@ -29,6 +36,8 @@ * The rest is standard, i.e. popcount_epi64 to perform popcount on * 64 bits and then reduce_add to sum the result. */ +#ifdef AVX512SUPP + uint64_t rank_512(const uint64_t *x, uint64_t count) { __m512i a = _mm512_maskz_set1_epi64((1ull << ((count >> 6))) - 1, std::numeric_limits::max()); @@ -36,20 +45,42 @@ uint64_t rank_512(const uint64_t *x, uint64_t count) { std::numeric_limits::max()); __m512i mask = _mm512_shldv_epi64(a, b, _mm512_set1_epi64(count % 64)); - __m512i res = _mm512_load_epi64(x); + __m512i res = _mm512_loadu_epi64(x); res = _mm512_and_epi64(res, mask); __m512i cnt = _mm512_popcnt_epi64(res); return _mm512_reduce_add_epi64(cnt); } +#else + +uint64_t rank_512(const uint64_t *x, uint64_t count) { + + uint64_t last_uint = count >> 6; + + uint64_t pop_val = 0; + + for (int i = 0; i < last_uint; i++) { + pop_val += std::popcount(x[i]); + } + + uint64_t final = x[last_uint] & first_bits_mask(count & 63); + + pop_val += std::popcount(final); + return pop_val; +} + +#endif + + /** * @brief Return position of @p rank 1 bit in @p x */ uint64_t select_64(uint64_t x, uint64_t rank) { - auto a = _pdep_u64(1ull << rank, x); - return _tzcnt_u64(a); + auto a = _pdep_u64(1ull << rank, x); + return _tzcnt_u64(a); } + /** * @brief Return position of @p rank 1 bit in @p x * @details Selecting within 64-bit word can be done @@ -67,14 +98,132 @@ uint64_t select_64(uint64_t x, uint64_t rank) { * It can also be used as an alternative for linear search but i don't * see a proper SIMD algorithm to make it faster. */ +#ifdef AVX512SUPP + uint64_t select_512(const uint64_t *x, uint64_t rank) { - __m512i res = _mm512_load_epi64(x); + __m512i res = _mm512_loadu_epi64(x); std::array counts; - _mm512_store_epi64(counts.data(), _mm512_popcnt_epi64(res)); + _mm512_storeu_epi64(counts.data(), _mm512_popcnt_epi64(res)); size_t i = 0; while (i < 8 && counts[i] <= rank) { rank -= counts[i++]; } return i * 64 + select_64(x[i], rank); -} \ No newline at end of file +} + +#else + +uint64_t select_512(const uint64_t* bits, uint64_t rank) { + size_t i = 0; + int popcount = std::popcount(bits[0]); + while (i < 7 && popcount <= rank) { + rank -= popcount; + popcount = std::popcount(bits[++i]); + } + return i * 64 + select_64(bits[i], rank); +} + +#endif + + +/** + * @brief Compare 4 64-bit numbers of @p x with @p y and + * return the length of the prefix where @p y is less then + */ +#ifdef AVX512SUPP + +uint16_t cmpl_pref_len_256(const uint64_t* x, uint64_t y) { + + auto y_4 = _mm256_set1_epi64x(y); + auto reg_256 = _mm256_loadu_epi64(x); + auto cmp = _mm256_cmpge_epu64_mask(reg_256, y_4); + + return _tzcnt_u16(cmp); +} + +#else + +uint16_t cmpl_pref_len_256(const uint64_t* x, uint64_t y) { + + auto y_4 = _mm256_set1_epi64x(y); + __m256i reg_256 = _mm256_loadu_si256(reinterpret_cast(x)); + + const __m256i offset = _mm256_set1_epi64x(0x8000000000000000ULL); + __m256i x_offset = _mm256_xor_si256(reg_256, offset); + __m256i y_offset = _mm256_xor_si256(y_4, offset); + auto mask = _mm256_movemask_epi8(_mm256_cmpgt_epi64(x_offset, _mm256_sub_epi64(y_offset, _mm256_set1_epi64x(1)))); + + return _tzcnt_u32(mask) >> 3; +} + +#endif + + +/** + * @brief Compare 8 64-bit numbers of @p x with @p y and + * return the length of the prefix where @p y is less then @p x + */ +#ifdef AVX512SUPP + +uint16_t cmpl_pref_len_512(const uint64_t* x, uint64_t y) { + + auto y_8 = _mm512_set1_epi64(y); + auto reg_512 = _mm512_loadu_epi64(x); + auto cmp = _mm512_cmpge_epu64_mask(reg_512, y_8); + + return _tzcnt_u16(cmp); +} + +#else + +uint16_t cmpl_pref_len_512(const uint64_t* x, uint64_t y) { + + uint16_t len = cmpl_pref_len_256(x, y); + + if (len < 4) + return len; + + return len + cmpl_pref_len_256(x + 4, y); +} + +#endif + + +/** + * @brief Compare 32 16-bit numbers of @p x with @p y and + * return the count of numbers where @p x is less then @p y + */ +#ifdef AVX512SUPP + +uint16_t cmpl_count_512(const uint16_t* x, uint64_t y) { + + auto y_32 = _mm512_set1_epi16(y); + auto reg_512 = _mm512_loadu_epi16(x); + auto cmp = _mm512_cmplt_epu16_mask(reg_512, y_32); + return std::popcount(cmp); +} + +#else + +uint16_t cmpl_count_512(const uint16_t* x, uint64_t y) { + + auto y_16 = _mm256_set1_epi16(y); + __m256i reg_256 = _mm256_loadu_si256(reinterpret_cast(x)); + + const __m256i offset = _mm256_set1_epi16(0x8000); + __m256i x_offset = _mm256_xor_si256(reg_256, offset); + __m256i y_offset = _mm256_xor_si256(y_16, offset); + uint32_t mask = _mm256_movemask_epi8(_mm256_cmpgt_epi16(y_offset, x_offset)); + + uint16_t count = std::popcount(mask) >> 1; + + reg_256 = _mm256_loadu_si256(reinterpret_cast(x + 16)); + + x_offset = _mm256_xor_si256(reg_256, offset); + mask = _mm256_movemask_epi8(_mm256_cmpgt_epi16(y_offset, x_offset)); + + return count + (std::popcount(mask) >> 1); +} + +#endif \ No newline at end of file diff --git a/include/bitvector.h b/include/bitvector.h index 6d742d8..1472339 100644 --- a/include/bitvector.h +++ b/include/bitvector.h @@ -1,491 +1,478 @@ -#pragma once - -#include -#include -#include -#include -#include -#include - -#include "bits.h" - -namespace pixie { - -/** - * @brief Non-interleaved non owning version of the bit vector class with rank - * and select operations. - * - * @details - * Implementation follows ideas from - * - * {1} "SPIDER: Improved Succint Rank and Select Performance" - * Matthew D. Laws, Jocelyn Bliven, Kit Conklin, Elyes Laalai, Samuel McCauley, - * Zach S. Sturdevant - * https://github.com/williams-cs/spider - * - * {2} "Engineering compact data structures for rank and select queries on - * bit vectors" Kurpicz F. - * https://github.com/pasta-toolbox/bit_vector - * - * Essentially it is a 2 level structure of - * - Super blocks of 2^16 bits - * - Basic blocks of 512 bits - * - Super block ranks of 64 bits resulting in ~ 0.98% overhead - * - Basic block ranks of 16 bits resulting in ~ 3.125% overhead - * - 64 bit Select samples of each 16384 bit ~ 0.39% overhead - * - * Rank is 2 table lookups + SIMD popcount in a 512 block. - * - * Select is - * - Initial super block guess by precomputed sampling - * - SIMD supported linear scan to find superblock - * --- For random data we expect the scan to quickly find the block - * --- As 8 64-bit ranks fits into a cache line we probably expect a single - * --- cache miss to find a superblock. SIMD allows to perform a single scan - * --- can of 8 ranks in 3-4 ops. - * - SIMD supported linear scan to find basicblock - * --- Super block consits of 128 basic blocks 16 bits rank each - * --- This fits into 4x512-bit cache lines. The best way seems to - * --- be just to linearly scan linearly, AVX-512 allowes to scan - * --- the whole cache line at once. - * - * Currently I didn't compare it against {1} or {2}, performance is compatible - * and is bottlenecked by cache misses. Usage of AVX-512 makes it limited for - * current generation processors. - * - * This implementation does not involve interleaving, it is notable that - * interleaving makes linear scans harder. I'd suggest that if interleaving is - * to be involved it's probably super-basic rank interleaving like in {2}. We - * lose ability to effectively perform linear scan on superblocks as we already - * expect to have good guesses from select samples. Interleanign original data - * with basic block ranks like in {1} some special guessing technique should be - * involved (like {1} does). - */ -class BitVector { -private: - constexpr static size_t kWordSize = 64; - constexpr static size_t kSuperBlockRankIntSize = 64; - constexpr static size_t kBasicBlockRankIntSize = 16; - constexpr static size_t kBasicBlockSize = 512; - constexpr static size_t kSuperBlockSize = 65536; - constexpr static size_t kWordsPerBlock = 8; - constexpr static size_t kSelectSampleFrequency = 16384; - constexpr static size_t kBlocksPerSuperBlock = 128; - - std::vector super_block_rank; - std::vector basic_block_rank; - std::vector select_samples; - const size_t num_bits_; - const size_t padded_size_; - size_t max_rank; - - std::span bits; - - inline uint64_t first_bits_mask(size_t num) const { - return num >= 64 ? UINT64_MAX : ((1llu << num) - 1); - } - - /** - * Precompute rank for fast queries - */ - void build_rank() { - size_t num_superblocks = 8 + (padded_size_ - 1) / kSuperBlockSize; - // Add more blocks to ease SIMD processing - // num_basicblocks to fully cover superblock, i.e. 128 - // This reduces branching in select - num_superblocks = ((num_superblocks + 7) / 8) * 8; - size_t num_basicblocks = num_superblocks * kBlocksPerSuperBlock; - super_block_rank.resize(num_superblocks); - basic_block_rank.resize(num_basicblocks); - - uint64_t super_block_sum = 0; - uint16_t basic_block_sum = 0; - - for (size_t i = 0; i / kBasicBlockSize < basic_block_rank.size(); - i += kWordSize) { - if (i % kSuperBlockSize == 0) { - super_block_sum += basic_block_sum; - super_block_rank[i / kSuperBlockSize] = super_block_sum; - basic_block_sum = 0; - } - if (i % kBasicBlockSize == 0) { - basic_block_rank[i / kBasicBlockSize] = basic_block_sum; - } - if (i / kWordSize < bits.size()) { - basic_block_sum += std::popcount(bits[i / kWordSize]); - } - } - max_rank = super_block_sum + basic_block_sum; - } - - /** - * @brief Calculate select samples - */ - void build_select() { - uint64_t milestone = kSelectSampleFrequency; - uint64_t rank = 0; - select_samples.emplace_back(0); - for (size_t i = 0; i < bits.size(); ++i) { - auto ones = std::popcount(bits[i]); - if (rank + ones >= milestone) { - auto pos = select_64(bits[i], milestone - rank - 1); - // TODO: try including global rank into select samples to save - // a cache miss on global rank scan - select_samples.emplace_back((64 * i + pos) / kSuperBlockSize); - milestone += kSelectSampleFrequency; - } - rank += ones; - } - } - - /** - * @brief first step of the select operation - */ - uint64_t find_superblock(uint64_t rank) const { - uint64_t left = select_samples[rank / kSelectSampleFrequency]; - - auto rank_8 = _mm512_set1_epi64(rank); - - while (left + 7 < super_block_rank.size()) { - auto reg_512 = _mm512_loadu_epi64(&super_block_rank[left]); - auto cmp = _mm512_cmpge_epu64_mask(reg_512, rank_8); - if (cmp) { - return left + _tzcnt_u16(cmp) - 1; - } - left += 8; - } - if (left + 3 < super_block_rank.size()) { - auto reg_256 = _mm256_loadu_epi64(&super_block_rank[left]); - auto cmp = _mm256_cmpge_epu64_mask(reg_256, *(__m256i *)(&rank_8)); - if (cmp) { - return left + _tzcnt_u16(cmp) - 1; - } - left += 4; - } - while (left < super_block_rank.size() && super_block_rank[left] < rank) - left++; - return left - 1; - } - - /** - * @brief SIMD-optimized linear scan - * @details - * Processes 32 16-bit entries at once (full cache line), so there is at most - * 4 iterations. - */ - uint64_t find_basicblock(uint16_t local_rank, uint64_t s_block) const { - auto rank_32 = _mm512_set1_epi16(local_rank); - auto pos = 0; - - for (size_t i = 0; i < 4; ++i) { - auto batch = - _mm512_loadu_epi16(&basic_block_rank[128 * s_block + 32 * i]); - auto cmp = _mm512_cmplt_epu16_mask(batch, rank_32); - auto count = std::popcount(cmp); - if (count < 32) { - return kBlocksPerSuperBlock * s_block + pos + count - 1; - } - pos += 32; - } - return kBlocksPerSuperBlock * s_block + pos - 1; - } - - /** - * @brief Interpolation search with SIMD optimization - * @details - * Similar to find_basicblock but initial guess is based on linear - * interpolation, for random data it should make initial guess correct - * most of the times, we start from the 32 wide block with interpolation - * guess at the center, if we see that select result lie in lower blocks - * we backoff to find_basicblock - */ - uint64_t find_basicblock_is(uint16_t local_rank, uint64_t s_block) const { - auto lower = super_block_rank[s_block]; - auto upper = super_block_rank[s_block + 1]; - - uint64_t pos = kBlocksPerSuperBlock * local_rank / (upper - lower); - pos = pos + 16 < 32 ? 0 : (pos - 16); - pos = pos > 96 ? 96 : pos; - auto rank_32 = _mm512_set1_epi16(local_rank); - while (pos < 96) { - auto batch = _mm512_loadu_epi16(&basic_block_rank[128 * s_block + pos]); - auto cmp = _mm512_cmplt_epu16_mask(batch, rank_32); - auto count = std::popcount(cmp); - if (count == 0) { - return find_basicblock(local_rank, s_block); - } - if (count < 32) { - return kBlocksPerSuperBlock * s_block + pos + count - 1; - } - pos += 32; - } - pos = 96; - auto batch = _mm512_loadu_epi16(&basic_block_rank[128 * s_block + pos]); - auto cmp = _mm512_cmplt_epu16_mask(batch, rank_32); - auto count = std::popcount(cmp); - if (count == 0) { - return find_basicblock(local_rank, s_block); - } - return kBlocksPerSuperBlock * s_block + pos + count - 1; - } - -public: - /** - * Constructor from an external array of uint64_t. - */ - explicit BitVector(std::span bit_vector, size_t num_bits) - : num_bits_(std::min(num_bits, bit_vector.size() * kWordSize)), - padded_size_(((num_bits_ + kWordSize - 1) / kWordSize) * kWordSize), - bits(bit_vector) { - build_rank(); - build_select(); - } - - /** - * Returns the number of bits in the vector - */ - size_t size() const { return num_bits_; } - - /** - * Get the bit at the specified position - */ - int operator[](size_t pos) const { - size_t word_idx = pos / kWordSize; - size_t bit_off = pos % kWordSize; - - return (bits[word_idx] >> bit_off) & 1; - } - - /** - * Rank operation: count the number of 1-bits up to position pos (exclusive) - * rank_1(pos) = number of 1s in positions [0...pos-1] - */ - uint64_t rank(size_t pos) const { - uint64_t b_block = pos / kBasicBlockSize; - uint64_t s_block = pos / kSuperBlockSize; - // Precomputed rank - uint64_t result = super_block_rank[s_block] + basic_block_rank[b_block]; - // Basic block tail - result += rank_512(&bits[b_block * kWordsPerBlock], - pos - (b_block * kBasicBlockSize)); - return result; - } - - uint64_t select(size_t rank) const { - if (rank > max_rank) [[unlikely]] { - return num_bits_; - } - if (rank == 0) [[unlikely]] { - return 0; - } - uint64_t s_block = find_superblock(rank); - rank -= super_block_rank[s_block]; - auto pos = find_basicblock_is(rank, s_block); - rank -= basic_block_rank[pos]; - pos *= kWordsPerBlock; - - // Final search - if (pos + kWordsPerBlock - 1 < kWordsPerBlock) [[unlikely]] { - size_t ones = std::popcount(bits[pos]); - while (pos < bits.size() && ones < rank) { - rank -= ones; - ones = std::popcount(bits[++pos]); - } - return kWordSize * pos + select_64(bits[pos], rank - 1); - } - return kWordSize * pos + select_512(&bits[pos], rank - 1); - } - - /** - * Convert the bit vector to a binary string, for debug purpose only - */ - std::string to_string() const { - std::string result; - result.reserve(num_bits_); - - for (size_t i = 0; i < num_bits_; i++) { - result.push_back(operator[](i) ? '1' : '0'); - } - - return result; - } -}; - -/** - * Iterleaved owning version of the bit vector class with rank and - * select operations. Interleaving in SDS is the optimization - * that stores auxiliary index and original in the same array which - * reduces the average number of cache misses. - * - * If we can afford copying original data then interleaved version - * is prefered over regular one. - * - * Implementation follows the paper - * - * "SPIDER: Improved Succint Rank and Select Performance" - * Matthew D. Laws, Jocelyn Bliven, Kit Conklin, Elyes Laalai, Samuel McCauley, - * Zach S. Sturdevant - * - */ -class BitVectorInterleaved { -private: - constexpr static size_t kWordSize = 64; - constexpr static size_t kSuperBlockRankIntSize = 64; - constexpr static size_t kBasicBlockRankIntSize = 16; - /** - * 496 bits data + 16 bit local rank - */ - constexpr static size_t kBasicBlockSize = 496; - /** - * 63488 = 496 * 128, so position of superblock can be obtained - * from the position of basic block by dividing on 128 or - * right shift on 7 bits which is cheaper then performing another - * division. - */ - constexpr static size_t kSuperBlockSize = 63488; - constexpr static size_t kBlocksPerSuperBlock = 128; - constexpr static size_t kWordsPerBlock = 8; - - const size_t num_bits_; - std::vector bits_interleaved; - std::vector super_block_rank; - - class BitReader { - size_t iterator_64_ = 0; - size_t offset_size_ = 0; - size_t offset_bits_ = 0; - std::span bits_; - - public: - BitReader(std::span bits) : bits_(bits) {} - uint64_t ReadBits64(size_t num_bits) { - if (num_bits > 64) { - num_bits = 64; - } - uint64_t result = offset_bits_ & first_bits_mask(num_bits); - if (offset_size_ >= num_bits) { - offset_bits_ >>= num_bits; - offset_size_ -= num_bits; - return result; - } - uint64_t next = bits_[iterator_64_++]; - result ^= (next & first_bits_mask(num_bits - offset_size_)) - << offset_size_; - offset_bits_ = (num_bits - offset_size_ == 64) - ? 0 - : next >> (num_bits - offset_size_); - offset_size_ = 64 - (num_bits - offset_size_); - return result; - } - }; - -public: - /** - * Constructor from an external array of uint64_t. - */ - explicit BitVectorInterleaved(std::span bit_vector, - size_t num_bits) - : num_bits_(std::min(num_bits, bit_vector.size() * kWordSize)) { - build_rank_interleaved(bit_vector, num_bits); - } - - static inline uint64_t first_bits_mask(size_t num) { - return num >= 64 ? UINT64_MAX : ((1llu << num) - 1); - } - - /** - * Returns the number of bits in the vector - */ - size_t size() const { return num_bits_; } - - /** - * Get the bit at the specified position - */ - int operator[](size_t pos) const { - size_t block_id = pos / kBasicBlockSize; - size_t block_bit = pos - block_id * kBasicBlockSize; - size_t word_id = block_id * kWordsPerBlock + block_bit / kWordSize; - size_t word_bit = block_bit % kWordSize; - kWordSize; - - return (bits_interleaved[word_id] >> word_bit) & 1; - } - - /** - * Precompute rank for fast queries - */ - void build_rank_interleaved(std::span bits, size_t num_bits) { - size_t num_superblocks = 1 + (num_bits_ - 1) / kSuperBlockSize; - super_block_rank.resize(num_superblocks); - size_t num_basicblocks = 1 + (num_bits_ - 1) / kBasicBlockSize; - bits_interleaved.resize(num_basicblocks * (512 / kWordSize)); - - uint64_t super_block_sum = 0; - uint16_t basic_block_sum = 0; - auto bit_reader = BitReader(bits); - - for (size_t i = 0; i * kBasicBlockSize < num_bits; ++i) { - if (i % (kSuperBlockSize / kBasicBlockSize) == 0) { - super_block_sum += basic_block_sum; - super_block_rank[i / (kSuperBlockSize / kBasicBlockSize)] = - super_block_sum; - basic_block_sum = 0; - } - bits_interleaved[i * (kWordsPerBlock) + 7] = - static_cast(basic_block_sum) << 48; - - for (size_t j = 0; j < 7 && kWordSize * (i + j) < num_bits; ++j) { - bits_interleaved[i * (kWordsPerBlock) + j] = bit_reader.ReadBits64( - std::min(64ull, num_bits - i * kBasicBlockSize + j * kWordSize)); - basic_block_sum += - std::popcount(bits_interleaved[i * (kWordsPerBlock) + j]); - } - if ((i + 7) * kWordSize < num_bits) { - auto v = bit_reader.ReadBits64( - std::min(48ull, num_bits - (i * kBasicBlockSize + 7 * kWordSize))); - bits_interleaved[i * (kWordsPerBlock) + 7] ^= v; - basic_block_sum += std::popcount(v); - } - } - } - - /** - * Rank operation: count the number of 1-bits up to position pos (exclusive) - * rank_1(pos) = number of 1s in positions [0...pos-1] - */ - uint64_t rank(size_t pos) const { - // Multiplication/devisions - uint64_t b_block = pos / kBasicBlockSize; - uint64_t s_block = b_block / kBlocksPerSuperBlock; - uint64_t b_block_pos = b_block * kWordsPerBlock; - // Super block rank - uint64_t result = super_block_rank[s_block]; - /** - * Ok, so here's quite the important factor to load 512-bit region - * at &bits_interleaved[b_block_pos], we store local rank as 16 last - * bits of it. Prefetch should guarantee but seems like there is no - * need for it. - */ - // __builtin_prefetch(&bits_interleaved[b_block_pos]); - result += rank_512(&bits_interleaved[b_block_pos], - pos - (b_block * kBasicBlockSize)); - result += bits_interleaved[b_block_pos + 7] >> 48; - return result; - } - - /** - * Convert the bit vector to a binary string, for debug purpose only - */ - std::string to_string() const { - std::string result; - result.reserve(num_bits_); - - for (size_t i = 0; i < num_bits_; i++) { - result.push_back(operator[](i) ? '1' : '0'); - } - - return result; - } -}; - -} // namespace pixie +#pragma once + +#include +#include +#include +#include +#include +#include + +#include "bits.h" + +namespace pixie { + +/** + * @brief Non-interleaved non owning version of the bit vector class with rank + * and select operations. + * + * @details + * Implementation follows ideas from + * + * {1} "SPIDER: Improved Succint Rank and Select Performance" + * Matthew D. Laws, Jocelyn Bliven, Kit Conklin, Elyes Laalai, Samuel McCauley, + * Zach S. Sturdevant + * https://github.com/williams-cs/spider + * + * {2} "Engineering compact data structures for rank and select queries on + * bit vectors" Kurpicz F. + * https://github.com/pasta-toolbox/bit_vector + * + * Essentially it is a 2 level structure of + * - Super blocks of 2^16 bits + * - Basic blocks of 512 bits + * - Super block ranks of 64 bits resulting in ~ 0.98% overhead + * - Basic block ranks of 16 bits resulting in ~ 3.125% overhead + * - 64 bit Select samples of each 16384 bit ~ 0.39% overhead + * + * Rank is 2 table lookups + SIMD popcount in a 512 block. + * + * Select is + * - Initial super block guess by precomputed sampling + * - SIMD supported linear scan to find superblock + * --- For random data we expect the scan to quickly find the block + * --- As 8 64-bit ranks fits into a cache line we probably expect a single + * --- cache miss to find a superblock. SIMD allows to perform a single scan + * --- can of 8 ranks in 3-4 ops. + * - SIMD supported linear scan to find basicblock + * --- Super block consits of 128 basic blocks 16 bits rank each + * --- This fits into 4x512-bit cache lines. The best way seems to + * --- be just to linearly scan linearly, AVX-512 allowes to scan + * --- the whole cache line at once. + * + * Currently I didn't compare it against {1} or {2}, performance is compatible + * and is bottlenecked by cache misses. Usage of AVX-512 makes it limited for + * current generation processors. + * + * This implementation does not involve interleaving, it is notable that + * interleaving makes linear scans harder. I'd suggest that if interleaving is + * to be involved it's probably super-basic rank interleaving like in {2}. We + * lose ability to effectively perform linear scan on superblocks as we already + * expect to have good guesses from select samples. Interleanign original data + * with basic block ranks like in {1} some special guessing technique should be + * involved (like {1} does). + */ +class BitVector { +private: + constexpr static size_t kWordSize = 64; + constexpr static size_t kSuperBlockRankIntSize = 64; + constexpr static size_t kBasicBlockRankIntSize = 16; + constexpr static size_t kBasicBlockSize = 512; + constexpr static size_t kSuperBlockSize = 65536; + constexpr static size_t kWordsPerBlock = 8; + constexpr static size_t kSelectSampleFrequency = 16384; + constexpr static size_t kBlocksPerSuperBlock = 128; + + std::vector super_block_rank; + std::vector basic_block_rank; + std::vector select_samples; + const size_t num_bits_; + const size_t padded_size_; + size_t max_rank; + + std::span bits; + + inline uint64_t first_bits_mask(size_t num) const { + return num >= 64 ? UINT64_MAX : ((1llu << num) - 1); + } + + /** + * Precompute rank for fast queries + */ + void build_rank() { + size_t num_superblocks = 8 + (padded_size_ - 1) / kSuperBlockSize; + // Add more blocks to ease SIMD processing + // num_basicblocks to fully cover superblock, i.e. 128 + // This reduces branching in select + num_superblocks = ((num_superblocks + 7) / 8) * 8; + size_t num_basicblocks = num_superblocks * kBlocksPerSuperBlock; + super_block_rank.resize(num_superblocks); + basic_block_rank.resize(num_basicblocks); + + uint64_t super_block_sum = 0; + uint16_t basic_block_sum = 0; + + for (size_t i = 0; i / kBasicBlockSize < basic_block_rank.size(); + i += kWordSize) { + if (i % kSuperBlockSize == 0) { + super_block_sum += basic_block_sum; + super_block_rank[i / kSuperBlockSize] = super_block_sum; + basic_block_sum = 0; + } + if (i % kBasicBlockSize == 0) { + basic_block_rank[i / kBasicBlockSize] = basic_block_sum; + } + if (i / kWordSize < bits.size()) { + basic_block_sum += std::popcount(bits[i / kWordSize]); + } + } + max_rank = super_block_sum + basic_block_sum; + } + + /** + * @brief Calculate select samples + */ + void build_select() { + uint64_t milestone = kSelectSampleFrequency; + uint64_t rank = 0; + select_samples.emplace_back(0); + for (size_t i = 0; i < bits.size(); ++i) { + auto ones = std::popcount(bits[i]); + if (rank + ones >= milestone) { + auto pos = select_64(bits[i], milestone - rank - 1); + // TODO: try including global rank into select samples to save + // a cache miss on global rank scan + select_samples.emplace_back((64 * i + pos) / kSuperBlockSize); + milestone += kSelectSampleFrequency; + } + rank += ones; + } + } + + /** + * @brief first step of the select operation + */ + uint64_t find_superblock(uint64_t rank) const { + uint64_t left = select_samples[rank / kSelectSampleFrequency]; + + while (left + 7 < super_block_rank.size()) { + auto len = cmpl_pref_len_512(&super_block_rank[left], rank); + if (len < 8) { + return left + len - 1; + } + left += 8; + } + if (left + 3 < super_block_rank.size()) { + auto len = cmpl_pref_len_256(&super_block_rank[left], rank); + if (len < 4) { + return left + len - 1; + } + left += 4; + } + while (left < super_block_rank.size() && super_block_rank[left] < rank) + left++; + return left - 1; + } + + /** + * @brief SIMD-optimized linear scan + * @details + * Processes 32 16-bit entries at once (full cache line), so there is at most + * 4 iterations. + */ + uint64_t find_basicblock(uint16_t local_rank, uint64_t s_block) const { + auto pos = 0; + + for (size_t i = 0; i < 4; ++i) { + auto count = cmpl_count_512(&basic_block_rank[128 * s_block + 32 * i], local_rank); + if (count < 32) { + return kBlocksPerSuperBlock * s_block + pos + count - 1; + } + pos += 32; + } + return kBlocksPerSuperBlock * s_block + pos - 1; + } + + /** + * @brief Interpolation search with SIMD optimization + * @details + * Similar to find_basicblock but initial guess is based on linear + * interpolation, for random data it should make initial guess correct + * most of the times, we start from the 32 wide block with interpolation + * guess at the center, if we see that select result lie in lower blocks + * we backoff to find_basicblock + */ + uint64_t find_basicblock_is(uint16_t local_rank, uint64_t s_block) const { + auto lower = super_block_rank[s_block]; + auto upper = super_block_rank[s_block + 1]; + + uint64_t pos = kBlocksPerSuperBlock * local_rank / (upper - lower); + pos = pos + 16 < 32 ? 0 : (pos - 16); + pos = pos > 96 ? 96 : pos; + while (pos < 96) { + auto count = cmpl_count_512(&basic_block_rank[128 * s_block + pos], local_rank); + if (count == 0) { + return find_basicblock(local_rank, s_block); + } + if (count < 32) { + return kBlocksPerSuperBlock * s_block + pos + count - 1; + } + pos += 32; + } + pos = 96; + auto count = cmpl_count_512(&basic_block_rank[128 * s_block + pos], local_rank); + if (count == 0) { + return find_basicblock(local_rank, s_block); + } + return kBlocksPerSuperBlock * s_block + pos + count - 1; + } + +public: + /** + * Constructor from an external array of uint64_t. + */ + explicit BitVector(std::span bit_vector, size_t num_bits) + : num_bits_(std::min(num_bits, bit_vector.size() * kWordSize)), + padded_size_(((num_bits_ + kWordSize - 1) / kWordSize) * kWordSize), + bits(bit_vector) { + build_rank(); + build_select(); + } + + /** + * Returns the number of bits in the vector + */ + size_t size() const { return num_bits_; } + + /** + * Get the bit at the specified position + */ + int operator[](size_t pos) const { + size_t word_idx = pos / kWordSize; + size_t bit_off = pos % kWordSize; + + return (bits[word_idx] >> bit_off) & 1; + } + + /** + * Rank operation: count the number of 1-bits up to position pos (exclusive) + * rank_1(pos) = number of 1s in positions [0...pos-1] + */ + uint64_t rank(size_t pos) const { + uint64_t b_block = pos / kBasicBlockSize; + uint64_t s_block = pos / kSuperBlockSize; + // Precomputed rank + uint64_t result = super_block_rank[s_block] + basic_block_rank[b_block]; + // Basic block tail + result += rank_512(&bits[b_block * kWordsPerBlock], + pos - (b_block * kBasicBlockSize)); + return result; + } + + uint64_t select(size_t rank) const { + if (rank > max_rank) [[unlikely]] { + return num_bits_; + } + if (rank == 0) [[unlikely]] { + return 0; + } + uint64_t s_block = find_superblock(rank); + rank -= super_block_rank[s_block]; + auto pos = find_basicblock_is(rank, s_block); + rank -= basic_block_rank[pos]; + pos *= kWordsPerBlock; + + // Final search + if (pos + kWordsPerBlock - 1 < kWordsPerBlock) [[unlikely]] { + size_t ones = std::popcount(bits[pos]); + while (pos < bits.size() && ones < rank) { + rank -= ones; + ones = std::popcount(bits[++pos]); + } + return kWordSize * pos + select_64(bits[pos], rank - 1); + } + return kWordSize * pos + select_512(&bits[pos], rank - 1); + } + + /** + * Convert the bit vector to a binary string, for debug purpose only + */ + std::string to_string() const { + std::string result; + result.reserve(num_bits_); + + for (size_t i = 0; i < num_bits_; i++) { + result.push_back(operator[](i) ? '1' : '0'); + } + + return result; + } +}; + +/** + * Iterleaved owning version of the bit vector class with rank and + * select operations. Interleaving in SDS is the optimization + * that stores auxiliary index and original in the same array which + * reduces the average number of cache misses. + * + * If we can afford copying original data then interleaved version + * is prefered over regular one. + * + * Implementation follows the paper + * + * "SPIDER: Improved Succint Rank and Select Performance" + * Matthew D. Laws, Jocelyn Bliven, Kit Conklin, Elyes Laalai, Samuel McCauley, + * Zach S. Sturdevant + * + */ +class BitVectorInterleaved { +private: + constexpr static size_t kWordSize = 64; + constexpr static size_t kSuperBlockRankIntSize = 64; + constexpr static size_t kBasicBlockRankIntSize = 16; + /** + * 496 bits data + 16 bit local rank + */ + constexpr static size_t kBasicBlockSize = 496; + /** + * 63488 = 496 * 128, so position of superblock can be obtained + * from the position of basic block by dividing on 128 or + * right shift on 7 bits which is cheaper then performing another + * division. + */ + constexpr static size_t kSuperBlockSize = 63488; + constexpr static size_t kBlocksPerSuperBlock = 128; + constexpr static size_t kWordsPerBlock = 8; + + const size_t num_bits_; + std::vector bits_interleaved; + std::vector super_block_rank; + + class BitReader { + size_t iterator_64_ = 0; + size_t offset_size_ = 0; + size_t offset_bits_ = 0; + std::span bits_; + + public: + BitReader(std::span bits) : bits_(bits) {} + uint64_t ReadBits64(size_t num_bits) { + if (num_bits > 64) { + num_bits = 64; + } + uint64_t result = offset_bits_ & first_bits_mask(num_bits); + if (offset_size_ >= num_bits) { + offset_bits_ >>= num_bits; + offset_size_ -= num_bits; + return result; + } + uint64_t next = bits_[iterator_64_++]; + result ^= (next & first_bits_mask(num_bits - offset_size_)) + << offset_size_; + offset_bits_ = (num_bits - offset_size_ == 64) + ? 0 + : next >> (num_bits - offset_size_); + offset_size_ = 64 - (num_bits - offset_size_); + return result; + } + }; + +public: + /** + * Constructor from an external array of uint64_t. + */ + explicit BitVectorInterleaved(std::span bit_vector, + size_t num_bits) + : num_bits_(std::min(num_bits, bit_vector.size() * kWordSize)) { + build_rank_interleaved(bit_vector, num_bits); + } + + static inline uint64_t first_bits_mask(size_t num) { + return num >= 64 ? UINT64_MAX : ((1llu << num) - 1); + } + + /** + * Returns the number of bits in the vector + */ + size_t size() const { return num_bits_; } + + /** + * Get the bit at the specified position + */ + int operator[](size_t pos) const { + size_t block_id = pos / kBasicBlockSize; + size_t block_bit = pos - block_id * kBasicBlockSize; + size_t word_id = block_id * kWordsPerBlock + block_bit / kWordSize; + size_t word_bit = block_bit % kWordSize; + kWordSize; + + return (bits_interleaved[word_id] >> word_bit) & 1; + } + + /** + * Precompute rank for fast queries + */ + void build_rank_interleaved(std::span bits, size_t num_bits) { + size_t num_superblocks = 1 + (num_bits_ - 1) / kSuperBlockSize; + super_block_rank.resize(num_superblocks); + size_t num_basicblocks = 1 + (num_bits_ - 1) / kBasicBlockSize; + bits_interleaved.resize(num_basicblocks * (512 / kWordSize)); + + uint64_t super_block_sum = 0; + uint16_t basic_block_sum = 0; + auto bit_reader = BitReader(bits); + + for (size_t i = 0; i * kBasicBlockSize < num_bits; ++i) { + if (i % (kSuperBlockSize / kBasicBlockSize) == 0) { + super_block_sum += basic_block_sum; + super_block_rank[i / (kSuperBlockSize / kBasicBlockSize)] = + super_block_sum; + basic_block_sum = 0; + } + bits_interleaved[i * (kWordsPerBlock) + 7] = + static_cast(basic_block_sum) << 48; + + for (size_t j = 0; j < 7 && kWordSize * (i + j) < num_bits; ++j) { + bits_interleaved[i * (kWordsPerBlock) + j] = bit_reader.ReadBits64( + std::min(64ull, num_bits - i * kBasicBlockSize + j * kWordSize)); + basic_block_sum += + std::popcount(bits_interleaved[i * (kWordsPerBlock) + j]); + } + if ((i + 7) * kWordSize < num_bits) { + auto v = bit_reader.ReadBits64( + std::min(48ull, num_bits - (i * kBasicBlockSize + 7 * kWordSize))); + bits_interleaved[i * (kWordsPerBlock) + 7] ^= v; + basic_block_sum += std::popcount(v); + } + } + } + + /** + * Rank operation: count the number of 1-bits up to position pos (exclusive) + * rank_1(pos) = number of 1s in positions [0...pos-1] + */ + uint64_t rank(size_t pos) const { + // Multiplication/devisions + uint64_t b_block = pos / kBasicBlockSize; + uint64_t s_block = b_block / kBlocksPerSuperBlock; + uint64_t b_block_pos = b_block * kWordsPerBlock; + // Super block rank + uint64_t result = super_block_rank[s_block]; + /** + * Ok, so here's quite the important factor to load 512-bit region + * at &bits_interleaved[b_block_pos], we store local rank as 16 last + * bits of it. Prefetch should guarantee but seems like there is no + * need for it. + */ + // __builtin_prefetch(&bits_interleaved[b_block_pos]); + result += rank_512(&bits_interleaved[b_block_pos], + pos - (b_block * kBasicBlockSize)); + result += bits_interleaved[b_block_pos + 7] >> 48; + return result; + } + + /** + * Convert the bit vector to a binary string, for debug purpose only + */ + std::string to_string() const { + std::string result; + result.reserve(num_bits_); + + for (size_t i = 0; i < num_bits_; i++) { + result.push_back(operator[](i) ? '1' : '0'); + } + + return result; + } +}; + +} // namespace pixie diff --git a/src/unittests.cpp b/src/unittests.cpp index 53f16e2..33f337a 100644 --- a/src/unittests.cpp +++ b/src/unittests.cpp @@ -1,468 +1,520 @@ -#include -#include - -#include "bits.h" -#include "bitvector.h" -#include - -using pixie::BitVector; -using pixie::BitVectorInterleaved; - -TEST(Rank512, Ones) { - std::array a{std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max()}; - for (size_t i = 0; i < 512; ++i) { - auto p = rank_512(a.data(), i); - EXPECT_EQ(p, i); - } -} - -TEST(Rank512, Random) { - std::array a{std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max()}; - - std::mt19937_64 rng(42); - for (size_t t = 0; t < 1000; ++t) { - for (size_t i = 0; i < 8; ++i) { - a[i] = rng(); - } - size_t rank = 0; - for (size_t i = 0; i <= 512; ++i) { - auto p = rank_512(a.data(), i); - ASSERT_EQ(p, rank); - rank += 1 & (a[i >> 6] >> (i & 63)); - } - } -} - -TEST(Select64, Ones) { - uint64_t x = std::numeric_limits::max(); - for (size_t i = 0; i < 64; ++i) { - auto p = select_64(x, i); - EXPECT_EQ(p, i); - } -} - -TEST(Select64, Random) { - uint64_t a; - - std::mt19937_64 rng(42); - for (size_t t = 0; t < 1000; ++t) { - a = rng(); - size_t rank = 0; - for (size_t i = 0; i < 64; ++i) { - if (1 & (a >> i)) { - auto p = select_64(a, rank++); - ASSERT_EQ(p, i); - } - } - } -} - -TEST(Select512, Ones) { - std::array a{std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max()}; - for (size_t i = 0; i < 512; ++i) { - auto p = select_512(a.data(), i); - EXPECT_EQ(p, i); - } -} - -TEST(Select512, Random) { - std::array a{std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max()}; - - std::mt19937_64 rng(42); - for (size_t t = 0; t < 1000; ++t) { - for (auto &x : a) { - x = rng(); - } - size_t rank = 0; - for (size_t i = 0; i < 512; ++i) { - if (1 & (a[i >> 6] >> (i & 63))) { - auto p = select_512(a.data(), rank++); - ASSERT_EQ(p, i); - } - } - } -} - -TEST(Select512, RankCompativility) { - std::array a{std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max(), - std::numeric_limits::max()}; - - std::mt19937_64 rng(42); - for (size_t t = 0; t < 1000; ++t) { - for (auto &x : a) { - x = rng(); - } - size_t rank = 0; - for (size_t i = 0; i < 512; ++i) { - if (1 & (a[i >> 6] >> (i & 63))) { - auto p = select_512(a.data(), rank++); - ASSERT_EQ(p, i); - auto r = rank_512(a.data(), p); - ASSERT_EQ(r + 1, rank); - r = rank_512(a.data(), p + 1); - ASSERT_EQ(r, rank); - } - } - } -} - -TEST(BitVectorTest, Basic) { - std::vector bits = {0b101010101}; - BitVector bv(bits, 9); - EXPECT_EQ(bv.size(), 9); - EXPECT_EQ(bv[0], 1); - EXPECT_EQ(bv[1], 0); - EXPECT_EQ(bv[2], 1); - EXPECT_EQ(bv[3], 0); - EXPECT_EQ(bv[4], 1); - EXPECT_EQ(bv[5], 0); - EXPECT_EQ(bv[6], 1); - EXPECT_EQ(bv[7], 0); - EXPECT_EQ(bv[8], 1); -} - -TEST(BitVectorTest, MultipleWords) { - std::vector bits = {0xAAAAAAAAAAAAAAAA, 0x5555555555555555}; - BitVector bv(bits, 128); - EXPECT_EQ(bv.size(), 128); - - // First word: 0xAAAAAAAAAAAAAAAA (alternating 1s and 0s, starting with 0) - for (size_t i = 0; i < 64; i++) { - EXPECT_EQ(bv[i], (i % 2 == 1)); - } - - // Second word: 0x5555555555555555 (alternating 0s and 1s, starting with 1) - for (size_t i = 64; i < 128; i++) { - EXPECT_EQ(bv[i], (i % 2 == 0)); - } -} - -TEST(BitVectorTest, ToString) { - std::vector bits = {0b10101010101}; - BitVector bv(bits, 11); - EXPECT_EQ(bv.to_string(), "10101010101"); -} - -TEST(BitVectorTest, RankBasic) { - std::vector bits = {0b10110}; - BitVector bv(bits, 5); - - EXPECT_EQ(bv.rank(0), 0); // No bits - EXPECT_EQ(bv.rank(1), 0); // 0 - EXPECT_EQ(bv.rank(2), 1); // 10 - EXPECT_EQ(bv.rank(3), 2); // 110 - EXPECT_EQ(bv.rank(4), 2); // 0110 - EXPECT_EQ(bv.rank(5), 3); // 10110 -} - -TEST(BitVectorTest, RankWithZeros) { - std::vector bits = {0}; - BitVector bv(bits, 5); - - for (size_t i = 0; i <= 5; i++) { - EXPECT_EQ(bv.rank(i), 0); - } -} - -TEST(BitVectorTest, SelectBasic) { - std::vector bits = {0b1100010110010110}; - BitVector bv(bits, 5); - - EXPECT_EQ(bv.select(1), 1); - EXPECT_EQ(bv.select(2), 2); - EXPECT_EQ(bv.select(3), 4); - EXPECT_EQ(bv.select(4), 7); - EXPECT_EQ(bv.select(5), 8); - EXPECT_EQ(bv.select(6), 10); - EXPECT_EQ(bv.select(7), 14); - EXPECT_EQ(bv.select(8), 15); -} - -TEST(BitVectorTest, MainRankTest) { - std::mt19937_64 rng(42); - std::vector bits(65536 * 32); - for (size_t i = 0; i < 65536 * 32; i++) { - bits[i] = rng(); - } - - BitVector bv(bits, 65536 * 32 * 64); - uint64_t rank = 0; - for (size_t i = 0; i < bv.size(); ++i) { - ASSERT_EQ(rank, bv.rank(i)); - rank += bv[i]; - } -} - -TEST(BitVectorTest, MainSelectTest) { - std::mt19937_64 rng(42); - std::vector bits(65536 * 32); - for (size_t i = 0; i < 65536 * 32; i++) { - bits[i] = rng(); - } - - BitVector bv(bits, 65536 * 32 * 64); - uint64_t rank = 0; - - for (size_t i = 0; i < bv.size(); ++i) { - if (bv[i]) { - ASSERT_EQ(bv.select(++rank), i); - ASSERT_EQ(bv.rank(i), rank - 1); - ASSERT_EQ(bv.rank(i + 1), rank); - } - } -} - -TEST(BitVectorTest, BenchmarkSelectTest) { - for (size_t n = 4; n <= (1ull << 34); n <<= 2) { - std::mt19937_64 rng(42); - - std::vector bits(1 + n / 64); - for (auto &x : bits) { - x = rng(); - } - pixie::BitVector bv(bits, n); - - auto max_rank = bv.rank(bv.size()); - // size_t rank = 0; - auto r = bv.rank(bv.size()); - auto s = bv.select(r); - - // for (size_t i = 0; i < 10000 && r - i > 0; ++i) { - // s = bv.select(r - i); - // } - - for (size_t i = 0; i < 20000000; ++i) { - uint64_t rank = rng() % max_rank; - bv.select(rank); - } - } -} - -TEST(BitVectorInterleavedTest, AtTest) { - std::mt19937_64 rng(42); - /* - * TODO: need to check edge cases, at least - * non-multiple of 64 size - */ - std::vector bits(65536 * 32); - for (size_t i = 0; i < 65536 * 32; i++) { - bits[i] = rng(); - } - - BitVectorInterleaved bv(bits, 65536 * 32 * 64); - uint64_t rank = 0; - for (size_t i = 0; i < 65536 * 32; ++i) { - for (size_t j = 0; j < 64; ++j) { - ASSERT_EQ((bits[i] >> j) & 1, bv[i * 64 + j]); - } - } -} - -TEST(BitVectorInterleavedTest, MainTest) { - std::mt19937_64 rng(42); - std::vector bits(65536 * 32); - for (size_t i = 0; i < 65536 * 32; i++) { - bits[i] = rng(); - } - - BitVectorInterleaved bv(bits, 65536 * 32 * 64); - uint64_t rank = 0; - for (size_t i = 0; i < bv.size(); ++i) { - ASSERT_EQ(rank, bv.rank(i)); - rank += bv[i]; - } -} - -// TEST(BitVectorTest, SelectBasic) { -// std::vector bits = {0b10110}; -// BitVector bv(bits, 5); - -// EXPECT_EQ(bv.select(0), 0); // Position of the 1st 1-bit -// EXPECT_EQ(bv.select(1), 2); // Position of the 2nd 1-bit -// EXPECT_EQ(bv.select(2), 3); // Position of the 3rd 1-bit - -// EXPECT_THROW(bv.select(3), std::out_of_range); // Only 3 bits set -// } - -// TEST(BitVectorTest, SelectWithZeros) { -// std::vector bits = {0}; -// BitVector bv(bits, 5); - -// EXPECT_THROW(bv.select(0), std::out_of_range); // No bits set -// } - -// TEST(BitVectorTest, SelectWithAllOnes) { -// std::vector bits = {0b11111}; -// BitVector bv(bits, 5); - -// for (size_t i = 0; i < 5; i++) { -// EXPECT_EQ(bv.select(i), i); // i-th 1-bit is at position i -// } -// } - -// TEST(BitVectorTest, SelectMultipleWords) { -// // Create a bitvector with bits set at every third position -// std::vector bits(10, 0); // 10 words, 640 bits total -// size_t num_bits = 10 * 64; - -// for (size_t i = 0; i < num_bits; i++) { -// if (i % 3 == 0) { -// size_t word_idx = i / 64; -// size_t bit_offset = i % 64; -// bits[word_idx] |= (1ULL << bit_offset); -// } -// } - -// BitVector bv(bits, num_bits); - -// // Test select for positions where bits are set -// for (size_t i = 0; i < num_bits / 3; i++) { -// EXPECT_EQ(bv.select(i), i * 3); -// } -// } - -// TEST(BitVectorTest, CrossBlockBoundaries) { -// // Create a bitvector with specific bits set to test block boundaries -// std::vector bits(32, 0); // 32 words = 2048 bits -// size_t num_bits = 32 * 64; - -// // Set bits at basic block boundaries (multiples of 512 bits) -// for (size_t i = 0; i < num_bits; i += 512) { -// size_t word_idx = i / 64; -// size_t bit_offset = i % 64; -// bits[word_idx] |= (1ULL << bit_offset); -// } - -// // Set bits at the end of each basic block (511, 1023, etc.) -// for (size_t i = 511; i < num_bits; i += 512) { -// size_t word_idx = i / 64; -// size_t bit_offset = i % 64; -// bits[word_idx] |= (1ULL << bit_offset); -// } - -// BitVector bv(bits, num_bits); - -// // Test rank at basic block boundaries -// for (size_t i = 0; i <= 4; i++) { -// EXPECT_EQ(bv.rank(i * 512), i * 2); -// } - -// // Test select for these specific bits -// for (size_t i = 0; i < 8; i++) { -// if (i % 2 == 0) { -// EXPECT_EQ(bv.select(i), (i / 2) * 512); -// } else { -// EXPECT_EQ(bv.select(i), (i / 2) * 512 + 511); -// } -// } -// } - -// TEST(BitVectorTest, RandomPattern) { -// // Create a bitvector with a random pattern to test rank and select -// std::vector bits(16, 0); // 16 words = 1024 bits -// size_t num_bits = 1000; // Use 1000 bits (not the full 1024) - -// std::mt19937 rng(42); // Fixed seed for reproducibility - -// // Set random bits -// for (size_t i = 0; i < num_bits; i++) { -// if (rng() % 2 == 1) { // ~50% chance of setting each bit -// size_t word_idx = i / 64; -// size_t bit_offset = i % 64; -// bits[word_idx] |= (1ULL << bit_offset); -// } -// } - -// BitVector bv(bits, num_bits); - -// // Count the number of set bits up to each position -// std::vector expected_ranks(num_bits + 1, 0); -// for (size_t i = 0; i < num_bits; i++) { -// expected_ranks[i+1] = expected_ranks[i] + (bv[i] ? 1 : 0); -// } - -// // Test rank for all positions -// for (size_t i = 0; i <= num_bits; i++) { -// EXPECT_EQ(bv.rank(i), expected_ranks[i]); -// } - -// // Test select for all set bits -// size_t total_ones = expected_ranks[num_bits]; -// std::vector one_positions; -// for (size_t i = 0; i < num_bits; i++) { -// if (bv[i]) { -// one_positions.push_back(i); -// } -// } - -// for (size_t i = 0; i < total_ones; i++) { -// EXPECT_EQ(bv.select(i), one_positions[i]); -// } -// } - -// TEST(BitVectorTest, LargeBitVector) { -// // Test with a larger bitvector to ensure the two-level structure works -// correctly size_t num_words = 1024; // 1024 words = 65536 bits (matches -// kSuperBlockSize) std::vector bits(num_words, 0); size_t -// num_bits = num_words * 64; - -// // Set every 100th bit -// for (size_t i = 0; i < num_bits; i += 100) { -// size_t word_idx = i / 64; -// size_t bit_offset = i % 64; -// bits[word_idx] |= (1ULL << bit_offset); -// } - -// BitVector bv(bits, num_bits); - -// // Test rank at multiples of 100 -// for (size_t i = 0; i <= 650; i++) { -// EXPECT_EQ(bv.rank(i * 100), i); -// } - -// // Test select -// for (size_t i = 0; i < 650; i++) { -// EXPECT_EQ(bv.select(i), i * 100); -// } -// } - -// // Tests for our BitVector implementation -// TEST(BitVectorTest2, EmptyConstructor) { -// std::vector bits = {0b101010101}; -// BitVector bv(bits, 8); -// EXPECT_EQ(bv.size(), 8); -// } +#include +#include + +#include "bits.h" +#include "bitvector.h" +#include + +using pixie::BitVector; +using pixie::BitVectorInterleaved; + +TEST(Rank512, Ones) { + std::array a{std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max()}; + for (size_t i = 0; i < 512; ++i) { + auto p = rank_512(a.data(), i); + EXPECT_EQ(p, i); + } +} + +TEST(Rank512, Random) { + std::array a{std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max()}; + + std::mt19937_64 rng(42); + for (size_t t = 0; t < 1000; ++t) { + for (size_t i = 0; i < 8; ++i) { + a[i] = rng(); + } + size_t rank = 0; + for (size_t i = 0; i <= 512; ++i) { + auto p = rank_512(a.data(), i); + ASSERT_EQ(p, rank); + rank += 1 & (a[i >> 6] >> (i & 63)); + } + } +} + +TEST(Select64, Ones) { + uint64_t x = std::numeric_limits::max(); + for (size_t i = 0; i < 64; ++i) { + auto p = select_64(x, i); + EXPECT_EQ(p, i); + } +} + +TEST(Select64, Random) { + uint64_t a; + + std::mt19937_64 rng(42); + for (size_t t = 0; t < 1000; ++t) { + a = rng(); + size_t rank = 0; + for (size_t i = 0; i < 64; ++i) { + if (1 & (a >> i)) { + auto p = select_64(a, rank++); + ASSERT_EQ(p, i); + } + } + } +} + +TEST(Select512, Ones) { + std::array a{std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max()}; + for (size_t i = 0; i < 512; ++i) { + auto p = select_512(a.data(), i); + EXPECT_EQ(p, i); + } +} + +TEST(Select512, Random) { + std::array a{std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max()}; + + std::mt19937_64 rng(42); + for (size_t t = 0; t < 1000; ++t) { + for (auto &x : a) { + x = rng(); + } + size_t rank = 0; + for (size_t i = 0; i < 512; ++i) { + if (1 & (a[i >> 6] >> (i & 63))) { + auto p = select_512(a.data(), rank++); + ASSERT_EQ(p, i); + } + } + } +} + +TEST(Select512, RankCompativility) { + std::array a{std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max()}; + + std::mt19937_64 rng(42); + for (size_t t = 0; t < 1000; ++t) { + for (auto &x : a) { + x = rng(); + } + size_t rank = 0; + for (size_t i = 0; i < 512; ++i) { + if (1 & (a[i >> 6] >> (i & 63))) { + auto p = select_512(a.data(), rank++); + ASSERT_EQ(p, i); + auto r = rank_512(a.data(), p); + ASSERT_EQ(r + 1, rank); + r = rank_512(a.data(), p + 1); + ASSERT_EQ(r, rank); + } + } + } +} + +TEST(CMPLPREFLEN256, Random) { + std::vector x(4); + std::mt19937_64 rng(42); + for (size_t i = 0; i < 1000; i++) { + uint64_t y = rng(); + uint16_t cnt = 0; + bool fl = 1; + for (size_t j = 0; j < 4;j++) { + x[j] = rng(); + fl &= x[j] < y; + cnt += fl; + } + if (cnt < 4) + EXPECT_EQ(cmpl_pref_len_256(x.data(), y), cnt); + else + EXPECT_GE(cmpl_pref_len_256(x.data(), y), cnt); + } +} + +TEST(CMPLPREFLEN512, Random) { + std::vector x(8); + std::mt19937_64 rng(42); + for (size_t i = 0; i < 1000; i++) { + uint64_t y = rng(); + uint16_t cnt = 0; + bool fl = 1; + for (size_t j = 0; j < 8;j++) { + x[j] = rng(); + fl &= x[j] < y; + cnt += fl; + } + if (cnt < 8) + EXPECT_EQ(cmpl_pref_len_512(x.data(), y), cnt); + else + EXPECT_GE(cmpl_pref_len_512(x.data(), y), cnt); + } +} + +TEST(CMPLCOUNT512, Random) { + std::vector x(32); + std::mt19937 rng(42); + for (size_t i = 0; i < 1000; i++) { + uint16_t y = rng(); + uint16_t cnt = 0; + for (size_t j = 0; j < 32;j++) { + x[j] = rng(); + cnt += x[j] < y; + } + EXPECT_EQ(cmpl_count_512(x.data(), y), cnt); + } +} + +TEST(BitVectorTest, Basic) { + std::vector bits = {0b101010101}; + BitVector bv(bits, 9); + EXPECT_EQ(bv.size(), 9); + EXPECT_EQ(bv[0], 1); + EXPECT_EQ(bv[1], 0); + EXPECT_EQ(bv[2], 1); + EXPECT_EQ(bv[3], 0); + EXPECT_EQ(bv[4], 1); + EXPECT_EQ(bv[5], 0); + EXPECT_EQ(bv[6], 1); + EXPECT_EQ(bv[7], 0); + EXPECT_EQ(bv[8], 1); +} + +TEST(BitVectorTest, MultipleWords) { + std::vector bits = {0xAAAAAAAAAAAAAAAA, 0x5555555555555555}; + BitVector bv(bits, 128); + EXPECT_EQ(bv.size(), 128); + + // First word: 0xAAAAAAAAAAAAAAAA (alternating 1s and 0s, starting with 0) + for (size_t i = 0; i < 64; i++) { + EXPECT_EQ(bv[i], (i % 2 == 1)); + } + + // Second word: 0x5555555555555555 (alternating 0s and 1s, starting with 1) + for (size_t i = 64; i < 128; i++) { + EXPECT_EQ(bv[i], (i % 2 == 0)); + } +} + +TEST(BitVectorTest, ToString) { + std::vector bits = {0b10101010101}; + BitVector bv(bits, 11); + EXPECT_EQ(bv.to_string(), "10101010101"); +} + +TEST(BitVectorTest, RankBasic) { + std::vector bits = {0b10110}; + BitVector bv(bits, 5); + + EXPECT_EQ(bv.rank(0), 0); // No bits + EXPECT_EQ(bv.rank(1), 0); // 0 + EXPECT_EQ(bv.rank(2), 1); // 10 + EXPECT_EQ(bv.rank(3), 2); // 110 + EXPECT_EQ(bv.rank(4), 2); // 0110 + EXPECT_EQ(bv.rank(5), 3); // 10110 +} + +TEST(BitVectorTest, RankWithZeros) { + std::vector bits = {0}; + BitVector bv(bits, 5); + + for (size_t i = 0; i <= 5; i++) { + EXPECT_EQ(bv.rank(i), 0); + } +} + +TEST(BitVectorTest, SelectBasic) { + std::vector bits = {0b1100010110010110}; + BitVector bv(bits, 5); + + EXPECT_EQ(bv.select(1), 1); + EXPECT_EQ(bv.select(2), 2); + EXPECT_EQ(bv.select(3), 4); + EXPECT_EQ(bv.select(4), 7); + EXPECT_EQ(bv.select(5), 8); + EXPECT_EQ(bv.select(6), 10); + EXPECT_EQ(bv.select(7), 14); + EXPECT_EQ(bv.select(8), 15); +} + +TEST(BitVectorTest, MainRankTest) { + std::mt19937_64 rng(42); + std::vector bits(65536 * 32); + for (size_t i = 0; i < 65536 * 32; i++) { + bits[i] = rng(); + } + + BitVector bv(bits, 65536 * 32 * 64); + uint64_t rank = 0; + for (size_t i = 0; i < bv.size(); ++i) { + ASSERT_EQ(rank, bv.rank(i)); + rank += bv[i]; + } +} + +TEST(BitVectorTest, MainSelectTest) { + std::mt19937_64 rng(42); + std::vector bits(65536 * 32); + for (size_t i = 0; i < 65536 * 32; i++) { + bits[i] = rng(); + } + + BitVector bv(bits, 65536 * 32 * 64); + uint64_t rank = 0; + + for (size_t i = 0; i < bv.size(); ++i) { + if (bv[i]) { + ASSERT_EQ(bv.select(++rank), i); + ASSERT_EQ(bv.rank(i), rank - 1); + ASSERT_EQ(bv.rank(i + 1), rank); + } + } +} + +TEST(BitVectorTest, BenchmarkSelectTest) { + for (size_t n = 4; n <= (1ull << 34); n <<= 2) { + std::mt19937_64 rng(42); + + std::vector bits(1 + n / 64); + for (auto &x : bits) { + x = rng(); + } + pixie::BitVector bv(bits, n); + + auto max_rank = bv.rank(bv.size()); + // size_t rank = 0; + auto r = bv.rank(bv.size()); + auto s = bv.select(r); + + // for (size_t i = 0; i < 10000 && r - i > 0; ++i) { + // s = bv.select(r - i); + // } + + for (size_t i = 0; i < 20000000; ++i) { + uint64_t rank = rng() % max_rank; + bv.select(rank); + } + } +} + +TEST(BitVectorInterleavedTest, AtTest) { + std::mt19937_64 rng(42); + /* + * TODO: need to check edge cases, at least + * non-multiple of 64 size + */ + std::vector bits(65536 * 32); + for (size_t i = 0; i < 65536 * 32; i++) { + bits[i] = rng(); + } + + BitVectorInterleaved bv(bits, 65536 * 32 * 64); + uint64_t rank = 0; + for (size_t i = 0; i < 65536 * 32; ++i) { + for (size_t j = 0; j < 64; ++j) { + ASSERT_EQ((bits[i] >> j) & 1, bv[i * 64 + j]); + } + } +} + +TEST(BitVectorInterleavedTest, MainTest) { + std::mt19937_64 rng(42); + std::vector bits(65536 * 32); + for (size_t i = 0; i < 65536 * 32; i++) { + bits[i] = rng(); + } + + BitVectorInterleaved bv(bits, 65536 * 32 * 64); + uint64_t rank = 0; + for (size_t i = 0; i < bv.size(); ++i) { + ASSERT_EQ(rank, bv.rank(i)); + rank += bv[i]; + } +} + +// TEST(BitVectorTest, SelectBasic) { +// std::vector bits = {0b10110}; +// BitVector bv(bits, 5); + +// EXPECT_EQ(bv.select(0), 0); // Position of the 1st 1-bit +// EXPECT_EQ(bv.select(1), 2); // Position of the 2nd 1-bit +// EXPECT_EQ(bv.select(2), 3); // Position of the 3rd 1-bit + +// EXPECT_THROW(bv.select(3), std::out_of_range); // Only 3 bits set +// } + +// TEST(BitVectorTest, SelectWithZeros) { +// std::vector bits = {0}; +// BitVector bv(bits, 5); + +// EXPECT_THROW(bv.select(0), std::out_of_range); // No bits set +// } + +// TEST(BitVectorTest, SelectWithAllOnes) { +// std::vector bits = {0b11111}; +// BitVector bv(bits, 5); + +// for (size_t i = 0; i < 5; i++) { +// EXPECT_EQ(bv.select(i), i); // i-th 1-bit is at position i +// } +// } + +// TEST(BitVectorTest, SelectMultipleWords) { +// // Create a bitvector with bits set at every third position +// std::vector bits(10, 0); // 10 words, 640 bits total +// size_t num_bits = 10 * 64; + +// for (size_t i = 0; i < num_bits; i++) { +// if (i % 3 == 0) { +// size_t word_idx = i / 64; +// size_t bit_offset = i % 64; +// bits[word_idx] |= (1ULL << bit_offset); +// } +// } + +// BitVector bv(bits, num_bits); + +// // Test select for positions where bits are set +// for (size_t i = 0; i < num_bits / 3; i++) { +// EXPECT_EQ(bv.select(i), i * 3); +// } +// } + +// TEST(BitVectorTest, CrossBlockBoundaries) { +// // Create a bitvector with specific bits set to test block boundaries +// std::vector bits(32, 0); // 32 words = 2048 bits +// size_t num_bits = 32 * 64; + +// // Set bits at basic block boundaries (multiples of 512 bits) +// for (size_t i = 0; i < num_bits; i += 512) { +// size_t word_idx = i / 64; +// size_t bit_offset = i % 64; +// bits[word_idx] |= (1ULL << bit_offset); +// } + +// // Set bits at the end of each basic block (511, 1023, etc.) +// for (size_t i = 511; i < num_bits; i += 512) { +// size_t word_idx = i / 64; +// size_t bit_offset = i % 64; +// bits[word_idx] |= (1ULL << bit_offset); +// } + +// BitVector bv(bits, num_bits); + +// // Test rank at basic block boundaries +// for (size_t i = 0; i <= 4; i++) { +// EXPECT_EQ(bv.rank(i * 512), i * 2); +// } + +// // Test select for these specific bits +// for (size_t i = 0; i < 8; i++) { +// if (i % 2 == 0) { +// EXPECT_EQ(bv.select(i), (i / 2) * 512); +// } else { +// EXPECT_EQ(bv.select(i), (i / 2) * 512 + 511); +// } +// } +// } + +// TEST(BitVectorTest, RandomPattern) { +// // Create a bitvector with a random pattern to test rank and select +// std::vector bits(16, 0); // 16 words = 1024 bits +// size_t num_bits = 1000; // Use 1000 bits (not the full 1024) + +// std::mt19937 rng(42); // Fixed seed for reproducibility + +// // Set random bits +// for (size_t i = 0; i < num_bits; i++) { +// if (rng() % 2 == 1) { // ~50% chance of setting each bit +// size_t word_idx = i / 64; +// size_t bit_offset = i % 64; +// bits[word_idx] |= (1ULL << bit_offset); +// } +// } + +// BitVector bv(bits, num_bits); + +// // Count the number of set bits up to each position +// std::vector expected_ranks(num_bits + 1, 0); +// for (size_t i = 0; i < num_bits; i++) { +// expected_ranks[i+1] = expected_ranks[i] + (bv[i] ? 1 : 0); +// } + +// // Test rank for all positions +// for (size_t i = 0; i <= num_bits; i++) { +// EXPECT_EQ(bv.rank(i), expected_ranks[i]); +// } + +// // Test select for all set bits +// size_t total_ones = expected_ranks[num_bits]; +// std::vector one_positions; +// for (size_t i = 0; i < num_bits; i++) { +// if (bv[i]) { +// one_positions.push_back(i); +// } +// } + +// for (size_t i = 0; i < total_ones; i++) { +// EXPECT_EQ(bv.select(i), one_positions[i]); +// } +// } + +// TEST(BitVectorTest, LargeBitVector) { +// // Test with a larger bitvector to ensure the two-level structure works +// correctly size_t num_words = 1024; // 1024 words = 65536 bits (matches +// kSuperBlockSize) std::vector bits(num_words, 0); size_t +// num_bits = num_words * 64; + +// // Set every 100th bit +// for (size_t i = 0; i < num_bits; i += 100) { +// size_t word_idx = i / 64; +// size_t bit_offset = i % 64; +// bits[word_idx] |= (1ULL << bit_offset); +// } + +// BitVector bv(bits, num_bits); + +// // Test rank at multiples of 100 +// for (size_t i = 0; i <= 650; i++) { +// EXPECT_EQ(bv.rank(i * 100), i); +// } + +// // Test select +// for (size_t i = 0; i < 650; i++) { +// EXPECT_EQ(bv.select(i), i * 100); +// } +// } + +// // Tests for our BitVector implementation +// TEST(BitVectorTest2, EmptyConstructor) { +// std::vector bits = {0b101010101}; +// BitVector bv(bits, 8); +// EXPECT_EQ(bv.size(), 8); +// }