From 79c3782def1dfb6b340119b3bdbfaccb8efe20d1 Mon Sep 17 00:00:00 2001 From: Denis Koronchik Date: Mon, 11 Aug 2014 18:25:50 +0300 Subject: [PATCH] [3party] Add succinct library --- 3party/succinct/.gitignore | 24 + 3party/succinct/CMakeLists.txt | 68 ++ 3party/succinct/LICENSE | 13 + 3party/succinct/README.md | 72 ++ 3party/succinct/bit_vector.hpp | 511 +++++++++++++ 3party/succinct/bp_vector.cpp | 709 ++++++++++++++++++ 3party/succinct/bp_vector.hpp | 101 +++ 3party/succinct/broadword.hpp | 185 +++++ 3party/succinct/cartesian_tree.hpp | 182 +++++ 3party/succinct/darray.hpp | 162 ++++ 3party/succinct/darray64.hpp | 115 +++ 3party/succinct/elias_fano.hpp | 274 +++++++ .../succinct/elias_fano_compressed_list.hpp | 74 ++ 3party/succinct/elias_fano_list.hpp | 68 ++ 3party/succinct/forward_enumerator.hpp | 25 + 3party/succinct/gamma_bit_vector.hpp | 134 ++++ 3party/succinct/gamma_vector.hpp | 119 +++ 3party/succinct/intrinsics.hpp | 101 +++ 3party/succinct/mappable_vector.hpp | 124 +++ 3party/succinct/mapper.hpp | 280 +++++++ 3party/succinct/nibble_vector.hpp | 68 ++ 3party/succinct/perftest/.gitignore | 4 + 3party/succinct/perftest/CMakeLists.txt | 9 + .../succinct/perftest/perftest_bp_vector.cpp | 119 +++ .../perftest/perftest_bp_vector_rmq.cpp | 77 ++ .../perftest/perftest_cartesian_tree.cpp | 74 ++ 3party/succinct/perftest/perftest_common.hpp | 33 + .../succinct/perftest/perftest_elias_fano.cpp | 114 +++ 3party/succinct/rs_bit_vector.cpp | 73 ++ 3party/succinct/rs_bit_vector.hpp | 192 +++++ 3party/succinct/succinct_config.hpp.in | 16 + 3party/succinct/tables.hpp | 63 ++ 3party/succinct/test_bit_vector.cpp | 196 +++++ 3party/succinct/test_bp_vector.cpp | 92 +++ 3party/succinct/test_bp_vector_common.hpp | 68 ++ 3party/succinct/test_bp_vector_rmq.cpp | 103 +++ 3party/succinct/test_cartesian_tree.cpp | 121 +++ 3party/succinct/test_common.hpp | 20 + 3party/succinct/test_darray.cpp | 76 ++ 3party/succinct/test_elias_fano.cpp | 71 ++ .../test_elias_fano_compressed_list.cpp | 29 + 3party/succinct/test_gamma_bit_vector.cpp | 62 ++ 3party/succinct/test_gamma_vector.cpp | 62 ++ 3party/succinct/test_mapper.cpp | 77 ++ 3party/succinct/test_rank_select_common.hpp | 123 +++ 3party/succinct/test_rs_bit_vector.cpp | 59 ++ 3party/succinct/test_topk_vector.cpp | 110 +++ 3party/succinct/topk_vector.hpp | 188 +++++ 3party/succinct/util.hpp | 273 +++++++ 3party/succinct/vbyte.hpp | 41 + 50 files changed, 5954 insertions(+) create mode 100644 3party/succinct/.gitignore create mode 100644 3party/succinct/CMakeLists.txt create mode 100644 3party/succinct/LICENSE create mode 100644 3party/succinct/README.md create mode 100644 3party/succinct/bit_vector.hpp create mode 100644 3party/succinct/bp_vector.cpp create mode 100644 3party/succinct/bp_vector.hpp create mode 100644 3party/succinct/broadword.hpp create mode 100644 3party/succinct/cartesian_tree.hpp create mode 100644 3party/succinct/darray.hpp create mode 100644 3party/succinct/darray64.hpp create mode 100644 3party/succinct/elias_fano.hpp create mode 100644 3party/succinct/elias_fano_compressed_list.hpp create mode 100644 3party/succinct/elias_fano_list.hpp create mode 100644 3party/succinct/forward_enumerator.hpp create mode 100644 3party/succinct/gamma_bit_vector.hpp create mode 100644 3party/succinct/gamma_vector.hpp create mode 100644 3party/succinct/intrinsics.hpp create mode 100644 3party/succinct/mappable_vector.hpp create mode 100644 3party/succinct/mapper.hpp create mode 100644 3party/succinct/nibble_vector.hpp create mode 100644 3party/succinct/perftest/.gitignore create mode 100644 3party/succinct/perftest/CMakeLists.txt create mode 100644 3party/succinct/perftest/perftest_bp_vector.cpp create mode 100644 3party/succinct/perftest/perftest_bp_vector_rmq.cpp create mode 100644 3party/succinct/perftest/perftest_cartesian_tree.cpp create mode 100644 3party/succinct/perftest/perftest_common.hpp create mode 100644 3party/succinct/perftest/perftest_elias_fano.cpp create mode 100644 3party/succinct/rs_bit_vector.cpp create mode 100644 3party/succinct/rs_bit_vector.hpp create mode 100644 3party/succinct/succinct_config.hpp.in create mode 100644 3party/succinct/tables.hpp create mode 100644 3party/succinct/test_bit_vector.cpp create mode 100644 3party/succinct/test_bp_vector.cpp create mode 100644 3party/succinct/test_bp_vector_common.hpp create mode 100644 3party/succinct/test_bp_vector_rmq.cpp create mode 100644 3party/succinct/test_cartesian_tree.cpp create mode 100644 3party/succinct/test_common.hpp create mode 100644 3party/succinct/test_darray.cpp create mode 100644 3party/succinct/test_elias_fano.cpp create mode 100644 3party/succinct/test_elias_fano_compressed_list.cpp create mode 100644 3party/succinct/test_gamma_bit_vector.cpp create mode 100644 3party/succinct/test_gamma_vector.cpp create mode 100644 3party/succinct/test_mapper.cpp create mode 100644 3party/succinct/test_rank_select_common.hpp create mode 100644 3party/succinct/test_rs_bit_vector.cpp create mode 100644 3party/succinct/test_topk_vector.cpp create mode 100644 3party/succinct/topk_vector.hpp create mode 100644 3party/succinct/util.hpp create mode 100644 3party/succinct/vbyte.hpp diff --git a/3party/succinct/.gitignore b/3party/succinct/.gitignore new file mode 100644 index 0000000000..eb61b823be --- /dev/null +++ b/3party/succinct/.gitignore @@ -0,0 +1,24 @@ +*~ +libsuccinct.* + +test_bit_vector +test_bp_vector +test_bp_vector_rmq +test_elias_fano +test_mapper +test_rs_bit_vector +test_darray +test_elias_fano_compressed_list +test_gamma_bit_vector +test_gamma_vector +test_cartesian_tree +test_topk_vector + + +# cmake +Makefile +CMakeCache.txt +*.cmake +Testing/ +CMakeFiles/ +succinct_config.hpp diff --git a/3party/succinct/CMakeLists.txt b/3party/succinct/CMakeLists.txt new file mode 100644 index 0000000000..efdcee594f --- /dev/null +++ b/3party/succinct/CMakeLists.txt @@ -0,0 +1,68 @@ +cmake_minimum_required(VERSION 2.6) +project(SUCCINCT) + +configure_file( + ${SUCCINCT_SOURCE_DIR}/succinct_config.hpp.in + ${SUCCINCT_SOURCE_DIR}/succinct_config.hpp) + +option(SUCCINCT_USE_LIBCXX + "Use libc++ with Clang instead of libstdc++ (must be same as that used to compile Boost)" + OFF) +option(SUCCINCT_USE_INTRINSICS + "Use a set of intrinsics available on all x86-64 architectures" + ON) +option(SUCCINCT_USE_POPCNT + "Use popcount intrinsic. Available on x86-64 since SSE4.2." + OFF) + +if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang") + if (SUCCINCT_USE_LIBCXX) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -stdlib=libc++") + endif () +endif () + +if (SUCCINCT_USE_POPCNT) + if (UNIX) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -msse4.2") + endif () + # XXX(ot): what to do for MSVC? +endif () + + +# XXX(ot): enable this on all compilers +if (UNIX) + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wno-missing-braces") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wconversion") +endif () + + +find_package(Boost 1.42.0 COMPONENTS + unit_test_framework iostreams system filesystem REQUIRED) +include_directories(${Boost_INCLUDE_DIRS}) +link_directories (${Boost_LIBRARY_DIRS}) + +include_directories(${PROJECT_SOURCE_DIR}) + +set(SUCCINCT_SOURCES + rs_bit_vector.cpp + bp_vector.cpp + ) + +add_library(succinct STATIC ${SUCCINCT_SOURCES}) + +add_subdirectory(perftest) + +# make and run tests only if library is compiled stand-alone +if (CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME) + enable_testing() + file(GLOB SUCCINCT_TEST_SOURCES test_*.cpp) + foreach(TEST_SRC ${SUCCINCT_TEST_SOURCES}) + get_filename_component (TEST_SRC_NAME ${TEST_SRC} NAME_WE) + add_executable(${TEST_SRC_NAME} ${TEST_SRC}) + target_link_libraries(${TEST_SRC_NAME} + succinct + ${Boost_LIBRARIES} + ) + add_test(${TEST_SRC_NAME} ${TEST_SRC_NAME}) + endforeach(TEST_SRC) +endif () diff --git a/3party/succinct/LICENSE b/3party/succinct/LICENSE new file mode 100644 index 0000000000..77eb33208f --- /dev/null +++ b/3party/succinct/LICENSE @@ -0,0 +1,13 @@ +Copyright 2011 Giuseppe Ottaviano + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. diff --git a/3party/succinct/README.md b/3party/succinct/README.md new file mode 100644 index 0000000000..1e53178e5c --- /dev/null +++ b/3party/succinct/README.md @@ -0,0 +1,72 @@ +succinct +======== + +This library contains the implementation of some succinct data +structures. It is rather undocumented now, but better documentation is +under way. On the other hand, the code is quite extensively +unit-tested. + +The library is meant to be imported as a git submodule in other +projects and then included as a CMake subdirectory. See the unit +tests, and the [semi_index](https://github.com/ot/semi_index) and +[path_decomposed_tries](https://github.com/ot/path_decomposed_tries) +projects for examples. + +How to build the code +--------------------- + +### Dependencies ### + +The following dependencies have to be installed to compile the library. + +* CMake >= 2.6, for the build system +* Boost >= 1.42 + +### Supported systems ### + +The library is developed and tested mainly on Linux and Mac OS X, and +it has been tested also on Windows 7. + +The code is designed for 64-bit architectures. It has been tested on +32-bit Linux as well, but it is significantly slower. To compile the +library on 32-bit architectures it is necessary to disable intrinsics +support, passing -DSUCCINCT_USE_INTRINSICS=OFF to cmake. + +### Building on Unix ### + +The project uses CMake. To build it on Unix systems it should be +sufficient to do the following: + + $ cmake . + $ make + +It is also advised to perform a `make test`, which runs the unit +tests. + +### Builing on Mac OS X ### + +Same instructions for Unix apply, with one exception: the library must +be compiled with the same standard library used to compile Boost. So, +if libc++ was used with Clang, the following command must be used: + + $ cmake . -DSUCCINCT_USE_LIBCXX=ON + + +### Building on Windows ### + +On Windows, Boost and zlib are not installed in default locations, so +it is necessary to set some environment variables to allow the build +system to find them. + +* For Boost `BOOST_ROOT` must be set to the directory which contains + the `boost` include directory. +* The directories that contain the Boost must be added to `PATH` so + that the executables find them + +Once the env variables are set, the quickest way to build the code is +by using NMake (instead of the default Visual Studio). Run the +following commands in a Visual Studio x64 Command Prompt: + + $ cmake -G "NMake Makefiles" . + $ nmake + $ nmake test diff --git a/3party/succinct/bit_vector.hpp b/3party/succinct/bit_vector.hpp new file mode 100644 index 0000000000..2f8bc4ab8e --- /dev/null +++ b/3party/succinct/bit_vector.hpp @@ -0,0 +1,511 @@ +#pragma once + +#include + +#include + +#include "mappable_vector.hpp" +#include "broadword.hpp" +#include "util.hpp" + +namespace succinct { + + namespace detail { + inline size_t words_for(uint64_t n) + { + return util::ceil_div(n, 64); + } + } + + class bit_vector; + + class bit_vector_builder : boost::noncopyable { + public: + + typedef std::vector bits_type; + + bit_vector_builder(uint64_t size = 0, bool init = 0) + : m_size(size) + { + m_bits.resize(detail::words_for(size), uint64_t(-init)); + if (size) { + m_cur_word = &m_bits.back(); + // clear padding bits + if (init && size % 64) { + *m_cur_word >>= 64 - (size % 64); + } + } + } + + void reserve(uint64_t size) { + m_bits.reserve(detail::words_for(size)); + } + + inline void push_back(bool b) { + uint64_t pos_in_word = m_size % 64; + if (pos_in_word == 0) { + m_bits.push_back(0); + m_cur_word = &m_bits.back(); + } + *m_cur_word |= (uint64_t)b << pos_in_word; + ++m_size; + } + + inline void set(uint64_t pos, bool b) { + uint64_t word = pos / 64; + uint64_t pos_in_word = pos % 64; + + m_bits[word] &= ~(uint64_t(1) << pos_in_word); + m_bits[word] |= uint64_t(b) << pos_in_word; + } + + inline void set_bits(uint64_t pos, uint64_t bits, size_t len) + { + assert(pos + len <= size()); + // check there are no spurious bits + assert(len == 64 || (bits >> len) == 0); + if (!len) return; + uint64_t mask = (len == 64) ? uint64_t(-1) : ((uint64_t(1) << len) - 1); + uint64_t word = pos / 64; + uint64_t pos_in_word = pos % 64; + + m_bits[word] &= ~(mask << pos_in_word); + m_bits[word] |= bits << pos_in_word; + + uint64_t stored = 64 - pos_in_word; + if (stored < len) { + m_bits[word + 1] &= ~(mask >> stored); + m_bits[word + 1] |= bits >> stored; + } + } + + inline void append_bits(uint64_t bits, size_t len) + { + // check there are no spurious bits + assert(len == 64 || (bits >> len) == 0); + if (!len) return; + uint64_t pos_in_word = m_size % 64; + m_size += len; + if (pos_in_word == 0) { + m_bits.push_back(bits); + } else { + *m_cur_word |= bits << pos_in_word; + if (len > 64 - pos_in_word) { + m_bits.push_back(bits >> (64 - pos_in_word)); + } + } + m_cur_word = &m_bits.back(); + } + + inline void zero_extend(uint64_t n) { + m_size += n; + uint64_t needed = detail::words_for(m_size) - m_bits.size(); + if (needed) { + m_bits.insert(m_bits.end(), needed, 0); + m_cur_word = &m_bits.back(); + } + } + + inline void one_extend(uint64_t n) + { + while (n >= 64) { + append_bits(uint64_t(-1), 64); + n -= 64; + } + if (n) { + append_bits(uint64_t(-1) >> (64 - n), n); + } + } + + void append(bit_vector_builder const& rhs) + { + if (!rhs.size()) return; + + uint64_t pos = m_bits.size(); + uint64_t shift = size() % 64; + m_size = size() + rhs.size(); + m_bits.resize(detail::words_for(m_size)); + + if (shift == 0) { // word-aligned, easy case + std::copy(rhs.m_bits.begin(), rhs.m_bits.end(), + m_bits.begin() + ptrdiff_t(pos)); + } else { + uint64_t* cur_word = &m_bits.front() + pos - 1; + for (size_t i = 0; i < rhs.m_bits.size() - 1; ++i) { + uint64_t w = rhs.m_bits[i]; + *cur_word |= w << shift; + *++cur_word = w >> (64 - shift); + } + *cur_word |= rhs.m_bits.back() << shift; + if (cur_word < &m_bits.back()) { + *++cur_word = rhs.m_bits.back() >> (64 - shift); + } + } + m_cur_word = &m_bits.back(); + } + + // reverse in place + void reverse() + { + uint64_t shift = 64 - (size() % 64); + + uint64_t remainder = 0; + for (size_t i = 0; i < m_bits.size(); ++i) { + uint64_t cur_word; + if (shift != 64) { // this should be hoisted out + cur_word = remainder | (m_bits[i] << shift); + remainder = m_bits[i] >> (64 - shift); + } else { + cur_word = m_bits[i]; + } + m_bits[i] = broadword::reverse_bits(cur_word); + } + assert(remainder == 0); + std::reverse(m_bits.begin(), m_bits.end()); + } + + bits_type& move_bits() { + assert(detail::words_for(m_size) == m_bits.size()); + return m_bits; + } + + uint64_t size() const { + return m_size; + } + + void swap(bit_vector_builder& other) + { + m_bits.swap(other.m_bits); + std::swap(m_size, other.m_size); + std::swap(m_cur_word, other.m_cur_word); + } + + private: + bits_type m_bits; + uint64_t m_size; + uint64_t* m_cur_word; + }; + + class bit_vector { + public: + bit_vector() + : m_size(0) + {} + + template + bit_vector(Range const& from) { + std::vector bits; + const uint64_t first_mask = uint64_t(1); + uint64_t mask = first_mask; + uint64_t cur_val = 0; + m_size = 0; + for (typename boost::range_const_iterator::type iter = boost::begin(from); + iter != boost::end(from); + ++iter) { + if (*iter) { + cur_val |= mask; + } + mask <<= 1; + m_size += 1; + if (!mask) { + bits.push_back(cur_val); + mask = first_mask; + cur_val = 0; + } + } + if (mask != first_mask) { + bits.push_back(cur_val); + } + m_bits.steal(bits); + } + + bit_vector(bit_vector_builder* from) { + m_size = from->size(); + m_bits.steal(from->move_bits()); + } + + template + void map(Visitor& visit) { + visit + (m_size, "m_size") + (m_bits, "m_bits"); + } + + void swap(bit_vector& other) { + std::swap(other.m_size, m_size); + other.m_bits.swap(m_bits); + } + + inline size_t size() const { + return m_size; + } + + inline bool operator[](uint64_t pos) const { + assert(pos < m_size); + uint64_t block = pos / 64; + assert(block < m_bits.size()); + uint64_t shift = pos % 64; + return (m_bits[block] >> shift) & 1; + } + + inline uint64_t get_bits(uint64_t pos, uint64_t len) const { + assert(pos + len <= size()); + if (!len) { + return 0; + } + uint64_t block = pos / 64; + uint64_t shift = pos % 64; + uint64_t mask = -(len == 64) | ((1ULL << len) - 1); + if (shift + len <= 64) { + return m_bits[block] >> shift & mask; + } else { + return (m_bits[block] >> shift) | (m_bits[block + 1] << (64 - shift) & mask); + } + } + + // same as get_bits(pos, 64) but it can extend further size(), padding with zeros + inline uint64_t get_word(uint64_t pos) const + { + assert(pos < size()); + uint64_t block = pos / 64; + uint64_t shift = pos % 64; + uint64_t word = m_bits[block] >> shift; + if (shift && block + 1 < m_bits.size()) { + word |= m_bits[block + 1] << (64 - shift); + } + return word; + } + + // unsafe and fast version of get_word, it retrieves at least 56 bits + inline uint64_t get_word56(uint64_t pos) const + { + // XXX check endianness? + const char* ptr = reinterpret_cast(m_bits.data()); + return *(reinterpret_cast(ptr + pos / 8)) >> (pos % 8); + } + + inline uint64_t predecessor0(uint64_t pos) const { + assert(pos < m_size); + uint64_t block = pos / 64; + uint64_t shift = 64 - pos % 64 - 1; + uint64_t word = ~m_bits[block]; + word = (word << shift) >> shift; + + unsigned long ret; + while (!broadword::msb(word, ret)) { + assert(block); + word = ~m_bits[--block]; + }; + return block * 64 + ret; + } + + inline uint64_t successor0(uint64_t pos) const { + assert(pos < m_size); + uint64_t block = pos / 64; + uint64_t shift = pos % 64; + uint64_t word = (~m_bits[block] >> shift) << shift; + + unsigned long ret; + while (!broadword::lsb(word, ret)) { + ++block; + assert(block < m_bits.size()); + word = ~m_bits[block]; + }; + return block * 64 + ret; + } + + inline uint64_t predecessor1(uint64_t pos) const { + assert(pos < m_size); + uint64_t block = pos / 64; + uint64_t shift = 64 - pos % 64 - 1; + uint64_t word = m_bits[block]; + word = (word << shift) >> shift; + + unsigned long ret; + while (!broadword::msb(word, ret)) { + assert(block); + word = m_bits[--block]; + }; + return block * 64 + ret; + } + + inline uint64_t successor1(uint64_t pos) const { + assert(pos < m_size); + uint64_t block = pos / 64; + uint64_t shift = pos % 64; + uint64_t word = (m_bits[block] >> shift) << shift; + + unsigned long ret; + while (!broadword::lsb(word, ret)) { + ++block; + assert(block < m_bits.size()); + word = m_bits[block]; + }; + return block * 64 + ret; + } + + mapper::mappable_vector const& data() const + { + return m_bits; + } + + struct enumerator { + enumerator() + : m_bv(0) + , m_pos(uint64_t(-1)) + {} + + enumerator(bit_vector const& bv, size_t pos) + : m_bv(&bv) + , m_pos(pos) + , m_buf(0) + , m_avail(0) + { + m_bv->data().prefetch(m_pos / 64); + } + + inline bool next() + { + if (!m_avail) fill_buf(); + bool b = m_buf & 1; + m_buf >>= 1; + m_avail -= 1; + m_pos += 1; + return b; + } + + inline uint64_t take(size_t l) + { + if (m_avail < l) fill_buf(); + uint64_t val; + if (l != 64) { + val = m_buf & ((uint64_t(1) << l) - 1); + m_buf >>= l; + } else { + val = m_buf; + } + m_avail -= l; + m_pos += l; + return val; + } + + inline uint64_t skip_zeros() + { + uint64_t zs = 0; + // XXX the loop may be optimized by aligning access + while (!m_buf) { + m_pos += m_avail; + zs += m_avail; + m_avail = 0; + fill_buf(); + } + + uint64_t l = broadword::lsb(m_buf); + m_buf >>= l; + m_buf >>= 1; + m_avail -= l + 1; + m_pos += l + 1; + return zs + l; + } + + inline uint64_t position() const + { + return m_pos; + } + + private: + + inline void fill_buf() + { + m_buf = m_bv->get_word(m_pos); + m_avail = 64; + } + + bit_vector const* m_bv; + size_t m_pos; + uint64_t m_buf; + size_t m_avail; + }; + + struct unary_enumerator { + unary_enumerator() + : m_data(0) + , m_position(0) + , m_buf(0) + {} + + unary_enumerator(bit_vector const& bv, uint64_t pos) + { + m_data = bv.data().data(); + m_position = pos; + m_buf = m_data[pos / 64]; + // clear low bits + m_buf &= uint64_t(-1) << (pos % 64); + } + + uint64_t position() const + { + return m_position; + } + + uint64_t next() + { + unsigned long pos_in_word; + uint64_t buf = m_buf; + while (!broadword::lsb(buf, pos_in_word)) { + m_position += 64; + buf = m_data[m_position / 64]; + } + + m_buf = buf & (buf - 1); // clear LSB + m_position = (m_position & ~uint64_t(63)) + pos_in_word; + return m_position; + } + + // skip to the k-th one after the current position + void skip(uint64_t k) + { + uint64_t skipped = 0; + uint64_t buf = m_buf; + uint64_t w = 0; + while (skipped + (w = broadword::popcount(buf)) <= k) { + skipped += w; + m_position += 64; + buf = m_data[m_position / 64]; + } + assert(buf); + uint64_t pos_in_word = broadword::select_in_word(buf, k - skipped); + m_buf = buf & (uint64_t(-1) << pos_in_word); + m_position = (m_position & ~uint64_t(63)) + pos_in_word; + } + + // skip to the k-th zero after the current position + void skip0(uint64_t k) + { + uint64_t skipped = 0; + uint64_t pos_in_word = m_position % 64; + uint64_t buf = ~m_buf & (uint64_t(-1) << pos_in_word); + uint64_t w = 0; + while (skipped + (w = broadword::popcount(buf)) <= k) { + skipped += w; + m_position += 64; + buf = ~m_data[m_position / 64]; + } + assert(buf); + pos_in_word = broadword::select_in_word(buf, k - skipped); + m_buf = ~buf & (uint64_t(-1) << pos_in_word); + m_position = (m_position & ~uint64_t(63)) + pos_in_word; + } + + private: + uint64_t const* m_data; + uint64_t m_position; + uint64_t m_buf; + }; + + protected: + size_t m_size; + mapper::mappable_vector m_bits; + }; + +} diff --git a/3party/succinct/bp_vector.cpp b/3party/succinct/bp_vector.cpp new file mode 100644 index 0000000000..79df0f47f7 --- /dev/null +++ b/3party/succinct/bp_vector.cpp @@ -0,0 +1,709 @@ +#include "bp_vector.hpp" +#include "util.hpp" + +namespace succinct { + + namespace { + + // XXX(ot): remove useless tables + + class excess_tables + { + public: + excess_tables() { + for (int c = 0; c < 256; ++c) { + for (uint8_t i = 0; i < 9; ++i) { + m_fwd_pos[c][i] = 0; + m_bwd_pos[c][i] = 0; + } + // populate m_fwd_pos, m_fwd_min, and m_fwd_min_idx + int excess = 0; + m_fwd_min[c] = 0; + m_fwd_min_idx[c] = 0; + + for (char i = 0; i < 8; ++i) { + if ((c >> i) & 1) { // opening + ++excess; + } else { // closing + --excess; + if (excess < 0 && + m_fwd_pos[c][-excess] == 0) { // not already found + m_fwd_pos[c][-excess] = uint8_t(i + 1); + } + } + + if (-excess > m_fwd_min[c]) { + m_fwd_min[c] = uint8_t(-excess); + m_fwd_min_idx[c] = uint8_t(i + 1); + } + } + m_fwd_exc[c] = (char)excess; + + // populate m_bwd_pos and m_bwd_min + excess = 0; + m_bwd_min[c] = 0; + + for (uint8_t i = 0; i < 8; ++i) { + if ((c << i) & 128) { // opening + ++excess; + if (excess > 0 && + m_bwd_pos[c][(uint8_t)excess] == 0) { // not already found + m_bwd_pos[c][(uint8_t)excess] = uint8_t(i + 1); + } + } else { // closing + --excess; + } + + m_bwd_min[c] = uint8_t(std::max(excess, (int)m_bwd_min[c])); + } + } + } + + char m_fwd_exc[256]; + + uint8_t m_fwd_pos[256][9]; + uint8_t m_bwd_pos[256][9]; + + uint8_t m_bwd_min[256]; + uint8_t m_fwd_min[256]; + + uint8_t m_fwd_min_idx[256]; + }; + + const static excess_tables tables; + + inline bool find_close_in_word(uint64_t word, uint64_t byte_counts, bp_vector::excess_t cur_exc, uint64_t& ret) + { + assert(cur_exc > 0 && cur_exc <= 64); + const uint64_t cum_exc_step_8 = (uint64_t(cur_exc) + ((2 * byte_counts - 8 * broadword::ones_step_8) << 8)) * broadword::ones_step_8; + + uint64_t min_exc_step_8 = 0; + for (size_t i = 0; i < 8; ++i) { + size_t shift = i * 8; + min_exc_step_8 |= ((uint64_t)(tables.m_fwd_min[(word >> shift) & 0xFF])) << shift; + } + + const uint64_t has_result = broadword::leq_step_8(cum_exc_step_8, min_exc_step_8); + + unsigned long shift; + if (broadword::lsb(has_result, shift)) { + uint8_t bit_pos = tables.m_fwd_pos[(word >> shift) & 0xFF][(cum_exc_step_8 >> shift) & 0xFF]; + assert(bit_pos > 0); + ret = shift + bit_pos - 1; + return true; + } + return false; + } + + inline bool find_open_in_word(uint64_t word, uint64_t byte_counts, bp_vector::excess_t cur_exc, uint64_t& ret) { + assert(cur_exc > 0 && cur_exc <= 64); + const uint64_t rev_byte_counts = broadword::reverse_bytes(byte_counts); + const uint64_t cum_exc_step_8 = (uint64_t(cur_exc) - ((2 * rev_byte_counts - 8 * broadword::ones_step_8) << 8)) * broadword::ones_step_8; + + uint64_t max_exc_step_8 = 0; + for (size_t i = 0; i < 8; ++i) { + size_t shift = i * 8; + max_exc_step_8 |= ((uint64_t)(tables.m_bwd_min[(word >> (64 - shift - 8)) & 0xFF])) << shift; + } + + const uint64_t has_result = broadword::leq_step_8(cum_exc_step_8, max_exc_step_8); + + unsigned long shift; + if (broadword::lsb(has_result, shift)) { + uint8_t bit_pos = tables.m_bwd_pos[(word >> (64 - shift - 8)) & 0xFF][(cum_exc_step_8 >> shift) & 0xFF]; + assert(bit_pos > 0); + ret = 64 - (shift + bit_pos); + return true; + } + return false; + } + + inline void + excess_rmq_in_word(uint64_t word, bp_vector::excess_t& exc, uint64_t word_start, + bp_vector::excess_t& min_exc, uint64_t& min_exc_idx) + { + bp_vector::excess_t min_byte_exc = min_exc; + uint64_t min_byte_idx = 0; + + for (size_t i = 0; i < 8; ++i) { + size_t shift = i * 8; + size_t byte = (word >> shift) & 0xFF; + // m_fwd_min is negated + bp_vector::excess_t cur_min = exc - tables.m_fwd_min[byte]; + + min_byte_idx = (cur_min < min_byte_exc) ? i : min_byte_idx; + min_byte_exc = (cur_min < min_byte_exc) ? cur_min : min_byte_exc; + + exc += tables.m_fwd_exc[byte]; + } + + if (min_byte_exc < min_exc) { + min_exc = min_byte_exc; + uint64_t shift = min_byte_idx * 8; + min_exc_idx = word_start + shift + tables.m_fwd_min_idx[(word >> shift) & 0xFF]; + } + } + } + + inline bool bp_vector::find_close_in_block(uint64_t block_offset, bp_vector::excess_t excess, uint64_t start, uint64_t& ret) const { + if (excess > excess_t((bp_block_size - start) * 64)) { + return false; + } + assert(excess > 0); + for (uint64_t sub_block_offset = start; sub_block_offset < bp_block_size; ++sub_block_offset) { + uint64_t sub_block = block_offset + sub_block_offset; + uint64_t word = m_bits[sub_block]; + uint64_t byte_counts = broadword::byte_counts(word); + assert(excess > 0); + if (excess <= 64) { + if (find_close_in_word(word, byte_counts, excess, ret)) { + ret += sub_block * 64; + return true; + } + } + excess += static_cast(2 * broadword::bytes_sum(byte_counts) - 64); + } + return false; + } + + uint64_t bp_vector::find_close(uint64_t pos) const + { + assert((*this)[pos]); // check there is an opening parenthesis in pos + uint64_t ret = -1U; + // Search in current word + uint64_t word_pos = (pos + 1) / 64; + uint64_t shift = (pos + 1) % 64; + uint64_t shifted_word = m_bits[word_pos] >> shift; + // Pad with "open" + uint64_t padded_word = shifted_word | (-!!shift & (~0ULL << (64 - shift))); + uint64_t byte_counts = broadword::byte_counts(padded_word); + + excess_t word_exc = 1; + if (find_close_in_word(padded_word, byte_counts, word_exc, ret)) { + ret += pos + 1; + return ret; + } + + // Otherwise search in the local block + uint64_t block = word_pos / bp_block_size; + uint64_t block_offset = block * bp_block_size; + uint64_t sub_block = word_pos % bp_block_size; + uint64_t local_rank = broadword::bytes_sum(byte_counts) - shift; // subtract back the padding + excess_t local_excess = static_cast((2 * local_rank) - (64 - shift)); + if (find_close_in_block(block_offset, local_excess + 1, sub_block + 1, ret)) { + return ret; + } + + // Otherwise, find the first appropriate block + excess_t pos_excess = excess(pos); + uint64_t found_block = search_min_tree<1>(block + 1, pos_excess); + uint64_t found_block_offset = found_block * bp_block_size; + excess_t found_block_excess = get_block_excess(found_block); + + // Search in the found block + bool found = find_close_in_block(found_block_offset, found_block_excess - pos_excess, 0, ret); + assert(found); (void)found; + return ret; + } + + inline bool bp_vector::find_open_in_block(uint64_t block_offset, bp_vector::excess_t excess, uint64_t start, uint64_t& ret) const { + if (excess > excess_t(start * 64)) { + return false; + } + assert(excess >= 0); + + for (uint64_t sub_block_offset = start - 1; sub_block_offset + 1 > 0; --sub_block_offset) { + assert(excess > 0); + uint64_t sub_block = block_offset + sub_block_offset; + uint64_t word = m_bits[sub_block]; + uint64_t byte_counts = broadword::byte_counts(word); + if (excess <= 64) { + if (find_open_in_word(word, byte_counts, excess, ret)) { + ret += sub_block * 64; + return true; + } + } + excess -= static_cast(2 * broadword::bytes_sum(byte_counts) - 64); + } + return false; + } + + uint64_t bp_vector::find_open(uint64_t pos) const + { + assert(pos); + uint64_t ret = -1U; + // Search in current word + uint64_t word_pos = (pos / 64); + uint64_t len = pos % 64; + // Rest is padded with "close" + uint64_t shifted_word = -!!len & (m_bits[word_pos] << (64 - len)); + uint64_t byte_counts = broadword::byte_counts(shifted_word); + + excess_t word_exc = 1; + if (find_open_in_word(shifted_word, byte_counts, word_exc, ret)) { + ret += pos - 64; + return ret; + } + + // Otherwise search in the local block + uint64_t block = word_pos / bp_block_size; + uint64_t block_offset = block * bp_block_size; + uint64_t sub_block = word_pos % bp_block_size; + uint64_t local_rank = broadword::bytes_sum(byte_counts); // no need to subtract the padding + excess_t local_excess = -static_cast((2 * local_rank) - len); + if (find_open_in_block(block_offset, local_excess + 1, sub_block, ret)) { + return ret; + } + + // Otherwise, find the first appropriate block + excess_t pos_excess = excess(pos) - 1; + uint64_t found_block = search_min_tree<0>(block - 1, pos_excess); + uint64_t found_block_offset = found_block * bp_block_size; + // Since search is backwards, have to add the current block + excess_t found_block_excess = get_block_excess(found_block + 1); + + // Search in the found block + bool found = find_open_in_block(found_block_offset, found_block_excess - pos_excess, bp_block_size, ret); + assert(found); (void)found; + return ret; + } + + template + inline bool bp_vector::search_block_in_superblock(uint64_t block, excess_t excess, size_t& found_block) const + { + size_t superblock = block / superblock_size; + excess_t superblock_excess = get_block_excess(superblock * superblock_size); + if (direction) { + for (size_t cur_block = block; + cur_block < std::min((superblock + 1) * superblock_size, (size_t)m_block_excess_min.size()); + ++cur_block) { + if (excess >= superblock_excess + m_block_excess_min[cur_block]) { + found_block = cur_block; + return true; + } + } + } else { + for (size_t cur_block = block; + cur_block + 1 >= (superblock * superblock_size) + 1; + --cur_block) { + if (excess >= superblock_excess + m_block_excess_min[cur_block]) { + found_block = cur_block; + return true; + } + } + } + + return false; + } + + inline bp_vector::excess_t bp_vector::get_block_excess(uint64_t block) const { + uint64_t sub_block_idx = block * bp_block_size; + uint64_t block_pos = sub_block_idx * 64; + excess_t excess = static_cast(2 * sub_block_rank(sub_block_idx) - block_pos); + assert(excess >= 0); + return excess; + } + + inline bool bp_vector::in_node_range(uint64_t node, excess_t excess) const { + assert(m_superblock_excess_min[node] != excess_t(size())); + return excess >= m_superblock_excess_min[node]; + } + + template + inline uint64_t bp_vector::search_min_tree(uint64_t block, excess_t excess) const + { + size_t found_block = -1U; + if (search_block_in_superblock(block, excess, found_block)) { + return found_block; + } + + size_t cur_superblock = block / superblock_size; + size_t cur_node = m_internal_nodes + cur_superblock; + while (true) { + assert(cur_node); + bool going_back = (cur_node & 1) == direction; + if (!going_back) { + size_t next_node = direction ? (cur_node + 1) : (cur_node - 1); + if (in_node_range(next_node, excess)) { + cur_node = next_node; + break; + } + } + cur_node /= 2; + } + + assert(cur_node); + + while (cur_node < m_internal_nodes) { + uint64_t next_node = cur_node * 2 + (1 - direction); + if (in_node_range(next_node, excess)) { + cur_node = next_node; + continue; + } + + next_node = direction ? (next_node + 1) : (next_node - 1); + // if it is not one child, it must be the other + assert(in_node_range(next_node, excess)); + cur_node = next_node; + } + + size_t next_superblock = cur_node - m_internal_nodes; + bool ret = search_block_in_superblock(next_superblock * superblock_size + (1 - direction) * (superblock_size - 1), + excess, found_block); + assert(ret); (void)ret; + + return found_block; + } + + + bp_vector::excess_t + bp_vector::excess(uint64_t pos) const + { + return static_cast(2 * rank(pos) - pos); + } + + void + bp_vector::excess_rmq_in_block(uint64_t start, uint64_t end, + bp_vector::excess_t& exc, + bp_vector::excess_t& min_exc, + uint64_t& min_exc_idx) const + { + assert(start <= end); + if (start == end) return; + + assert((start / bp_block_size) == ((end - 1) / bp_block_size)); + for (size_t w = start; w < end; ++w) { + excess_rmq_in_word(m_bits[w], exc, w * 64, + min_exc, min_exc_idx); + } + } + + void + bp_vector::excess_rmq_in_superblock(uint64_t block_start, uint64_t block_end, + bp_vector::excess_t& block_min_exc, + uint64_t& block_min_idx) const + { + assert(block_start <= block_end); + if (block_start == block_end) return; + + uint64_t superblock = block_start / superblock_size; + + assert(superblock == ((block_end - 1) / superblock_size)); + excess_t superblock_excess = get_block_excess(superblock * superblock_size); + + for (uint64_t block = block_start; block < block_end; ++block) { + if (superblock_excess + m_block_excess_min[block] < block_min_exc) { + block_min_exc = superblock_excess + m_block_excess_min[block]; + block_min_idx = block; + } + } + } + + + void + bp_vector::find_min_superblock(uint64_t superblock_start, uint64_t superblock_end, + bp_vector::excess_t& superblock_min_exc, + uint64_t& superblock_min_idx) const { + + if (superblock_start == superblock_end) return; + + uint64_t cur_node = m_internal_nodes + superblock_start; + uint64_t rightmost_span = superblock_start; + + excess_t node_min_exc = m_superblock_excess_min[cur_node]; + uint64_t node_min_idx = cur_node; + + // code below assumes that there is at least one right-turn in + // the node-root-node path, so we must handle this case + // separately + if (superblock_end - superblock_start == 1) { + superblock_min_exc = node_min_exc; + superblock_min_idx = superblock_start; + return; + } + + // go up the tree until we find the lowest node that spans the + // whole superblock range + size_t h = 0; + while (true) { + assert(cur_node); + + if ((cur_node & 1) == 0) { // is a left child + // add right subtree to candidate superblocks + uint64_t right_sibling = cur_node + 1; + rightmost_span += uint64_t(1) << h; + + if (rightmost_span < superblock_end && + m_superblock_excess_min[right_sibling] < node_min_exc) { + node_min_exc = m_superblock_excess_min[right_sibling]; + node_min_idx = right_sibling; + } + + if (rightmost_span >= superblock_end - 1) { + cur_node += 1; + break; + } + } + + cur_node /= 2; // parent + h += 1; + } + + assert(cur_node); + + // go down until we reach superblock_end + while (rightmost_span > superblock_end - 1) { + assert(cur_node < m_superblock_excess_min.size()); + assert(h > 0); + + h -= 1; + uint64_t left_child = cur_node * 2; + uint64_t right_child_span = uint64_t(1) << h; + if ((rightmost_span - right_child_span) >= (superblock_end - 1)) { + // go to left child + rightmost_span -= right_child_span; + cur_node = left_child; + } else { + // go to right child and add left subtree to candidate + // subblocks + if (m_superblock_excess_min[left_child] < node_min_exc) { + node_min_exc = m_superblock_excess_min[left_child]; + node_min_idx = left_child; + } + cur_node = left_child + 1; + } + } + + // check last left-turn + if (rightmost_span < superblock_end && + m_superblock_excess_min[cur_node] < node_min_exc) { + node_min_exc = m_superblock_excess_min[cur_node]; + node_min_idx = cur_node; + } + + assert(rightmost_span == superblock_end - 1); + + // now reach the minimum leaf in the found subtree (cur_node), + // which is entirely contained in the range + if (node_min_exc < superblock_min_exc) { + cur_node = node_min_idx; + while (cur_node < m_internal_nodes) { + cur_node *= 2; + // remember that past-the-end nodes are filled with size() + if (m_superblock_excess_min[cur_node + 1] < + m_superblock_excess_min[cur_node]) { + cur_node += 1; + } + } + + assert(m_superblock_excess_min[cur_node] == node_min_exc); + superblock_min_exc = node_min_exc; + superblock_min_idx = cur_node - m_internal_nodes; + + assert(superblock_min_idx >= superblock_start); + assert(superblock_min_idx < superblock_end); + } + } + + uint64_t bp_vector::excess_rmq(uint64_t a, uint64_t b, excess_t& min_exc) const + { + assert(a <= b); + + excess_t cur_exc = excess(a); + min_exc = cur_exc; + uint64_t min_exc_idx = a; + + if (a == b) { + return min_exc_idx; + } + + uint64_t range_len = b - a; + + uint64_t word_a_idx = a / 64; + uint64_t word_b_idx = (b - 1) / 64; + + // search in word_a + uint64_t shift_a = a % 64; + uint64_t shifted_word_a = m_bits[word_a_idx] >> shift_a; + uint64_t subword_len_a = std::min(64 - shift_a, range_len); + + uint64_t padded_word_a = + (subword_len_a == 64) + ? shifted_word_a + : (shifted_word_a | (~0ULL << subword_len_a)); + + excess_rmq_in_word(padded_word_a, cur_exc, a, + min_exc, min_exc_idx); + + if (word_a_idx == word_b_idx) { + // single word + return min_exc_idx; + } + + uint64_t block_a = word_a_idx / bp_block_size; + uint64_t block_b = word_b_idx / bp_block_size; + + cur_exc -= 64 - excess_t(subword_len_a); // remove padding + + if (block_a == block_b) { + // same block + excess_rmq_in_block(word_a_idx + 1, word_b_idx, + cur_exc, min_exc, min_exc_idx); + + } else { + // search in partial block of word_a + excess_rmq_in_block(word_a_idx + 1, (block_a + 1) * bp_block_size, + cur_exc, min_exc, min_exc_idx); + + // search in blocks + excess_t block_min_exc = min_exc; + uint64_t block_min_idx = -1U; + + uint64_t superblock_a = (block_a + 1) / superblock_size; + uint64_t superblock_b = block_b / superblock_size; + + if (superblock_a == superblock_b) { + // same superblock + excess_rmq_in_superblock(block_a + 1, block_b, + block_min_exc, block_min_idx); + } else { + // partial superblock of a + excess_rmq_in_superblock(block_a + 1, + (superblock_a + 1) * superblock_size, + block_min_exc, + block_min_idx); + + // search min superblock in the min tree + excess_t superblock_min_exc = min_exc; + uint64_t superblock_min_idx = -1U; + find_min_superblock(superblock_a + 1, superblock_b, + superblock_min_exc, superblock_min_idx); + + if (superblock_min_exc < min_exc) { + excess_rmq_in_superblock(superblock_min_idx * superblock_size, + (superblock_min_idx + 1) * superblock_size, + block_min_exc, + block_min_idx); + } + + // partial superblock of b + excess_rmq_in_superblock(superblock_b * superblock_size, + block_b, + block_min_exc, + block_min_idx); + } + + if (block_min_exc < min_exc) { + cur_exc = get_block_excess(block_min_idx); + excess_rmq_in_block(block_min_idx * bp_block_size, + (block_min_idx + 1) * bp_block_size, + cur_exc, min_exc, min_exc_idx); + assert(min_exc == block_min_exc); + } + + // search in partial block of word_b + cur_exc = get_block_excess(block_b); + excess_rmq_in_block(block_b * bp_block_size, word_b_idx, + cur_exc, min_exc, min_exc_idx); + } + + // search in word_b + uint64_t word_b = m_bits[word_b_idx]; + uint64_t offset_b = b % 64; + uint64_t padded_word_b = + (offset_b == 0) + ? word_b + : (word_b | (~0ULL << offset_b)); + + excess_rmq_in_word(padded_word_b, cur_exc, word_b_idx * 64, + min_exc, min_exc_idx); + + assert(min_exc_idx >= a); + assert(min_exc == excess(min_exc_idx)); + + return min_exc_idx; + } + + + void bp_vector::build_min_tree() + { + if (!size()) return; + + std::vector block_excess_min; + excess_t cur_block_min = 0, cur_superblock_excess = 0; + for (uint64_t sub_block = 0; sub_block < m_bits.size(); ++sub_block) { + if (sub_block % bp_block_size == 0) { + if (sub_block % (bp_block_size * superblock_size) == 0) { + cur_superblock_excess = 0; + } + if (sub_block) { + assert(cur_block_min >= std::numeric_limits::min()); + assert(cur_block_min <= std::numeric_limits::max()); + block_excess_min.push_back((block_min_excess_t)cur_block_min); + cur_block_min = cur_superblock_excess; + } + } + uint64_t word = m_bits[sub_block]; + uint64_t mask = 1ULL; + // for last block stop at bit boundary + uint64_t n_bits = + (sub_block == m_bits.size() - 1 && size() % 64) + ? size() % 64 + : 64; + // XXX(ot) use tables.m_fwd_{min,max} + for (uint64_t i = 0; i < n_bits; ++i) { + cur_superblock_excess += (word & mask) ? 1 : -1; + cur_block_min = std::min(cur_block_min, cur_superblock_excess); + mask <<= 1; + } + } + // Flush last block mins + assert(cur_block_min >= std::numeric_limits::min()); + assert(cur_block_min <= std::numeric_limits::max()); + block_excess_min.push_back((block_min_excess_t)cur_block_min); + + size_t n_blocks = util::ceil_div(data().size(), bp_block_size); + assert(n_blocks == block_excess_min.size()); + + size_t n_superblocks = (n_blocks + superblock_size - 1) / superblock_size; + + size_t n_complete_leaves = 1; + while (n_complete_leaves < n_superblocks) n_complete_leaves <<= 1; // XXX(ot): I'm sure this can be done with broadword::msb... + // n_complete_leaves is the smallest power of 2 >= n_superblocks + m_internal_nodes = n_complete_leaves; + size_t treesize = m_internal_nodes + n_superblocks; + + std::vector superblock_excess_min(treesize); + + // Fill in the leaves of the tree + for (size_t superblock = 0; superblock < n_superblocks; ++superblock) { + excess_t cur_super_min = static_cast(size()); + excess_t superblock_excess = get_block_excess(superblock * superblock_size); + + for (size_t block = superblock * superblock_size; + block < std::min((superblock + 1) * superblock_size, n_blocks); + ++block) { + cur_super_min = std::min(cur_super_min, superblock_excess + block_excess_min[block]); + } + assert(cur_super_min >= 0 && cur_super_min < excess_t(size())); + + superblock_excess_min[m_internal_nodes + superblock] = cur_super_min; + } + + // fill in the internal nodes with past-the-boundary values + // (they will also serve as sentinels in debug) + for (size_t node = 0; node < m_internal_nodes; ++node) { + superblock_excess_min[node] = static_cast(size()); + } + + // Fill bottom-up the other layers: each node updates the parent + for (size_t node = treesize - 1; node > 1; --node) { + size_t parent = node / 2; + superblock_excess_min[parent] = std::min(superblock_excess_min[parent], // same node + superblock_excess_min[node]); + } + + m_block_excess_min.steal(block_excess_min); + m_superblock_excess_min.steal(superblock_excess_min); + } +} diff --git a/3party/succinct/bp_vector.hpp b/3party/succinct/bp_vector.hpp new file mode 100644 index 0000000000..55e35fb77f --- /dev/null +++ b/3party/succinct/bp_vector.hpp @@ -0,0 +1,101 @@ +#pragma once + +#include +#include +#include + +#include + +#include "rs_bit_vector.hpp" + +namespace succinct { + + class bp_vector : public rs_bit_vector { + public: + bp_vector() + : rs_bit_vector() + {} + + template + bp_vector(Range const& from, + bool with_select_hints = false, + bool with_select0_hints = false) + : rs_bit_vector(from, with_select_hints, with_select0_hints) + { + build_min_tree(); + } + + template + void map(Visitor& visit) { + rs_bit_vector::map(visit); + visit + (m_internal_nodes, "m_internal_nodes") + (m_block_excess_min, "m_block_excess_min") + (m_superblock_excess_min, "m_superblock_excess_min") + ; + } + + void swap(bp_vector& other) { + rs_bit_vector::swap(other); + std::swap(m_internal_nodes, other.m_internal_nodes); + m_block_excess_min.swap(other.m_block_excess_min); + m_superblock_excess_min.swap(other.m_superblock_excess_min); + } + + uint64_t find_open(uint64_t pos) const; + uint64_t find_close(uint64_t pos) const; + uint64_t enclose(uint64_t pos) const { + assert((*this)[pos]); + return find_open(pos); + } + + typedef int32_t excess_t; // Allow at most 2^31 depth of the tree + + excess_t excess(uint64_t pos) const; + uint64_t excess_rmq(uint64_t a, uint64_t b, excess_t& min_exc) const; + inline uint64_t excess_rmq(uint64_t a, uint64_t b) const { + excess_t foo; + return excess_rmq(a, b, foo); + } + + + protected: + + static const size_t bp_block_size = 4; // to increase confusion, bp block_size is not necessarily rs_bit_vector block_size + static const size_t superblock_size = 32; // number of blocks in superblock + + typedef int16_t block_min_excess_t; // superblock must be at most 2^15 - 1 bits + + bool find_close_in_block(uint64_t pos, excess_t excess, + uint64_t max_sub_blocks, uint64_t& ret) const; + bool find_open_in_block(uint64_t pos, excess_t excess, + uint64_t max_sub_blocks, uint64_t& ret) const; + + void excess_rmq_in_block(uint64_t start, uint64_t end, + bp_vector::excess_t& exc, + bp_vector::excess_t& min_exc, + uint64_t& min_exc_idx) const; + void excess_rmq_in_superblock(uint64_t block_start, uint64_t block_end, + bp_vector::excess_t& block_min_exc, + uint64_t& block_min_idx) const; + void find_min_superblock(uint64_t superblock_start, uint64_t superblock_end, + bp_vector::excess_t& superblock_min_exc, + uint64_t& superblock_min_idx) const; + + + inline excess_t get_block_excess(uint64_t block) const; + inline bool in_node_range(uint64_t node, excess_t excess) const; + + template + inline bool search_block_in_superblock(uint64_t block, excess_t excess, size_t& found_block) const; + + template + inline uint64_t search_min_tree(uint64_t block, excess_t excess) const; + + void build_min_tree(); + + uint64_t m_internal_nodes; + mapper::mappable_vector m_block_excess_min; + mapper::mappable_vector m_superblock_excess_min; + }; +} diff --git a/3party/succinct/broadword.hpp b/3party/succinct/broadword.hpp new file mode 100644 index 0000000000..9e62317201 --- /dev/null +++ b/3party/succinct/broadword.hpp @@ -0,0 +1,185 @@ +#pragma once + +#include +#include "intrinsics.hpp" +#include "tables.hpp" + +namespace succinct { namespace broadword { + + static const uint64_t ones_step_4 = 0x1111111111111111ULL; + static const uint64_t ones_step_8 = 0x0101010101010101ULL; + static const uint64_t ones_step_9 = 1ULL << 0 | 1ULL << 9 | 1ULL << 18 | 1ULL << 27 | 1ULL << 36 | 1ULL << 45 | 1ULL << 54; + static const uint64_t msbs_step_8 = 0x80ULL * ones_step_8; + static const uint64_t msbs_step_9 = 0x100ULL * ones_step_9; + static const uint64_t incr_step_8 = 0x80ULL << 56 | 0x40ULL << 48 | 0x20ULL << 40 | 0x10ULL << 32 | 0x8ULL << 24 | 0x4ULL << 16 | 0x2ULL << 8 | 0x1; + static const uint64_t inv_count_step_9 = 1ULL << 54 | 2ULL << 45 | 3ULL << 36 | 4ULL << 27 | 5ULL << 18 | 6ULL << 9 | 7ULL; + + static const uint64_t magic_mask_1 = 0x5555555555555555ULL; + static const uint64_t magic_mask_2 = 0x3333333333333333ULL; + static const uint64_t magic_mask_3 = 0x0F0F0F0F0F0F0F0FULL; + static const uint64_t magic_mask_4 = 0x00FF00FF00FF00FFULL; + static const uint64_t magic_mask_5 = 0x0000FFFF0000FFFFULL; + static const uint64_t magic_mask_6 = 0x00000000FFFFFFFFULL; + + inline uint64_t leq_step_8(uint64_t x, uint64_t y) + { + return ((((y | msbs_step_8) - (x & ~msbs_step_8)) ^ (x ^ y)) & msbs_step_8) >> 7; + } + + inline uint64_t uleq_step_8(uint64_t x, uint64_t y) + { + return (((((y | msbs_step_8) - (x & ~msbs_step_8)) ^ (x ^ y)) ^ (x & ~y)) & msbs_step_8) >> 7; + } + + inline uint64_t zcompare_step_8(uint64_t x) + { + return ((x | ((x | msbs_step_8) - ones_step_8)) & msbs_step_8) >> 7; + } + + inline uint64_t uleq_step_9(uint64_t x, uint64_t y) + { + return (((((y | msbs_step_9) - (x & ~msbs_step_9)) | (x ^ y)) ^ (x & ~y)) & msbs_step_9 ) >> 8; + } + + inline uint64_t byte_counts(uint64_t x) + { + x = x - ((x & 0xa * ones_step_4) >> 1); + x = (x & 3 * ones_step_4) + ((x >> 2) & 3 * ones_step_4); + x = (x + (x >> 4)) & 0x0f * ones_step_8; + return x; + } + + inline uint64_t bytes_sum(uint64_t x) + { + return x * ones_step_8 >> 56; + } + + inline uint64_t popcount(uint64_t x) + { +#if SUCCINCT_USE_POPCNT + return intrinsics::popcount(x); +#else + return bytes_sum(byte_counts(x)); +#endif + } + + inline uint64_t reverse_bytes(uint64_t x) + { +#if SUCCINCT_USE_INTRINSICS + return intrinsics::byteswap64(x); +#else + x = ((x >> 8) & magic_mask_4) | ((x & magic_mask_4) << 8); + x = ((x >> 16) & magic_mask_5) | ((x & magic_mask_5) << 16); + x = ((x >> 32) ) | ((x ) << 32); + return x; +#endif + } + + inline uint64_t reverse_bits(uint64_t x) + { + x = ((x >> 1) & magic_mask_1) | ((x & magic_mask_1) << 1); + x = ((x >> 2) & magic_mask_2) | ((x & magic_mask_2) << 2); + x = ((x >> 4) & magic_mask_3) | ((x & magic_mask_3) << 4); + return reverse_bytes(x); + } + + inline uint64_t select_in_word(const uint64_t x, const uint64_t k) + { + assert(k < popcount(x)); + + uint64_t byte_sums = byte_counts(x) * ones_step_8; + + const uint64_t k_step_8 = k * ones_step_8; + const uint64_t geq_k_step_8 = (((k_step_8 | msbs_step_8) - byte_sums) & msbs_step_8); +#if SUCCINCT_USE_POPCNT + const uint64_t place = intrinsics::popcount(geq_k_step_8) * 8; +#else + const uint64_t place = ((geq_k_step_8 >> 7) * ones_step_8 >> 53) & ~uint64_t(0x7); +#endif + const uint64_t byte_rank = k - (((byte_sums << 8 ) >> place) & uint64_t(0xFF)); + return place + tables::select_in_byte[((x >> place) & 0xFF ) | (byte_rank << 8)]; + } + + inline uint64_t same_msb(uint64_t x, uint64_t y) + { + return (x ^ y) <= (x & y); + } + + namespace detail { + // Adapted from LSB of Chess Programming Wiki + static const uint8_t debruijn64_mapping[64] = { + 63, 0, 58, 1, 59, 47, 53, 2, + 60, 39, 48, 27, 54, 33, 42, 3, + 61, 51, 37, 40, 49, 18, 28, 20, + 55, 30, 34, 11, 43, 14, 22, 4, + 62, 57, 46, 52, 38, 26, 32, 41, + 50, 36, 17, 19, 29, 10, 13, 21, + 56, 45, 25, 31, 35, 16, 9, 12, + 44, 24, 15, 8, 23, 7, 6, 5 + }; + static const uint64_t debruijn64 = 0x07EDD5E59A4E28C2ULL; + } + + // return the position of the single bit set in the word x + inline uint8_t bit_position(uint64_t x) + { + assert(popcount(x) == 1); + return detail::debruijn64_mapping + [(x * detail::debruijn64) >> 58]; + } + + inline uint8_t msb(uint64_t x, unsigned long& ret) + { +#if SUCCINCT_USE_INTRINSICS + return intrinsics::bsr64(&ret, x); +#else + if (!x) { + return false; + } + + // right-saturate the word + x |= x >> 1; + x |= x >> 2; + x |= x >> 4; + x |= x >> 8; + x |= x >> 16; + x |= x >> 32; + + // isolate the MSB + x ^= x >> 1; + ret = bit_position(x); + + return true; +#endif + } + + inline uint8_t msb(uint64_t x) + { + assert(x); + unsigned long ret = -1U; + msb(x, ret); + return (uint8_t)ret; + } + + inline uint8_t lsb(uint64_t x, unsigned long& ret) + { +#if SUCCINCT_USE_INTRINSICS + return intrinsics::bsf64(&ret, x); +#else + if (!x) { + return false; + } + ret = bit_position(x & -int64_t(x)); + return true; +#endif + } + + inline uint8_t lsb(uint64_t x) + { + assert(x); + unsigned long ret = -1U; + lsb(x, ret); + return (uint8_t)ret; + } + +}} diff --git a/3party/succinct/cartesian_tree.hpp b/3party/succinct/cartesian_tree.hpp new file mode 100644 index 0000000000..8d6214255a --- /dev/null +++ b/3party/succinct/cartesian_tree.hpp @@ -0,0 +1,182 @@ +#pragma once + +#include + +#include + +#include "bp_vector.hpp" +#include "util.hpp" + +namespace succinct { + + // This class implements a cartesian-tree-based RMQ data + // structure, using the 2d-Min-Heap DFUDS representation described + // in "Space-Efficient Preprocessing Schemes for Range Minimum + // Queries on Static Arrays", Johannes Fischer and Volker Heun, + // SIAM J. Comput., 40(2), 465–492. + + // We made a few variations: + // + // - The rmq() operation in the paper checks whether x is parent + // of w - 1, which can be written as select0(x - 1) < + // find_open(w - 1). We use instead the fact that the excess + // between x and w (both excluded) is strictly greater than the + // excess of w, so the formula above holds iff excess(select0(x + // - 1) + 1) <= excess(w). This is faster because a select0 is + // faster than find_open+rank0. + // + // - The construction is done in reverse order so that the input + // array can be traversed left-to-right. This involves + // re-mapping all the indices at query time. Since the array is + // reversed, in ties the leftmost element wins + // + // - Our data structures have 0-based indices, so the operations + // are slightly different from those in the paper + + class cartesian_tree : boost::noncopyable { + public: + + template + class builder { + public: + builder(uint64_t expected_size = 0) + { + if (expected_size) { + m_bp.reserve(2 * expected_size + 2); + } + } + + template + void push_back(T const& val, Comparator const& comp) + { + m_bp.push_back(0); + + while (!m_stack.empty() + && comp(val, m_stack.back())) { // val < m_stack.back() + m_stack.pop_back(); + m_bp.push_back(1); + } + + m_stack.push_back(val); + } + + bit_vector_builder& finalize() + { + // super-root + m_bp.push_back(0); + while (!m_stack.empty()) { + m_stack.pop_back(); + m_bp.push_back(1); + } + m_bp.push_back(1); + + m_bp.reverse(); + return m_bp; + } + + friend class cartesian_tree; + private: + std::vector m_stack; + bit_vector_builder m_bp; + }; + + cartesian_tree() {} + + template + cartesian_tree(builder* b) + { + bp_vector(&b->finalize(), false, true).swap(m_bp); + } + + template + cartesian_tree(Range const& v) + { + build_from_range(v, std::less::type>()); + } + + template + cartesian_tree(Range const& v, Comparator const& comp) + { + build_from_range(v, comp); + } + + // NOTE: this is RMQ in the interval [a, b], b inclusive + // XXX(ot): maybe change this to [a, b), for consistency with + // the rest of the library? + uint64_t rmq(uint64_t a, uint64_t b) const + { + typedef bp_vector::excess_t excess_t; + + assert(a <= b); + if (a == b) return a; + + uint64_t n = size(); + + uint64_t t = m_bp.select0(n - b - 1); + excess_t exc_t = excess_t(t - 2 * (n - b - 1)); + assert(exc_t - 1 == m_bp.excess(t + 1)); + + uint64_t x = m_bp.select0(n - b); + uint64_t y = m_bp.select0(n - a); + + excess_t exc_w; + uint64_t w = m_bp.excess_rmq(x, y, exc_w); + uint64_t rank0_w = (w - uint64_t(exc_w)) / 2; + assert(m_bp[w - 1] == 0); + + uint64_t ret; + if (exc_w >= exc_t - 1) { + ret = b; + } else { + ret = n - rank0_w; + } + + assert(ret >= a); + assert(ret <= b); + return ret; + } + + bp_vector const& get_bp() const + { + return m_bp; + } + + uint64_t size() const + { + return m_bp.size() / 2 - 1; + } + + template + void map(Visitor& visit) + { + visit + (m_bp, "m_bp"); + } + + void swap(cartesian_tree& other) + { + other.m_bp.swap(m_bp); + } + + protected: + + template + void build_from_range(Range const& v, Comparator const& comp) + { + typedef typename + boost::range_value::type value_type; + typedef typename + boost::range_iterator::type iter_type; + + builder b; + for (iter_type it = boost::begin(v); it != boost::end(v); ++it) { + b.push_back(*it, comp); + } + cartesian_tree(&b).swap(*this); + } + + + bp_vector m_bp; + }; + +} diff --git a/3party/succinct/darray.hpp b/3party/succinct/darray.hpp new file mode 100644 index 0000000000..f098c2e811 --- /dev/null +++ b/3party/succinct/darray.hpp @@ -0,0 +1,162 @@ +#pragma once + +#include "bit_vector.hpp" + +namespace succinct { + + namespace detail { + + template + class darray { + public: + + darray() + : m_positions() + {} + + darray(bit_vector const& bv) + : m_positions() + { + mapper::mappable_vector const& data = bv.data(); + + std::vector cur_block_positions; + std::vector block_inventory; + std::vector subblock_inventory; + std::vector overflow_positions; + + for (size_t word_idx = 0; word_idx < data.size(); ++word_idx) { + size_t cur_pos = word_idx * 64; + uint64_t cur_word = WordGetter()(data, word_idx); + unsigned long l; + while (broadword::lsb(cur_word, l)) { + cur_pos += l; + cur_word >>= l; + if (cur_pos >= bv.size()) break; + + cur_block_positions.push_back(cur_pos); + + if (cur_block_positions.size() == block_size) { + flush_cur_block(cur_block_positions, block_inventory, subblock_inventory, overflow_positions); + } + + // can't do >>= l + 1, can be 64 + cur_word >>= 1; + cur_pos += 1; + m_positions += 1; + } + } + if (cur_block_positions.size()) { + flush_cur_block(cur_block_positions, block_inventory, subblock_inventory, overflow_positions); + } + + m_block_inventory.steal(block_inventory); + m_subblock_inventory.steal(subblock_inventory); + m_overflow_positions.steal(overflow_positions); + } + + template + void map(Visitor& visit) { + visit + (m_positions, "m_positions") + (m_block_inventory, "m_block_inventory") + (m_subblock_inventory, "m_subblock_inventory") + (m_overflow_positions, "m_overflow_positions") + ; + } + + void swap(darray& other) { + std::swap(other.m_positions, m_positions); + m_block_inventory.swap(other.m_block_inventory); + m_subblock_inventory.swap(other.m_subblock_inventory); + m_overflow_positions.swap(other.m_overflow_positions); + } + + inline uint64_t select(bit_vector const& bv, uint64_t idx) const + { + assert(idx < num_positions()); + uint64_t block = idx / block_size; + int64_t block_pos = m_block_inventory[block]; + if (block_pos < 0) { + uint64_t overflow_pos = uint64_t(-block_pos - 1); + return m_overflow_positions[overflow_pos + (idx % block_size)]; + } + + size_t subblock = idx / subblock_size; + size_t start_pos = uint64_t(block_pos) + m_subblock_inventory[subblock]; + size_t reminder = idx % subblock_size; + mapper::mappable_vector const& data = bv.data(); + + if (!reminder) { + return start_pos; + } else { + size_t word_idx = start_pos / 64; + size_t word_shift = start_pos % 64; + uint64_t word = WordGetter()(data, word_idx) & (uint64_t(-1) << word_shift); + + while (true) { + size_t popcnt = broadword::popcount(word); + if (reminder < popcnt) break; + reminder -= popcnt; + word = WordGetter()(data, ++word_idx); + } + + return 64 * word_idx + broadword::select_in_word(word, reminder); + } + } + + inline uint64_t num_positions() const { + return m_positions; + } + + protected: + + static void flush_cur_block(std::vector& cur_block_positions, std::vector& block_inventory, + std::vector& subblock_inventory, std::vector& overflow_positions) + { + if (cur_block_positions.back() - cur_block_positions.front() < max_in_block_distance) { + block_inventory.push_back(int64_t(cur_block_positions.front())); + for (size_t i = 0; i < cur_block_positions.size(); i += subblock_size) { + subblock_inventory.push_back(uint16_t(cur_block_positions[i] - cur_block_positions.front())); + } + } else { + block_inventory.push_back(-int64_t(overflow_positions.size()) - 1); + for (size_t i = 0; i < cur_block_positions.size(); ++i) { + overflow_positions.push_back(cur_block_positions[i]); + } + for (size_t i = 0; i < cur_block_positions.size(); i += subblock_size) { + subblock_inventory.push_back(uint16_t(-1)); + } + } + cur_block_positions.clear(); + } + + + + static const size_t block_size = 1024; + static const size_t subblock_size = 32; + static const size_t max_in_block_distance = 1 << 16; + + size_t m_positions; + mapper::mappable_vector m_block_inventory; + mapper::mappable_vector m_subblock_inventory; + mapper::mappable_vector m_overflow_positions; + }; + + struct identity_getter + { + uint64_t operator()(mapper::mappable_vector const& data, size_t idx) const { + return data[idx]; + } + }; + + struct negating_getter + { + uint64_t operator()(mapper::mappable_vector const& data, size_t idx) const { + return ~data[idx]; + } + }; + } + + typedef detail::darray darray1; + typedef detail::darray darray0; +} diff --git a/3party/succinct/darray64.hpp b/3party/succinct/darray64.hpp new file mode 100644 index 0000000000..fdd4c3104b --- /dev/null +++ b/3party/succinct/darray64.hpp @@ -0,0 +1,115 @@ +#pragma once + +#include "bit_vector.hpp" +#include "broadword.hpp" + +namespace succinct { + + struct darray64 + { + struct builder { + builder() + : n_ones(0) + {} + + void append1(size_t skip0 = 0) + { + bits.append_bits(0, skip0); + bits.push_back(1); + + if (n_ones % block_size == 0) { + block_inventory.push_back(bits.size() - 1); + } + if (n_ones % subblock_size == 0) { + subblock_inventory.push_back(uint16_t(bits.size() - 1 - block_inventory[n_ones / block_size])); + } + + n_ones += 1; + } + + size_t n_ones; + bit_vector_builder bits; + std::vector block_inventory; + std::vector subblock_inventory; + }; + + darray64() + : m_num_ones(0) + {} + + darray64(builder* b) + { + m_num_ones = b->n_ones; + bit_vector(&b->bits).swap(m_bits); + m_block_inventory.steal(b->block_inventory); + m_subblock_inventory.steal(b->subblock_inventory); + } + + void swap(darray64& other) + { + std::swap(m_num_ones, other.m_num_ones); + m_bits.swap(other.m_bits); + m_block_inventory.swap(other.m_block_inventory); + m_subblock_inventory.swap(other.m_subblock_inventory); + } + + template + void map(Visitor& visit) { + visit + (m_num_ones, "m_num_ones") + (m_bits, "m_bits") + (m_block_inventory, "m_block_inventory") + (m_subblock_inventory, "m_subblock_inventory") + ; + } + + size_t num_ones() const + { + return m_num_ones; + } + + bit_vector const& bits() const + { + return m_bits; + } + + size_t select(size_t idx) const + { + assert(idx < num_ones()); + size_t block = idx / block_size; + size_t block_pos = m_block_inventory[block]; + + size_t subblock = idx / subblock_size; + size_t start_pos = block_pos + m_subblock_inventory[subblock]; + size_t reminder = idx % subblock_size; + + if (!reminder) { + return start_pos; + } else { + size_t word_idx = start_pos / 64; + size_t word_shift = start_pos % 64; + uint64_t word = m_bits.data()[word_idx] & (uint64_t(-1) << word_shift); + + while (true) { + size_t popcnt = broadword::popcount(word); + if (reminder < popcnt) break; + reminder -= popcnt; + word = m_bits.data()[++word_idx]; + } + + return 64 * word_idx + broadword::select_in_word(word, reminder); + } + } + + protected: + + static const size_t block_size = 1024; // 64 * block_size must fit in an uint16_t (64 is the max sparsity of bits) + static const size_t subblock_size = 64; + + size_t m_num_ones; + bit_vector m_bits; + mapper::mappable_vector m_block_inventory; + mapper::mappable_vector m_subblock_inventory; + + }; +} diff --git a/3party/succinct/elias_fano.hpp b/3party/succinct/elias_fano.hpp new file mode 100644 index 0000000000..4991fcc9b8 --- /dev/null +++ b/3party/succinct/elias_fano.hpp @@ -0,0 +1,274 @@ +#pragma once + +#include "bit_vector.hpp" +#include "darray.hpp" + +namespace succinct { + + class elias_fano { + public: + elias_fano() + : m_size(0) + {} + + struct elias_fano_builder { + elias_fano_builder(uint64_t n, uint64_t m) + : m_n(n) + , m_m(m) + , m_pos(0) + , m_last(0) + , m_l(uint8_t((m && n / m) ? broadword::msb(n / m) : 0)) + , m_high_bits((m + 1) + (n >> m_l) + 1) + { + assert(m_l < 64); // for the correctness of low_mask + m_low_bits.reserve(m * m_l); + } + + inline void push_back(uint64_t i) { + assert(i >= m_last && i <= m_n); + m_last = i; + uint64_t low_mask = (1ULL << m_l) - 1; + + if (m_l) { + m_low_bits.append_bits(i & low_mask, m_l); + } + m_high_bits.set((i >> m_l) + m_pos, 1); + ++m_pos; + assert(m_pos <= m_m); (void)m_m; + } + + friend class elias_fano; + private: + uint64_t m_n; + uint64_t m_m; + uint64_t m_pos; + uint64_t m_last; + uint8_t m_l; + bit_vector_builder m_high_bits; + bit_vector_builder m_low_bits; + }; + + + elias_fano(bit_vector_builder* bvb, bool with_rank_index = true) + { + bit_vector_builder::bits_type& bits = bvb->move_bits(); + uint64_t n = bvb->size(); + + uint64_t m = 0; + for (size_t i = 0; i < bits.size(); ++i) { + m += broadword::popcount(bits[i]); + } + + bit_vector bv(bvb); + elias_fano_builder builder(n, m); + + uint64_t i = 0; + for (uint64_t pos = 0; pos < m; ++pos) { + i = bv.successor1(i); + builder.push_back(i); + ++i; + } + + build(builder, with_rank_index); + } + + elias_fano(elias_fano_builder* builder, bool with_rank_index = true) + { + build(*builder, with_rank_index); + } + + template + void map(Visitor& visit) { + visit + (m_size, "m_size") + (m_high_bits, "m_high_bits") + (m_high_bits_d1, "m_high_bits_d1") + (m_high_bits_d0, "m_high_bits_d0") + (m_low_bits, "m_low_bits") + (m_l, "m_l") + ; + } + + void swap(elias_fano& other) { + std::swap(other.m_size, m_size); + other.m_high_bits.swap(m_high_bits); + other.m_high_bits_d1.swap(m_high_bits_d1); + other.m_high_bits_d0.swap(m_high_bits_d0); + other.m_low_bits.swap(m_low_bits); + std::swap(other.m_l, m_l); + } + + inline uint64_t size() const { + return m_size; + } + + inline uint64_t num_ones() const { + return m_high_bits_d1.num_positions(); + } + + inline bool operator[](uint64_t pos) const { + assert(pos < size()); + assert(m_high_bits_d0.num_positions()); // needs rank index + uint64_t h_rank = pos >> m_l; + uint64_t h_pos = m_high_bits_d0.select(m_high_bits, h_rank); + uint64_t rank = h_pos - h_rank; + uint64_t l_pos = pos & ((1ULL << m_l) - 1); + + while (h_pos > 0 + && m_high_bits[h_pos - 1]) { + --rank; + --h_pos; + uint64_t cur_low_bits = m_low_bits.get_bits(rank * m_l, m_l); + if (cur_low_bits == l_pos) { + return true; + } else if (cur_low_bits < l_pos) { + return false; + } + } + + return false; + } + + inline uint64_t select(uint64_t n) const { + return + ((m_high_bits_d1.select(m_high_bits, n) - n) << m_l) + | m_low_bits.get_bits(n * m_l, m_l); + } + + inline uint64_t rank(uint64_t pos) const { + assert(pos <= m_size); + assert(m_high_bits_d0.num_positions()); // needs rank index + if (pos == size()) { + return num_ones(); + } + + uint64_t h_rank = pos >> m_l; + uint64_t h_pos = m_high_bits_d0.select(m_high_bits, h_rank); + uint64_t rank = h_pos - h_rank; + uint64_t l_pos = pos & ((1ULL << m_l) - 1); + + while (h_pos > 0 + && m_high_bits[h_pos - 1] + && m_low_bits.get_bits((rank - 1) * m_l, m_l) >= l_pos) { + --rank; + --h_pos; + } + + return rank; + } + + inline uint64_t predecessor1(uint64_t pos) const { + return select(rank(pos + 1) - 1); + } + + inline uint64_t successor1(uint64_t pos) const { + return select(rank(pos)); + } + + + // Equivalent to select(n) - select(n - 1) (and select(0) for n = 0) + // Involves a linear search for predecessor in high bits. + // Efficient only if there are no large gaps in high bits + // XXX(ot): could make this adaptive + inline uint64_t delta(uint64_t n) const { + uint64_t high_val = m_high_bits_d1.select(m_high_bits, n); + uint64_t low_val = m_low_bits.get_bits(n * m_l, m_l); + if (n) { + return + // need a + here instead of an | for carry + ((high_val - m_high_bits.predecessor1(high_val - 1) - 1) << m_l) + + low_val - m_low_bits.get_bits((n - 1) * m_l, m_l); + } else { + return + ((high_val - n) << m_l) + | low_val; + } + } + + + // same as delta() + inline std::pair select_range(uint64_t n) const + { + assert(n + 1 < num_ones()); + uint64_t high_val_b = m_high_bits_d1.select(m_high_bits, n); + uint64_t low_val_b = m_low_bits.get_bits(n * m_l, m_l); + uint64_t high_val_e = m_high_bits.successor1(high_val_b + 1); + uint64_t low_val_e = m_low_bits.get_bits((n + 1) * m_l, m_l); + return std::make_pair(((high_val_b - n) << m_l) | low_val_b, + ((high_val_e - n - 1) << m_l) | low_val_e); + } + + struct select_enumerator { + + select_enumerator(elias_fano const& ef, uint64_t i) + : m_ef(&ef) + , m_i(i) + , m_l(ef.m_l) + { + m_low_mask = (uint64_t(1) << m_l) - 1; + m_low_buf = 0; + if (m_l) { + m_chunks_in_word = 64 / m_l; + m_chunks_avail = 0; + } else { + m_chunks_in_word = 0; + m_chunks_avail = m_ef->num_ones(); + } + + if (!m_ef->num_ones()) return; + uint64_t pos = m_ef->m_high_bits_d1.select(m_ef->m_high_bits, m_i); + m_high_enum = bit_vector::unary_enumerator(m_ef->m_high_bits, pos); + assert(m_l < 64); + } + + uint64_t next() { + if (!m_chunks_avail--) { + m_low_buf = m_ef->m_low_bits.get_word(m_i * m_l); + m_chunks_avail = m_chunks_in_word - 1; + } + + uint64_t high = m_high_enum.next(); + assert(high == m_ef->m_high_bits_d1.select(m_ef->m_high_bits, m_i)); + uint64_t low = m_low_buf & m_low_mask; + uint64_t ret = + ((high - m_i) << m_l) + | low; + m_i += 1; + m_low_buf >>= m_l; + + return ret; + } + + private: + + elias_fano const* m_ef; + uint64_t m_i; + uint64_t m_l; + bit_vector::unary_enumerator m_high_enum; + uint64_t m_low_buf; + uint64_t m_low_mask; + uint64_t m_chunks_in_word; + uint64_t m_chunks_avail; + }; + + protected: + void build(elias_fano_builder& builder, bool with_rank_index) { + m_size = builder.m_n; + m_l = builder.m_l; + bit_vector(&builder.m_high_bits).swap(m_high_bits); + darray1(m_high_bits).swap(m_high_bits_d1); + if (with_rank_index) { + darray0(m_high_bits).swap(m_high_bits_d0); + } + bit_vector(&builder.m_low_bits).swap(m_low_bits); + } + + uint64_t m_size; + bit_vector m_high_bits; + darray1 m_high_bits_d1; + darray0 m_high_bits_d0; + bit_vector m_low_bits; + uint8_t m_l; + }; + +} diff --git a/3party/succinct/elias_fano_compressed_list.hpp b/3party/succinct/elias_fano_compressed_list.hpp new file mode 100644 index 0000000000..51d1848710 --- /dev/null +++ b/3party/succinct/elias_fano_compressed_list.hpp @@ -0,0 +1,74 @@ +#pragma once + +#include "elias_fano.hpp" + +namespace succinct { + + struct elias_fano_compressed_list + { + typedef uint64_t value_type; + + elias_fano_compressed_list() {} + + template + elias_fano_compressed_list(Range const& ints) + { + typedef typename boost::range_const_iterator::type iterator_t; + + size_t s = 0; + size_t n = 0; + for (iterator_t iter = boost::begin(ints); + iter != boost::end(ints); + ++iter) { + s += broadword::msb(*iter + 1); + n += 1; + } + + elias_fano::elias_fano_builder ef_builder(s + 1, n + 1); + bit_vector_builder bits_builder; + + ef_builder.push_back(bits_builder.size()); + for (iterator_t iter = boost::begin(ints); + iter != boost::end(ints); + ++iter) { + size_t val = *iter + 1; + size_t l = broadword::msb(val); + bits_builder.append_bits(val ^ (uint64_t(1) << l), l); + ef_builder.push_back(bits_builder.size()); + } + elias_fano(&ef_builder, false).swap(m_ef); + bit_vector(&bits_builder).swap(m_bits); + } + + value_type operator[](size_t idx) const + { + std::pair r = m_ef.select_range(idx); + size_t l = r.second - r.first; + return ((uint64_t(1) << l) | m_bits.get_bits(r.first, l)) - 1; + } + + size_t size() const + { + return m_ef.num_ones() - 1; + } + + void swap(elias_fano_compressed_list& other) + { + m_ef.swap(other.m_ef); + m_bits.swap(other.m_bits); + } + + template + void map(Visitor& visit) { + visit + (m_ef, "m_ef") + (m_bits, "m_bits") + ; + } + + private: + elias_fano m_ef; + bit_vector m_bits; + }; + +} diff --git a/3party/succinct/elias_fano_list.hpp b/3party/succinct/elias_fano_list.hpp new file mode 100644 index 0000000000..e8191c806b --- /dev/null +++ b/3party/succinct/elias_fano_list.hpp @@ -0,0 +1,68 @@ +#pragma once + +#include "elias_fano.hpp" + +namespace succinct { + + struct elias_fano_list + { + typedef uint64_t value_type; + + elias_fano_list() {} + + template + elias_fano_list(Range const& ints) + { + typedef typename boost::range_const_iterator::type iterator_t; + + size_t s = 0; + size_t n = 0; + for (iterator_t iter = boost::begin(ints); + iter != boost::end(ints); + ++iter) { + s += *iter; + n += 1; + } + + elias_fano::elias_fano_builder ef_builder(s + 1, n); + size_t cur_base = 0; + for (iterator_t iter = boost::begin(ints); + iter != boost::end(ints); + ++iter) { + cur_base += *iter; + ef_builder.push_back(cur_base); + } + elias_fano(&ef_builder, false).swap(m_ef); + } + + value_type operator[](size_t idx) const + { + return m_ef.delta(idx); + } + + size_t size() const + { + return m_ef.num_ones(); + } + + size_t sum() const { + return m_ef.size() - 1; + } + + void swap(elias_fano_list& other) + { + m_ef.swap(other.m_ef); + } + + template + void map(Visitor& visit) { + visit + (m_ef, "m_ef") + ; + } + + private: + elias_fano m_ef; + }; + +} diff --git a/3party/succinct/forward_enumerator.hpp b/3party/succinct/forward_enumerator.hpp new file mode 100644 index 0000000000..a0e1ba0252 --- /dev/null +++ b/3party/succinct/forward_enumerator.hpp @@ -0,0 +1,25 @@ +#pragma once + +namespace succinct { + + template + struct forward_enumerator + { + typedef typename Container::value_type value_type; + + forward_enumerator(Container const& c, size_t idx = 0) + : m_c(&c) + , m_idx(idx) + {} + + value_type next() + { + return (*m_c)[m_idx++]; + } + + private: + Container const* m_c; + size_t m_idx; + }; + +} diff --git a/3party/succinct/gamma_bit_vector.hpp b/3party/succinct/gamma_bit_vector.hpp new file mode 100644 index 0000000000..89d5fb66cc --- /dev/null +++ b/3party/succinct/gamma_bit_vector.hpp @@ -0,0 +1,134 @@ +#pragma once + +#include "broadword.hpp" +#include "forward_enumerator.hpp" +#include "darray64.hpp" + +namespace succinct { + + // Compressed random-access vector to store unsigned integers + // using gamma codes. This implementation optimizes for integers + // whose representation is at least one bit long. It can be used, + // for example, to represent signed integers (with uniform sign + // distribution) by putting the sign in the LSB. For generic + // unsigned integers, use gamma_vector + + struct gamma_bit_vector + { + typedef uint64_t value_type; + + gamma_bit_vector() {} + + template + gamma_bit_vector(Range const& vals) + { + darray64::builder high_bits; + bit_vector_builder low_bits; + + high_bits.append1(); + + typedef typename boost::range_const_iterator::type iterator_t; + for (iterator_t iter = boost::begin(vals); + iter != boost::end(vals); + ++iter) { + const uint64_t val = *iter + 2; // increment the second bit + + uint8_t l = broadword::msb(val); + + assert(l > 0); + high_bits.append1(l - 1); + low_bits.append_bits(val ^ (uint64_t(1) << l), l); + } + + darray64(&high_bits).swap(m_high_bits); + bit_vector(&low_bits).swap(m_low_bits); + } + + value_type operator[](size_t idx) const + { + size_t pos = m_high_bits.select(idx); + size_t l; // ignored + return retrieve_value(pos, l); + } + + size_t size() const + { + return m_high_bits.num_ones() - 1; + } + + void swap(gamma_bit_vector& other) + { + m_high_bits.swap(other.m_high_bits); + m_low_bits.swap(other.m_low_bits); + } + + template + void map(Visitor& visit) { + visit + (m_high_bits, "m_high_bits") + (m_low_bits, "m_low_bits") + ; + } + + private: + + value_type retrieve_value(size_t pos, size_t& l) const + { + assert(m_high_bits.bits()[pos] == 1); + l = broadword::lsb(m_high_bits.bits().get_word(pos + 1)); + uint64_t chunk = m_low_bits.get_bits(pos, l + 1); // bit . val + uint64_t val = ((uint64_t(1) << (l + 1)) | chunk) - 2; + return val; + } + + friend struct forward_enumerator; + + darray64 m_high_bits; + bit_vector m_low_bits; + }; + + template <> + struct forward_enumerator + { + typedef gamma_bit_vector::value_type value_type; + + forward_enumerator(gamma_bit_vector const& c, size_t idx = 0) + : m_c(&c) + , m_idx(idx) + , m_pos(0) + { + if (idx < m_c->size()) { + m_pos = m_c->m_high_bits.select(idx); + m_high_bits_enumerator = + bit_vector::unary_enumerator(m_c->m_high_bits.bits(), m_pos + 1); + m_low_bits_enumerator = bit_vector::enumerator(m_c->m_low_bits, m_pos); + } + } + + void skip(size_t k) + { + // XXX actually implement this + while (k--) next(); + } + + value_type next() + { + assert(m_idx <= m_c->size()); + size_t next_pos = m_high_bits_enumerator.next(); + size_t l = next_pos - m_pos - 1; + m_pos = next_pos; + uint64_t chunk = m_low_bits_enumerator.take(l + 1); + uint64_t val = (chunk | (uint64_t(1) << (l + 1))) - 2; + m_idx += 1; + return val; + } + + private: + gamma_bit_vector const* m_c; + size_t m_idx; + size_t m_pos; + + bit_vector::unary_enumerator m_high_bits_enumerator; + bit_vector::enumerator m_low_bits_enumerator; + }; +} diff --git a/3party/succinct/gamma_vector.hpp b/3party/succinct/gamma_vector.hpp new file mode 100644 index 0000000000..9da81149b8 --- /dev/null +++ b/3party/succinct/gamma_vector.hpp @@ -0,0 +1,119 @@ +#pragma once + +#include "broadword.hpp" +#include "forward_enumerator.hpp" +#include "darray64.hpp" + +namespace succinct { + + // Compressed random-access vector to store unsigned integers + // using gamma codes. + struct gamma_vector + { + typedef uint64_t value_type; + + gamma_vector() {} + + template + gamma_vector(Range const& ints) + { + darray64::builder high_bits; + bit_vector_builder low_bits; + + high_bits.append1(); + + typedef typename boost::range_const_iterator::type iterator_t; + for (iterator_t iter = boost::begin(ints); + iter != boost::end(ints); + ++iter) { + const value_type val = *iter + 1; + + uint8_t l = broadword::msb(val); + + low_bits.append_bits(val ^ (uint64_t(1) << l), l); + high_bits.append1(l); + } + + darray64(&high_bits).swap(m_high_bits); + bit_vector(&low_bits).swap(m_low_bits); + } + + value_type operator[](size_t idx) const + { + size_t pos = m_high_bits.select(idx); + size_t l; // ignored + return retrieve_value(idx, pos, l); + } + + size_t size() const + { + return m_high_bits.num_ones() - 1; + } + + void swap(gamma_vector& other) + { + m_high_bits.swap(other.m_high_bits); + m_low_bits.swap(other.m_low_bits); + } + + template + void map(Visitor& visit) { + visit + (m_high_bits, "m_high_bits") + (m_low_bits, "m_low_bits") + ; + } + + private: + + value_type retrieve_value(size_t idx, size_t pos, size_t& l) const + { + assert(m_high_bits.bits()[pos] == 1); + l = broadword::lsb(m_high_bits.bits().get_word(pos + 1)); + return ((uint64_t(1) << l) | m_low_bits.get_bits(pos - idx, l)) - 1; + } + + friend struct forward_enumerator; + darray64 m_high_bits; + bit_vector m_low_bits; + }; + + template <> + struct forward_enumerator + { + typedef gamma_vector::value_type value_type; + + forward_enumerator(gamma_vector const& c, size_t idx = 0) + : m_c(&c) + , m_idx(idx) + , m_pos(0) + { + if (idx < m_c->size()) { + m_pos = m_c->m_high_bits.select(idx); + m_high_bits_enumerator = + bit_vector::unary_enumerator(m_c->m_high_bits.bits(), m_pos + 1); + m_low_bits_enumerator = bit_vector::enumerator(m_c->m_low_bits, m_pos - m_idx); + } + } + + value_type next() + { + assert(m_idx <= m_c->size()); + size_t next_pos = m_high_bits_enumerator.next(); + size_t l = next_pos - m_pos - 1; + m_pos = next_pos; + uint64_t chunk = m_low_bits_enumerator.take(l); + uint64_t val = (chunk | (uint64_t(1) << (l))) - 1; + m_idx += 1; + return val; + } + + private: + gamma_vector const* m_c; + size_t m_idx; + size_t m_pos; + + bit_vector::unary_enumerator m_high_bits_enumerator; + bit_vector::enumerator m_low_bits_enumerator; + }; +} diff --git a/3party/succinct/intrinsics.hpp b/3party/succinct/intrinsics.hpp new file mode 100644 index 0000000000..90f731697d --- /dev/null +++ b/3party/succinct/intrinsics.hpp @@ -0,0 +1,101 @@ +#pragma once + +#include +#include "succinct_config.hpp" + +#if SUCCINCT_USE_INTRINSICS +#include + +#if defined(__GNUC__) || defined(__clang__) +# define __INTRIN_INLINE inline __attribute__((__always_inline__)) +#elif defined(_MSC_VER) +# define __INTRIN_INLINE inline __forceinline +#else +# define __INTRIN_INLINE inline +#endif + +#endif + +#if SUCCINCT_USE_POPCNT +# if !SUCCINCT_USE_INTRINSICS +# error "Intrinsics support needed for popcnt" +# endif +#include +#endif + + + +namespace succinct { namespace intrinsics { + + +#if SUCCINCT_USE_INTRINSICS + + __INTRIN_INLINE uint64_t byteswap64(uint64_t value) + { +#if defined(__GNUC__) || defined(__clang__) + return __builtin_bswap64(value); +#elif defined(_MSC_VER) + return _byteswap_uint64(value); +#else +# error Unsupported platform +#endif + } + + __INTRIN_INLINE bool bsf64(unsigned long* const index, const uint64_t mask) + { +#if defined(__GNUC__) || defined(__clang__) + if (mask) { + *index = (unsigned long)__builtin_ctzll(mask); + return true; + } else { + return false; + } +#elif defined(_MSC_VER) + return _BitScanForward64(index, mask) != 0; +#else +# error Unsupported platform +#endif + } + + __INTRIN_INLINE bool bsr64(unsigned long* const index, const uint64_t mask) + { +#if defined(__GNUC__) || defined(__clang__) + if (mask) { + *index = (unsigned long)(63 - __builtin_clzll(mask)); + return true; + } else { + return false; + } +#elif defined(_MSC_VER) + return _BitScanReverse64(index, mask) != 0; +#else +# error Unsupported platform +#endif + } + + template + __INTRIN_INLINE void prefetch(T const* ptr) + { + _mm_prefetch((const char*)ptr, _MM_HINT_T0); + } + +#else /* SUCCINCT_USE_INTRINSICS */ + + template + inline void prefetch(T const* /* ptr */) + { + /* do nothing */ + } + +#endif /* SUCCINCT_USE_INTRINSICS */ + +#if SUCCINCT_USE_POPCNT + + __INTRIN_INLINE uint64_t popcount(uint64_t x) + { + return uint64_t(_mm_popcnt_u64(x)); + } + +#endif /* SUCCINCT_USE_POPCNT */ + +}} diff --git a/3party/succinct/mappable_vector.hpp b/3party/succinct/mappable_vector.hpp new file mode 100644 index 0000000000..0c1c9a76ea --- /dev/null +++ b/3party/succinct/mappable_vector.hpp @@ -0,0 +1,124 @@ +#pragma once + +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include "intrinsics.hpp" + +namespace succinct { namespace mapper { + + namespace detail { + class freeze_visitor; + class map_visitor; + class sizeof_visitor; + } + + typedef boost::function deleter_t; + + template // T must be a POD + class mappable_vector : boost::noncopyable { + public: + typedef T value_type; + typedef const T* iterator; + typedef const T* const_iterator; + + mappable_vector() + : m_data(0) + , m_size(0) + , m_deleter() + {} + + template + mappable_vector(Range const& from) + : m_data(0) + , m_size(0) + { + size_t size = boost::size(from); + T* data = new T[size]; + m_deleter = boost::lambda::bind(boost::lambda::delete_array(), data); + + std::copy(boost::begin(from), + boost::end(from), + data); + m_data = data; + m_size = size; + } + + ~mappable_vector() { + if (m_deleter) { + m_deleter(); + } + } + + void swap(mappable_vector& other) { + using std::swap; + swap(m_data, other.m_data); + swap(m_size, other.m_size); + swap(m_deleter, other.m_deleter); + } + + void clear() { + mappable_vector().swap(*this); + } + + void steal(std::vector& vec) { + clear(); + m_size = vec.size(); + if (m_size) { + std::vector* new_vec = new std::vector; + new_vec->swap(vec); + m_deleter = boost::lambda::bind(boost::lambda::delete_ptr(), new_vec); + m_data = &(*new_vec)[0]; + } + } + + template + void assign(Range const& from) { + clear(); + mappable_vector(from).swap(*this); + } + + uint64_t size() const { + return m_size; + } + + inline const_iterator begin() const { + return m_data; + } + + inline const_iterator end() const { + return m_data + m_size; + } + + inline T const& operator[](uint64_t i) const { + assert(i < m_size); + return m_data[i]; + } + + inline T const* data() const { + return m_data; + } + + inline void prefetch(size_t i) const { + succinct::intrinsics::prefetch(m_data + i); + } + + friend class detail::freeze_visitor; + friend class detail::map_visitor; + friend class detail::sizeof_visitor; + + protected: + const T* m_data; + uint64_t m_size; + deleter_t m_deleter; + }; + +}} diff --git a/3party/succinct/mapper.hpp b/3party/succinct/mapper.hpp new file mode 100644 index 0000000000..381e1992ef --- /dev/null +++ b/3party/succinct/mapper.hpp @@ -0,0 +1,280 @@ +#pragma once + +#include +#include + +#include +#include +#include +#include +#include + +#include "mappable_vector.hpp" + +namespace succinct { namespace mapper { + + struct freeze_flags { + // enum { + // }; + }; + + struct map_flags { + enum { + warmup = 1 + }; + }; + + struct size_node; + typedef boost::shared_ptr size_node_ptr; + + struct size_node + { + size_node() + : size(0) + {} + + std::string name; + size_t size; + std::vector children; + + void dump(std::ostream& os = std::cerr, size_t depth = 0) { + os << std::string(depth * 4, ' ') + << name << ": " + << size << '\n'; + for (size_t i = 0; i < children.size(); ++i) { + children[i]->dump(os, depth + 1); + } + } + }; + + namespace detail { + class freeze_visitor : boost::noncopyable { + public: + freeze_visitor(std::ofstream& fout, uint64_t flags) + : m_fout(fout) + , m_flags(flags) + , m_written(0) + { + // Save freezing flags + m_fout.write(reinterpret_cast(&m_flags), sizeof(m_flags)); + m_written += sizeof(m_flags); + } + + template + typename boost::disable_if, freeze_visitor&>::type + operator()(T& val, const char* /* friendly_name */) { + val.map(*this); + return *this; + } + + template + typename boost::enable_if, freeze_visitor&>::type + operator()(T& val, const char* /* friendly_name */) { + m_fout.write(reinterpret_cast(&val), sizeof(T)); + m_written += sizeof(T); + return *this; + } + + template + freeze_visitor& + operator()(mappable_vector& vec, const char* /* friendly_name */) { + (*this)(vec.m_size, "size"); + + size_t n_bytes = static_cast(vec.m_size * sizeof(T)); + m_fout.write(reinterpret_cast(vec.m_data), long(n_bytes)); + m_written += n_bytes; + + return *this; + } + + size_t written() const { + return m_written; + } + + protected: + std::ofstream& m_fout; + const uint64_t m_flags; + uint64_t m_written; + }; + + class map_visitor : boost::noncopyable { + public: + map_visitor(const char* base_address, uint64_t flags) + : m_base(base_address) + , m_cur(m_base) + , m_flags(flags) + { + m_freeze_flags = *reinterpret_cast(m_cur); + m_cur += sizeof(m_freeze_flags); + } + + template + typename boost::disable_if, map_visitor&>::type + operator()(T& val, const char* /* friendly_name */) { + val.map(*this); + return *this; + } + + template + typename boost::enable_if, map_visitor&>::type + operator()(T& val, const char* /* friendly_name */) { + val = *reinterpret_cast(m_cur); + m_cur += sizeof(T); + return *this; + } + + template + map_visitor& + operator()(mappable_vector& vec, const char* /* friendly_name */) { + vec.clear(); + (*this)(vec.m_size, "size"); + + vec.m_data = reinterpret_cast(m_cur); + size_t bytes = vec.m_size * sizeof(T); + + if (m_flags & map_flags::warmup) { + T foo; + volatile T* bar = &foo; + for (size_t i = 0; i < vec.m_size; ++i) { + *bar = vec.m_data[i]; + } + } + + m_cur += bytes; + return *this; + } + + size_t bytes_read() const { + return size_t(m_cur - m_base); + } + + protected: + const char* m_base; + const char* m_cur; + const uint64_t m_flags; + uint64_t m_freeze_flags; + }; + + class sizeof_visitor : boost::noncopyable { + public: + sizeof_visitor(bool with_tree = false) + : m_size(0) + { + if (with_tree) { + m_cur_size_node = boost::make_shared(); + } + } + + template + typename boost::disable_if, sizeof_visitor&>::type + operator()(T& val, const char* friendly_name) { + size_t checkpoint = m_size; + size_node_ptr parent_node; + if (m_cur_size_node) { + parent_node = m_cur_size_node; + m_cur_size_node = make_node(friendly_name); + } + + val.map(*this); + + if (m_cur_size_node) { + m_cur_size_node->size = m_size - checkpoint; + m_cur_size_node = parent_node; + } + return *this; + } + + template + typename boost::enable_if, sizeof_visitor&>::type + operator()(T& /* val */, const char* /* friendly_name */) { + // don't track PODs in the size tree (they are constant sized) + m_size += sizeof(T); + return *this; + } + + template + sizeof_visitor& + operator()(mappable_vector& vec, const char* friendly_name) { + size_t checkpoint = m_size; + (*this)(vec.m_size, "size"); + m_size += static_cast(vec.m_size * sizeof(T)); + + if (m_cur_size_node) { + make_node(friendly_name)->size = m_size - checkpoint; + } + + return *this; + } + + size_t size() const { + return m_size; + } + + size_node_ptr size_tree() const { + assert(m_cur_size_node); + return m_cur_size_node; + } + + protected: + + size_node_ptr make_node(const char* name) + { + size_node_ptr node = boost::make_shared(); + m_cur_size_node->children.push_back(node); + node->name = name; + return node; + } + + size_t m_size; + size_node_ptr m_cur_size_node; + }; + + } + + template + size_t freeze(T& val, std::ofstream& fout, uint64_t flags = 0, const char* friendly_name = "") + { + detail::freeze_visitor freezer(fout, flags); + freezer(val, friendly_name); + return freezer.written(); + } + + template + size_t freeze(T& val, const char* filename, uint64_t flags = 0, const char* friendly_name = "") + { + std::ofstream fout(filename, std::ios::binary); + return freeze(val, fout, flags, friendly_name); + } + + template + size_t map(T& val, const char* base_address, uint64_t flags = 0, const char* friendly_name = "") + { + detail::map_visitor mapper(base_address, flags); + mapper(val, friendly_name); + return mapper.bytes_read(); + } + + template + size_t map(T& val, boost::iostreams::mapped_file_source const& m, uint64_t flags = 0, const char* friendly_name = "") + { + return map(val, m.data(), flags, friendly_name); + } + + template + size_t size_of(T& val) + { + detail::sizeof_visitor sizer; + sizer(val, ""); + return sizer.size(); + } + + template + size_node_ptr size_tree_of(T& val, const char* friendly_name = "") + { + detail::sizeof_visitor sizer(true); + sizer(val, friendly_name); + assert(sizer.size_tree()->children.size()); + return sizer.size_tree()->children[0]; + } + +}} diff --git a/3party/succinct/nibble_vector.hpp b/3party/succinct/nibble_vector.hpp new file mode 100644 index 0000000000..81758b7ab2 --- /dev/null +++ b/3party/succinct/nibble_vector.hpp @@ -0,0 +1,68 @@ +#pragma once + +#include + +#include + +#include "succinct/mappable_vector.hpp" + +namespace succinct { + + class nibble_vector { + public: + nibble_vector() + : m_size(0) + {} + + template + nibble_vector(Range const& from) + : m_size(0) + { + std::vector nibbles; + bool parity = 0; + uint8_t cur_byte = 0; + for (typename boost::range_const_iterator::type iter = boost::begin(from); + iter != boost::end(from); + ++iter) { + assert(*iter < 16); + cur_byte |= *iter << (parity * 4); + parity = !parity; + if (!parity) { + nibbles.push_back(cur_byte); + cur_byte = 0; + } + ++m_size; + } + if (parity) { + nibbles.push_back(cur_byte); + } + m_nibbles.steal(nibbles); + } + + template + void map(Visitor& visit) { + visit + (m_size, "m_size") + (m_nibbles, "m_nibbles"); + } + + void swap(nibble_vector& other) { + std::swap(other.m_size, m_size); + other.m_nibbles.swap(m_nibbles); + } + + size_t size() const { + return m_size; + } + + uint8_t operator[](uint64_t pos) const { + assert(pos < m_size); + return (m_nibbles[pos / 2] >> ((pos % 2) * 4)) & 0x0F; + } + + protected: + size_t m_size; + mapper::mappable_vector m_nibbles; + }; + +} diff --git a/3party/succinct/perftest/.gitignore b/3party/succinct/perftest/.gitignore new file mode 100644 index 0000000000..47a0e27a29 --- /dev/null +++ b/3party/succinct/perftest/.gitignore @@ -0,0 +1,4 @@ +perftest_bp_vector +perftest_bp_vector_rmq +perftest_cartesian_tree +perftest_elias_fano diff --git a/3party/succinct/perftest/CMakeLists.txt b/3party/succinct/perftest/CMakeLists.txt new file mode 100644 index 0000000000..8240a1bca3 --- /dev/null +++ b/3party/succinct/perftest/CMakeLists.txt @@ -0,0 +1,9 @@ +file(GLOB SUCCINCT_TEST_SOURCES perftest_*.cpp) +foreach(TEST_SRC ${SUCCINCT_TEST_SOURCES}) + get_filename_component (TEST_SRC_NAME ${TEST_SRC} NAME_WE) + add_executable(${TEST_SRC_NAME} ${TEST_SRC}) + target_link_libraries(${TEST_SRC_NAME} + succinct + ${Boost_LIBRARIES} + ) +endforeach(TEST_SRC) diff --git a/3party/succinct/perftest/perftest_bp_vector.cpp b/3party/succinct/perftest/perftest_bp_vector.cpp new file mode 100644 index 0000000000..da56aee540 --- /dev/null +++ b/3party/succinct/perftest/perftest_bp_vector.cpp @@ -0,0 +1,119 @@ +#include +#include + +#include + +#include "util.hpp" +#include "test_bp_vector_common.hpp" + +#include "bp_vector.hpp" +#include "mapper.hpp" + +#include "perftest_common.hpp" + +// this generic trait enables easy comparisons with other BP +// implementations + +struct succinct_bp_vector_traits +{ + typedef succinct::bit_vector_builder builder_type; + typedef succinct::bp_vector bp_vector_type; + + static inline void build(builder_type& builder, bp_vector_type& bp) + { + bp_vector_type(&builder, true, false).swap(bp); + } + + static inline std::string log_header() + { + return std::string("SUCCINCT"); + } + + static inline double bits_per_bp(bp_vector_type& vec) + { + return double(succinct::mapper::size_of(vec)) + * 8.0 / double(vec.size()); + } + +}; + +template +double time_visit(BpVector const& bp, size_t sample_size = 1000000) +{ + std::vector random_bits; + for (size_t i = 0; i < sample_size; ++i) { + random_bits.push_back(rand() > (RAND_MAX / 2)); + } + + volatile size_t foo = 0; // to prevent the compiler to optimize away the loop + + size_t find_close_performed = 0; + size_t steps_done = 0; + double elapsed; + SUCCINCT_TIMEIT(elapsed) { + while (steps_done < sample_size) { + size_t cur_node = 1; // root + + while (bp[cur_node] && steps_done < sample_size) { + if (random_bits[steps_done++]) { + size_t next_node = bp.find_close(cur_node); + cur_node = next_node + 1; + find_close_performed += 1; + } else { + cur_node += 1; + } + } + foo = cur_node; + } + } + + (void)foo; // silence warning + return elapsed / double(find_close_performed); +} + +template +void build_random_binary_tree(typename BpVectorTraits::bp_vector_type& bp, size_t size) +{ + typename BpVectorTraits::builder_type builder; + succinct::random_binary_tree(builder, size); + BpVectorTraits::build(builder, bp); +} + +template +void bp_benchmark(size_t runs) +{ + srand(42); // make everything deterministic + static const size_t sample_size = 10000000; + + std::cout << BpVectorTraits::log_header() << std::endl; + std::cout << "log_height" "\t" "find_close_us" "\t" "bits_per_bp" << std::endl; + + for (size_t ln = 10; ln <= 28; ln += 2) { + size_t n = 1 << ln; + double elapsed = 0; + double bits_per_bp = 0; + for (size_t run = 0; run < runs; ++run) { + typename BpVectorTraits::bp_vector_type bp; + build_random_binary_tree(bp, n); + elapsed += time_visit(bp, sample_size); + bits_per_bp += BpVectorTraits::bits_per_bp(bp); + } + std::cout << ln + << "\t" << elapsed / double(runs) + << "\t" << bits_per_bp / double(runs) + << std::endl; + } +} + +int main(int argc, char** argv) +{ + size_t runs; + + if (argc == 2) { + runs = boost::lexical_cast(argv[1]); + } else { + runs = 1; + } + + bp_benchmark(runs); +} diff --git a/3party/succinct/perftest/perftest_bp_vector_rmq.cpp b/3party/succinct/perftest/perftest_bp_vector_rmq.cpp new file mode 100644 index 0000000000..49ec80d0a0 --- /dev/null +++ b/3party/succinct/perftest/perftest_bp_vector_rmq.cpp @@ -0,0 +1,77 @@ +#include +#include + +#include + +#include "util.hpp" +#include "test_bp_vector_common.hpp" + +#include "bp_vector.hpp" + +#include "perftest_common.hpp" + +double time_avg_rmq(succinct::bp_vector const& bp, size_t sample_size = 1000000) +{ + typedef std::pair range_pair; + std::vector pairs_sample; + for (size_t i = 0; i < sample_size; ++i) { + uint64_t a = uint64_t(rand()) % bp.size(); + uint64_t b = a + (uint64_t(rand()) % (bp.size() - a)); + pairs_sample.push_back(range_pair(a, b)); + } + + volatile uint64_t foo; // to prevent the compiler to optimize away the loop + + size_t rmq_performed = 0; + double elapsed; + SUCCINCT_TIMEIT(elapsed) { + for (size_t i = 0; i < pairs_sample.size(); ++i) { + range_pair r = pairs_sample[i]; + foo = bp.excess_rmq(r.first, r.second); + rmq_performed += 1; + } + } + + (void)foo; // silence warning + return elapsed / double(rmq_performed); +} + +void build_random_binary_tree(succinct::bp_vector& bp, size_t size) +{ + succinct::bit_vector_builder builder; + succinct::random_binary_tree(builder, size); + succinct::bp_vector(&builder, true, false).swap(bp); +} + +void rmq_benchmark(size_t runs) +{ + srand(42); // make everything deterministic + static const size_t sample_size = 10000000; + + std::cout << "SUCCINCT_EXCESS_RMQ" << std::endl; + std::cout << "log_height" "\t" "excess_rmq_us" << std::endl; + + for (size_t ln = 10; ln <= 28; ln += 2) { + size_t n = 1 << ln; + double elapsed = 0; + for (size_t run = 0; run < runs; ++run) { + succinct::bp_vector bp; + build_random_binary_tree(bp, n); + elapsed += time_avg_rmq(bp, sample_size); + } + std::cout << ln << "\t" << elapsed / double(runs) << std::endl; + } +} + +int main(int argc, char** argv) +{ + size_t runs; + + if (argc == 2) { + runs = boost::lexical_cast(argv[1]); + } else { + runs = 1; + } + + rmq_benchmark(runs); +} diff --git a/3party/succinct/perftest/perftest_cartesian_tree.cpp b/3party/succinct/perftest/perftest_cartesian_tree.cpp new file mode 100644 index 0000000000..0b9f5bfd77 --- /dev/null +++ b/3party/succinct/perftest/perftest_cartesian_tree.cpp @@ -0,0 +1,74 @@ +#include +#include + +#include + +#include "util.hpp" +#include "test_bp_vector_common.hpp" + +#include "cartesian_tree.hpp" + +#include "perftest_common.hpp" + +double time_avg_rmq(succinct::cartesian_tree const& tree, size_t sample_size = 1000000) +{ + typedef std::pair range_pair; + std::vector pairs_sample; + for (size_t i = 0; i < sample_size; ++i) { + uint64_t a = uint64_t(rand()) % tree.size(); + uint64_t b = a + (uint64_t(rand()) % (tree.size() - a)); + pairs_sample.push_back(range_pair(a, b)); + } + + volatile uint64_t foo; // to prevent the compiler to optimize away the loop + + size_t rmq_performed = 0; + double elapsed; + SUCCINCT_TIMEIT(elapsed) { + for (size_t i = 0; i < pairs_sample.size(); ++i) { + range_pair r = pairs_sample[i]; + foo = tree.rmq(r.first, r.second); + rmq_performed += 1; + } + } + + (void)foo; // silence warning + return elapsed / double(rmq_performed); +} + +void rmq_benchmark(size_t runs) +{ + srand(42); // make everything deterministic + static const size_t sample_size = 10000000; + + std::cout << "SUCCINCT_CARTESIAN_TREE_RMQ" << std::endl; + std::cout << "log_height" "\t" "excess_rmq_us" << std::endl; + + for (size_t ln = 10; ln <= 28; ln += 2) { + size_t n = 1 << ln; + double elapsed = 0; + for (size_t run = 0; run < runs; ++run) { + std::vector v(n); + for (size_t i = 0; i < v.size(); ++i) { + v[i] = uint64_t(rand()) % 1024; + } + + succinct::cartesian_tree tree(v); + elapsed += time_avg_rmq(tree, sample_size); + } + std::cout << ln << "\t" << elapsed / double(runs) << std::endl; + } +} + +int main(int argc, char** argv) +{ + size_t runs; + + if (argc == 2) { + runs = boost::lexical_cast(argv[1]); + } else { + runs = 1; + } + + rmq_benchmark(runs); +} diff --git a/3party/succinct/perftest/perftest_common.hpp b/3party/succinct/perftest/perftest_common.hpp new file mode 100644 index 0000000000..f07b5d8bd0 --- /dev/null +++ b/3party/succinct/perftest/perftest_common.hpp @@ -0,0 +1,33 @@ +#pragma once + +#include + +namespace succinct { +namespace detail { + + struct timer { + timer() + : m_tick(boost::posix_time::microsec_clock::universal_time()) + , m_done(false) + {} + + bool done() { return m_done; } + + void report(double& elapsed) { + elapsed = (double)(boost::posix_time::microsec_clock::universal_time() - m_tick).total_microseconds(); + m_done = true; + } + + const std::string m_msg; + boost::posix_time::ptime m_tick; + bool m_done; + }; + +} +} + +#define SUCCINCT_TIMEIT(elapsed) \ + for (::succinct::detail::timer SUCCINCT_TIMEIT_timer; \ + !SUCCINCT_TIMEIT_timer.done(); \ + SUCCINCT_TIMEIT_timer.report(elapsed)) \ + /**/ diff --git a/3party/succinct/perftest/perftest_elias_fano.cpp b/3party/succinct/perftest/perftest_elias_fano.cpp new file mode 100644 index 0000000000..d892610b50 --- /dev/null +++ b/3party/succinct/perftest/perftest_elias_fano.cpp @@ -0,0 +1,114 @@ +#include +#include + +#include +#include +#include +#include + +#include "util.hpp" + +#include "elias_fano.hpp" +#include "mapper.hpp" + +#include "perftest_common.hpp" + +struct monotone_generator +{ + monotone_generator(uint64_t m, uint8_t bits, unsigned int seed) + : m_gen(seed) + , m_bits(bits) + { + m_stack.push_back(state_t(0, m, 0)); + } + + uint64_t next() + { + uint64_t cur_word, cur_m; + uint8_t cur_depth; + + assert(m_stack.size()); + boost::tie(cur_word, cur_m, cur_depth) = m_stack.back(); + m_stack.pop_back(); + + while (cur_depth < m_bits) { + boost::random::uniform_int_distribution dist(0, cur_m); + uint64_t left_m = dist(m_gen); + uint64_t right_m = cur_m - left_m; + + // push left and right children, if present + if (right_m > 0) { + m_stack.push_back(state_t(cur_word | (uint64_t(1) << (m_bits - cur_depth - 1)), + right_m, cur_depth + 1)); + } + if (left_m > 0) { + m_stack.push_back(state_t(cur_word, left_m, cur_depth + 1)); + + } + + // pop next child in visit + boost::tie(cur_word, cur_m, cur_depth) = m_stack.back(); + m_stack.pop_back(); + } + + if (cur_m > 1) { + // push back the current leaf, with cur_m decreased by one + m_stack.push_back(state_t(cur_word, cur_m - 1, cur_depth)); + } + + return cur_word; + } + + bool done() const + { + return m_stack.empty(); + } + +private: + typedef boost::tuple state_t; + std::vector m_stack; + boost::random::mt19937 m_gen; + uint8_t m_bits; +}; + +void ef_enumeration_benchmark(uint64_t m, uint8_t bits) +{ + succinct::elias_fano::elias_fano_builder bvb(uint64_t(1) << bits, m); + monotone_generator mgen(m, bits, 37); + for (size_t i = 0; i < m; ++i) { + bvb.push_back(mgen.next()); + } + assert(mgen.done()); + + succinct::elias_fano ef(&bvb); + succinct::mapper::size_tree_of(ef)->dump(); + + + double elapsed; + uint64_t foo = 0; + SUCCINCT_TIMEIT(elapsed) { + succinct::elias_fano::select_enumerator it(ef, 0); + for (size_t i = 0; i < m; ++i) { + foo ^= it.next(); + } + } + volatile uint64_t vfoo = foo; + (void)vfoo; // silence warning + + std::cerr << "Elapsed: " << elapsed / 1000 << " msec\n" + << double(m) / elapsed << " Mcodes/s" << std::endl; +} + +int main(int argc, char** argv) +{ + if (argc != 3) { + std::cerr << "Invalid arguments" << std::endl; + std::terminate(); + } + size_t m = boost::lexical_cast(argv[1]); + uint8_t bits = uint8_t(boost::lexical_cast(argv[2])); + + ef_enumeration_benchmark(m, bits); +} diff --git a/3party/succinct/rs_bit_vector.cpp b/3party/succinct/rs_bit_vector.cpp new file mode 100644 index 0000000000..6b47017ccd --- /dev/null +++ b/3party/succinct/rs_bit_vector.cpp @@ -0,0 +1,73 @@ +#include "rs_bit_vector.hpp" + +namespace succinct { + + void rs_bit_vector::build_indices(bool with_select_hints, bool with_select0_hints) + { + { + using broadword::popcount; + std::vector block_rank_pairs; + uint64_t next_rank = 0; + uint64_t cur_subrank = 0; + uint64_t subranks = 0; + block_rank_pairs.push_back(0); + for (uint64_t i = 0; i < m_bits.size(); ++i) { + uint64_t word_pop = popcount(m_bits[i]); + uint64_t shift = i % block_size; + if (shift) { + subranks <<= 9; + subranks |= cur_subrank; + } + next_rank += word_pop; + cur_subrank += word_pop; + + if (shift == block_size - 1) { + block_rank_pairs.push_back(subranks); + block_rank_pairs.push_back(next_rank); + subranks = 0; + cur_subrank = 0; + } + } + uint64_t left = block_size - m_bits.size() % block_size; + for (uint64_t i = 0; i < left; ++i) { + subranks <<= 9; + subranks |= cur_subrank; + } + block_rank_pairs.push_back(subranks); + + if (m_bits.size() % block_size) { + block_rank_pairs.push_back(next_rank); + block_rank_pairs.push_back(0); + } + + m_block_rank_pairs.steal(block_rank_pairs); + } + + if (with_select_hints) { + std::vector select_hints; + uint64_t cur_ones_threshold = select_ones_per_hint; + for (uint64_t i = 0; i < num_blocks(); ++i) { + if (block_rank(i + 1) > cur_ones_threshold) { + select_hints.push_back(i); + cur_ones_threshold += select_ones_per_hint; + } + } + select_hints.push_back(num_blocks()); + m_select_hints.steal(select_hints); + } + + if (with_select0_hints) { + std::vector select0_hints; + uint64_t cur_zeros_threshold = select_zeros_per_hint; + for (uint64_t i = 0; i < num_blocks(); ++i) { + if (block_rank0(i + 1) > cur_zeros_threshold) { + select0_hints.push_back(i); + cur_zeros_threshold += select_zeros_per_hint; + } + } + select0_hints.push_back(num_blocks()); + m_select0_hints.steal(select0_hints); + } + } + +} diff --git a/3party/succinct/rs_bit_vector.hpp b/3party/succinct/rs_bit_vector.hpp new file mode 100644 index 0000000000..1197ffb100 --- /dev/null +++ b/3party/succinct/rs_bit_vector.hpp @@ -0,0 +1,192 @@ +#pragma once + +#include +#include + +#include "bit_vector.hpp" +#include "broadword.hpp" + +namespace succinct { + + class rs_bit_vector : public bit_vector { + public: + rs_bit_vector() + : bit_vector() + {} + + template + rs_bit_vector(Range const& from, + bool with_select_hints = false, + bool with_select0_hints = false) + : bit_vector(from) + { + build_indices(with_select_hints, with_select0_hints); + } + + template + void map(Visitor& visit) { + bit_vector::map(visit); + visit + (m_block_rank_pairs, "m_block_rank_pairs") + (m_select_hints, "m_select_hints") + (m_select0_hints, "m_select0_hints") + ; + } + + void swap(rs_bit_vector& other) { + bit_vector::swap(other); + m_block_rank_pairs.swap(other.m_block_rank_pairs); + m_select_hints.swap(other.m_select_hints); + m_select0_hints.swap(other.m_select0_hints); + } + + inline uint64_t num_ones() const { + return *(m_block_rank_pairs.end() - 2); + } + + inline uint64_t num_zeros() const { + return size() - num_ones(); + } + + inline uint64_t rank(uint64_t pos) const { + assert(pos <= size()); + if (pos == size()) { + return num_ones(); + } + + uint64_t sub_block = pos / 64; + uint64_t r = sub_block_rank(sub_block); + uint64_t sub_left = pos % 64; + if (sub_left) { + r += broadword::popcount(m_bits[sub_block] << (64 - sub_left)); + } + return r; + } + + inline uint64_t rank0(uint64_t pos) const { + return pos - rank(pos); + } + + inline uint64_t select(uint64_t n) const { + using broadword::popcount; + using broadword::select_in_word; + assert(n < num_ones()); + uint64_t a = 0; + uint64_t b = num_blocks(); + if (m_select_hints.size()) { + uint64_t chunk = n / select_ones_per_hint; + if (chunk != 0) { + a = m_select_hints[chunk - 1]; + } + b = m_select_hints[chunk] + 1; + } + + uint64_t block = 0; + while (b - a > 1) { + uint64_t mid = a + (b - a) / 2; + uint64_t x = block_rank(mid); + if (x <= n) { + a = mid; + } else { + b = mid; + } + } + block = a; + + assert(block < num_blocks()); + uint64_t block_offset = block * block_size; + uint64_t cur_rank = block_rank(block); + assert(cur_rank <= n); + + + uint64_t rank_in_block_parallel = (n - cur_rank) * broadword::ones_step_9; + uint64_t sub_ranks = sub_block_ranks(block); + uint64_t sub_block_offset = broadword::uleq_step_9(sub_ranks, rank_in_block_parallel) * broadword::ones_step_9 >> 54 & 0x7; + cur_rank += sub_ranks >> (7 - sub_block_offset) * 9 & 0x1FF; + assert(cur_rank <= n); + + uint64_t word_offset = block_offset + sub_block_offset; + return word_offset * 64 + select_in_word(m_bits[word_offset], n - cur_rank); + } + + // TODO(ot): share code between select and select0 + inline uint64_t select0(uint64_t n) const { + using broadword::popcount; + using broadword::select_in_word; + assert(n < num_zeros()); + uint64_t a = 0; + uint64_t b = num_blocks(); + if (m_select0_hints.size()) { + uint64_t chunk = n / select_zeros_per_hint; + if (chunk != 0) { + a = m_select0_hints[chunk - 1]; + } + b = m_select0_hints[chunk] + 1; + } + + uint64_t block = 0; + while (b - a > 1) { + uint64_t mid = a + (b - a) / 2; + uint64_t x = block_rank0(mid); + if (x <= n) { + a = mid; + } else { + b = mid; + } + } + block = a; + + assert(block < num_blocks()); + uint64_t block_offset = block * block_size; + uint64_t cur_rank0 = block_rank0(block); + assert(cur_rank0 <= n); + + uint64_t rank_in_block_parallel = (n - cur_rank0) * broadword::ones_step_9; + uint64_t sub_ranks = 64 * broadword::inv_count_step_9 - sub_block_ranks(block); + uint64_t sub_block_offset = broadword::uleq_step_9(sub_ranks, rank_in_block_parallel) * broadword::ones_step_9 >> 54 & 0x7; + cur_rank0 += sub_ranks >> (7 - sub_block_offset) * 9 & 0x1FF; + assert(cur_rank0 <= n); + + uint64_t word_offset = block_offset + sub_block_offset; + return word_offset * 64 + select_in_word(~m_bits[word_offset], n - cur_rank0); + } + + protected: + + inline uint64_t num_blocks() const { + return m_block_rank_pairs.size() / 2 - 1; + } + + inline uint64_t block_rank(uint64_t block) const { + return m_block_rank_pairs[block * 2]; + } + + inline uint64_t sub_block_rank(uint64_t sub_block) const { + uint64_t r = 0; + uint64_t block = sub_block / block_size; + r += block_rank(block); + uint64_t left = sub_block % block_size; + r += sub_block_ranks(block) >> ((7 - left) * 9) & 0x1FF; + return r; + } + + inline uint64_t sub_block_ranks(uint64_t block) const { + return m_block_rank_pairs[block * 2 + 1]; + } + + inline uint64_t block_rank0(uint64_t block) const { + return block * block_size * 64 - m_block_rank_pairs[block * 2]; + } + + void build_indices(bool with_select_hints, bool with_select0_hints); + + static const uint64_t block_size = 8; // in 64bit words + static const uint64_t select_ones_per_hint = 64 * block_size * 2; // must be > block_size * 64 + static const uint64_t select_zeros_per_hint = select_ones_per_hint; + + typedef mapper::mappable_vector uint64_vec; + uint64_vec m_block_rank_pairs; + uint64_vec m_select_hints; + uint64_vec m_select0_hints; + }; +} diff --git a/3party/succinct/succinct_config.hpp.in b/3party/succinct/succinct_config.hpp.in new file mode 100644 index 0000000000..a7f6b83ed7 --- /dev/null +++ b/3party/succinct/succinct_config.hpp.in @@ -0,0 +1,16 @@ +#pragma once + +#cmakedefine SUCCINCT_USE_LIBCXX 1 +#ifndef SUCCINCT_USE_LIBCXX +# define SUCCINCT_USE_LIBCXX 0 +#endif + +#cmakedefine SUCCINCT_USE_INTRINSICS 1 +#ifndef SUCCINCT_USE_INTRINSICS +# define SUCCINCT_USE_INTRINSICS 0 +#endif + +#cmakedefine SUCCINCT_USE_POPCNT 1 +#ifndef SUCCINCT_USE_POPCNT +# define SUCCINCT_USE_POPCNT 0 +#endif diff --git a/3party/succinct/tables.hpp b/3party/succinct/tables.hpp new file mode 100644 index 0000000000..83fdc481ea --- /dev/null +++ b/3party/succinct/tables.hpp @@ -0,0 +1,63 @@ +#pragma once + +#include + +namespace succinct { namespace tables { + + const uint8_t select_in_byte[2048] = { + 8, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5, 0, 1, 0, 2, 0, 1, + 0, 3, 0, 1, 0, 2, 0, 1, 0, 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 6, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, + 1, 0, 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4, 0, 1, 0, 2, + 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 7, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, + 2, 0, 1, 0, 5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 6, 0, 1, + 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5, 0, 1, 0, 2, 0, 1, 0, 3, 0, + 1, 0, 2, 0, 1, 0, 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 8, 8, 8, 1, 8, 2, 2, 1, 8, 3, 3, 1, 3, 2, 2, 1, 8, + 4, 4, 1, 4, 2, 2, 1, 4, 3, 3, 1, 3, 2, 2, 1, 8, 5, 5, 1, 5, 2, 2, 1, 5, 3, 3, 1, 3, 2, 2, 1, 5, 4, 4, 1, 4, 2, 2, 1, + 4, 3, 3, 1, 3, 2, 2, 1, 8, 6, 6, 1, 6, 2, 2, 1, 6, 3, 3, 1, 3, 2, 2, 1, 6, 4, 4, 1, 4, 2, 2, 1, 4, 3, 3, 1, 3, 2, 2, + 1, 6, 5, 5, 1, 5, 2, 2, 1, 5, 3, 3, 1, 3, 2, 2, 1, 5, 4, 4, 1, 4, 2, 2, 1, 4, 3, 3, 1, 3, 2, 2, 1, 8, 7, 7, 1, 7, 2, + 2, 1, 7, 3, 3, 1, 3, 2, 2, 1, 7, 4, 4, 1, 4, 2, 2, 1, 4, 3, 3, 1, 3, 2, 2, 1, 7, 5, 5, 1, 5, 2, 2, 1, 5, 3, 3, 1, 3, + 2, 2, 1, 5, 4, 4, 1, 4, 2, 2, 1, 4, 3, 3, 1, 3, 2, 2, 1, 7, 6, 6, 1, 6, 2, 2, 1, 6, 3, 3, 1, 3, 2, 2, 1, 6, 4, 4, 1, + 4, 2, 2, 1, 4, 3, 3, 1, 3, 2, 2, 1, 6, 5, 5, 1, 5, 2, 2, 1, 5, 3, 3, 1, 3, 2, 2, 1, 5, 4, 4, 1, 4, 2, 2, 1, 4, 3, 3, + 1, 3, 2, 2, 1, 8, 8, 8, 8, 8, 8, 8, 2, 8, 8, 8, 3, 8, 3, 3, 2, 8, 8, 8, 4, 8, 4, 4, 2, 8, 4, 4, 3, 4, 3, 3, 2, 8, 8, + 8, 5, 8, 5, 5, 2, 8, 5, 5, 3, 5, 3, 3, 2, 8, 5, 5, 4, 5, 4, 4, 2, 5, 4, 4, 3, 4, 3, 3, 2, 8, 8, 8, 6, 8, 6, 6, 2, 8, + 6, 6, 3, 6, 3, 3, 2, 8, 6, 6, 4, 6, 4, 4, 2, 6, 4, 4, 3, 4, 3, 3, 2, 8, 6, 6, 5, 6, 5, 5, 2, 6, 5, 5, 3, 5, 3, 3, 2, + 6, 5, 5, 4, 5, 4, 4, 2, 5, 4, 4, 3, 4, 3, 3, 2, 8, 8, 8, 7, 8, 7, 7, 2, 8, 7, 7, 3, 7, 3, 3, 2, 8, 7, 7, 4, 7, 4, 4, + 2, 7, 4, 4, 3, 4, 3, 3, 2, 8, 7, 7, 5, 7, 5, 5, 2, 7, 5, 5, 3, 5, 3, 3, 2, 7, 5, 5, 4, 5, 4, 4, 2, 5, 4, 4, 3, 4, 3, + 3, 2, 8, 7, 7, 6, 7, 6, 6, 2, 7, 6, 6, 3, 6, 3, 3, 2, 7, 6, 6, 4, 6, 4, 4, 2, 6, 4, 4, 3, 4, 3, 3, 2, 7, 6, 6, 5, 6, + 5, 5, 2, 6, 5, 5, 3, 5, 3, 3, 2, 6, 5, 5, 4, 5, 4, 4, 2, 5, 4, 4, 3, 4, 3, 3, 2, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 3, 8, 8, 8, 8, 8, 8, 8, 4, 8, 8, 8, 4, 8, 4, 4, 3, 8, 8, 8, 8, 8, 8, 8, 5, 8, 8, 8, 5, 8, 5, 5, 3, 8, 8, 8, + 5, 8, 5, 5, 4, 8, 5, 5, 4, 5, 4, 4, 3, 8, 8, 8, 8, 8, 8, 8, 6, 8, 8, 8, 6, 8, 6, 6, 3, 8, 8, 8, 6, 8, 6, 6, 4, 8, 6, + 6, 4, 6, 4, 4, 3, 8, 8, 8, 6, 8, 6, 6, 5, 8, 6, 6, 5, 6, 5, 5, 3, 8, 6, 6, 5, 6, 5, 5, 4, 6, 5, 5, 4, 5, 4, 4, 3, 8, + 8, 8, 8, 8, 8, 8, 7, 8, 8, 8, 7, 8, 7, 7, 3, 8, 8, 8, 7, 8, 7, 7, 4, 8, 7, 7, 4, 7, 4, 4, 3, 8, 8, 8, 7, 8, 7, 7, 5, + 8, 7, 7, 5, 7, 5, 5, 3, 8, 7, 7, 5, 7, 5, 5, 4, 7, 5, 5, 4, 5, 4, 4, 3, 8, 8, 8, 7, 8, 7, 7, 6, 8, 7, 7, 6, 7, 6, 6, + 3, 8, 7, 7, 6, 7, 6, 6, 4, 7, 6, 6, 4, 6, 4, 4, 3, 8, 7, 7, 6, 7, 6, 6, 5, 7, 6, 6, 5, 6, 5, 5, 3, 7, 6, 6, 5, 6, 5, + 5, 4, 6, 5, 5, 4, 5, 4, 4, 3, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 4, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 5, 8, 8, 8, 8, 8, 8, 8, 5, 8, 8, 8, 5, 8, 5, 5, 4, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 6, 8, 8, 8, 8, 8, 8, 8, 6, 8, 8, 8, 6, 8, 6, 6, 4, 8, 8, 8, 8, 8, 8, 8, 6, 8, 8, 8, + 6, 8, 6, 6, 5, 8, 8, 8, 6, 8, 6, 6, 5, 8, 6, 6, 5, 6, 5, 5, 4, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 7, 8, 8, + 8, 8, 8, 8, 8, 7, 8, 8, 8, 7, 8, 7, 7, 4, 8, 8, 8, 8, 8, 8, 8, 7, 8, 8, 8, 7, 8, 7, 7, 5, 8, 8, 8, 7, 8, 7, 7, 5, 8, + 7, 7, 5, 7, 5, 5, 4, 8, 8, 8, 8, 8, 8, 8, 7, 8, 8, 8, 7, 8, 7, 7, 6, 8, 8, 8, 7, 8, 7, 7, 6, 8, 7, 7, 6, 7, 6, 6, 4, + 8, 8, 8, 7, 8, 7, 7, 6, 8, 7, 7, 6, 7, 6, 6, 5, 8, 7, 7, 6, 7, 6, 6, 5, 7, 6, 6, 5, 6, 5, 5, 4, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 5, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 6, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 6, 8, 8, 8, 8, 8, 8, 8, 6, 8, 8, 8, 6, + 8, 6, 6, 5, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 7, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 7, 8, 8, 8, 8, 8, 8, 8, 7, 8, 8, 8, 7, 8, 7, 7, 5, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 7, 8, 8, 8, 8, 8, 8, 8, 7, 8, 8, 8, 7, 8, 7, 7, 6, 8, 8, 8, 8, 8, 8, 8, 7, 8, 8, 8, 7, 8, 7, 7, 6, 8, + 8, 8, 7, 8, 7, 7, 6, 8, 7, 7, 6, 7, 6, 6, 5, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 6, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 7, 8, 8, 8, 8, 8, 8, 8, 7, 8, 8, 8, 7, 8, 7, 7, 6, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, + 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 7 + }; + +}} diff --git a/3party/succinct/test_bit_vector.cpp b/3party/succinct/test_bit_vector.cpp new file mode 100644 index 0000000000..eaa08ef63b --- /dev/null +++ b/3party/succinct/test_bit_vector.cpp @@ -0,0 +1,196 @@ +#define BOOST_TEST_MODULE bit_vector +#include "test_common.hpp" +#include "test_rank_select_common.hpp" + +#include +#include + +#include "mapper.hpp" +#include "bit_vector.hpp" + +BOOST_AUTO_TEST_CASE(bit_vector) +{ + srand(42); + + std::vector v = random_bit_vector(); + + { + succinct::bit_vector_builder bvb; + for (size_t i = 0; i < v.size(); ++i) { + bvb.push_back(v[i]); + } + + succinct::bit_vector bitmap(&bvb); + test_equal_bits(v, bitmap, "Random bits (push_back)"); + } + + { + succinct::bit_vector_builder bvb(v.size()); + for (size_t i = 0; i < v.size(); ++i) { + bvb.set(i, v[i]); + } + bvb.push_back(0); + v.push_back(0); + bvb.push_back(1); + v.push_back(1); + + succinct::bit_vector bitmap(&bvb); + test_equal_bits(v, bitmap, "Random bits (set)"); + } + + uint64_t ints[] = {uint64_t(-1), uint64_t(1) << 63, 1, 1, 1, 3, 5, 7, 0xFFF, 0xF0F, 1, 0xFFFFFF, 0x123456, uint64_t(1) << 63, uint64_t(-1)}; + { + succinct::bit_vector_builder bvb; + BOOST_FOREACH(uint64_t i, ints) { + uint64_t len = succinct::broadword::msb(i) + 1; + bvb.append_bits(i, len); + } + succinct::bit_vector bitmap(&bvb); + uint64_t pos = 0; + BOOST_FOREACH(uint64_t i, ints) { + uint64_t len = succinct::broadword::msb(i) + 1; + BOOST_REQUIRE_EQUAL(i, bitmap.get_bits(pos, len)); + pos += len; + } + } + + { + using succinct::broadword::msb; + std::vector positions(1); + BOOST_FOREACH(uint64_t i, ints) { + positions.push_back(positions.back() + msb(i) + 1); + } + + succinct::bit_vector_builder bvb(positions.back()); + + for (size_t i = 0; i < positions.size() - 1; ++i) { + uint64_t v = ints[i]; + uint64_t len = positions[i + 1] - positions[i]; + bvb.set_bits(positions[i], v, len); + } + + succinct::bit_vector bitmap(&bvb); + for (size_t i = 0; i < positions.size() - 1; ++i) { + uint64_t v = ints[i]; + uint64_t len = positions[i + 1] - positions[i]; + BOOST_REQUIRE_EQUAL(v, bitmap.get_bits(positions[i], len)); + } + } +} + +BOOST_AUTO_TEST_CASE(bit_vector_enumerator) +{ + srand(42); + std::vector v = random_bit_vector(); + succinct::bit_vector bitmap(v); + + size_t i = 0; + size_t pos = 0; + + succinct::bit_vector::enumerator e(bitmap, pos); + while (pos < bitmap.size()) { + bool next = e.next(); + MY_REQUIRE_EQUAL(next, v[pos], "pos = " << pos << " i = " << i); + pos += 1; + + pos += size_t(rand()) % (bitmap.size() - pos + 1); + e = succinct::bit_vector::enumerator(bitmap, pos); + i += 1; + } +} + +BOOST_AUTO_TEST_CASE(bit_vector_unary_enumerator) +{ + srand(42); + uint64_t n = 20000; + std::vector v = random_bit_vector(n); + + // punch some long gaps in v + for (size_t g = 0; g < n / 1000; ++g) { + ssize_t l = std::min(ssize_t(rand() % 256), ssize_t(v.size() - g)); + std::fill(v.begin(), v.begin() + l, 0); + } + + succinct::bit_vector bitmap(v); + + std::vector ones; + for (size_t i = 0; i < v.size(); ++i) { + if (bitmap[i]) { + ones.push_back(i); + } + } + + { + succinct::bit_vector::unary_enumerator e(bitmap, 0); + + for (size_t r = 0; r < ones.size(); ++r) { + uint64_t pos = e.next(); + MY_REQUIRE_EQUAL(ones[r], pos, + "r = " << r); + } + } + + { + succinct::bit_vector::unary_enumerator e(bitmap, 0); + + for (size_t r = 0; r < ones.size(); ++r) { + for (size_t k = 0; k < std::min(size_t(256), size_t(ones.size() - r)); ++k) { + succinct::bit_vector::unary_enumerator ee(e); + ee.skip(k); + uint64_t pos = ee.next(); + MY_REQUIRE_EQUAL(ones[r + k], pos, + "r = " << r << " k = " << k); + } + e.next(); + } + } + + { + succinct::bit_vector::unary_enumerator e(bitmap, 0); + + for (size_t pos = 0; pos < v.size(); ++pos) { + uint64_t skip = 0; + for (size_t d = 0; d < std::min(size_t(256), size_t(v.size() - pos)); ++d) { + if (v[pos + d] == 0) { + succinct::bit_vector::unary_enumerator ee(bitmap, pos); + ee.skip0(skip); + + uint64_t expected_pos = pos + d; + for (; !v[expected_pos] && expected_pos < v.size(); ++expected_pos); + if (!v[expected_pos]) break; + uint64_t pos = ee.next(); + MY_REQUIRE_EQUAL(expected_pos, pos, + "pos = " << pos << " skip = " << skip); + + skip += 1; + } + } + } + } +} + +void test_bvb_reverse(size_t n) +{ + std::vector v = random_bit_vector(n); + succinct::bit_vector_builder bvb; + for (size_t i = 0; i < v.size(); ++i) { + bvb.push_back(v[i]); + } + + std::reverse(v.begin(), v.end()); + bvb.reverse(); + + succinct::bit_vector bitmap(&bvb); + test_equal_bits(v, bitmap, "In-place reverse"); +} + +BOOST_AUTO_TEST_CASE(bvb_reverse) +{ + srand(42); + + test_bvb_reverse(0); + test_bvb_reverse(63); + test_bvb_reverse(64); + test_bvb_reverse(1000); + test_bvb_reverse(1024); +} diff --git a/3party/succinct/test_bp_vector.cpp b/3party/succinct/test_bp_vector.cpp new file mode 100644 index 0000000000..d63cf8c9d6 --- /dev/null +++ b/3party/succinct/test_bp_vector.cpp @@ -0,0 +1,92 @@ +#define BOOST_TEST_MODULE bp_vector +#include "test_common.hpp" + +#include +#include + +#include "mapper.hpp" +#include "bp_vector.hpp" +#include "test_bp_vector_common.hpp" + +template +void test_parentheses(std::vector const& v, BPVector const& bitmap, std::string test_name) +{ + std::stack stack; + std::vector open(v.size()); + std::vector close(v.size()); + std::vector enclose(v.size(), uint64_t(-1)); + + for (size_t i = 0; i < v.size(); ++i) { + if (v[i]) { // opening + if (!stack.empty()) { + enclose[i] = stack.top(); + } + stack.push(i); + } else { // closing + BOOST_REQUIRE(!stack.empty()); // this is more a test on the test + size_t opening = stack.top(); + stack.pop(); + close[opening] = i; + open[i] = opening; + + } + } + BOOST_REQUIRE_EQUAL(0U, stack.size()); // ditto as above + + for (size_t i = 0; i < bitmap.size(); ++i) { + if (v[i]) { // opening + if (enclose[i] != uint64_t(-1)) { + MY_REQUIRE_EQUAL(enclose[i], bitmap.enclose(i), + "enclose (" << test_name << "): i = " << i); + } + MY_REQUIRE_EQUAL(close[i], bitmap.find_close(i), + "find_close (" << test_name << "): i = " << i); + } else { // closing + MY_REQUIRE_EQUAL(open[i], bitmap.find_open(i), + "find_open (" << test_name << "): i = " << i); + } + } +} + +BOOST_AUTO_TEST_CASE(bp_vector) +{ + srand(42); + + { + std::vector v; + succinct::bp_vector bitmap(v); + test_parentheses(v, bitmap, "Empty vector"); + } + + { + std::vector v; + succinct::random_bp(v, 100000); + succinct::bp_vector bitmap(v); + test_parentheses(v, bitmap, "Random parentheses"); + } + + { + size_t sizes[] = {2, 4, 512, 514, 8190, 8192, 8194, 16384, 16386, 100000}; + for (size_t i = 0; i < sizeof(sizes) / sizeof(sizes[0]); ++i) { + std::vector v; + succinct::random_binary_tree(v, sizes[i]); + succinct::bp_vector bitmap(v); + test_parentheses(v, bitmap, "Random binary tree"); + } + } + + { + size_t sizes[] = {2, 4, 512, 514, 8190, 8192, 8194, 16384, 16386, 32768, 32770}; + size_t iterations[] = {1, 2, 3}; + for (size_t s = 0; s < sizeof(sizes) / sizeof(sizes[0]); ++s) { + for (size_t r = 0; r < sizeof(iterations) / sizeof(iterations[0]); ++r) { + std::vector v; + for (size_t i = 0; i < iterations[r]; ++i) { + succinct::bp_path(v, sizes[s]); + } + succinct::bp_vector bitmap(v); + test_parentheses(v, bitmap, "Nested parentheses"); + } + } + } +} diff --git a/3party/succinct/test_bp_vector_common.hpp b/3party/succinct/test_bp_vector_common.hpp new file mode 100644 index 0000000000..6ceb0a59e8 --- /dev/null +++ b/3party/succinct/test_bp_vector_common.hpp @@ -0,0 +1,68 @@ +#pragma once + +namespace succinct { + + namespace detail { + template + void random_binary_tree_helper(BitVectorBuilder& builder, size_t size) + { + assert((size & 1) == 1); // binary trees can only have an odd number of nodes (internal + leaves) + if (size == 1) { + builder.push_back(0); // can only be a leaf + return; + } + + builder.push_back(1); + size_t left_subtree_size = 2 * (size_t(rand()) % (size - 1) / 2) + 1; + assert(left_subtree_size >= 1); + size_t right_subtree_size = size - 1 - left_subtree_size; + assert(right_subtree_size >= 1); + assert(left_subtree_size + right_subtree_size + 1 == size); + + random_binary_tree_helper(builder, left_subtree_size); + random_binary_tree_helper(builder, right_subtree_size); + } + } + + + template + void random_binary_tree(BitVectorBuilder& builder, size_t size) + { + assert((size & 1) == 0 && size >= 2); + + builder.push_back(1); // fake root + detail::random_binary_tree_helper(builder, size - 1); + } + + template + void random_bp(BitVectorBuilder& builder, size_t size_est) + { + int excess = 0; + for (size_t i = 0; i < size_est; ++i) { + bool val = rand() > (RAND_MAX / 2); + if (excess <= 1 && !val) { + val = 1; + } + excess += (val ? 1 : -1); + builder.push_back(val); + } + + for (size_t i = 0; i < size_t(excess); ++i) { + builder.push_back(0); // close all parentheses + } + } + + template + void bp_path(BitVectorBuilder& builder, size_t size) + { + assert((size & 1) == 0); + for (size_t i = 0; i < size / 2; ++i) { + builder.push_back(1); + } + for (size_t i = 0; i < size / 2; ++i) { + builder.push_back(0); + } + } + + +} diff --git a/3party/succinct/test_bp_vector_rmq.cpp b/3party/succinct/test_bp_vector_rmq.cpp new file mode 100644 index 0000000000..1f45999d9c --- /dev/null +++ b/3party/succinct/test_bp_vector_rmq.cpp @@ -0,0 +1,103 @@ +#define BOOST_TEST_MODULE bp_vector_rmq +#include "test_common.hpp" + +#include +#include + +#include "mapper.hpp" +#include "bp_vector.hpp" +#include "test_bp_vector_common.hpp" + + +template +void test_rmq(std::vector const& v, BPVector const& bitmap, std::string test_name) +{ + // XXX should we test this case? + if (v.empty()) return; + + // test all values from a to v.size() for a in specific locations + // plus a few random + + std::vector tests; + tests.push_back(1); + tests.push_back(8); + tests.push_back(64); + tests.push_back(8192); + tests.push_back(v.size()); + for (size_t t = 0; t < 10; ++t) { + tests.push_back(size_t(rand()) % v.size()); + } + + for(size_t t = 0; t < tests.size(); ++t) { + uint64_t a = tests[t]; + if (a > v.size()) continue; + + typename BPVector::enumerator bp_it(bitmap, a); + typename BPVector::excess_t cur_exc = bitmap.excess(a); + typename BPVector::excess_t min_exc = cur_exc, found_min_exc; + uint64_t min_idx = a; + + BOOST_REQUIRE_EQUAL(min_idx, bitmap.excess_rmq(a, a, found_min_exc)); + + for (uint64_t b = a + 1; b < v.size(); ++b) { + cur_exc += bp_it.next() ? 1 : -1; + if (cur_exc < min_exc) { + min_exc = cur_exc; + min_idx = b; + + assert(min_exc == bitmap.excess(min_idx)); + } + + MY_REQUIRE_EQUAL(min_idx, bitmap.excess_rmq(a, b, found_min_exc), + "excess_rmq (" << test_name << "):" + << " a = " << a + << " b = " << b + << " min_exc = " << min_exc + << " found_min_exc = " << found_min_exc + ); + } + } +} + +BOOST_AUTO_TEST_CASE(bp_vector_rmq) +{ + srand(42); + + { + std::vector v; + succinct::bp_vector bitmap(v); + test_rmq(v, bitmap, "Empty vector"); + } + + { + std::vector v; + succinct::random_bp(v, 100000); + succinct::bp_vector bitmap(v); + test_rmq(v, bitmap, "Random parentheses"); + } + + { + size_t sizes[] = {2, 4, 512, 514, 8190, 8192, 8194, 16384, 16386, 100000}; + for (size_t i = 0; i < sizeof(sizes) / sizeof(sizes[0]); ++i) { + std::vector v; + succinct::random_binary_tree(v, sizes[i]); + succinct::bp_vector bitmap(v); + test_rmq(v, bitmap, "Random binary tree"); + } + } + + { + size_t sizes[] = {2, 4, 512, 514, 8190, 8192, 8194, 16384, 16386, 32768, 32770}; + size_t iterations[] = {1, 2, 3}; + for (size_t s = 0; s < sizeof(sizes) / sizeof(sizes[0]); ++s) { + for (size_t r = 0; r < sizeof(iterations) / sizeof(iterations[0]); ++r) { + std::vector v; + for (size_t i = 0; i < iterations[r]; ++i) { + succinct::bp_path(v, sizes[s]); + } + succinct::bp_vector bitmap(v); + test_rmq(v, bitmap, "Nested parentheses"); + } + } + } +} diff --git a/3party/succinct/test_cartesian_tree.cpp b/3party/succinct/test_cartesian_tree.cpp new file mode 100644 index 0000000000..8d66814e45 --- /dev/null +++ b/3party/succinct/test_cartesian_tree.cpp @@ -0,0 +1,121 @@ +#define BOOST_TEST_MODULE cartesian_tree +#include "test_common.hpp" + +#include +#include + +#include "mapper.hpp" +#include "cartesian_tree.hpp" + +typedef uint64_t value_type; + +// XXX test (de)serialization + +template +void test_rmq(std::vector const& v, succinct::cartesian_tree const& tree, + Comparator const& comp, std::string test_name) +{ + BOOST_REQUIRE_EQUAL(v.size(), tree.size()); + + if (v.empty()) return; + + std::vector tests; + // A few special cases + tests.push_back(0); + tests.push_back(1); + // This is the global minimum of the vector + tests.push_back(uint64_t(std::min_element(v.begin(), v.end(), comp) - v.begin())); + + // Plus some random... + for (size_t t = 0; t < 10; ++t) { + tests.push_back(size_t(rand()) % v.size()); + } + + for(size_t t = 0; t < tests.size(); ++t) { + uint64_t a = tests[t]; + if (a > v.size()) continue; + + uint64_t min_idx = a; + value_type cur_min = v[a]; + + BOOST_REQUIRE_EQUAL(min_idx, tree.rmq(a, a)); + + for (uint64_t b = a + 1; b < v.size(); ++b) { + if (comp(v[b], cur_min)) { + cur_min = v[b]; + min_idx = b; + } + + uint64_t found_idx = tree.rmq(a, b); + + MY_REQUIRE_EQUAL(min_idx, found_idx, + "rmq (" << test_name << "):" + << " a = " << a + << " b = " << b + << " min = " << cur_min + << " found_min = " << v[found_idx] + ); + } + } +} + +BOOST_AUTO_TEST_CASE(cartesian_tree) +{ + srand(42); + + { + std::vector v; + succinct::cartesian_tree t(v); + test_rmq(v, t, std::less(), "Empty vector"); + } + + { + std::vector v(20000); + for (size_t i = 0; i < v.size(); ++i) { + v[i] = i; + } + + { + succinct::cartesian_tree t(v); + test_rmq(v, t, std::less(), "Increasing values"); + } + { + succinct::cartesian_tree t(v, std::greater()); + test_rmq(v, t, std::greater(), "Decreasing values"); + } + } + + { + std::vector v(20000); + for (size_t i = 0; i < v.size(); ++i) { + if (i < v.size() / 2) { + v[i] = i; + } else { + v[i] = v.size() - i; + } + } + + { + succinct::cartesian_tree t(v); + test_rmq(v, t, std::less(), "Convex values"); + } + + { + succinct::cartesian_tree t(v, std::greater()); + test_rmq(v, t, std::greater(), "Concave values"); + } + } + + { + size_t sizes[] = {2, 4, 512, 514, 8190, 8192, 8194, 16384, 16386, 100000}; + for (size_t i = 0; i < sizeof(sizes) / sizeof(sizes[0]); ++i) { + std::vector v(sizes[i]); + for (size_t i = 0; i < v.size(); ++i) { + v[i] = size_t(rand()) % 1024; + } + + succinct::cartesian_tree t(v); + test_rmq(v, t, std::less(), "Random values"); + } + } +} diff --git a/3party/succinct/test_common.hpp b/3party/succinct/test_common.hpp new file mode 100644 index 0000000000..39227221b4 --- /dev/null +++ b/3party/succinct/test_common.hpp @@ -0,0 +1,20 @@ +#pragma once + +#define BOOST_TEST_DYN_LINK +#include + +#include +#include +#include + +#define MY_REQUIRE_EQUAL(A, B, MSG) \ + BOOST_REQUIRE_MESSAGE((A) == (B), BOOST_PP_STRINGIZE(A) << " == " << BOOST_PP_STRINGIZE(B) << " [" << A << " != " << B << "] " << MSG) + +inline std::vector random_bit_vector(size_t n = 10000, double density = 0.5) +{ + std::vector v; + for (size_t i = 0; i < n; ++i) { + v.push_back(rand() < (RAND_MAX * density)); + } + return v; +} diff --git a/3party/succinct/test_darray.cpp b/3party/succinct/test_darray.cpp new file mode 100644 index 0000000000..519dd44544 --- /dev/null +++ b/3party/succinct/test_darray.cpp @@ -0,0 +1,76 @@ +#define BOOST_TEST_MODULE darray +#include "test_common.hpp" +#include "test_rank_select_common.hpp" + +#include +#include + +#include "mapper.hpp" +#include "darray.hpp" + +void test_darray(std::vector const& v, const char* test_name) +{ + succinct::bit_vector bv(v); + succinct::darray1 d1(bv); + succinct::darray0 d0(bv); + + size_t cur_rank = 0; + size_t cur_rank0 = 0; + for (size_t i = 0; i < v.size(); ++i) { + if (v[i]) { + MY_REQUIRE_EQUAL(i, d1.select(bv, cur_rank), + "select (" << test_name << "): cur_rank = " << cur_rank << ", i = " << i << ", v[i] = " << v[i]); + cur_rank += 1; + } else { + MY_REQUIRE_EQUAL(i, d0.select(bv, cur_rank0), + "select0 (" << test_name << "): cur_rank0 = " << cur_rank0 << ", i = " << i << ", v[i] = " << v[i]); + cur_rank0 += 1; + } + } + + BOOST_REQUIRE_EQUAL(cur_rank, d1.num_positions()); + BOOST_REQUIRE_EQUAL(cur_rank0, d0.num_positions()); +} + +BOOST_AUTO_TEST_CASE(darray) +{ + srand(42); + size_t N = 10000; + + { + // Random bitmap + std::vector v = random_bit_vector(N); + test_darray(v, "random"); + } + + { + // Empty bitmap + std::vector v; + test_darray(v, "empty"); + } + + { + // Only one value + std::vector v(N); + v[37] = 1; + test_darray(v, "singleton"); + } + + { + // Full bitmap + std::vector v(N, 1); + test_darray(v, "full"); + } + + { + // Very sparse random bitmap + size_t bigN = (1 << 16) * 4; + std::vector v(bigN); + size_t cur_pos = 0; + while(cur_pos < bigN) { + v[cur_pos] = 1; + cur_pos += rand() % 1024; + } + test_darray(v, "sparse"); + } +} diff --git a/3party/succinct/test_elias_fano.cpp b/3party/succinct/test_elias_fano.cpp new file mode 100644 index 0000000000..8578ac838b --- /dev/null +++ b/3party/succinct/test_elias_fano.cpp @@ -0,0 +1,71 @@ +#define BOOST_TEST_MODULE elias_fano +#include "test_common.hpp" +#include "test_rank_select_common.hpp" + +#include +#include + +#include "mapper.hpp" +#include "elias_fano.hpp" + +BOOST_AUTO_TEST_CASE(elias_fano) +{ + srand(42); + size_t N = 10000; + + { + // Random bitmap + for (size_t d = 1; d < 8; ++d) { + double density = 1.0 / (1 << d); + std::vector v = random_bit_vector(N, density); + + succinct::bit_vector_builder bvb; + for (size_t i = 0; i < v.size(); ++i) { + bvb.push_back(v[i]); + } + + succinct::elias_fano bitmap(&bvb); + test_equal_bits(v, bitmap, "Random bitmap"); + test_rank_select1(v, bitmap, "Random bitmap"); + test_delta(bitmap, "Random bitmap"); + test_select_enumeration(v, bitmap, "Random bitmap"); + } + } + + { + // Empty bitmap + succinct::bit_vector_builder bvb(N); + succinct::elias_fano bitmap(&bvb); + BOOST_REQUIRE_EQUAL(0U, bitmap.num_ones()); + test_equal_bits(std::vector(N), bitmap, "Empty bitmap"); + test_select_enumeration(std::vector(N), bitmap, "Empty bitmap"); + } + + { + // Only one value + std::vector v(N); + succinct::bit_vector_builder bvb(N); + bvb.set(37, 1); + v[37] = 1; + succinct::elias_fano bitmap(&bvb); + test_equal_bits(v, bitmap, "Only one value"); + test_rank_select1(v, bitmap, "Only one value"); + test_delta(bitmap, "Only one value"); + test_select_enumeration(v, bitmap, "Only one value"); + BOOST_REQUIRE_EQUAL(1U, bitmap.num_ones()); + } + + { + // Full bitmap + std::vector v(N, 1); + succinct::bit_vector_builder bvb; + for (size_t i = 0; i < N; ++i) { + bvb.push_back(1); + } + succinct::elias_fano bitmap(&bvb); + test_equal_bits(v, bitmap, "Full bitmap"); + test_rank_select1(v, bitmap, "Full bitmap"); + test_delta(bitmap, "Full bitmap"); + test_select_enumeration(v, bitmap, "Full bitmap"); + } +} diff --git a/3party/succinct/test_elias_fano_compressed_list.cpp b/3party/succinct/test_elias_fano_compressed_list.cpp new file mode 100644 index 0000000000..014b9a298c --- /dev/null +++ b/3party/succinct/test_elias_fano_compressed_list.cpp @@ -0,0 +1,29 @@ +#define BOOST_TEST_MODULE elias_fano_compressed_list +#include "test_common.hpp" + +#include + +#include "elias_fano_compressed_list.hpp" + +BOOST_AUTO_TEST_CASE(elias_fano_compressed_list) +{ + srand(42); + const size_t test_size = 12345; + + std::vector v; + + for (size_t i = 0; i < test_size; ++i) { + if (rand() < (RAND_MAX / 3)) { + v.push_back(0); + } else { + v.push_back(size_t(rand())); + } + } + + succinct::elias_fano_compressed_list vv(v); + + BOOST_REQUIRE_EQUAL(v.size(), vv.size()); + for (size_t i = 0; i < v.size(); ++i) { + MY_REQUIRE_EQUAL(v[i], vv[i], "i = " << i); + } +} diff --git a/3party/succinct/test_gamma_bit_vector.cpp b/3party/succinct/test_gamma_bit_vector.cpp new file mode 100644 index 0000000000..3239460ce3 --- /dev/null +++ b/3party/succinct/test_gamma_bit_vector.cpp @@ -0,0 +1,62 @@ +#define BOOST_TEST_MODULE gamma_bit_vector +#include "test_common.hpp" + +#include + +#include "gamma_bit_vector.hpp" + +typedef std::vector std_vector_type; + +std_vector_type random_vector(size_t test_size) +{ + std_vector_type v; + + for (size_t i = 0; i < test_size; ++i) { + bool b = uint64_t(rand()) & 1; + if (rand() < (RAND_MAX / 3)) { + v.push_back(b); + } else { + v.push_back((uint64_t(rand()) << 1) | b); + } + } + + return v; +} + +BOOST_AUTO_TEST_CASE(gamma_bit_vector) +{ + srand(42); + const size_t test_size = 12345; + std_vector_type v = random_vector(test_size); + + succinct::gamma_bit_vector vv(v); + + BOOST_REQUIRE_EQUAL(v.size(), vv.size()); + for (size_t i = 0; i < v.size(); ++i) { + MY_REQUIRE_EQUAL(v[i], vv[i], "i = " << i); + } +} + +BOOST_AUTO_TEST_CASE(gamma_bit_enumerator) +{ + srand(42); + const size_t test_size = 12345; + std_vector_type v = random_vector(test_size); + + succinct::gamma_bit_vector vv(v); + + size_t i = 0; + size_t pos = 0; + + succinct::forward_enumerator e(vv, pos); + while (pos < vv.size()) { + succinct::gamma_bit_vector::value_type next = e.next(); + MY_REQUIRE_EQUAL(next, v[pos], "pos = " << pos << " i = " << i); + pos += 1; + + size_t step = uint64_t(rand()) % (vv.size() - pos + 1); + pos += step; + e = succinct::forward_enumerator(vv, pos); + i += 1; + } +} diff --git a/3party/succinct/test_gamma_vector.cpp b/3party/succinct/test_gamma_vector.cpp new file mode 100644 index 0000000000..097f07eebb --- /dev/null +++ b/3party/succinct/test_gamma_vector.cpp @@ -0,0 +1,62 @@ +#define BOOST_TEST_MODULE gamma_vector +#include "test_common.hpp" + +#include + +#include "gamma_vector.hpp" + +BOOST_AUTO_TEST_CASE(gamma_vector) +{ + srand(42); + const size_t test_size = 12345; + + std::vector v; + + for (size_t i = 0; i < test_size; ++i) { + if (rand() < (RAND_MAX / 3)) { + v.push_back(0); + } else { + v.push_back(uint64_t(rand())); + } + } + + succinct::gamma_vector vv(v); + + BOOST_REQUIRE_EQUAL(v.size(), vv.size()); + for (size_t i = 0; i < v.size(); ++i) { + MY_REQUIRE_EQUAL(v[i], vv[i], "i = " << i); + } +} + +BOOST_AUTO_TEST_CASE(gamma_enumerator) +{ + srand(42); + const size_t test_size = 12345; + + std::vector v; + + for (size_t i = 0; i < test_size; ++i) { + if (rand() < (RAND_MAX / 3)) { + v.push_back(0); + } else { + v.push_back(uint64_t(rand())); + } + } + + succinct::gamma_vector vv(v); + + size_t i = 0; + size_t pos = 0; + + succinct::forward_enumerator e(vv, pos); + while (pos < vv.size()) { + uint64_t next = e.next(); + MY_REQUIRE_EQUAL(next, v[pos], "pos = " << pos << " i = " << i); + pos += 1; + + size_t step = uint64_t(rand()) % (vv.size() - pos + 1); + pos += step; + e = succinct::forward_enumerator(vv, pos); + i += 1; + } +} diff --git a/3party/succinct/test_mapper.cpp b/3party/succinct/test_mapper.cpp new file mode 100644 index 0000000000..46840ae456 --- /dev/null +++ b/3party/succinct/test_mapper.cpp @@ -0,0 +1,77 @@ +#define BOOST_TEST_MODULE mapper +#include "test_common.hpp" + +#include + +#include "mapper.hpp" + +BOOST_AUTO_TEST_CASE(basic_map) +{ + succinct::mapper::mappable_vector vec; + BOOST_REQUIRE_EQUAL(vec.size(), 0U); + + int nums[] = {1, 2, 3, 4}; + vec.assign(nums); + + BOOST_REQUIRE_EQUAL(4U, vec.size()); + BOOST_REQUIRE_EQUAL(1, vec[0]); + BOOST_REQUIRE_EQUAL(4, vec[3]); + + succinct::mapper::freeze(vec, "temp.bin"); + + { + succinct::mapper::mappable_vector mapped_vec; + boost::iostreams::mapped_file_source m("temp.bin"); + succinct::mapper::map(mapped_vec, m); + BOOST_REQUIRE_EQUAL(vec.size(), mapped_vec.size()); + BOOST_REQUIRE(std::equal(vec.begin(), vec.end(), mapped_vec.begin())); + } + + boost::filesystem::remove("temp.bin"); +} + +class complex_struct { +public: + complex_struct() + : m_a(0) + {} + + void init() { + m_a = 42; + uint32_t b[] = {1, 2}; + m_b.assign(b); + } + + template + void map(Visitor& visit) { + visit + (m_a, "m_a") + (m_b, "m_b") + ; + } + + uint64_t m_a; + succinct::mapper::mappable_vector m_b; +}; + +BOOST_AUTO_TEST_CASE(complex_struct_map) +{ + complex_struct s; + s.init(); + succinct::mapper::freeze(s, "temp.bin"); + + BOOST_REQUIRE_EQUAL(24, succinct::mapper::size_of(s)); + + complex_struct mapped_s; + BOOST_REQUIRE_EQUAL(0, mapped_s.m_a); + BOOST_REQUIRE_EQUAL(0U, mapped_s.m_b.size()); + + { + boost::iostreams::mapped_file_source m("temp.bin"); + succinct::mapper::map(mapped_s, m); + BOOST_REQUIRE_EQUAL(s.m_a, mapped_s.m_a); + BOOST_REQUIRE_EQUAL(s.m_b.size(), mapped_s.m_b.size()); + } + + boost::filesystem::remove("temp.bin"); +} diff --git a/3party/succinct/test_rank_select_common.hpp b/3party/succinct/test_rank_select_common.hpp new file mode 100644 index 0000000000..f96318f1d0 --- /dev/null +++ b/3party/succinct/test_rank_select_common.hpp @@ -0,0 +1,123 @@ +#pragma once + +#include "test_common.hpp" + +template +inline void test_equal_bits(std::vector const& v, Vector const& bitmap, const char* test_name) +{ + BOOST_REQUIRE_EQUAL(v.size(), bitmap.size()); + for (size_t i = 0; i < v.size(); ++i) { + MY_REQUIRE_EQUAL((bool)v[i], bitmap[i], + "operator[] (" << test_name << "): i=" << i); + } +} + +template +void test_rank_select0(std::vector const& v, Vector const& bitmap, const char* test_name) +{ + uint64_t cur_rank0 = 0; + uint64_t last_zero = uint64_t(-1); + + for (size_t i = 0; i < v.size(); ++i) { + MY_REQUIRE_EQUAL(cur_rank0, bitmap.rank0(i), + "rank0 (" << test_name << "): cur_rank0 = " << cur_rank0 << ", i = " << i << ", v[i] = " << v[i]); + if (!v[i]) { + last_zero = i; + MY_REQUIRE_EQUAL(last_zero, bitmap.select0(cur_rank0), + "select0 (" << test_name << "): cur_rank0 = " << cur_rank0 << ", i = " << i << ", v[i] = " << v[i] << ", last_zero = " << last_zero); + ++cur_rank0; + } + if (last_zero != uint64_t(-1)) { + MY_REQUIRE_EQUAL(last_zero, bitmap.predecessor0(i), + "predecessor0 (" << test_name << "): last_zero = " << last_zero <<", i = " << i << ",v[i] = " << v[i]); + } + } + + last_zero = uint64_t(-1); + for (size_t i = v.size() - 1; i + 1 > 0; --i) { + if (!v[i]) { + last_zero = i; + } + + if (last_zero != uint64_t(-1)) { + MY_REQUIRE_EQUAL(last_zero, bitmap.successor0(i), + "successor0 (" << test_name << "): last_zero = " << last_zero <<", i = " << i << ",v[i] = " << v[i]); + } + } +} + +template +void test_rank_select1(std::vector const& v, Vector const& bitmap, const char* test_name) +{ + uint64_t cur_rank = 0; + uint64_t last_one = uint64_t(-1); + + for (size_t i = 0; i < v.size(); ++i) { + MY_REQUIRE_EQUAL(cur_rank, bitmap.rank(i), + "rank (" << test_name << "): cur_rank = " << cur_rank << ", i = " << i << ", v[i] = " << v[i]); + + if (v[i]) { + last_one = i; + MY_REQUIRE_EQUAL(last_one, bitmap.select(cur_rank), + "select (" << test_name << "): cur_rank = " << cur_rank << ", i = " << i << ", v[i] = " << v[i] << ", last_one = " << last_one); + ++cur_rank; + } + + if (last_one != uint64_t(-1)) { + MY_REQUIRE_EQUAL(last_one, bitmap.predecessor1(i), + "predecessor1 (" << test_name << "): last_one = " << last_one <<", i = " << i << ",v[i] = " << v[i]); + } + } + + last_one = uint64_t(-1); + for (size_t i = v.size() - 1; i + 1 > 0; --i) { + if (v[i]) { + last_one = i; + } + + if (last_one != uint64_t(-1)) { + MY_REQUIRE_EQUAL(last_one, bitmap.successor1(i), + "successor1 (" << test_name << "): last_one = " << last_one <<", i = " << i << ",v[i] = " << v[i]); + } + } +} + +template +void test_rank_select(std::vector const& v, Vector const& bitmap, const char* test_name) +{ + test_rank_select0(v, bitmap, test_name); + test_rank_select1(v, bitmap, test_name); +} + +template +void test_delta(Vector const& bitmap, const char* test_name) +{ + for (size_t i = 0; i < bitmap.num_ones(); ++i) { + if (i) { + MY_REQUIRE_EQUAL(bitmap.select(i) - bitmap.select(i - 1), + bitmap.delta(i), + "delta (" << test_name << "), i = " << i); + } else { + MY_REQUIRE_EQUAL(bitmap.select(i), + bitmap.delta(i), + "delta (" << test_name << "), i = " << i); + } + } +} + +template +void test_select_enumeration(std::vector const& v, Vector const& bitmap, const char* test_name) +{ + // XXX test other starting points + typename Vector::select_enumerator it(bitmap, 0); + + for (size_t i = 0; i < v.size(); ++i) { + if (v[i]) { + uint64_t res = it.next(); + MY_REQUIRE_EQUAL(i, + res, + "select_iterator next (" << test_name << "), i = " << i + << ", n = " << bitmap.size() << ", m = " << bitmap.num_ones()); + } + } +} diff --git a/3party/succinct/test_rs_bit_vector.cpp b/3party/succinct/test_rs_bit_vector.cpp new file mode 100644 index 0000000000..fc4c854e29 --- /dev/null +++ b/3party/succinct/test_rs_bit_vector.cpp @@ -0,0 +1,59 @@ +#define BOOST_TEST_MODULE rs_bit_vector +#include "test_common.hpp" +#include "test_rank_select_common.hpp" + +#include +#include + +#include "mapper.hpp" +#include "rs_bit_vector.hpp" + +BOOST_AUTO_TEST_CASE(rs_bit_vector) +{ + srand(42); + + // empty vector + std::vector v; + succinct::rs_bit_vector bitmap; + + succinct::rs_bit_vector(v).swap(bitmap); + BOOST_REQUIRE_EQUAL(v.size(), bitmap.size()); + succinct::rs_bit_vector(v, true).swap(bitmap); + BOOST_REQUIRE_EQUAL(v.size(), bitmap.size()); + + // random vector + v = random_bit_vector(); + + succinct::rs_bit_vector(v).swap(bitmap); + BOOST_REQUIRE_EQUAL(v.size(), bitmap.size()); + test_equal_bits(v, bitmap, "RS - Uniform bits"); + test_rank_select(v, bitmap, "Uniform bits"); + + succinct::rs_bit_vector(v, true, true).swap(bitmap); + test_rank_select(v, bitmap, "Uniform bits - with hints"); + + v.resize(10000); + v[9999] = 1; + v[9000] = 1; + succinct::rs_bit_vector(v).swap(bitmap); + + BOOST_REQUIRE_EQUAL(v.size(), bitmap.size()); + test_rank_select(v, bitmap, "Long runs of 0"); + succinct::rs_bit_vector(v, true, true).swap(bitmap); + test_rank_select(v, bitmap, "Long runs of 0 - with hints"); + + // corner cases + v.clear(); + v.resize(10000); + v[0] = 1; + v[511] = 1; + v[512] = 1; + v[1024] = 1; + v[2112] = 1; + succinct::rs_bit_vector(v).swap(bitmap); + + BOOST_REQUIRE_EQUAL(v.size(), bitmap.size()); + test_rank_select(v, bitmap, "Corner cases"); + succinct::rs_bit_vector(v, true).swap(bitmap); + test_rank_select(v, bitmap, "Corner cases - with hints"); +} diff --git a/3party/succinct/test_topk_vector.cpp b/3party/succinct/test_topk_vector.cpp new file mode 100644 index 0000000000..677737b598 --- /dev/null +++ b/3party/succinct/test_topk_vector.cpp @@ -0,0 +1,110 @@ +#define BOOST_TEST_MODULE topk_vector +#include "test_common.hpp" + +#include +#include +#include + +#include "mapper.hpp" +#include "mappable_vector.hpp" +#include "topk_vector.hpp" +#include "elias_fano_compressed_list.hpp" + +typedef uint64_t value_type; + +// XXX test (de)serialization + +struct value_index_comparator { + template + bool operator()(Tuple const& a, Tuple const& b) const + { + using boost::get; + // lexicographic, decreasing on value and increasing + // on index + return (get<0>(a) > get<0>(b) || + (get<0>(a) == get<0>(b) && + get<1>(a) < get<1>(b))); + } +}; + +template +void test_topk(std::vector const& v, TopKVector const& topkv, std::string /* test_name */) +{ + BOOST_REQUIRE_EQUAL(v.size(), topkv.size()); + + if (v.empty()) return; + + // test random pairs + const size_t sample_size = 100; + typedef std::pair range_pair; + std::vector pairs_sample; + for (size_t i = 0; i < sample_size; ++i) { + uint64_t a = size_t(rand()) % v.size(); + uint64_t b = a + (size_t(rand()) % (v.size() - a)); + pairs_sample.push_back(range_pair(a, b)); + } + + typedef typename TopKVector::entry_type entry_type; + + size_t k = 10; + + for (size_t i = 0; i < pairs_sample.size(); ++i) { + range_pair r = pairs_sample[i]; + uint64_t a = r.first, b = r.second; + + std::vector expected; + for (uint64_t i = a; i <= b; ++i) { + expected.push_back(entry_type(v[i], i)); + } + std::sort(expected.begin(), expected.end(), value_index_comparator()); // XXX + expected.resize(std::min(expected.size(), k)); + + std::vector found = topkv.topk(a, b, k); + + BOOST_REQUIRE_EQUAL_COLLECTIONS(expected.begin(), expected.end(), + found.begin(), found.end()); + } +} + +BOOST_AUTO_TEST_CASE(topk_vector) +{ + srand(42); + + //typedef succinct::topk_vector > topk_type; + typedef succinct::topk_vector topk_type; + + { + std::vector v; + topk_type t(v); + test_topk(v, t, "Empty vector"); + } + + { + std::vector v(20000); + for (size_t i = 0; i < v.size(); ++i) { + if (i < v.size() / 2) { + v[i] = i; + } else { + v[i] = v.size() - i; + } + } + + { + topk_type t(v); + test_topk(v, t, "Convex values"); + } + } + + { + size_t sizes[] = {2, 4, 512, 514, 8190, 8192, 8194, 16384, 16386, 100000}; + for (size_t i = 0; i < sizeof(sizes) / sizeof(sizes[0]); ++i) { + std::vector v(sizes[i]); + for (size_t i = 0; i < v.size(); ++i) { + v[i] = size_t(rand()) % 1024; + } + + topk_type t(v); + test_topk(v, t, "Random values"); + } + } +} diff --git a/3party/succinct/topk_vector.hpp b/3party/succinct/topk_vector.hpp new file mode 100644 index 0000000000..c64895f414 --- /dev/null +++ b/3party/succinct/topk_vector.hpp @@ -0,0 +1,188 @@ +#pragma once + +#include +#include + +#include +#include + +#include "cartesian_tree.hpp" + +namespace succinct { + + // XXX(ot): implement arbitrary comparator + template + class topk_vector : boost::noncopyable { + public: + typedef Vector vector_type; + typedef typename vector_type::value_type value_type; + typedef boost::tuple entry_type; + typedef std::vector entry_vector_type; + + topk_vector() + {} + + template + topk_vector(Range const& v) + { + cartesian_tree(v, std::greater::type>()) + .swap(m_cartesian_tree); + vector_type(v).swap(m_v); + } + + value_type const + operator[](uint64_t idx) const + { + return m_v[idx]; + } + + uint64_t size() const + { + return m_v.size(); + } + + class enumerator + { + public: + enumerator() + : m_topkv(0) + {} + + bool next() + { + using boost::tie; + if (m_q.empty()) return false; + + value_type cur_mid_val; + uint64_t cur_mid, cur_a, cur_b; + + std::pop_heap(m_q.begin(), m_q.end(), value_index_comparator()); + tie(cur_mid_val, cur_mid, cur_a, cur_b) = m_q.back(); + m_q.pop_back(); + + m_cur = entry_type(cur_mid_val, cur_mid); + + if (cur_mid != cur_a) { + uint64_t m = m_topkv->m_cartesian_tree.rmq(cur_a, cur_mid - 1); + m_q.push_back(queue_element_type(m_topkv->m_v[m], m, cur_a, cur_mid - 1)); + std::push_heap(m_q.begin(), m_q.end(), value_index_comparator()); + } + + if (cur_mid != cur_b) { + uint64_t m = m_topkv->m_cartesian_tree.rmq(cur_mid + 1, cur_b); + m_q.push_back(queue_element_type(m_topkv->m_v[m], m, cur_mid + 1, cur_b)); + std::push_heap(m_q.begin(), m_q.end(), value_index_comparator()); + } + + return true; + } + + entry_type const& value() const + { + return m_cur; + } + + friend class topk_vector; + + void swap(enumerator& other) + { + using std::swap; + swap(m_topkv, other.m_topkv); + swap(m_q, other.m_q); + swap(m_cur, other.m_cur); + } + + private: + + void set(topk_vector const* topkv, uint64_t a, uint64_t b) + { + assert(a <= b); + clear(); + m_topkv = topkv; + + uint64_t m = m_topkv->m_cartesian_tree.rmq(a, b); + m_q.push_back(queue_element_type(m_topkv->m_v[m], m, a, b)); + } + + typedef boost::tuple queue_element_type; + + struct value_index_comparator { + template + bool operator()(Tuple const& a, Tuple const& b) const + { + using boost::get; + // lexicographic, increasing on value and decreasing + // on index + return (get<0>(a) < get<0>(b) || + (get<0>(a) == get<0>(b) && + get<1>(a) > get<1>(b))); + } + }; + + public: + void clear() + { + m_topkv = 0; + m_q.clear(); + } + + private: + topk_vector const* m_topkv; + std::vector m_q; + entry_type m_cur; + }; + + // NOTE this is b inclusive + // XXX switch to [a, b) ? + void get_topk_enumerator(uint64_t a, uint64_t b, enumerator& ret) const + { + ret.set(this, a, b); + } + + enumerator get_topk_enumerator(uint64_t a, uint64_t b) const + { + enumerator ret; + get_topk_enumerator(a, b, ret); + return ret; + } + + entry_vector_type + topk(uint64_t a, uint64_t b, size_t k) const + { + entry_vector_type ret(std::min(size_t(b - a + 1), k)); + enumerator it = get_topk_enumerator(a, b); + + bool hasnext; + for (size_t i = 0; i < ret.size(); ++i) { + hasnext = it.next(); + assert(hasnext); (void)hasnext; + ret[i] = it.value(); + } + + assert(ret.size() == k || !it.next()); + + return ret; + } + + + template + void map(Visitor& visit) + { + visit + (m_v, "m_v") + (m_cartesian_tree, "m_cartesian_tree"); + } + + void swap(topk_vector& other) + { + other.m_v.swap(m_v); + other.m_cartesian_tree.swap(m_cartesian_tree); + } + + protected: + + vector_type m_v; + cartesian_tree m_cartesian_tree; + }; + +} diff --git a/3party/succinct/util.hpp b/3party/succinct/util.hpp new file mode 100644 index 0000000000..c8c079fcf4 --- /dev/null +++ b/3party/succinct/util.hpp @@ -0,0 +1,273 @@ +#pragma once + +#include +#include +#include +#include + +#include + +#include +#include + +namespace succinct { namespace util { + + inline void trim_newline_chars(std::string& s) + { + size_t l = s.size(); + while (l && (s[l-1] == '\r' || + s[l-1] == '\n')) { + --l; + } + s.resize(l); + } + + // this is considerably faster than std::getline + inline bool fast_getline(std::string& line, FILE* input = stdin, bool trim_newline = false) + { + line.clear(); + static const size_t max_buffer = 65536; + char buffer[max_buffer]; + bool done = false; + while (!done) { + if (!fgets(buffer, max_buffer, input)) { + if (!line.size()) { + return false; + } else { + done = true; + } + } + line += buffer; + if (*line.rbegin() == '\n') { + done = true; + } + } + if (trim_newline) { + trim_newline_chars(line); + } + return true; + } + + class line_iterator + : public boost::iterator_facade + { + public: + line_iterator() + : m_file(0) + {} + + explicit line_iterator(FILE* input, bool trim_newline = false) + : m_file(input) + , m_trim_newline(trim_newline) + {} + + private: + friend class boost::iterator_core_access; + + void increment() { + assert(m_file); + if (!fast_getline(m_line, m_file, m_trim_newline)) { + m_file = 0; + } + } + + bool equal(line_iterator const& other) const + { + return this->m_file == other.m_file; + } + + std::string const& dereference() const { + return m_line; + } + + std::string m_line; + FILE* m_file; + bool m_trim_newline; + }; + + typedef std::pair line_range_t; + + inline line_range_t lines(FILE* ifs, bool trim_newline = false) { + return std::make_pair(line_iterator(ifs, trim_newline), line_iterator()); + } + + struct auto_file { + + auto_file(const char* name, const char* mode = "rb") + : m_file(0) + { + m_file = fopen(name, mode); + if(!m_file) { + std::string msg("Unable to open file '"); + msg += name; + msg += "'."; + throw std::invalid_argument(msg); + + } + } + + ~auto_file() + { + if(m_file) { + fclose(m_file); + } + } + + FILE* get() + { + return m_file; + } + + private: + auto_file(); + auto_file( const auto_file & ); + auto_file & operator=( const auto_file & ); + + FILE * m_file; + }; + + typedef std::pair char_range; + + struct identity_adaptor + { + char_range operator()(char_range s) const + { + return s; + } + }; + + struct stl_string_adaptor + { + char_range operator()(std::string const& s) const + { + const uint8_t* buf = reinterpret_cast(s.c_str()); + const uint8_t* end = buf + s.size() + 1; // add the null terminator + return char_range(buf, end); + } + }; + + class buffer_line_iterator + : public boost::iterator_facade + { + public: + buffer_line_iterator() + : m_buffer(0) + , m_end(0) + , m_cur_pos(0) + {} + + buffer_line_iterator(const char* buffer, size_t size) + : m_buffer(buffer) + , m_end(buffer + size) + , m_cur_pos(buffer) + { + increment(); + } + + private: + friend class boost::iterator_core_access; + + void increment() { + assert(m_cur_pos); + if (m_cur_pos >= m_end) { + m_cur_pos = 0; + return; + } + const char* begin = m_cur_pos; + while (m_cur_pos < m_end && *m_cur_pos != '\n') { + ++m_cur_pos; + } + const char* end = m_cur_pos; + ++m_cur_pos; // skip the newline + + if (begin != end && *(end - 1) == '\r') { + --end; + } + m_cur_value = std::string(begin, size_t(end - begin)); + } + + bool equal(buffer_line_iterator const& other) const + { + return m_cur_pos == other.m_cur_pos; } + + std::string const& dereference() const + { + assert(m_cur_pos); + return m_cur_value; + } + + const char* m_buffer; + const char* m_end; + const char* m_cur_pos; + std::string m_cur_value; + }; + + struct mmap_lines + { + typedef buffer_line_iterator iterator; + typedef buffer_line_iterator const_iterator; + + mmap_lines(std::string filename) + : m_map(filename) + {} + + const_iterator begin() const + { + return const_iterator(m_map.data(), m_map.size()); + } + + const_iterator end() const + { + return const_iterator(); + } + + private: + boost::iostreams::mapped_file_source m_map; + }; + + struct input_error : std::invalid_argument + { + input_error(std::string const& what) + : invalid_argument(what) + {} + }; + + template + inline void dispose(T& t) + { + T().swap(t); + } + + inline uint64_t int2nat(int64_t x) + { + if (x < 0) { + return uint64_t(-2 * x - 1); + } else { + return uint64_t(2 * x); + } + } + + inline int64_t nat2int(uint64_t n) + { + if (n % 2) { + return -int64_t((n + 1) / 2); + } else { + return int64_t(n / 2); + } + } + + template + inline IntType1 ceil_div(IntType1 dividend, IntType2 divisor) + { + // XXX(ot): put some static check that IntType1 >= IntType2 + IntType1 d = IntType1(divisor); + return IntType1(dividend + d - 1) / d; + } + +}} diff --git a/3party/succinct/vbyte.hpp b/3party/succinct/vbyte.hpp new file mode 100644 index 0000000000..8f08e10daa --- /dev/null +++ b/3party/succinct/vbyte.hpp @@ -0,0 +1,41 @@ +#pragma once + +#include "broadword.hpp" + +namespace succinct { + + inline size_t vbyte_size(size_t val) + { + unsigned long bits; + if (!broadword::msb(val, bits)) bits = 0; + return util::ceil_div(bits + 1, 7); + } + + template + inline size_t append_vbyte(Vector& v, size_t val) + { + size_t chunks = vbyte_size(val); + for (size_t b = chunks - 1; b + 1 > 0; --b) { + uint8_t chunk = (val >> (b * 7)) & 0x7F; + chunk |= b ? 0x80 : 0; + v.push_back(chunk); + } + return chunks; + } + + template + inline size_t decode_vbyte(Vector const& v, size_t offset, size_t& val) + { + size_t pos = offset; + val = 0; + uint8_t chunk; + do { + chunk = v[pos++]; + val <<= 7; + val |= chunk & 0x7F; + } while (chunk & 0x80); + + return pos - offset; + } + +}