From ac65326242af579d6e1a7bd71730f1c78c8bde9b Mon Sep 17 00:00:00 2001 From: David Robillard Date: Sun, 19 Aug 2018 18:22:26 +0200 Subject: Reorganize headers to make a clear public/private distinction --- chilbert/BitVecIndex.hpp | 49 ---- chilbert/BitVecIterator.hpp | 98 -------- chilbert/BitVecMask.hpp | 72 ------ chilbert/BoundedBitVec.hpp | 25 +- chilbert/DynamicBitVec.hpp | 16 +- chilbert/GrayCodeRank.hpp | 180 ------------- chilbert/Hilbert.hpp | 102 -------- chilbert/Hilbert.ipp | 501 ------------------------------------- chilbert/MultiBitVec.hpp | 379 ---------------------------- chilbert/Operations.hpp | 166 ------------ chilbert/Operators.hpp | 94 ------- chilbert/SmallBitVec.hpp | 13 +- chilbert/StaticBitVec.hpp | 20 +- chilbert/Traits.hpp | 37 --- chilbert/chilbert.hpp | 102 ++++++++ chilbert/chilbert.ipp | 501 +++++++++++++++++++++++++++++++++++++ chilbert/detail/BitVecIndex.hpp | 51 ++++ chilbert/detail/BitVecIterator.hpp | 100 ++++++++ chilbert/detail/BitVecMask.hpp | 74 ++++++ chilbert/detail/MultiBitVec.hpp | 383 ++++++++++++++++++++++++++++ chilbert/detail/gray_code_rank.hpp | 182 ++++++++++++++ chilbert/detail/operations.hpp | 168 +++++++++++++ chilbert/detail/traits.hpp | 39 +++ chilbert/operators.hpp | 96 +++++++ test/test_bitvec.cpp | 10 +- test/test_gray_code_rank.cpp | 19 +- test/test_hilbert.cpp | 2 +- 27 files changed, 1759 insertions(+), 1720 deletions(-) delete mode 100644 chilbert/BitVecIndex.hpp delete mode 100644 chilbert/BitVecIterator.hpp delete mode 100644 chilbert/BitVecMask.hpp delete mode 100644 chilbert/GrayCodeRank.hpp delete mode 100644 chilbert/Hilbert.hpp delete mode 100644 chilbert/Hilbert.ipp delete mode 100644 chilbert/MultiBitVec.hpp delete mode 100644 chilbert/Operations.hpp delete mode 100644 chilbert/Operators.hpp delete mode 100644 chilbert/Traits.hpp create mode 100644 chilbert/chilbert.hpp create mode 100644 chilbert/chilbert.ipp create mode 100644 chilbert/detail/BitVecIndex.hpp create mode 100644 chilbert/detail/BitVecIterator.hpp create mode 100644 chilbert/detail/BitVecMask.hpp create mode 100644 chilbert/detail/MultiBitVec.hpp create mode 100644 chilbert/detail/gray_code_rank.hpp create mode 100644 chilbert/detail/operations.hpp create mode 100644 chilbert/detail/traits.hpp create mode 100644 chilbert/operators.hpp diff --git a/chilbert/BitVecIndex.hpp b/chilbert/BitVecIndex.hpp deleted file mode 100644 index 4038105..0000000 --- a/chilbert/BitVecIndex.hpp +++ /dev/null @@ -1,49 +0,0 @@ -/* - Copyright (C) 2018 David Robillard - Copyright (C) 2006-2007 Chris Hamilton - - This program is free software: you can redistribute it and/or modify it under - the terms of the GNU General Public License as published by the Free Software - Foundation, either version 2 of the License, or (at your option) any later - version. - - This program is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - FOR A PARTICULAR PURPOSE. See the GNU General Public License for more - details. - - You should have received a copy of the GNU General Public License along with - this program. If not, see . -*/ - -#ifndef CHILBERT_BITVECINDEX_HPP -#define CHILBERT_BITVECINDEX_HPP - -#include -#include -#include - -namespace chilbert { - -/// Index into a multi-rack bit vector -template -struct BitVecIndex -{ - using Rack = typename BitVec::Rack; - - static constexpr size_t bits_per_rack = sizeof(Rack) * CHAR_BIT; - - BitVecIndex(const size_t bits) - : rack{bits / bits_per_rack} - , bit{bits - rack * bits_per_rack} - { - assert(bit < bits_per_rack); - } - - size_t rack; - size_t bit; -}; - -} // namespace chilbert - -#endif diff --git a/chilbert/BitVecIterator.hpp b/chilbert/BitVecIterator.hpp deleted file mode 100644 index bbafd42..0000000 --- a/chilbert/BitVecIterator.hpp +++ /dev/null @@ -1,98 +0,0 @@ -/* - Copyright (C) 2018 David Robillard - Copyright (C) 2006-2007 Chris Hamilton - - This program is free software: you can redistribute it and/or modify it under - the terms of the GNU General Public License as published by the Free Software - Foundation, either version 2 of the License, or (at your option) any later - version. - - This program is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - FOR A PARTICULAR PURPOSE. See the GNU General Public License for more - details. - - You should have received a copy of the GNU General Public License along with - this program. If not, see . -*/ - -#ifndef CHILBERT_BITVECITERATOR_HPP -#define CHILBERT_BITVECITERATOR_HPP - -#include "chilbert/BitVecMask.hpp" - -#include - -namespace chilbert { - -template -class BitVecIteratorBase : public BitVecMask -{ -public: - using Mask = typename BitVec::Mask; - - BitVecIteratorBase& operator++() - { - Mask::operator++(); - return *this; - } - - BitVecIteratorBase& operator--() - { - Mask::operator--(); - return *this; - } - - bool operator==(const BitVecIteratorBase& rhs) const - { - return m_vec == rhs.m_vec && Mask::operator==(rhs); - } - - bool operator!=(const BitVecIteratorBase& rhs) const - { - return !operator==(rhs); - } - - bool operator*() const { return m_vec->test(*this); } - -protected: - BitVecIteratorBase(BitVec& vec, const size_t index) - : Mask{index} - , m_vec{&vec} - { - } - - BitVec* m_vec; -}; - -template -class BitVecIterator : public BitVecIteratorBase -{ -public: - void set() { this->m_vec->set(*this); } - void reset() { this->m_vec->reset(*this); } - -private: - friend BitVec; - - BitVecIterator(BitVec& vec, const size_t index) - : BitVecIteratorBase{vec, index} - { - } -}; - -template -class ConstBitVecIterator : public BitVecIteratorBase -{ -private: - friend BitVec; - - ConstBitVecIterator(const BitVec& vec, const size_t index) - : BitVecIteratorBase{vec, index} - { - } -}; - -} // namespace chilbert - -#endif diff --git a/chilbert/BitVecMask.hpp b/chilbert/BitVecMask.hpp deleted file mode 100644 index af98894..0000000 --- a/chilbert/BitVecMask.hpp +++ /dev/null @@ -1,72 +0,0 @@ -/* - Copyright (C) 2018 David Robillard - Copyright (C) 2006-2007 Chris Hamilton - - This program is free software: you can redistribute it and/or modify it under - the terms of the GNU General Public License as published by the Free Software - Foundation, either version 2 of the License, or (at your option) any later - version. - - This program is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - FOR A PARTICULAR PURPOSE. See the GNU General Public License for more - details. - - You should have received a copy of the GNU General Public License along with - this program. If not, see . -*/ - -#ifndef CHILBERT_BITVECMASK_HPP -#define CHILBERT_BITVECMASK_HPP - -#include -#include - -namespace chilbert { - -/** Mask for a bit that can be incremented like an index. - * - * This enables fast iteration while setting or resetting bits since it is - * internally stored as a mask which avoids repeated index math. - */ -template -struct BitVecMask -{ - static constexpr size_t bits_per_rack = sizeof(Rack) * CHAR_BIT; - - BitVecMask(const size_t index) - : rack{index / bits_per_rack} - , mask{Rack{1} << (index - rack * bits_per_rack)} - { - } - - void operator++() - { - if ((mask <<= 1) == 0) { - mask = 1; - ++rack; - } - } - - void operator--() - { - if ((mask >>= 1) == 0) { - mask = Rack{1} << (bits_per_rack - 1); - --rack; - } - } - - bool operator==(const BitVecMask& rhs) const - { - return rack == rhs.rack && mask == rhs.mask; - } - - bool operator!=(const BitVecMask& rhs) const { return !operator==(rhs); } - - size_t rack; - Rack mask; -}; - -} // namespace chilbert - -#endif diff --git a/chilbert/BoundedBitVec.hpp b/chilbert/BoundedBitVec.hpp index 79912cc..fa414ca 100644 --- a/chilbert/BoundedBitVec.hpp +++ b/chilbert/BoundedBitVec.hpp @@ -19,11 +19,11 @@ #ifndef CHILBERT_BOUNDEDBITVEC_HPP #define CHILBERT_BOUNDEDBITVEC_HPP -#include "chilbert/BitVecIndex.hpp" -#include "chilbert/BitVecIterator.hpp" -#include "chilbert/BitVecMask.hpp" -#include "chilbert/MultiBitVec.hpp" -#include "chilbert/Operations.hpp" +#include "chilbert/detail/BitVecIndex.hpp" +#include "chilbert/detail/BitVecIterator.hpp" +#include "chilbert/detail/BitVecMask.hpp" +#include "chilbert/detail/MultiBitVec.hpp" +#include "chilbert/detail/operations.hpp" #include #include @@ -42,12 +42,12 @@ namespace chilbert { * @tparam MaxN Maximum number of bits. */ template -class BoundedBitVec : public MultiBitVec> +class BoundedBitVec : public detail::MultiBitVec> { public: - using Rack = typename MultiBitVec>::Rack; + using Rack = typename detail::MultiBitVec>::Rack; - using MultiBitVec>::bits_per_rack; + using detail::MultiBitVec>::bits_per_rack; BoundedBitVec() = default; @@ -95,6 +95,8 @@ private: size_t m_size; }; +namespace detail { + template struct is_bitvec> { @@ -105,16 +107,19 @@ template void gray_code(BoundedBitVec& value) { - gray_code(static_cast>&>(value)); + gray_code(static_cast>&>(value)); } template void gray_code_inv(BoundedBitVec& value) { - gray_code_inv(static_cast>&>(value)); + gray_code_inv( + static_cast>&>(value)); } +} // namespace detail + } // namespace chilbert #endif diff --git a/chilbert/DynamicBitVec.hpp b/chilbert/DynamicBitVec.hpp index 7372629..722b689 100644 --- a/chilbert/DynamicBitVec.hpp +++ b/chilbert/DynamicBitVec.hpp @@ -19,11 +19,11 @@ #ifndef CHILBERT_DYNAMICBITVEC_HPP #define CHILBERT_DYNAMICBITVEC_HPP -#include "chilbert/BitVecIndex.hpp" -#include "chilbert/BitVecIterator.hpp" -#include "chilbert/BitVecMask.hpp" -#include "chilbert/MultiBitVec.hpp" -#include "chilbert/Operations.hpp" +#include "chilbert/detail/BitVecIndex.hpp" +#include "chilbert/detail/BitVecIterator.hpp" +#include "chilbert/detail/BitVecMask.hpp" +#include "chilbert/detail/MultiBitVec.hpp" +#include "chilbert/detail/operations.hpp" #include #include @@ -38,7 +38,7 @@ namespace chilbert { * This uses dynamic allocation internally and can be constructed with any * size, assuming sufficient memory is available. */ -class DynamicBitVec : public MultiBitVec +class DynamicBitVec : public detail::MultiBitVec { public: struct RacksDeleter @@ -128,6 +128,8 @@ private: size_t m_size; }; +namespace detail { + template <> struct is_bitvec { @@ -148,6 +150,8 @@ gray_code_inv(DynamicBitVec& value) gray_code_inv(static_cast&>(value)); } +} // namespace detail + } // namespace chilbert #endif diff --git a/chilbert/GrayCodeRank.hpp b/chilbert/GrayCodeRank.hpp deleted file mode 100644 index d7fbc77..0000000 --- a/chilbert/GrayCodeRank.hpp +++ /dev/null @@ -1,180 +0,0 @@ -/* - Copyright (C) 2018 David Robillard - Copyright (C) 2006-2007 Chris Hamilton - - This program is free software: you can redistribute it and/or modify it under - the terms of the GNU General Public License as published by the Free Software - Foundation, either version 2 of the License, or (at your option) any later - version. - - This program is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - FOR A PARTICULAR PURPOSE. See the GNU General Public License for more - details. - - You should have received a copy of the GNU General Public License along with - this program. If not, see . -*/ - -#ifndef CHILBERT_GRAYCODERANK_HPP -#define CHILBERT_GRAYCODERANK_HPP - -#include -#include -#include - -namespace chilbert { - -// This is the bulk of the cost in calculating -// a Compact Hilbert Index. It compresses a previously -// calculated index when provided the rotation -// at each level of precision. -template -inline void -compact_index(const size_t* const ms, - const size_t* const ds, - const size_t n, - const size_t m, - H& h, - HC& hc) -{ - assert(h.size() >= n * m); - assert(hc.size() >= std::accumulate(ms, ms + n, size_t(0))); - - reset_bits(hc); - - auto hm = h.mask(0); - auto hcm = hc.mask(0); - - // Run through the levels of precision - for (size_t i = 0; i < m; i++) { - // Run through the dimensions - size_t j = ds[i]; - do { - // This dimension contributes a bit? - if (ms[j] > i) { - if (h.test(hm)) { - hc.set(hcm); - } - ++hcm; - } - - if (++j == n) { - j = 0; - } - ++hm; - } while (j != ds[i]); - } -} - -template -inline void -gray_code_rank(const I& mask, const I& gi, const size_t n, I& r) -{ - assert(mask.size() == n); - assert(gi.size() == n); - assert(r.size() == n); - - r.reset(); - - auto mi = mask.begin(); - auto gii = gi.begin(); - auto ri = r.begin(); - - for (size_t i = 0; i < n; ++i, ++mi, ++gii) { - if (*mi) { - if (*gii) { - ri.set(); - } - ++ri; - } - } -} - -template -inline void -gray_code_rank_inv(const I& mask, - const I& ptrn, - const I& r, - const size_t n, - const size_t b, - I& g, - I& gi) -{ - assert(mask.size() == n); - assert(ptrn.size() == n); - assert(r.size() == n); - assert(g.size() == n); - assert(gi.size() == n); - - g.reset(); - gi.reset(); - - auto m = mask.mask(n - 1); - auto ri = r.begin(b - 1); - - typename I::Rack gi0 = 0; - typename I::Rack gi1 = 0; - typename I::Rack g0 = 0; - - for (size_t i = 0; i < n; ++i) { - if (mask.test(m)) { // Unconstrained bit - gi1 = gi0; - gi0 = *ri; - g0 = gi0 ^ gi1; - if (gi0) { - gi.set(m); - } - if (g0) { - g.set(m); - } - --ri; - } else { // Constrained bit - g0 = (ptrn.test(m) > 0); - gi1 = gi0; - gi0 = g0 ^ gi1; - if (gi0) { - gi.set(m); - } - if (g0) { - g.set(m); - } - } - - --m; - } -} - -template -inline void -extract_mask(const size_t* const ms, - const size_t n, - const size_t d, - const size_t i, - I& mask, - size_t& b) -{ - assert(0 <= d && d < n); - assert(mask.size() == n); - - mask.reset(); - b = 0; - - auto mi = mask.begin(); - size_t j = d; // #D j = (d==n-1) ? 0 : d+1; - do { - if (ms[j] > i) { - mi.set(); - ++b; - } - - ++mi; - if (++j == n) { - j = 0; - } - } while (j != d); -} - -} // namespace chilbert - -#endif diff --git a/chilbert/Hilbert.hpp b/chilbert/Hilbert.hpp deleted file mode 100644 index 8db5b1c..0000000 --- a/chilbert/Hilbert.hpp +++ /dev/null @@ -1,102 +0,0 @@ -/* - Copyright (C) 2018 David Robillard - Copyright (C) 2006-2007 Chris Hamilton - - This program is free software: you can redistribute it and/or modify it under - the terms of the GNU General Public License as published by the Free Software - Foundation, either version 2 of the License, or (at your option) any later - version. - - This program is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - FOR A PARTICULAR PURPOSE. See the GNU General Public License for more - details. - - You should have received a copy of the GNU General Public License along with - this program. If not, see . -*/ - -#ifndef CHILBERT_HILBERT_HPP -#define CHILBERT_HILBERT_HPP - -#include - -namespace chilbert { - -/** Map the point `p` to a Hilbert Index. - * - * @tparam P Type used to represent a value in a point dimension. - * @tparam H Hilbert Index. - * - * @param p Point with `n` dimensions. - * @param m Precision of each dimension in bits. - * @param n Number of dimensions. - * @param[out] h Hilbert Index. - */ -template -inline void coords_to_index(const P* const p, - const size_t m, - const size_t n, - H& h); - -/** Map the Hilbert Index `p` to a point. - * - * @tparam P Type used to represent a value in a point dimension. - * @tparam H Hilbert Index. - * - * @param[out] p Point with `n` dimensions. - * @param m Precision of each dimension in bits. - * @param n Number of dimensions. - * @param h Hilbert Index. - */ -template -inline void index_to_coords(P* const p, - const size_t m, - const size_t n, - const H& h); - -/** Map the point `p` to a Compact Hilbert Index. - * - * @tparam P Type used to represent a value in a point dimension. - * @tparam H Compact Hilbert Index. - * - * @param p Point with `n` dimensions. - * @param ms Array of `n` precision values for each dimension in bits. - * @param n Number of dimensions. - * @param[out] hc Compact Hilbert Index. - * @param M Optional net precision (sum of `ms`), the size of `hc` in bits. - * @param m Optional largest precision in `m`. - */ -template -inline void coords_to_compact_index(const P* const p, - const size_t* const ms, - const size_t n, - H& hc, - const size_t M = 0, - const size_t m = 0); - -/** Map the Compact Hilbert Index `hc` to a point. - * - * @tparam P Type used to represent a value in a point dimension. - * @tparam H Compact Hilbert Index. - * - * @param[out] p Point with `n` dimensions. - * @param ms Array of `n` precision values for each dimension in bits. - * @param n Number of dimensions. - * @param hc Compact Hilbert Index. - * @param M Optional net precision (sum of `ms`), the size of `hc` in bits. - * @param m Optional largest precision in `m`. - */ -template -inline void compact_index_to_coords(P* const p, - const size_t* const ms, - const size_t n, - const H& hc, - const size_t M = 0, - const size_t m = 0); - -} // namespace chilbert - -#include "chilbert/Hilbert.ipp" - -#endif diff --git a/chilbert/Hilbert.ipp b/chilbert/Hilbert.ipp deleted file mode 100644 index 837375e..0000000 --- a/chilbert/Hilbert.ipp +++ /dev/null @@ -1,501 +0,0 @@ -/* - Copyright (C) 2018 David Robillard - Copyright (C) 2006-2007 Chris Hamilton - - This program is free software: you can redistribute it and/or modify it under - the terms of the GNU General Public License as published by the Free Software - Foundation, either version 2 of the License, or (at your option) any later - version. - - This program is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - FOR A PARTICULAR PURPOSE. See the GNU General Public License for more - details. - - You should have received a copy of the GNU General Public License along with - this program. If not, see . -*/ - -#ifndef CHILBERT_ALGORITHM_HPP -#define CHILBERT_ALGORITHM_HPP - -#include "chilbert/BoundedBitVec.hpp" -#include "chilbert/DynamicBitVec.hpp" -#include "chilbert/GrayCodeRank.hpp" -#include "chilbert/Hilbert.hpp" -#include "chilbert/SmallBitVec.hpp" -#include "chilbert/StaticBitVec.hpp" - -#include -#include -#include -#include - -namespace chilbert { -namespace detail { - -/** The dimension across which the Hilbert curve travels principally. - * - * D0 is d + 1, where d is the actual dimension you want to walk across. - * MUST HAVE 0 <= D0 < n. - */ -static constexpr size_t D0 = 1; - -template -size_t -num_bits(const T&, std::enable_if_t::value>* = nullptr) -{ - return sizeof(T) * CHAR_BIT; -} - -template -size_t -num_bits(const T& vec, - std::enable_if_t::value || - std::is_same::value>* = nullptr) -{ - return vec.size(); -} - -template -size_t -num_bits(const StaticBitVec&, void* = nullptr) -{ - return N; -} - -template -size_t -num_bits(const BoundedBitVec& vec, void* = nullptr) -{ - return vec.size(); -} - -/** Copy a range of bits from one field to the start of another. - * - * @param h Source bit field - * @param n Number of bits - * @param i Start bit index in source - * @param[out] w Destination bit field - */ -template -inline void -get_bits(const H& h, const size_t n, const size_t i, I& w) -{ - for (size_t j = 0; j < n; j++) { - set_bit(w, j, test_bit(h, i + j)); - } -} - -/** Set a range of bits in one field from the start of another. - * - * @param[out] h Destination bit field - * @param n Number of bits - * @param i Start bit index in destination - * @param w Source bit field - */ -template -inline void -set_bits(H& h, const size_t n, const size_t i, const I& w) -{ - for (size_t j = 0; j < n; j++) { - set_bit(h, i + j, test_bit(w, j)); - } -} - -/** Set the leading bits in `l` to bit `i` from the corresponding value in `p`. - * - * @param p Point. - * @param n Number of bits to set. - * @param i Bit position to copy from values in `p`. - * @param[out] l Output bits. - */ -template -inline void -get_location(const P* const p, const size_t n, const size_t i, I& l) -{ - for (size_t j = 0; j < n; ++j) { - set_bit(l, j, test_bit(p[j], i)); - } -} - -/** Set bit `i` in values in `p` to the corresponding leadings bits in `l`. - * - * @param[out] p Point. - * @param n Number of bits to set. - * @param i Bit position to set in values in `p`. - * @param l Input bits. - */ -template -inline void -set_location(P* const p, const size_t n, const size_t i, const I& l) -{ - for (size_t j = 0; j < n; j++) { - set_bit(p[j], i, test_bit(l, j)); - } -} - -// 'Transforms' a point. -template -inline void -transform(const I& e, const size_t d, const size_t n, I& a) -{ - assert(a.size() == n); - a ^= e; - a.rotr(d); -} - -// Inverse 'transforms' a point. -template -inline void -transform_inv(const I& e, const size_t d, const size_t n, I& a) -{ - assert(a.size() == n); - a.rotl(d); - a ^= e; -} - -// Update for method 1 (GrayCodeInv in the loop) -template -inline void -update1(const I& l, const I& t, const I& w, const size_t n, I& e, size_t& d) -{ - assert(0 <= d && d < n); - e = l; - e.flip(d); //#D d == n-1 ? 0 : d+1 ); - - // Update direction - d += 1 + t.find_first(); - if (d >= n) { - d -= n; - } - if (d >= n) { - d -= n; - } - assert(0 <= d && d < n); - - if (!w.test(0)) { - e.flip(d == 0 ? n - 1 : d - 1); //#D d ); - } -} - -// Update for method 2 (GrayCodeInv out of loop) -template -inline void -update2(const I& l, const I& t, const size_t n, I& e, size_t& d) -{ - assert(0 <= d && d < n); - e = l; - e.flip(d); //#D d == n-1 ? 0 : d+1 ); - - // Update direction - d += 1 + t.find_first(); - if (d >= n) { - d -= n; - } - if (d >= n) { - d -= n; - } - assert(0 <= d && d < n); -} - -template -inline void -coords_to_index(const P* const p, - const size_t m, - const size_t n, - H& h, - I&& scratch, - size_t* const ds = nullptr) -{ - I e{std::move(scratch)}; - I l{e}; - I t{e}; - I w{e}; - - // Initialize - e.reset(); - l.reset(); - reset_bits(h); - - // Work from MSB to LSB - size_t d = D0; - size_t ho = m * n; - for (intptr_t i = static_cast(m - 1); i >= 0; i--) { - if (ds) { - ds[i] = d; - } - - // Get corner of sub-hypercube where point lies. - get_location(p, n, static_cast(i), l); - - // Mirror and reflect the location. - // t = T_{(e,d)}(l) - t = l; - transform(e, d, n, t); - - w = t; - if (static_cast(i) < m - 1) { - w.flip(n - 1); - } - - // Concatenate to the index. - ho -= n; - set_bits(h, n, ho, w); - - // Update the entry point and direction. - update2(l, t, n, e, d); - } - - gray_code_inv(h); -} - -template -inline void -index_to_coords(P* p, const size_t m, const size_t n, const H& h, I&& scratch) -{ - I e{std::move(scratch)}; - I l{e}; - I t{e}; - I w{e}; - - // Initialize - e.reset(); - l.reset(); - for (size_t j = 0; j < n; j++) { - reset_bits(p[j]); - } - - // Work from MSB to LSB - size_t d = D0; - size_t ho = m * n; - for (intptr_t i = static_cast(m - 1); i >= 0; i--) { - // Get the Hilbert index bits - ho -= n; - get_bits(h, n, ho, w); - - // t = GrayCode(w) - t = w; - gray_code(t); - - // Reverse the transform - // l = T^{-1}_{(e,d)}(t) - l = t; - transform_inv(e, d, n, l); - - // Distribute these bits - // to the coordinates. - set_location(p, n, static_cast(i), l); - - // Update the entry point and direction. - update1(l, t, w, n, e, d); - } -} - -template -inline void -coords_to_compact_index(const P* const p, - const size_t* const ms, - const size_t n, - HC& hc, - I&& scratch, - size_t M = 0, - size_t m = 0) -{ - // Get total precision and max precision if not supplied - if (M == 0 || m == 0) { - M = m = 0; - for (size_t i = 0; i < n; i++) { - assert(num_bits(p[i]) >= ms[i]); - if (ms[i] > m) { - m = ms[i]; - } - M += ms[i]; - } - } - - const size_t mn = m * n; - - assert(num_bits(hc) >= M); - - // If we could avoid allocation altogether (ie: have a - // fixed buffer allocated on the stack) then this increases - // speed by a bit (4% when n=4, m=20) - size_t* const ds = new size_t[m]; - - if (mn > SmallBitVec::bits_per_rack) { - DynamicBitVec h(mn); - detail::coords_to_index( - p, m, n, h, std::move(scratch), ds); - compact_index(ms, ds, n, m, h, hc); - } else { - SmallBitVec h(mn); - detail::coords_to_index( - p, m, n, h, std::move(scratch), ds); - compact_index(ms, ds, n, m, h, hc); - } - - delete[] ds; -} - -template -inline void -compact_index_to_coords(P* const p, - const size_t* ms, - const size_t n, - const HC& hc, - I&& scratch, - size_t M, - size_t m) -{ - I e{std::move(scratch)}; - I l{e}; - I t{e}; - I w{e}; - I r{e}; - I mask{e}; - I ptrn{e}; - - // Get total precision and max precision - // if not supplied - if (M == 0 || m == 0) { - M = m = 0; - for (size_t i = 0; i < n; i++) { - if (ms[i] > m) { - m = ms[i]; - } - M += ms[i]; - } - } - - assert(num_bits(hc) >= M); - assert(num_bits(p[0]) >= m); - - // Initialize - e.reset(); - l.reset(); - for (size_t j = 0; j < n; j++) { - reset_bits(p[j]); - } - - // Work from MSB to LSB - size_t d = D0; - - for (intptr_t i = static_cast(m - 1); i >= 0; i--) { - // Get the mask and ptrn - size_t b = 0; - extract_mask(ms, n, d, static_cast(i), mask, b); - ptrn = e; - assert(ptrn.size() == n); - ptrn.rotr(d); - - // Get the Hilbert index bits - M -= b; - r.reset(); // GetBits doesn't do this - get_bits(hc, b, M, r); - - // w = GrayCodeRankInv(r) - // t = GrayCode(w) - gray_code_rank_inv(mask, ptrn, r, n, b, t, w); - - // Reverse the transform - // l = T^{-1}_{(e,d)}(t) - l = t; - transform_inv(e, d, n, l); - - // Distribute these bits - // to the coordinates. - set_location(p, n, static_cast(i), l); - - // Update the entry point and direction. - update1(l, t, w, n, e, d); - } -} - -} // namespace detail - -template -inline void -coords_to_index(const P* const p, const size_t m, const size_t n, H& h) -{ - assert(detail::num_bits(h) >= n * m); - assert(detail::num_bits(p[0]) >= m); - - if (n <= SmallBitVec::bits_per_rack) { - // Intermediate variables will fit in fixed width - detail::coords_to_index(p, m, n, h, SmallBitVec(n)); - } else { - // Otherwise, they must be DynamicBitVecs - detail::coords_to_index( - p, m, n, h, DynamicBitVec(n)); - } -} - -template -inline void -index_to_coords(P* const p, const size_t m, const size_t n, const H& h) -{ - assert(m > 0); - assert(n > 0); - assert(detail::num_bits(h) >= n * m); - assert(detail::num_bits(p[0]) >= m); - - if (n <= SmallBitVec::bits_per_rack) { - // Intermediate variables will fit in fixed width - detail::index_to_coords(p, m, n, h, SmallBitVec(n)); - } else { - // Otherwise, they must be DynamicBitVecs - detail::index_to_coords( - p, m, n, h, DynamicBitVec(n)); - } -} - -template -inline void -coords_to_compact_index(const P* const p, - const size_t* const ms, - size_t n, - HC& hc, - const size_t M, - const size_t m) -{ - assert(hc.size() >= std::accumulate(ms, ms + n, size_t(0))); - - if (n <= SmallBitVec::bits_per_rack) { - // Intermediate variables will fit in fixed width - detail::coords_to_compact_index( - p, ms, n, hc, SmallBitVec(n), M, m); - } else { - // Otherwise, they must be DynamicBitVecs - detail::coords_to_compact_index( - p, ms, n, hc, DynamicBitVec(n), M, m); - } -} - -template -inline void -compact_index_to_coords(P* const p, - const size_t* ms, - const size_t n, - const HC& hc, - const size_t M, - const size_t m) -{ - assert(hc.size() >= std::accumulate(ms, ms + n, size_t(0))); - - if (n <= SmallBitVec::bits_per_rack) { - // Intermediate variables will fit in fixed width - SmallBitVec scratch(n); - detail::compact_index_to_coords( - p, ms, n, hc, std::move(scratch), M, m); - } else { - // Otherwise, they must be DynamicBitVecs - DynamicBitVec scratch(n); - detail::compact_index_to_coords( - p, ms, n, hc, std::move(scratch), M, m); - } -} - -} // namespace chilbert - -#endif diff --git a/chilbert/MultiBitVec.hpp b/chilbert/MultiBitVec.hpp deleted file mode 100644 index 1401fd7..0000000 --- a/chilbert/MultiBitVec.hpp +++ /dev/null @@ -1,379 +0,0 @@ -/* - Copyright (C) 2018 David Robillard - Copyright (C) 2006-2007 Chris Hamilton - - This program is free software: you can redistribute it and/or modify it under - the terms of the GNU General Public License as published by the Free Software - Foundation, either version 2 of the License, or (at your option) any later - version. - - This program is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - FOR A PARTICULAR PURPOSE. See the GNU General Public License for more - details. - - You should have received a copy of the GNU General Public License along with - this program. If not, see . -*/ - -#ifndef CHILBERT_MULTIBITVEC_HPP -#define CHILBERT_MULTIBITVEC_HPP - -#include "chilbert/BitVecIterator.hpp" -#include "chilbert/Operations.hpp" - -#include -#include -#include -#include - -namespace chilbert { - -template -class MultiBitVec -{ -public: - using Rack = uintptr_t; - using Mask = BitVecMask; - - using iterator = BitVecIterator>; - using const_iterator = ConstBitVecIterator>; - - static constexpr size_t bits_per_rack = sizeof(Rack) * CHAR_BIT; - - /// Return the value of the bit covered by `mask` - bool test(const Mask mask) const { return rack(mask.rack) & mask.mask; } - - /// Return the value of the `index`th bit - bool test(const size_t index) const { return test(mask(index)); } - - /// Set all bits to one - Derived& set() - { - if (size()) { - memset(data(), 0xFF, data_size()); - self()->truncate(); - } - - return *self(); - } - - /// Set the bit covered by `mask` to 1 - Derived& set(const Mask mask) - { - rack(mask.rack) |= mask.mask; - return *self(); - } - - /// Set the `index`th bit to 1 - Derived& set(const size_t index) { return set(mask(index)); } - - /// Set the bit covered by `mask` to `value` - Derived& set(const Mask mask, const bool value) - { - auto& r = rack(mask.rack); - r ^= (-Rack{value} ^ r) & mask.mask; - return *self(); - } - - /// Set the `index`th bit to `value` - Derived& set(const size_t index, const bool value) - { - return set(mask(index), value); - } - - /// Set all bits to zero - Derived& reset() - { - memset(data(), 0, data_size()); - return *self(); - } - - /// Reset the bit covered by `mask` to 0 - Derived& reset(const Mask mask) - { - rack(mask.rack) &= ~mask.mask; - return *self(); - } - - /// Reset the `index`th bit to 0 - Derived& reset(const size_t index) { return reset(mask(index)); } - - /// Flip all bits (one's complement) - Derived& flip() - { - for (size_t i = 0; i < num_racks(); ++i) { - rack(i) = ~rack(i); - } - return *self(); - } - - /// Flip the value of the bit covered by `mask` - Derived& flip(const Mask mask) - { - rack(mask.rack) ^= mask.mask; - return *self(); - } - - /// Flip the value of the `index`th bit - Derived& flip(const size_t index) { return flip(mask(index)); } - - /// Clear any bits in storage outside the valid range if necessary - void truncate() - { - if (const auto pad = num_racks() * bits_per_rack - size()) { - rack(num_racks() - 1) &= ~Rack{0} >> pad; - } - } - - /// Right-rotate by `bits` positions - Derived& rotr(const size_t bits) - { - assert(bits <= size()); - if (bits == 0 || bits == size()) { - return *self(); - } - - Derived t1(*self()); - *self() >>= bits; - t1 <<= (size() - bits); - *self() |= t1; - - truncate(); - return *self(); - } - - /// Left-rotate by `bits` positions - Derived& rotl(const size_t bits) - { - assert(bits <= size()); - if (bits == 0 || bits == size()) { - return *self(); - } - - Derived t1(*self()); - *self() <<= bits; - t1 >>= (size() - bits); - *self() |= t1; - - truncate(); - return *self(); - } - - /// Return true iff all bits are zero - bool none() const - { - for (size_t i = 0; i < num_racks(); ++i) { - if (rack(i)) { - return false; - } - } - return true; - } - - /// Return 1 + the index of the first set bit, or 0 if there are none - size_t find_first() const - { - for (size_t i = 0; i < num_racks(); ++i) { - const int j = chilbert::find_first(rack(i)); - if (j) { - return (i * bits_per_rack) + static_cast(j); - } - } - return 0; - } - - /// Return the number of set bits - size_t count() const - { - size_t c = 0; - for (size_t i = 0; i < num_racks(); ++i) { - c += static_cast(pop_count(rack(i))); - } - return c; - } - - /// Return a mask that covers the bit with index `i` - Mask mask(const size_t i = 0) const - { - assert(i <= size()); - return Mask{i}; - } - - bool operator==(const Derived& vec) const - { - return (num_racks() == vec.num_racks() && - (num_racks() == 0 || - !memcmp(data(), vec.data(), num_racks() * sizeof(Rack)))); - } - - bool operator!=(const Derived& vec) const { return !(*this == vec); } - - bool operator<(const Derived& vec) const - { - assert(size() == vec.size()); - - for (size_t ri = 0; ri < num_racks(); ++ri) { - size_t i = num_racks() - ri - 1; - if (rack(i) < vec.rack(i)) { - return true; - } else if (rack(i) > vec.rack(i)) { - return false; - } - } - return false; - } - - Derived& operator&=(const Derived& vec) - { - for (size_t i = 0; i < std::min(num_racks(), vec.num_racks()); ++i) { - rack(i) &= vec.rack(i); - } - - return *self(); - } - - Derived& operator|=(const Derived& vec) - { - for (size_t i = 0; i < std::min(num_racks(), vec.num_racks()); ++i) { - rack(i) |= vec.rack(i); - } - - return *self(); - } - - Derived& operator^=(const Derived& vec) - { - for (size_t i = 0; i < std::min(num_racks(), vec.num_racks()); ++i) { - rack(i) ^= vec.rack(i); - } - - return *self(); - } - - Derived& operator<<=(const size_t bits) - { - if (bits == 0) { - return *self(); - } else if (bits >= size()) { - reset(); - return *self(); - } - - const Index index{bits}; - - if (index.rack > 0) { - // Shift entire racks - for (size_t i = num_racks() - 1; i >= index.rack; --i) { - rack(i) = rack(i - index.rack); - } - for (size_t i = 0; i < index.rack; ++i) { - rack(i) = 0; - } - } - - if (index.bit > 0) { - // Shift bits within racks - size_t bi = bits_per_rack - index.bit; - size_t i; - for (i = num_racks() - 1; i >= index.rack + 1; --i) { - rack(i) <<= index.bit; - rack(i) |= rack(i - 1) >> bi; - } - rack(i) <<= index.bit; - } - - return *self(); - } - - Derived& operator>>=(const size_t bits) - { - if (bits == 0) { - return *self(); - } else if (bits >= size()) { - reset(); - return *self(); - } - - const Index index{bits}; - - if (index.rack > 0) { - // Shift entire racks - size_t i; - for (i = 0; i < num_racks() - index.rack; ++i) { - rack(i) = rack(i + index.rack); - } - for (; i < num_racks(); ++i) { - rack(i) = 0; - } - } - - if (index.bit > 0) { - // Shift bits within racks - size_t bi = bits_per_rack - index.bit; - size_t i; - for (i = 0; i < num_racks() - index.rack - 1; ++i) { - rack(i) >>= index.bit; - rack(i) |= rack(i + 1) << bi; - } - rack(i) >>= index.bit; - } - - return *self(); - } - - auto begin(const size_t i = 0) { return iterator(*self(), i); } - auto end() { return iterator(*self(), size()); } - auto begin(const size_t i = 0) const { return const_iterator(*self(), i); } - auto end() const { return const_iterator(self(), size()); } - - const Rack& rack(const size_t index) const { return self()->rack(index); } - Rack& rack(const size_t index) { return self()->rack(index); } - Rack* data() { return self()->data(); } - const Rack* data() const { return self()->data(); } - size_t num_racks() const { return self()->num_racks(); } - size_t size() const { return self()->size(); } - size_t data_size() const { return self()->data_size(); } - -private: - using Index = BitVecIndex; - - Derived* self() { return static_cast(this); } - const Derived* self() const { return static_cast(this); } -}; - -template -void -gray_code(MultiBitVec& value) -{ - typename MultiBitVec::Rack s = 0; - - for (size_t ri = 0; ri < value.num_racks(); ++ri) { - const size_t i = value.num_racks() - ri - 1; - const auto t = value.rack(i) & 1; - gray_code(value.rack(i)); - value.rack(i) ^= (s << (MultiBitVec::bits_per_rack - 1)); - s = t; - } -} - -template -void -gray_code_inv(MultiBitVec& value) -{ - typename MultiBitVec::Rack s = 0; - - for (size_t ri = 0; ri < value.num_racks(); ++ri) { - const size_t i = value.num_racks() - ri - 1; - auto& rack = value.rack(i); - gray_code_inv(rack); - if (s) { - rack = ~rack; - } - s = rack & 1; - } -} - -} // namespace chilbert - -#endif diff --git a/chilbert/Operations.hpp b/chilbert/Operations.hpp deleted file mode 100644 index 24474cc..0000000 --- a/chilbert/Operations.hpp +++ /dev/null @@ -1,166 +0,0 @@ -/* - Copyright (C) 2018 David Robillard - Copyright (C) 2006-2007 Chris Hamilton - - This program is free software: you can redistribute it and/or modify it under - the terms of the GNU General Public License as published by the Free Software - Foundation, either version 2 of the License, or (at your option) any later - version. - - This program is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - FOR A PARTICULAR PURPOSE. See the GNU General Public License for more - details. - - You should have received a copy of the GNU General Public License along with - this program. If not, see . -*/ - -#ifndef CHILBERT_OPERATIONS_HPP -#define CHILBERT_OPERATIONS_HPP - -#include "chilbert/Operators.hpp" -#include "chilbert/Traits.hpp" - -#include -#include -#include -#include -#include - -namespace chilbert { - -/// Reset all bits in `field` -template -std::enable_if_t::value> -reset_bits(T& field) -{ - field = static_cast(0); -} - -/// Reset all bits in `field` -template -std::enable_if_t> -reset_bits(T& field) -{ - field.reset(); -} - -/// Return the `index`th bit in `field` -template -std::enable_if_t::value, bool> -test_bit(const T& field, const size_t index) -{ - assert(size_t(index) < sizeof(field) * CHAR_BIT); - return field & (T{1} << index); -} - -/// Return the `index`th bit in `field` -template -std::enable_if_t, bool> -test_bit(const T& field, const size_t index) -{ - return field.test(index); -} - -/// Set the `index`th bit in `field` -template -std::enable_if_t::value> -set_bit(T& field, const size_t index) -{ - assert(size_t(index) < sizeof(field) * CHAR_BIT); - field |= (T{1} << index); - assert(test_bit(field, index)); -} - -/// Set the `index`th bit in `field` to `value` -template -std::enable_if_t::value> -set_bit(T& field, const size_t index, const bool value) -{ - assert(size_t(index) < sizeof(field) * CHAR_BIT); - field ^= (-T{value} ^ field) & (T{1U} << index); - assert(test_bit(field, index) == value); -} - -/// Set the `index`th bit in `field` -template -std::enable_if_t> -set_bit(T& field, const size_t index) -{ - field.set(index); -} - -/// Set the `index`th bit in `field` to `value` -template -std::enable_if_t> -set_bit(T& field, const size_t index, const bool value) -{ - field.set(index, value); -} - -/// Return 1 + the index of the least significant 1-bit of `field`, or zero -template -int pop_count(const T& field); - -template <> -int -pop_count(const unsigned long& field) -{ - return __builtin_popcountl(field); -} - -template <> -int -pop_count(const unsigned long long& field) -{ - return __builtin_popcountll(field); -} - -/// Return 1 + the index of the least significant 1-bit of `field`, or zero -template -int find_first(const T field); - -template <> -int -find_first(const unsigned long field) -{ - return __builtin_ffsl(static_cast(field)); -} - -template <> -int -find_first(const unsigned long long field) -{ - return __builtin_ffsll(static_cast(field)); -} - -/// Calculates the Gray Code of `value` in place -template -std::enable_if_t> gray_code(T& value); - -/// Implementation of grayCode for any integral type -template -std::enable_if_t::value> -gray_code(T& value) -{ - value ^= (value >> 1); -} - -/// Calculates the inverse Gray Code of `value` in place -template -std::enable_if_t> gray_code_inv(T& value); - -/// Implementation of gray_code_inv for any integral type -template -std::enable_if_t::value> -gray_code_inv(T& value) -{ - for (T shift = 1; shift < sizeof(T) * CHAR_BIT; shift <<= 1) { - value ^= (value >> shift); - } -} - -} // namespace chilbert - -#endif diff --git a/chilbert/Operators.hpp b/chilbert/Operators.hpp deleted file mode 100644 index 683c792..0000000 --- a/chilbert/Operators.hpp +++ /dev/null @@ -1,94 +0,0 @@ -/* - Copyright (C) 2018 David Robillard - Copyright (C) 2006-2007 Chris Hamilton - - This program is free software: you can redistribute it and/or modify it under - the terms of the GNU General Public License as published by the Free Software - Foundation, either version 2 of the License, or (at your option) any later - version. - - This program is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - FOR A PARTICULAR PURPOSE. See the GNU General Public License for more - details. - - You should have received a copy of the GNU General Public License along with - this program. If not, see . -*/ - -#ifndef CHILBERT_OPERATORS_HPP -#define CHILBERT_OPERATORS_HPP - -#include "chilbert/Traits.hpp" - -#include -#include - -namespace chilbert { - -template -std::enable_if_t, T> operator&(const T& lhs, const T& rhs) -{ - T r{lhs}; - r &= rhs; - return r; -} - -template -std::enable_if_t, T> -operator|(const T& lhs, const T& rhs) -{ - T r{lhs}; - r |= rhs; - return r; -} - -template -std::enable_if_t, T> -operator^(const T& lhs, const T& rhs) -{ - T r{lhs}; - r ^= rhs; - return r; -} - -template -std::enable_if_t, T> -operator~(const T& vec) -{ - T r{vec}; - r.flip(); - return r; -} - -template -std::enable_if_t, T> -operator<<(const T& vec, const size_t bits) -{ - T r{vec}; - r <<= bits; - return r; -} - -template -std::enable_if_t, T> -operator>>(const T& vec, const size_t bits) -{ - T r{vec}; - r >>= bits; - return r; -} - -template -inline std::enable_if_t, std::ostream>& -operator<<(std::ostream& os, const T& vec) -{ - for (size_t i = 0; i < vec.size(); ++i) { - os << vec.test(vec.size() - i - 1); - } - return os; -} - -} // namespace chilbert - -#endif diff --git a/chilbert/SmallBitVec.hpp b/chilbert/SmallBitVec.hpp index 308d943..478cbc9 100644 --- a/chilbert/SmallBitVec.hpp +++ b/chilbert/SmallBitVec.hpp @@ -19,7 +19,7 @@ #ifndef CHILBERT_SMALLBITVEC_HPP #define CHILBERT_SMALLBITVEC_HPP -#include "chilbert/Operations.hpp" +#include "chilbert/detail/operations.hpp" #include #include @@ -233,11 +233,14 @@ public: /// Return 1 + the index of the first set bit, or 0 if there are none size_t find_first() const { - return static_cast(chilbert::find_first(m_rack)); + return static_cast(detail::find_first(m_rack)); } /// Return the number of set bits - size_t count() const { return static_cast(pop_count(m_rack)); } + size_t count() const + { + return static_cast(detail::pop_count(m_rack)); + } /// Flip all bits (one's complement) SmallBitVec& flip() @@ -347,6 +350,8 @@ private: size_t m_size{}; }; +namespace detail { + template <> struct is_bitvec { @@ -367,6 +372,8 @@ gray_code_inv(SmallBitVec& value) gray_code_inv(value.rack()); } +} // namespace detail + } // namespace chilbert #endif diff --git a/chilbert/StaticBitVec.hpp b/chilbert/StaticBitVec.hpp index 66d136a..9aff3ad 100644 --- a/chilbert/StaticBitVec.hpp +++ b/chilbert/StaticBitVec.hpp @@ -19,11 +19,11 @@ #ifndef CHILBERT_STATICBITVEC_HPP #define CHILBERT_STATICBITVEC_HPP -#include "chilbert/BitVecIndex.hpp" -#include "chilbert/BitVecIterator.hpp" -#include "chilbert/BitVecMask.hpp" -#include "chilbert/MultiBitVec.hpp" -#include "chilbert/Operations.hpp" +#include "chilbert/detail/BitVecIndex.hpp" +#include "chilbert/detail/BitVecIterator.hpp" +#include "chilbert/detail/BitVecMask.hpp" +#include "chilbert/detail/MultiBitVec.hpp" +#include "chilbert/detail/operations.hpp" #include #include @@ -43,12 +43,12 @@ namespace chilbert { * @tparam N Number of bits. */ template -class StaticBitVec : public MultiBitVec> +class StaticBitVec : public detail::MultiBitVec> { public: - using Rack = typename MultiBitVec>::Rack; + using Rack = typename detail::MultiBitVec>::Rack; - using MultiBitVec>::bits_per_rack; + using detail::MultiBitVec>::bits_per_rack; StaticBitVec() = default; @@ -91,6 +91,8 @@ private: std::array m_racks{}; }; +namespace detail { + template struct is_bitvec> { @@ -111,6 +113,8 @@ gray_code_inv(StaticBitVec& value) gray_code_inv(static_cast>&>(value)); } +} // namespace detail + } // namespace chilbert #endif diff --git a/chilbert/Traits.hpp b/chilbert/Traits.hpp deleted file mode 100644 index 0c7807d..0000000 --- a/chilbert/Traits.hpp +++ /dev/null @@ -1,37 +0,0 @@ -/* - Copyright (C) 2018 David Robillard - Copyright (C) 2006-2007 Chris Hamilton - - This program is free software: you can redistribute it and/or modify it under - the terms of the GNU General Public License as published by the Free Software - Foundation, either version 2 of the License, or (at your option) any later - version. - - This program is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS - FOR A PARTICULAR PURPOSE. See the GNU General Public License for more - details. - - You should have received a copy of the GNU General Public License along with - this program. If not, see . -*/ - -#ifndef CHILBERT_TRAITS_HPP -#define CHILBERT_TRAITS_HPP - -namespace chilbert { - -/// Member `value` is true iff T is a chilbert bitset -template -struct is_bitvec -{ - static constexpr bool value = false; -}; - -/// True iff T is a chilbert bitset -template -static constexpr bool is_bitvec_v = is_bitvec::value; - -} // namespace chilbert - -#endif diff --git a/chilbert/chilbert.hpp b/chilbert/chilbert.hpp new file mode 100644 index 0000000..2954cf3 --- /dev/null +++ b/chilbert/chilbert.hpp @@ -0,0 +1,102 @@ +/* + Copyright (C) 2018 David Robillard + Copyright (C) 2006-2007 Chris Hamilton + + This program is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License as published by the Free Software + Foundation, either version 2 of the License, or (at your option) any later + version. + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + details. + + You should have received a copy of the GNU General Public License along with + this program. If not, see . +*/ + +#ifndef CHILBERT_HILBERT_HPP +#define CHILBERT_HILBERT_HPP + +#include + +namespace chilbert { + +/** Map the point `p` to a Hilbert Index. + * + * @tparam P Type used to represent a value in a point dimension. + * @tparam H Hilbert Index. + * + * @param p Point with `n` dimensions. + * @param m Precision of each dimension in bits. + * @param n Number of dimensions. + * @param[out] h Hilbert Index. + */ +template +inline void coords_to_index(const P* const p, + const size_t m, + const size_t n, + H& h); + +/** Map the Hilbert Index `p` to a point. + * + * @tparam P Type used to represent a value in a point dimension. + * @tparam H Hilbert Index. + * + * @param[out] p Point with `n` dimensions. + * @param m Precision of each dimension in bits. + * @param n Number of dimensions. + * @param h Hilbert Index. + */ +template +inline void index_to_coords(P* const p, + const size_t m, + const size_t n, + const H& h); + +/** Map the point `p` to a Compact Hilbert Index. + * + * @tparam P Type used to represent a value in a point dimension. + * @tparam H Compact Hilbert Index. + * + * @param p Point with `n` dimensions. + * @param ms Array of `n` precision values for each dimension in bits. + * @param n Number of dimensions. + * @param[out] hc Compact Hilbert Index. + * @param M Optional net precision (sum of `ms`), the size of `hc` in bits. + * @param m Optional largest precision in `m`. + */ +template +inline void coords_to_compact_index(const P* const p, + const size_t* const ms, + const size_t n, + H& hc, + const size_t M = 0, + const size_t m = 0); + +/** Map the Compact Hilbert Index `hc` to a point. + * + * @tparam P Type used to represent a value in a point dimension. + * @tparam H Compact Hilbert Index. + * + * @param[out] p Point with `n` dimensions. + * @param ms Array of `n` precision values for each dimension in bits. + * @param n Number of dimensions. + * @param hc Compact Hilbert Index. + * @param M Optional net precision (sum of `ms`), the size of `hc` in bits. + * @param m Optional largest precision in `m`. + */ +template +inline void compact_index_to_coords(P* const p, + const size_t* const ms, + const size_t n, + const H& hc, + const size_t M = 0, + const size_t m = 0); + +} // namespace chilbert + +#include "chilbert/chilbert.ipp" + +#endif diff --git a/chilbert/chilbert.ipp b/chilbert/chilbert.ipp new file mode 100644 index 0000000..c4e8249 --- /dev/null +++ b/chilbert/chilbert.ipp @@ -0,0 +1,501 @@ +/* + Copyright (C) 2018 David Robillard + Copyright (C) 2006-2007 Chris Hamilton + + This program is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License as published by the Free Software + Foundation, either version 2 of the License, or (at your option) any later + version. + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + details. + + You should have received a copy of the GNU General Public License along with + this program. If not, see . +*/ + +#ifndef CHILBERT_ALGORITHM_HPP +#define CHILBERT_ALGORITHM_HPP + +#include "chilbert/BoundedBitVec.hpp" +#include "chilbert/DynamicBitVec.hpp" +#include "chilbert/SmallBitVec.hpp" +#include "chilbert/StaticBitVec.hpp" +#include "chilbert/chilbert.hpp" +#include "chilbert/detail/gray_code_rank.hpp" + +#include +#include +#include +#include + +namespace chilbert { +namespace detail { + +/** The dimension across which the Hilbert curve travels principally. + * + * D0 is d + 1, where d is the actual dimension you want to walk across. + * MUST HAVE 0 <= D0 < n. + */ +static constexpr size_t D0 = 1; + +template +size_t +num_bits(const T&, std::enable_if_t::value>* = nullptr) +{ + return sizeof(T) * CHAR_BIT; +} + +template +size_t +num_bits(const T& vec, + std::enable_if_t::value || + std::is_same::value>* = nullptr) +{ + return vec.size(); +} + +template +size_t +num_bits(const StaticBitVec&, void* = nullptr) +{ + return N; +} + +template +size_t +num_bits(const BoundedBitVec& vec, void* = nullptr) +{ + return vec.size(); +} + +/** Copy a range of bits from one field to the start of another. + * + * @param h Source bit field + * @param n Number of bits + * @param i Start bit index in source + * @param[out] w Destination bit field + */ +template +inline void +get_bits(const H& h, const size_t n, const size_t i, I& w) +{ + for (size_t j = 0; j < n; j++) { + set_bit(w, j, test_bit(h, i + j)); + } +} + +/** Set a range of bits in one field from the start of another. + * + * @param[out] h Destination bit field + * @param n Number of bits + * @param i Start bit index in destination + * @param w Source bit field + */ +template +inline void +set_bits(H& h, const size_t n, const size_t i, const I& w) +{ + for (size_t j = 0; j < n; j++) { + set_bit(h, i + j, test_bit(w, j)); + } +} + +/** Set the leading bits in `l` to bit `i` from the corresponding value in `p`. + * + * @param p Point. + * @param n Number of bits to set. + * @param i Bit position to copy from values in `p`. + * @param[out] l Output bits. + */ +template +inline void +get_location(const P* const p, const size_t n, const size_t i, I& l) +{ + for (size_t j = 0; j < n; ++j) { + set_bit(l, j, test_bit(p[j], i)); + } +} + +/** Set bit `i` in values in `p` to the corresponding leadings bits in `l`. + * + * @param[out] p Point. + * @param n Number of bits to set. + * @param i Bit position to set in values in `p`. + * @param l Input bits. + */ +template +inline void +set_location(P* const p, const size_t n, const size_t i, const I& l) +{ + for (size_t j = 0; j < n; j++) { + set_bit(p[j], i, test_bit(l, j)); + } +} + +// 'Transforms' a point. +template +inline void +transform(const I& e, const size_t d, const size_t n, I& a) +{ + assert(a.size() == n); + a ^= e; + a.rotr(d); +} + +// Inverse 'transforms' a point. +template +inline void +transform_inv(const I& e, const size_t d, const size_t n, I& a) +{ + assert(a.size() == n); + a.rotl(d); + a ^= e; +} + +// Update for method 1 (GrayCodeInv in the loop) +template +inline void +update1(const I& l, const I& t, const I& w, const size_t n, I& e, size_t& d) +{ + assert(0 <= d && d < n); + e = l; + e.flip(d); //#D d == n-1 ? 0 : d+1 ); + + // Update direction + d += 1 + t.find_first(); + if (d >= n) { + d -= n; + } + if (d >= n) { + d -= n; + } + assert(0 <= d && d < n); + + if (!w.test(0)) { + e.flip(d == 0 ? n - 1 : d - 1); //#D d ); + } +} + +// Update for method 2 (GrayCodeInv out of loop) +template +inline void +update2(const I& l, const I& t, const size_t n, I& e, size_t& d) +{ + assert(0 <= d && d < n); + e = l; + e.flip(d); //#D d == n-1 ? 0 : d+1 ); + + // Update direction + d += 1 + t.find_first(); + if (d >= n) { + d -= n; + } + if (d >= n) { + d -= n; + } + assert(0 <= d && d < n); +} + +template +inline void +coords_to_index(const P* const p, + const size_t m, + const size_t n, + H& h, + I&& scratch, + size_t* const ds = nullptr) +{ + I e{std::move(scratch)}; + I l{e}; + I t{e}; + I w{e}; + + // Initialize + e.reset(); + l.reset(); + reset_bits(h); + + // Work from MSB to LSB + size_t d = D0; + size_t ho = m * n; + for (intptr_t i = static_cast(m - 1); i >= 0; i--) { + if (ds) { + ds[i] = d; + } + + // Get corner of sub-hypercube where point lies. + get_location(p, n, static_cast(i), l); + + // Mirror and reflect the location. + // t = T_{(e,d)}(l) + t = l; + transform(e, d, n, t); + + w = t; + if (static_cast(i) < m - 1) { + w.flip(n - 1); + } + + // Concatenate to the index. + ho -= n; + set_bits(h, n, ho, w); + + // Update the entry point and direction. + update2(l, t, n, e, d); + } + + gray_code_inv(h); +} + +template +inline void +index_to_coords(P* p, const size_t m, const size_t n, const H& h, I&& scratch) +{ + I e{std::move(scratch)}; + I l{e}; + I t{e}; + I w{e}; + + // Initialize + e.reset(); + l.reset(); + for (size_t j = 0; j < n; j++) { + reset_bits(p[j]); + } + + // Work from MSB to LSB + size_t d = D0; + size_t ho = m * n; + for (intptr_t i = static_cast(m - 1); i >= 0; i--) { + // Get the Hilbert index bits + ho -= n; + get_bits(h, n, ho, w); + + // t = GrayCode(w) + t = w; + gray_code(t); + + // Reverse the transform + // l = T^{-1}_{(e,d)}(t) + l = t; + transform_inv(e, d, n, l); + + // Distribute these bits + // to the coordinates. + set_location(p, n, static_cast(i), l); + + // Update the entry point and direction. + update1(l, t, w, n, e, d); + } +} + +template +inline void +coords_to_compact_index(const P* const p, + const size_t* const ms, + const size_t n, + HC& hc, + I&& scratch, + size_t M = 0, + size_t m = 0) +{ + // Get total precision and max precision if not supplied + if (M == 0 || m == 0) { + M = m = 0; + for (size_t i = 0; i < n; i++) { + assert(num_bits(p[i]) >= ms[i]); + if (ms[i] > m) { + m = ms[i]; + } + M += ms[i]; + } + } + + const size_t mn = m * n; + + assert(num_bits(hc) >= M); + + // If we could avoid allocation altogether (ie: have a + // fixed buffer allocated on the stack) then this increases + // speed by a bit (4% when n=4, m=20) + size_t* const ds = new size_t[m]; + + if (mn > SmallBitVec::bits_per_rack) { + DynamicBitVec h(mn); + detail::coords_to_index( + p, m, n, h, std::move(scratch), ds); + compact_index(ms, ds, n, m, h, hc); + } else { + SmallBitVec h(mn); + detail::coords_to_index( + p, m, n, h, std::move(scratch), ds); + compact_index(ms, ds, n, m, h, hc); + } + + delete[] ds; +} + +template +inline void +compact_index_to_coords(P* const p, + const size_t* ms, + const size_t n, + const HC& hc, + I&& scratch, + size_t M, + size_t m) +{ + I e{std::move(scratch)}; + I l{e}; + I t{e}; + I w{e}; + I r{e}; + I mask{e}; + I ptrn{e}; + + // Get total precision and max precision + // if not supplied + if (M == 0 || m == 0) { + M = m = 0; + for (size_t i = 0; i < n; i++) { + if (ms[i] > m) { + m = ms[i]; + } + M += ms[i]; + } + } + + assert(num_bits(hc) >= M); + assert(num_bits(p[0]) >= m); + + // Initialize + e.reset(); + l.reset(); + for (size_t j = 0; j < n; j++) { + reset_bits(p[j]); + } + + // Work from MSB to LSB + size_t d = D0; + + for (intptr_t i = static_cast(m - 1); i >= 0; i--) { + // Get the mask and ptrn + size_t b = 0; + extract_mask(ms, n, d, static_cast(i), mask, b); + ptrn = e; + assert(ptrn.size() == n); + ptrn.rotr(d); + + // Get the Hilbert index bits + M -= b; + r.reset(); // GetBits doesn't do this + get_bits(hc, b, M, r); + + // w = GrayCodeRankInv(r) + // t = GrayCode(w) + gray_code_rank_inv(mask, ptrn, r, n, b, t, w); + + // Reverse the transform + // l = T^{-1}_{(e,d)}(t) + l = t; + transform_inv(e, d, n, l); + + // Distribute these bits + // to the coordinates. + set_location(p, n, static_cast(i), l); + + // Update the entry point and direction. + update1(l, t, w, n, e, d); + } +} + +} // namespace detail + +template +inline void +coords_to_index(const P* const p, const size_t m, const size_t n, H& h) +{ + assert(detail::num_bits(h) >= n * m); + assert(detail::num_bits(p[0]) >= m); + + if (n <= SmallBitVec::bits_per_rack) { + // Intermediate variables will fit in fixed width + detail::coords_to_index(p, m, n, h, SmallBitVec(n)); + } else { + // Otherwise, they must be DynamicBitVecs + detail::coords_to_index( + p, m, n, h, DynamicBitVec(n)); + } +} + +template +inline void +index_to_coords(P* const p, const size_t m, const size_t n, const H& h) +{ + assert(m > 0); + assert(n > 0); + assert(detail::num_bits(h) >= n * m); + assert(detail::num_bits(p[0]) >= m); + + if (n <= SmallBitVec::bits_per_rack) { + // Intermediate variables will fit in fixed width + detail::index_to_coords(p, m, n, h, SmallBitVec(n)); + } else { + // Otherwise, they must be DynamicBitVecs + detail::index_to_coords( + p, m, n, h, DynamicBitVec(n)); + } +} + +template +inline void +coords_to_compact_index(const P* const p, + const size_t* const ms, + size_t n, + HC& hc, + const size_t M, + const size_t m) +{ + assert(hc.size() >= std::accumulate(ms, ms + n, size_t(0))); + + if (n <= SmallBitVec::bits_per_rack) { + // Intermediate variables will fit in fixed width + detail::coords_to_compact_index( + p, ms, n, hc, SmallBitVec(n), M, m); + } else { + // Otherwise, they must be DynamicBitVecs + detail::coords_to_compact_index( + p, ms, n, hc, DynamicBitVec(n), M, m); + } +} + +template +inline void +compact_index_to_coords(P* const p, + const size_t* ms, + const size_t n, + const HC& hc, + const size_t M, + const size_t m) +{ + assert(hc.size() >= std::accumulate(ms, ms + n, size_t(0))); + + if (n <= SmallBitVec::bits_per_rack) { + // Intermediate variables will fit in fixed width + SmallBitVec scratch(n); + detail::compact_index_to_coords( + p, ms, n, hc, std::move(scratch), M, m); + } else { + // Otherwise, they must be DynamicBitVecs + DynamicBitVec scratch(n); + detail::compact_index_to_coords( + p, ms, n, hc, std::move(scratch), M, m); + } +} + +} // namespace chilbert + +#endif diff --git a/chilbert/detail/BitVecIndex.hpp b/chilbert/detail/BitVecIndex.hpp new file mode 100644 index 0000000..e7b385e --- /dev/null +++ b/chilbert/detail/BitVecIndex.hpp @@ -0,0 +1,51 @@ +/* + Copyright (C) 2018 David Robillard + Copyright (C) 2006-2007 Chris Hamilton + + This program is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License as published by the Free Software + Foundation, either version 2 of the License, or (at your option) any later + version. + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + details. + + You should have received a copy of the GNU General Public License along with + this program. If not, see . +*/ + +#ifndef CHILBERT_BITVECINDEX_HPP +#define CHILBERT_BITVECINDEX_HPP + +#include +#include +#include + +namespace chilbert { +namespace detail { + +/// Index into a multi-rack bit vector +template +struct BitVecIndex +{ + using Rack = typename BitVec::Rack; + + static constexpr size_t bits_per_rack = sizeof(Rack) * CHAR_BIT; + + BitVecIndex(const size_t bits) + : rack{bits / bits_per_rack} + , bit{bits - rack * bits_per_rack} + { + assert(bit < bits_per_rack); + } + + size_t rack; + size_t bit; +}; + +} // namespace detail +} // namespace chilbert + +#endif diff --git a/chilbert/detail/BitVecIterator.hpp b/chilbert/detail/BitVecIterator.hpp new file mode 100644 index 0000000..8902747 --- /dev/null +++ b/chilbert/detail/BitVecIterator.hpp @@ -0,0 +1,100 @@ +/* + Copyright (C) 2018 David Robillard + Copyright (C) 2006-2007 Chris Hamilton + + This program is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License as published by the Free Software + Foundation, either version 2 of the License, or (at your option) any later + version. + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + details. + + You should have received a copy of the GNU General Public License along with + this program. If not, see . +*/ + +#ifndef CHILBERT_BITVECITERATOR_HPP +#define CHILBERT_BITVECITERATOR_HPP + +#include "chilbert/detail/BitVecMask.hpp" + +#include + +namespace chilbert { +namespace detail { + +template +class BitVecIteratorBase : public BitVecMask +{ +public: + using Mask = typename BitVec::Mask; + + BitVecIteratorBase& operator++() + { + Mask::operator++(); + return *this; + } + + BitVecIteratorBase& operator--() + { + Mask::operator--(); + return *this; + } + + bool operator==(const BitVecIteratorBase& rhs) const + { + return m_vec == rhs.m_vec && Mask::operator==(rhs); + } + + bool operator!=(const BitVecIteratorBase& rhs) const + { + return !operator==(rhs); + } + + bool operator*() const { return m_vec->test(*this); } + +protected: + BitVecIteratorBase(BitVec& vec, const size_t index) + : Mask{index} + , m_vec{&vec} + { + } + + BitVec* m_vec; +}; + +template +class BitVecIterator : public BitVecIteratorBase +{ +public: + void set() { this->m_vec->set(*this); } + void reset() { this->m_vec->reset(*this); } + +private: + friend BitVec; + + BitVecIterator(BitVec& vec, const size_t index) + : BitVecIteratorBase{vec, index} + { + } +}; + +template +class ConstBitVecIterator : public BitVecIteratorBase +{ +private: + friend BitVec; + + ConstBitVecIterator(const BitVec& vec, const size_t index) + : BitVecIteratorBase{vec, index} + { + } +}; + +} // namespace detail +} // namespace chilbert + +#endif diff --git a/chilbert/detail/BitVecMask.hpp b/chilbert/detail/BitVecMask.hpp new file mode 100644 index 0000000..03eaf5f --- /dev/null +++ b/chilbert/detail/BitVecMask.hpp @@ -0,0 +1,74 @@ +/* + Copyright (C) 2018 David Robillard + Copyright (C) 2006-2007 Chris Hamilton + + This program is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License as published by the Free Software + Foundation, either version 2 of the License, or (at your option) any later + version. + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + details. + + You should have received a copy of the GNU General Public License along with + this program. If not, see . +*/ + +#ifndef CHILBERT_BITVECMASK_HPP +#define CHILBERT_BITVECMASK_HPP + +#include +#include + +namespace chilbert { +namespace detail { + +/** Mask for a bit that can be incremented like an index. + * + * This enables fast iteration while setting or resetting bits since it is + * internally stored as a mask which avoids repeated index math. + */ +template +struct BitVecMask +{ + static constexpr size_t bits_per_rack = sizeof(Rack) * CHAR_BIT; + + BitVecMask(const size_t index) + : rack{index / bits_per_rack} + , mask{Rack{1} << (index - rack * bits_per_rack)} + { + } + + void operator++() + { + if ((mask <<= 1) == 0) { + mask = 1; + ++rack; + } + } + + void operator--() + { + if ((mask >>= 1) == 0) { + mask = Rack{1} << (bits_per_rack - 1); + --rack; + } + } + + bool operator==(const BitVecMask& rhs) const + { + return rack == rhs.rack && mask == rhs.mask; + } + + bool operator!=(const BitVecMask& rhs) const { return !operator==(rhs); } + + size_t rack; + Rack mask; +}; + +} // namespace detail +} // namespace chilbert + +#endif diff --git a/chilbert/detail/MultiBitVec.hpp b/chilbert/detail/MultiBitVec.hpp new file mode 100644 index 0000000..4bf5085 --- /dev/null +++ b/chilbert/detail/MultiBitVec.hpp @@ -0,0 +1,383 @@ +/* + Copyright (C) 2018 David Robillard + Copyright (C) 2006-2007 Chris Hamilton + + This program is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License as published by the Free Software + Foundation, either version 2 of the License, or (at your option) any later + version. + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + details. + + You should have received a copy of the GNU General Public License along with + this program. If not, see . +*/ + +#ifndef CHILBERT_MULTIBITVEC_HPP +#define CHILBERT_MULTIBITVEC_HPP + +#include "chilbert/detail/BitVecIndex.hpp" +#include "chilbert/detail/BitVecIterator.hpp" +#include "chilbert/detail/BitVecMask.hpp" +#include "chilbert/detail/operations.hpp" + +#include +#include +#include +#include + +namespace chilbert { +namespace detail { + +template +class MultiBitVec +{ +public: + using Rack = uintptr_t; + using Mask = BitVecMask; + + using iterator = BitVecIterator>; + using const_iterator = ConstBitVecIterator>; + + static constexpr size_t bits_per_rack = sizeof(Rack) * CHAR_BIT; + + /// Return the value of the bit covered by `mask` + bool test(const Mask mask) const { return rack(mask.rack) & mask.mask; } + + /// Return the value of the `index`th bit + bool test(const size_t index) const { return test(mask(index)); } + + /// Set all bits to one + Derived& set() + { + if (size()) { + memset(data(), 0xFF, data_size()); + self()->truncate(); + } + + return *self(); + } + + /// Set the bit covered by `mask` to 1 + Derived& set(const Mask mask) + { + rack(mask.rack) |= mask.mask; + return *self(); + } + + /// Set the `index`th bit to 1 + Derived& set(const size_t index) { return set(mask(index)); } + + /// Set the bit covered by `mask` to `value` + Derived& set(const Mask mask, const bool value) + { + auto& r = rack(mask.rack); + r ^= (-Rack{value} ^ r) & mask.mask; + return *self(); + } + + /// Set the `index`th bit to `value` + Derived& set(const size_t index, const bool value) + { + return set(mask(index), value); + } + + /// Set all bits to zero + Derived& reset() + { + memset(data(), 0, data_size()); + return *self(); + } + + /// Reset the bit covered by `mask` to 0 + Derived& reset(const Mask mask) + { + rack(mask.rack) &= ~mask.mask; + return *self(); + } + + /// Reset the `index`th bit to 0 + Derived& reset(const size_t index) { return reset(mask(index)); } + + /// Flip all bits (one's complement) + Derived& flip() + { + for (size_t i = 0; i < num_racks(); ++i) { + rack(i) = ~rack(i); + } + return *self(); + } + + /// Flip the value of the bit covered by `mask` + Derived& flip(const Mask mask) + { + rack(mask.rack) ^= mask.mask; + return *self(); + } + + /// Flip the value of the `index`th bit + Derived& flip(const size_t index) { return flip(mask(index)); } + + /// Clear any bits in storage outside the valid range if necessary + void truncate() + { + if (const auto pad = num_racks() * bits_per_rack - size()) { + rack(num_racks() - 1) &= ~Rack{0} >> pad; + } + } + + /// Right-rotate by `bits` positions + Derived& rotr(const size_t bits) + { + assert(bits <= size()); + if (bits == 0 || bits == size()) { + return *self(); + } + + Derived t1(*self()); + *self() >>= bits; + t1 <<= (size() - bits); + *self() |= t1; + + truncate(); + return *self(); + } + + /// Left-rotate by `bits` positions + Derived& rotl(const size_t bits) + { + assert(bits <= size()); + if (bits == 0 || bits == size()) { + return *self(); + } + + Derived t1(*self()); + *self() <<= bits; + t1 >>= (size() - bits); + *self() |= t1; + + truncate(); + return *self(); + } + + /// Return true iff all bits are zero + bool none() const + { + for (size_t i = 0; i < num_racks(); ++i) { + if (rack(i)) { + return false; + } + } + return true; + } + + /// Return 1 + the index of the first set bit, or 0 if there are none + size_t find_first() const + { + for (size_t i = 0; i < num_racks(); ++i) { + const int j = chilbert::detail::find_first(rack(i)); + if (j) { + return (i * bits_per_rack) + static_cast(j); + } + } + return 0; + } + + /// Return the number of set bits + size_t count() const + { + size_t c = 0; + for (size_t i = 0; i < num_racks(); ++i) { + c += static_cast(pop_count(rack(i))); + } + return c; + } + + /// Return a mask that covers the bit with index `i` + Mask mask(const size_t i = 0) const + { + assert(i <= size()); + return Mask{i}; + } + + bool operator==(const Derived& vec) const + { + return (num_racks() == vec.num_racks() && + (num_racks() == 0 || + !memcmp(data(), vec.data(), num_racks() * sizeof(Rack)))); + } + + bool operator!=(const Derived& vec) const { return !(*this == vec); } + + bool operator<(const Derived& vec) const + { + assert(size() == vec.size()); + + for (size_t ri = 0; ri < num_racks(); ++ri) { + size_t i = num_racks() - ri - 1; + if (rack(i) < vec.rack(i)) { + return true; + } else if (rack(i) > vec.rack(i)) { + return false; + } + } + return false; + } + + Derived& operator&=(const Derived& vec) + { + for (size_t i = 0; i < std::min(num_racks(), vec.num_racks()); ++i) { + rack(i) &= vec.rack(i); + } + + return *self(); + } + + Derived& operator|=(const Derived& vec) + { + for (size_t i = 0; i < std::min(num_racks(), vec.num_racks()); ++i) { + rack(i) |= vec.rack(i); + } + + return *self(); + } + + Derived& operator^=(const Derived& vec) + { + for (size_t i = 0; i < std::min(num_racks(), vec.num_racks()); ++i) { + rack(i) ^= vec.rack(i); + } + + return *self(); + } + + Derived& operator<<=(const size_t bits) + { + if (bits == 0) { + return *self(); + } else if (bits >= size()) { + reset(); + return *self(); + } + + const Index index{bits}; + + if (index.rack > 0) { + // Shift entire racks + for (size_t i = num_racks() - 1; i >= index.rack; --i) { + rack(i) = rack(i - index.rack); + } + for (size_t i = 0; i < index.rack; ++i) { + rack(i) = 0; + } + } + + if (index.bit > 0) { + // Shift bits within racks + size_t bi = bits_per_rack - index.bit; + size_t i; + for (i = num_racks() - 1; i >= index.rack + 1; --i) { + rack(i) <<= index.bit; + rack(i) |= rack(i - 1) >> bi; + } + rack(i) <<= index.bit; + } + + return *self(); + } + + Derived& operator>>=(const size_t bits) + { + if (bits == 0) { + return *self(); + } else if (bits >= size()) { + reset(); + return *self(); + } + + const Index index{bits}; + + if (index.rack > 0) { + // Shift entire racks + size_t i; + for (i = 0; i < num_racks() - index.rack; ++i) { + rack(i) = rack(i + index.rack); + } + for (; i < num_racks(); ++i) { + rack(i) = 0; + } + } + + if (index.bit > 0) { + // Shift bits within racks + size_t bi = bits_per_rack - index.bit; + size_t i; + for (i = 0; i < num_racks() - index.rack - 1; ++i) { + rack(i) >>= index.bit; + rack(i) |= rack(i + 1) << bi; + } + rack(i) >>= index.bit; + } + + return *self(); + } + + auto begin(const size_t i = 0) { return iterator(*self(), i); } + auto end() { return iterator(*self(), size()); } + auto begin(const size_t i = 0) const { return const_iterator(*self(), i); } + auto end() const { return const_iterator(self(), size()); } + + const Rack& rack(const size_t index) const { return self()->rack(index); } + Rack& rack(const size_t index) { return self()->rack(index); } + Rack* data() { return self()->data(); } + const Rack* data() const { return self()->data(); } + size_t num_racks() const { return self()->num_racks(); } + size_t size() const { return self()->size(); } + size_t data_size() const { return self()->data_size(); } + +private: + using Index = detail::BitVecIndex; + + Derived* self() { return static_cast(this); } + const Derived* self() const { return static_cast(this); } +}; + +template +void +gray_code(MultiBitVec& value) +{ + typename MultiBitVec::Rack s = 0; + + for (size_t ri = 0; ri < value.num_racks(); ++ri) { + const size_t i = value.num_racks() - ri - 1; + const auto t = value.rack(i) & 1; + gray_code(value.rack(i)); + value.rack(i) ^= (s << (MultiBitVec::bits_per_rack - 1)); + s = t; + } +} + +template +void +gray_code_inv(MultiBitVec& value) +{ + typename MultiBitVec::Rack s = 0; + + for (size_t ri = 0; ri < value.num_racks(); ++ri) { + const size_t i = value.num_racks() - ri - 1; + auto& rack = value.rack(i); + gray_code_inv(rack); + if (s) { + rack = ~rack; + } + s = rack & 1; + } +} + +} // namespace detail +} // namespace chilbert + +#endif diff --git a/chilbert/detail/gray_code_rank.hpp b/chilbert/detail/gray_code_rank.hpp new file mode 100644 index 0000000..815519e --- /dev/null +++ b/chilbert/detail/gray_code_rank.hpp @@ -0,0 +1,182 @@ +/* + Copyright (C) 2018 David Robillard + Copyright (C) 2006-2007 Chris Hamilton + + This program is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License as published by the Free Software + Foundation, either version 2 of the License, or (at your option) any later + version. + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + details. + + You should have received a copy of the GNU General Public License along with + this program. If not, see . +*/ + +#ifndef CHILBERT_GRAYCODERANK_HPP +#define CHILBERT_GRAYCODERANK_HPP + +#include +#include +#include + +namespace chilbert { +namespace detail { + +// This is the bulk of the cost in calculating +// a Compact Hilbert Index. It compresses a previously +// calculated index when provided the rotation +// at each level of precision. +template +inline void +compact_index(const size_t* const ms, + const size_t* const ds, + const size_t n, + const size_t m, + H& h, + HC& hc) +{ + assert(h.size() >= n * m); + assert(hc.size() >= std::accumulate(ms, ms + n, size_t(0))); + + reset_bits(hc); + + auto hm = h.mask(0); + auto hcm = hc.mask(0); + + // Run through the levels of precision + for (size_t i = 0; i < m; i++) { + // Run through the dimensions + size_t j = ds[i]; + do { + // This dimension contributes a bit? + if (ms[j] > i) { + if (h.test(hm)) { + hc.set(hcm); + } + ++hcm; + } + + if (++j == n) { + j = 0; + } + ++hm; + } while (j != ds[i]); + } +} + +template +inline void +gray_code_rank(const I& mask, const I& gi, const size_t n, I& r) +{ + assert(mask.size() == n); + assert(gi.size() == n); + assert(r.size() == n); + + r.reset(); + + auto mi = mask.begin(); + auto gii = gi.begin(); + auto ri = r.begin(); + + for (size_t i = 0; i < n; ++i, ++mi, ++gii) { + if (*mi) { + if (*gii) { + ri.set(); + } + ++ri; + } + } +} + +template +inline void +gray_code_rank_inv(const I& mask, + const I& ptrn, + const I& r, + const size_t n, + const size_t b, + I& g, + I& gi) +{ + assert(mask.size() == n); + assert(ptrn.size() == n); + assert(r.size() == n); + assert(g.size() == n); + assert(gi.size() == n); + + g.reset(); + gi.reset(); + + auto m = mask.mask(n - 1); + auto ri = r.begin(b - 1); + + typename I::Rack gi0 = 0; + typename I::Rack gi1 = 0; + typename I::Rack g0 = 0; + + for (size_t i = 0; i < n; ++i) { + if (mask.test(m)) { // Unconstrained bit + gi1 = gi0; + gi0 = *ri; + g0 = gi0 ^ gi1; + if (gi0) { + gi.set(m); + } + if (g0) { + g.set(m); + } + --ri; + } else { // Constrained bit + g0 = (ptrn.test(m) > 0); + gi1 = gi0; + gi0 = g0 ^ gi1; + if (gi0) { + gi.set(m); + } + if (g0) { + g.set(m); + } + } + + --m; + } +} + +template +inline void +extract_mask(const size_t* const ms, + const size_t n, + const size_t d, + const size_t i, + I& mask, + size_t& b) +{ + assert(0 <= d && d < n); + assert(mask.size() == n); + + mask.reset(); + b = 0; + + auto mi = mask.begin(); + size_t j = d; // #D j = (d==n-1) ? 0 : d+1; + do { + if (ms[j] > i) { + mi.set(); + ++b; + } + + ++mi; + if (++j == n) { + j = 0; + } + } while (j != d); +} + +} // namespace detail +} // namespace chilbert + +#endif diff --git a/chilbert/detail/operations.hpp b/chilbert/detail/operations.hpp new file mode 100644 index 0000000..94635b4 --- /dev/null +++ b/chilbert/detail/operations.hpp @@ -0,0 +1,168 @@ +/* + Copyright (C) 2018 David Robillard + Copyright (C) 2006-2007 Chris Hamilton + + This program is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License as published by the Free Software + Foundation, either version 2 of the License, or (at your option) any later + version. + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + details. + + You should have received a copy of the GNU General Public License along with + this program. If not, see . +*/ + +#ifndef CHILBERT_OPERATIONS_HPP +#define CHILBERT_OPERATIONS_HPP + +#include "chilbert/detail/traits.hpp" +#include "chilbert/operators.hpp" + +#include +#include +#include +#include +#include + +namespace chilbert { +namespace detail { + +/// Reset all bits in `field` +template +std::enable_if_t::value> +reset_bits(T& field) +{ + field = static_cast(0); +} + +/// Reset all bits in `field` +template +std::enable_if_t> +reset_bits(T& field) +{ + field.reset(); +} + +/// Return the `index`th bit in `field` +template +std::enable_if_t::value, bool> +test_bit(const T& field, const size_t index) +{ + assert(size_t(index) < sizeof(field) * CHAR_BIT); + return field & (T{1} << index); +} + +/// Return the `index`th bit in `field` +template +std::enable_if_t, bool> +test_bit(const T& field, const size_t index) +{ + return field.test(index); +} + +/// Set the `index`th bit in `field` +template +std::enable_if_t::value> +set_bit(T& field, const size_t index) +{ + assert(size_t(index) < sizeof(field) * CHAR_BIT); + field |= (T{1} << index); + assert(test_bit(field, index)); +} + +/// Set the `index`th bit in `field` to `value` +template +std::enable_if_t::value> +set_bit(T& field, const size_t index, const bool value) +{ + assert(size_t(index) < sizeof(field) * CHAR_BIT); + field ^= (-T{value} ^ field) & (T{1U} << index); + assert(test_bit(field, index) == value); +} + +/// Set the `index`th bit in `field` +template +std::enable_if_t> +set_bit(T& field, const size_t index) +{ + field.set(index); +} + +/// Set the `index`th bit in `field` to `value` +template +std::enable_if_t> +set_bit(T& field, const size_t index, const bool value) +{ + field.set(index, value); +} + +/// Return 1 + the index of the least significant 1-bit of `field`, or zero +template +int pop_count(const T& field); + +template <> +int +pop_count(const unsigned long& field) +{ + return __builtin_popcountl(field); +} + +template <> +int +pop_count(const unsigned long long& field) +{ + return __builtin_popcountll(field); +} + +/// Return 1 + the index of the least significant 1-bit of `field`, or zero +template +int find_first(const T field); + +template <> +int +find_first(const unsigned long field) +{ + return __builtin_ffsl(static_cast(field)); +} + +template <> +int +find_first(const unsigned long long field) +{ + return __builtin_ffsll(static_cast(field)); +} + +/// Calculates the Gray Code of `value` in place +template +std::enable_if_t> gray_code(T& value); + +/// Implementation of grayCode for any integral type +template +std::enable_if_t::value> +gray_code(T& value) +{ + value ^= (value >> 1); +} + +/// Calculates the inverse Gray Code of `value` in place +template +std::enable_if_t> gray_code_inv(T& value); + +/// Implementation of gray_code_inv for any integral type +template +std::enable_if_t::value> +gray_code_inv(T& value) +{ + for (T shift = 1; shift < sizeof(T) * CHAR_BIT; shift <<= 1) { + value ^= (value >> shift); + } +} + +} // namespace detail +} // namespace chilbert + +#endif diff --git a/chilbert/detail/traits.hpp b/chilbert/detail/traits.hpp new file mode 100644 index 0000000..0cb0f03 --- /dev/null +++ b/chilbert/detail/traits.hpp @@ -0,0 +1,39 @@ +/* + Copyright (C) 2018 David Robillard + Copyright (C) 2006-2007 Chris Hamilton + + This program is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License as published by the Free Software + Foundation, either version 2 of the License, or (at your option) any later + version. + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + details. + + You should have received a copy of the GNU General Public License along with + this program. If not, see . +*/ + +#ifndef CHILBERT_TRAITS_HPP +#define CHILBERT_TRAITS_HPP + +namespace chilbert { +namespace detail { + +/// Member `value` is true iff T is a chilbert bitset +template +struct is_bitvec +{ + static constexpr bool value = false; +}; + +/// True iff T is a chilbert bitset +template +static constexpr bool is_bitvec_v = is_bitvec::value; + +} // namespace detail +} // namespace chilbert + +#endif diff --git a/chilbert/operators.hpp b/chilbert/operators.hpp new file mode 100644 index 0000000..efad7f4 --- /dev/null +++ b/chilbert/operators.hpp @@ -0,0 +1,96 @@ +/* + Copyright (C) 2018 David Robillard + Copyright (C) 2006-2007 Chris Hamilton + + This program is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License as published by the Free Software + Foundation, either version 2 of the License, or (at your option) any later + version. + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + details. + + You should have received a copy of the GNU General Public License along with + this program. If not, see . +*/ + +#ifndef CHILBERT_OPERATORS_HPP +#define CHILBERT_OPERATORS_HPP + +#include "chilbert/detail/traits.hpp" + +#include +#include + +namespace chilbert { + +using detail::is_bitvec_v; + +template +std::enable_if_t, T> operator&(const T& lhs, const T& rhs) +{ + T r{lhs}; + r &= rhs; + return r; +} + +template +std::enable_if_t, T> +operator|(const T& lhs, const T& rhs) +{ + T r{lhs}; + r |= rhs; + return r; +} + +template +std::enable_if_t, T> +operator^(const T& lhs, const T& rhs) +{ + T r{lhs}; + r ^= rhs; + return r; +} + +template +std::enable_if_t, T> +operator~(const T& vec) +{ + T r{vec}; + r.flip(); + return r; +} + +template +std::enable_if_t, T> +operator<<(const T& vec, const size_t bits) +{ + T r{vec}; + r <<= bits; + return r; +} + +template +std::enable_if_t, T> +operator>>(const T& vec, const size_t bits) +{ + T r{vec}; + r >>= bits; + return r; +} + +template +inline std::enable_if_t, std::ostream>& +operator<<(std::ostream& os, const T& vec) +{ + for (size_t i = 0; i < vec.size(); ++i) { + os << vec.test(vec.size() - i - 1); + } + return os; +} + +} // namespace chilbert + +#endif diff --git a/test/test_bitvec.cpp b/test/test_bitvec.cpp index 1596042..bd944b5 100644 --- a/test/test_bitvec.cpp +++ b/test/test_bitvec.cpp @@ -252,13 +252,13 @@ test_gray_code(Context& ctx) { const T v = make_random_bitvec(ctx); T r = v; - gray_code(r); + chilbert::detail::gray_code(r); if (N > 0) { assert(N == 1 || r == (v ^ (v >> 1))); T s = r; - gray_code_inv(s); + chilbert::detail::gray_code_inv(s); assert(s == v); } } @@ -271,11 +271,11 @@ test_comparison(Context&) T b = make_zero_bitvec(); for (size_t bit = 1; bit < N; ++bit) { - set_bit(a, bit, 1); + chilbert::detail::set_bit(a, bit, 1); for (size_t i = 0; i < bit; ++i) { - set_bit(a, i, rand() % 2 == 0); - set_bit(b, i, rand() % 2 == 0); + chilbert::detail::set_bit(a, i, rand() % 2 == 0); + chilbert::detail::set_bit(b, i, rand() % 2 == 0); } assert(b < a); diff --git a/test/test_gray_code_rank.cpp b/test/test_gray_code_rank.cpp index 0c06a33..9c903c1 100644 --- a/test/test_gray_code_rank.cpp +++ b/test/test_gray_code_rank.cpp @@ -21,9 +21,9 @@ #include "chilbert/BoundedBitVec.hpp" #include "chilbert/DynamicBitVec.hpp" -#include "chilbert/GrayCodeRank.hpp" #include "chilbert/SmallBitVec.hpp" #include "chilbert/StaticBitVec.hpp" +#include "chilbert/detail/gray_code_rank.hpp" #include #include @@ -39,7 +39,7 @@ get_mask(const std::array& ms, const size_t d, const size_t step) { T mask = make_zero_bitvec(); size_t b; - chilbert::extract_mask(ms.data(), D, d, step, mask, b); + chilbert::detail::extract_mask(ms.data(), D, d, step, mask, b); assert(b == mask.count()); @@ -75,21 +75,21 @@ test_gray_code_rank(Context& ctx) const auto b = make_random_bitvec(ctx); auto ga = a; - gray_code(ga); + chilbert::detail::gray_code(ga); auto gb = b; - gray_code(gb); + chilbert::detail::gray_code(gb); // Calculate gray code ranks auto ra = make_zero_bitvec(); - chilbert::gray_code_rank(mask, ga, D, ra); + chilbert::detail::gray_code_rank(mask, ga, D, ra); auto rb = make_zero_bitvec(); - chilbert::gray_code_rank(mask, gb, D, rb); + chilbert::detail::gray_code_rank(mask, gb, D, rb); // Ensure ranks have at most mask.count() bits auto max = make_zero_bitvec(); - set_bit(max, mask.count(), 1); + chilbert::detail::set_bit(max, mask.count(), 1); assert(ra < max); assert(rb < max); @@ -102,11 +102,12 @@ test_gray_code_rank(Context& ctx) const auto pat = ~mask; auto ga_out = make_zero_bitvec(); auto gag_out = make_zero_bitvec(); - gray_code_rank_inv(mask, pat, ra, D, mask.count(), gag_out, ga_out); + chilbert::detail::gray_code_rank_inv( + mask, pat, ra, D, mask.count(), gag_out, ga_out); assert((ga_out & mask) == (ga & mask)); auto gag_check = ga_out; - gray_code(gag_check); + chilbert::detail::gray_code(gag_check); assert(gag_check == gag_out); } } diff --git a/test/test_hilbert.cpp b/test/test_hilbert.cpp index a217dc2..307dbae 100644 --- a/test/test_hilbert.cpp +++ b/test/test_hilbert.cpp @@ -21,9 +21,9 @@ #include "chilbert/BoundedBitVec.hpp" #include "chilbert/DynamicBitVec.hpp" -#include "chilbert/Hilbert.hpp" #include "chilbert/SmallBitVec.hpp" #include "chilbert/StaticBitVec.hpp" +#include "chilbert/chilbert.hpp" #include -- cgit v1.2.1