BitStore.h

Go to the documentation of this file.

#pragma once
// SPDX-FileCopyrightText: 2025 Nessan Fitzmaurice <nzznfitz+gh@icloud.com>
// SPDX-License-Identifier: MIT

 
#include <gf2/BitRef.h>
#include <gf2/assert.h>
#include <gf2/Iterators.h>
#include <gf2/RNG.h>
#include <gf2/unsigned_word.h>
#include <gf2/store_word.h>
 
#include <bitset>
#include <concepts>
#include <format>
#include <iostream>
#include <limits>
#include <optional>
#include <ranges>
#include <string>
namespace gf2 {

template<typename Store>
class BitStore {
private:
    // CRTP helper method that returns a reference to the derived type.
    constexpr Store& store() { return static_cast<Store&>(*this); }
 
    // CRTP helper method that returns a const reference to the derived type.
    constexpr const Store& store() const { return static_cast<const Store&>(*this); }
 
public:
    using word_type = typename store_word<Store>::word_type;

    static constexpr u8 bits_per_word = BITS<word_type>;
 
    // A default constructor plus the rule of 5 to stop picky compilers from complaining ...
    BitStore() = default;
    virtual ~BitStore() = default;
    BitStore(const BitStore&) = default;
    BitStore& operator=(const BitStore&) = default;
    BitStore(BitStore&&) = default;
    BitStore& operator=(BitStore&&) = default;


    constexpr usize size() const { return store().size(); }

    constexpr usize words() const { return store().words(); }

    constexpr word_type word(usize i) const { return store().word(i); }

    constexpr void set_word(usize i, word_type value) { store().set_word(i, value); }

    constexpr const word_type* data() const { return store().data(); }

    constexpr word_type* data() { return store().data(); }

    constexpr u8 offset() const { return store().offset(); }


    constexpr bool get(usize i) const {
        gf2_debug_assert(i < size(), "Index {} is out of bounds for a bit store of length {}", i, size());
        auto [word_index, mask] = index_and_mask<word_type>(i);
        return word(word_index) & mask;
    }

    constexpr bool operator[](usize index) const { return get(index); }

    constexpr bool front() const {
        gf2_debug_assert(!is_empty(), "The store is empty!");
        return get(0);
    }

    constexpr bool back() const {
        gf2_debug_assert(!is_empty(), "The store is empty!");
        return get(size() - 1);
    }


    auto set(usize index, bool value = true) {
        gf2_debug_assert(index < size(), "Index {} is out of bounds for a bit store of length {}", index, size());
        auto [word_index, mask] = index_and_mask<word_type>(index);
        auto word_value = word(word_index);
        auto bit_value = (word_value & mask) != 0;
        if (bit_value != value) set_word(word_index, word_value ^ mask);
        return *this;
    }

    constexpr auto operator[](usize index) { return BitRef<Store>{this, index}; }

    auto flip(usize index) {
        gf2_debug_assert(index < size(), "Index {} is out of bounds for a bit store of length {}", index, size());
        auto [word_index, mask] = index_and_mask<word_type>(index);
        set_word(word_index, word(word_index) ^ mask);
        return *this;
    }

    constexpr auto swap(usize i0, usize i1) {
        gf2_debug_assert(i0 < size(), "index {} is out of bounds for a store of length {}", i0, size());
        gf2_debug_assert(i1 < size(), "index {} is out of bounds for a store of length {}", i1, size());
        if (i0 != i1) {
            auto [w0, mask0] = index_and_mask<word_type>(i0);
            auto [w1, mask1] = index_and_mask<word_type>(i1);
            auto word0 = word(w0);
            auto word1 = word(w1);
            auto val0 = (word0 & mask0) != 0;
            auto val1 = (word1 & mask1) != 0;
            if (val0 != val1) {
                if (w0 == w1) {
                    // Both bits are in the same word
                    set_word(w0, word0 ^ mask0 ^ mask1);
                } else {
                    // BitsIter are in different words
                    set_word(w0, word0 ^ mask0);
                    set_word(w1, word1 ^ mask1);
                }
            }
        }
        return *this;
    }


    constexpr bool is_empty() const { return size() == 0; }

    constexpr bool any() const {
        for (auto i = 0uz; i < words(); ++i)
            if (word(i) != 0) return true;
        return false;
    }

    constexpr bool all() const {
        auto num_words = words();
        if (num_words == 0) { return true; }
 
        // Check the fully occupied words ...
        for (auto i = 0uz; i < num_words - 1; ++i) {
            if (word(i) != MAX<word_type>) return false;
        }
 
        // Last word might not be fully occupied ...
        auto unused_bits = bits_per_word - (size() % bits_per_word);
        auto last_word_max = MAX<word_type> >> unused_bits;
        if (word(num_words - 1) != last_word_max) return false;
 
        return true;
    }

    constexpr bool none() const { return !any(); }


    auto set_all(bool value = true) {
        auto word_value = value ? MAX<word_type> : word_type{0};
        for (auto i = 0uz; i < words(); ++i) set_word(i, word_value);
        return *this;
    }

    auto flip_all() {
        for (auto i = 0uz; i < words(); ++i) set_word(i, static_cast<word_type>(~word(i)));
        return *this;
    }


    template<unsigned_word Src>
    auto copy(Src src) {
        constexpr auto src_bits = BITS<Src>;
        gf2_assert(size() == src_bits, "Lengths do not match: {} != {}.", size(), src_bits);
 
        if constexpr (src_bits <= bits_per_word) {
            // The source fits into a single one of our words
            set_word(0, static_cast<word_type>(src));
        } else {
            // The source spans some number of our words (which we assume is an integer number).
            constexpr auto num_words = src_bits / bits_per_word;
            for (auto i = 0uz; i < num_words; ++i) {
                auto shift = i * bits_per_word;
                auto word_value = static_cast<word_type>((src >> shift) & MAX<word_type>);
                set_word(i, word_value);
            }
        }
        return *this;
    }

    template<typename SrcStore>
    auto copy(const BitStore<SrcStore>& src) {
        // The sizes must match ...
        gf2_assert(size() == src.size(), "Lengths do not match: {} != {}.", size(), src.size());
 
        // What is the source word type (without any const qualifier)?
        using src_word_type = std::remove_const_t<typename store_word<SrcStore>::word_type>;
 
        // Numbers of bits per respective words
        constexpr auto word_bits = BITS<word_type>;
        constexpr auto src_word_bits = BITS<src_word_type>;
 
        // What we do next depends on the size of the respective words.
        if constexpr (word_bits == src_word_bits) {
            // Easiest case: Same sized words -- we can copy whole words directly ...
            for (auto i = 0uz; i < words(); ++i) {
                auto src_word = static_cast<word_type>(src.word(i));
                set_word(i, src_word);
            }
        } else if constexpr (word_bits > src_word_bits) {
            // Our word type is larger -- some number of source words fit into each one of ours.
            // We assume that our word size is an integer multiple of the source word size.
            constexpr auto ratio = word_bits / src_word_bits;
 
            auto src_words = src.words();
            for (auto i = 0uz; i < words(); ++i) {
 
                // Combine `ratio` source words into a single destination word being careful not to run off the end ...
                word_type dst_word = 0;
                for (auto j = 0uz; j < ratio; ++j) {
                    auto      src_index = i * ratio + j;
                    word_type src_word = 0;
                    if (src_index < src_words) src_word = static_cast<word_type>(src.word(src_index));
                    dst_word |= (src_word << (j * src_word_bits));
                }
                set_word(i, dst_word);
            }
        } else {
            // Our word size is smaller -- each source word becomes multiple destination words.
            // We assume that the source word size is an integer multiple of our word size.
            constexpr auto ratio = src_word_bits / word_bits;
 
            auto dst_words = words();
            for (auto src_index = 0uz; src_index < src.words(); ++src_index) {
                auto src_word = src.word(src_index);
 
                // Split the source word into `ratio` destination words ...
                for (auto j = 0uz; j < ratio; ++j) {
                    auto dst_index = src_index * ratio + j;
                    if (dst_index >= dst_words) break;
                    auto shift = j * word_bits;
                    auto dst_word = static_cast<word_type>((src_word >> shift) & MAX<word_type>);
                    set_word(dst_index, dst_word);
                }
            }
        }
        return *this;
    }

    template<usize N>
    auto copy(const std::bitset<N>& src) {
        // The sizes must match
        gf2_assert(size() == N, "Lengths do not match: {} != {}.", size(), N);
 
        // Doing this the dumb way as there is no standard access to the bitset's underlying store
        set_all(false);
        for (auto i = 0uz; i < N; ++i)
            if (src[i]) set(i);
        return *this;
    }

    auto fill(std::invocable<usize> auto f) {
        set_all(false);
        for (auto i = 0uz; i < size(); ++i)
            if (f(i) == true) set(i);
        return *this;
    }

    auto random_fill(double p = 0.5, u64 seed = 0) {
        // Keep a single static RNG per thread for all calls to this method, seeded with entropy on the first call.
        thread_local RNG rng;
 
        // Edge case handling ...
        if (p < 0) return set_all(false);
 
        // Scale p by 2^64 to remove floating point arithmetic from the main loop below.
        // If we determine p rounds to 1 then we can just set all elements to 1 and return early.
        p = p * 0x1p64 + 0.5;
        if (p >= 0x1p64) return set_all();
 
        // p does not round to 1 so we use a 64-bit URNG and check each draw against the 64-bit scaled p.
        auto scaled_p = static_cast<u64>(p);
 
        // If a seed was provided, set the RNG's seed to it. Otherwise, we carry on from where we left off.
        u64 old_seed = rng.seed();
        if (seed != 0) rng.set_seed(seed);
 
        set_all(false);
        for (auto i = 0uz; i < size(); ++i)
            if (rng.u64() < scaled_p) set(i);
 
        // Restore the old seed if necessary.
        if (seed != 0) rng.set_seed(old_seed);
        return *this;
    }


    constexpr usize count_ones() const {
        usize count = 0;
        for (auto i = 0uz; i < words(); ++i) count += static_cast<usize>(gf2::count_ones(word(i)));
        return count;
    }

    constexpr usize count_zeros() const { return size() - count_ones(); }

    constexpr usize leading_zeros() const {
        // Note: Even if the last word is not fully occupied, we know unused bits are set to 0.
        for (auto i = 0uz; i < words(); ++i) {
            auto w = word(i);
            if (w != 0) return i * bits_per_word + static_cast<usize>(gf2::trailing_zeros(w));
        }
        return size();
    }

    constexpr usize trailing_zeros() const {
        if (is_empty()) return 0;
 
        // The last word may not be fully occupied so we need to subtract the number of unused bits.
        auto unused_bits = bits_per_word - (size() % bits_per_word);
        auto num_words = words();
        auto i = num_words;
        while (i--) {
            auto w = word(i);
            if (w != 0) {
                auto whole_count = (num_words - i - 1) * bits_per_word;
                auto partial_count = static_cast<usize>(gf2::leading_zeros(w)) - unused_bits;
                return whole_count + partial_count;
            }
        }
        return size();
    }


    std::optional<usize> first_set() const {
        // Iterate forward looking for a word with a set bit and use the lowest of those ...
        // Remember that any unused bits in the final word are guaranteed to be unset.
        for (auto i = 0uz; i < words(); ++i) {
            if (auto loc = lowest_set_bit(word(i))) return i * bits_per_word + loc.value();
        }
        return {};
    }

    std::optional<usize> last_set() const {
        // Iterate backward looking for a word with a set bit and use the highest of those ...
        // Remember that any unused bits in the final word are guaranteed to be unset.
        for (auto i = words(); i--;) {
            if (auto loc = highest_set_bit(word(i))) return i * bits_per_word + loc.value();
        }
        return {};
    }

    std::optional<usize> next_set(usize index) const {
        // Start our search at index + 1.
        index = index + 1;
 
        // Perhaps we are off the end? (This also handles the case of an empty type).
        if (index >= size()) return {};
 
        // Where is that starting index located in the word store?
        auto [word_index, bit] = index_and_offset<word_type>(index);
 
        // Iterate forward looking for a word with a new set bit and use the lowest one ...
        for (auto i = word_index; i < words(); ++i) {
            auto w = word(i);
            if (i == word_index) {
                // First word -- turn all the bits before our starting bit into zeros so we don't see them as set.
                reset_bits(w, 0, bit);
            }
            if (auto loc = lowest_set_bit(w)) return i * bits_per_word + loc.value();
        }
        return {};
    }

    std::optional<usize> previous_set(usize index) const {
        // Edge case: If the store is empty or we are at the start, there are no previous set bits.
        if (is_empty() || index == 0) return {};
 
        // Silently fix large indices and also adjust the index down a slot.
        if (index > size()) index = size();
        index--;
 
        // Where is that starting index located in the word store?
        auto [word_index, bit] = index_and_offset<word_type>(index);
 
        // Iterate backwards looking for a word with a new set bit and use the highest one ...
        for (auto i = word_index + 1; i--;) {
            auto w = word(i);
            if (i == word_index) {
                // First word -- turn all higher bits after our starting bit into zeros so we don't see them as set.
                reset_bits(w, bit + 1, bits_per_word);
            }
            if (auto loc = highest_set_bit(w)) return i * bits_per_word + loc.value();
        }
        return {};
    }


    std::optional<usize> first_unset() const {
        // Iterate forward looking for a word with a unset bit and use the lowest of those ...
        for (auto i = 0uz; i < words(); ++i) {
            auto w = word(i);
            if (i == words() - 1) {
                // Final word may have some unused zero bits that we need to replace with ones.
                auto last_occupied_bit = bit_offset<word_type>(size() - 1);
                set_bits(w, last_occupied_bit + 1, bits_per_word);
            }
            if (auto loc = lowest_unset_bit(w)) return i * bits_per_word + loc.value();
        }
        return {};
    }

    std::optional<usize> last_unset() const {
        // Iterate backwards looking for a word with a unset bit and use the highest of those ...
        for (auto i = words(); i--;) {
            auto w = word(i);
            if (i == words() - 1) {
                // Final word may have some unused zero bits that we need to replace with ones.
                auto last_occupied_bit = bit_offset<word_type>(size() - 1);
                set_bits(w, last_occupied_bit + 1, bits_per_word);
            }
            if (auto loc = highest_unset_bit(w)) return i * bits_per_word + loc.value();
        }
        return {};
    }

    std::optional<usize> next_unset(usize index) const {
        // Start our search at index + 1.
        index = index + 1;
 
        // Perhaps we are off the end? (This also handles the case of an empty type).
        if (index >= size()) return {};
 
        // Where is that starting index located in the word store?
        auto [word_index, bit] = index_and_offset<word_type>(index);
 
        // Iterate forward looking for a word with a new unset bit and use the lowest of those ...
        for (auto i = word_index; i < words(); ++i) {
            auto w = word(i);
            if (i == word_index) {
                // First word -- turn all the bits before our starting bit into ones so we don't see them as unset.
                set_bits(w, 0, bit);
            }
            if (i == words() - 1) {
                // Final word may have some unused zero bits that we need to replace with ones.
                auto last_occupied_bit = bit_offset<word_type>(size() - 1);
                set_bits(w, last_occupied_bit + 1, bits_per_word);
            }
            if (auto loc = lowest_unset_bit(w)) return i * bits_per_word + loc.value();
        }
        return {};
    }

    std::optional<usize> previous_unset(usize index) const {
        // Edge case: If the store is empty or we are at the start, there are no previous set bits.
        if (is_empty() || index == 0) return {};
 
        // Silently fix large indices and also adjust the index down a slot.
        if (index > size()) index = size();
        index--;
 
        // Where is that starting index located in the word store?
        auto [word_index, bit] = index_and_offset<word_type>(index);
 
        // Iterate backwards looking for a word with a new unset bit and use the highest of those ...
        for (auto i = word_index + 1; i--;) {
            auto w = word(i);
            if (i == word_index) {
                // First word -- set all the bits after our starting bit into ones so we don't see them as unset.
                set_bits(w, bit + 1, bits_per_word);
            }
            if (i == words() - 1) {
                // Final word may have some unused zero bits that we need to replace with ones.
                auto last_occupied_bit = bit_offset<word_type>(size() - 1);
                set_bits(w, last_occupied_bit + 1, bits_per_word);
            }
            if (auto loc = highest_unset_bit(w)) return i * bits_per_word + loc.value();
        }
        return {};
    }


    constexpr auto bits() const { return BitsIter<Store, true>{this, 0}; }

    constexpr auto bits() { return BitsIter<Store, false>{this, 0}; }

    constexpr auto set_bit_indices() const { return SetBitsIter{this}; }

    constexpr auto unset_bit_indices() const { return UnsetBitsIter{this}; }

    constexpr auto store_words() const { return WordsIter{this}; }

    constexpr auto to_words() const { return std::ranges::to<std::vector>(store_words()); }


    constexpr auto span(usize begin, usize end) const {
        // Check that the range is correctly formed and non-empty.
        gf2_assert(begin < end, "bit-span [{}, {}) is invalid.", begin, end);
 
        // Check that the range does not extend beyond the end of the bit-store.
        gf2_assert(end <= size(), "bit-span end {} extends beyond the store end {}.", end, size());
 
        // Get the zero-based word index and bit offset in that word for the begin location.
        auto [index, bit_offset] = index_and_offset<word_type>(begin);
 
        // If the store is already a span over another store, we need to adjust those values.
        bit_offset += offset();
        if (bit_offset >= bits_per_word) {
            index += 1;
            bit_offset = bit_offset % bits_per_word;
        }
 
        // Return the bit-span over the underlying word data.
        return BitSpan<const word_type>(data() + index, bit_offset, end - begin);
    }

    constexpr auto span(usize begin, usize end) {
        // Check that the range is correctly formed and non-empty.
        gf2_assert(begin < end, "bit-span [{}, {}) is invalid.", begin, end);
 
        // Check that the range does not extend beyond the end of the bit-store.
        gf2_assert(end <= size(), "bit-span end {} extends beyond the store end {}.", end, size());
 
        // Get the zero-based word index and bit offset in that word for the begin location.
        auto [index, bit_offset] = index_and_offset<word_type>(begin);
 
        // If the store is already a span over another store, we need to adjust those values.
        bit_offset += offset();
        if (bit_offset >= bits_per_word) {
            index += 1;
            bit_offset = bit_offset % bits_per_word;
        }
 
        // Return the bit-span over the underlying word data.
        return BitSpan<word_type>(data() + index, bit_offset, end - begin);
    }


    constexpr auto sub(usize begin, usize end) const { return BitVec<word_type>::from(span(begin, end)); }

    constexpr void split_at(usize at, BitVec<word_type>& left, BitVec<word_type>& right) const {
        gf2_assert(at <= size(), "split point {} is beyond the end of the bit-vector", at);
        left.clear();
        right.clear();
        left.append(span(0, at));
        right.append(span(at, size()));
    }

    constexpr auto split_at(usize at) const {
        BitVec<word_type> left, right;
        split_at(at, left, right);
        return std::pair{left, right};
    }


    constexpr void riffle_into(BitVec<word_type>& dst) const {
        // Nothing to do if the bit-vector is empty or has just one element.
        // With two bits `ab` we fill `dst` with `a0b`.
        if (size() <= 1) {
            dst.copy(*this);
            return;
        }
 
        // Make sure `dst` is large enough to hold the riffled bits (a bit too big but we fix that below).
        dst.resize(2 * size());
        auto dst_words = dst.words();
 
        // Riffle each word in `this` into two adjacent words in `dst`.
        for (auto i = 0uz; i < words(); ++i) {
            auto [lo, hi] = riffle(word(i));
            dst.set_word(2 * i, lo);
 
            // Note that `hi` may be completely superfluous ...
            if (2 * i + 1 < dst_words) dst.set_word(2 * i + 1, hi);
        }
 
        // If this bit-store was say `abcde` then `dst` will now be `a0b0c0d0e0`. Pop the last 0.
        dst.pop();
    }

    constexpr auto riffled() const {
        BitVec<word_type> result;
        riffle_into(result);
        return result;
    }


    constexpr void operator<<=(usize shift) {
        // Edge cases where there is nothing to do.
        if (shift == 0 || size() == 0) return;
 
        // Perhaps we are shifting all the bits out and are left with all zeros.
        if (shift >= size()) {
            set_all(false);
            return;
        }
 
        // For larger shifts, we can efficiently shift by whole words first.
        usize word_shift = shift / bits_per_word;
        usize end_word = words() - word_shift;
 
        // Do the whole word shifts first, pushing in zero words to fill the empty slots.
        if (word_shift > 0) {
            // Shift whole words -- starting at the beginning of the store.
            for (auto i = 0uz; i < end_word; ++i) set_word(i, word(i + word_shift));
 
            // Fill in the high order words with zeros.
            for (auto i = end_word; i < words(); ++i) set_word(i, 0);
 
            // How many bits are left to shift?
            shift -= word_shift * bits_per_word;
        }
 
        // Perhaps there are some partial word shifts left to do.
        if (shift != 0) {
            // Do the "interior" words where the shift moves bits from one word to the next.
            if (end_word > 0) {
                auto shift_complement = bits_per_word - shift;
                for (auto i = 0uz; i < end_word - 1; ++i) {
                    auto lo = static_cast<word_type>(word(i) >> shift);
                    auto hi = static_cast<word_type>(word(i + 1) << shift_complement);
                    set_word(i, lo | hi);
                }
            }
 
            // Do the last word.
            auto value = word(end_word - 1);
            set_word(end_word - 1, static_cast<word_type>(value >> shift));
        }
    }

    constexpr void operator>>=(usize shift) {
        // Edge cases where there is nothing to do.
        if (shift == 0 || size() == 0) return;
 
        // Perhaps we are shifting all the bits out and are left with all zeros.
        if (shift >= size()) {
            set_all(false);
            return;
        }
 
        // For larger shifts, we can efficiently shift by whole words first.
        usize word_shift = shift / bits_per_word;
 
        // Do the whole word shifts first, pushing in zero words to fill the empty slots.
        if (word_shift > 0) {
            // Shift whole words -- starting at the end of the store.
            for (auto i = words(); i-- > word_shift;) set_word(i, word(i - word_shift));
 
            // Fill in the low order words with zeros.
            for (auto i = 0uz; i < word_shift; ++i) set_word(i, 0);
 
            // How many bits are left to shift?
            shift -= word_shift * bits_per_word;
        }
 
        // Perhaps there are some partial word shifts left to do.
        if (shift != 0) {
            // Do the "interior" words where the shift moves bits from one word to the next.
            auto shift_complement = bits_per_word - shift;
            for (auto i = words(); i-- > word_shift + 1;) {
                auto lo = static_cast<word_type>(word(i - 1) >> shift_complement);
                auto hi = static_cast<word_type>(word(i) << shift);
                set_word(i, lo | hi);
            }
        }
 
        // Do the "first" word.
        auto value = word(word_shift);
        set_word(word_shift, static_cast<word_type>(value << shift));
    }

    constexpr auto operator<<(usize shift) const {
        auto result = BitVec<word_type>::from(store());
        result <<= shift;
        return result;
    }

    constexpr auto operator>>(usize shift) const {
        auto result = BitVec<word_type>::from(store());
        result >>= shift;
        return result;
    }


    template<typename Rhs>
        requires store_words_match<Store, Rhs>
    void operator^=(const BitStore<Rhs>& rhs) {
        gf2_assert(size() == rhs.size(), "Lengths do not match: {} != {}.", size(), rhs.size());
        for (auto i = 0uz; i < words(); ++i) set_word(i, word(i) ^ rhs.word(i));
    }

    template<typename Rhs>
        requires store_words_match<Store, Rhs>
    void operator&=(const BitStore<Rhs>& rhs) {
        gf2_assert(size() == rhs.size(), "Lengths do not match: {} != {}.", size(), rhs.size());
        for (auto i = 0uz; i < words(); ++i) set_word(i, word(i) & rhs.word(i));
    }

    template<typename Rhs>
        requires store_words_match<Store, Rhs>
    void operator|=(const BitStore<Rhs>& rhs) {
        gf2_assert(size() == rhs.size(), "Lengths do not match: {} != {}.", size(), rhs.size());
        for (auto i = 0uz; i < words(); ++i) set_word(i, word(i) | rhs.word(i));
    }

    auto operator~() const {
        auto result = BitVec<word_type>::from(store());
        result.flip_all();
        return result;
    }

    template<typename Rhs>
        requires store_words_match<Store, Rhs>
    auto operator^(const BitStore<Rhs>& rhs) const {
        auto result = BitVec<word_type>::from(store());
        result ^= rhs;
        return result;
    }

    template<typename Rhs>
        requires store_words_match<Store, Rhs>
    auto operator&(const BitStore<Rhs>& rhs) const {
        auto result = BitVec<word_type>::from(store());
        result &= rhs;
        return result;
    }

    template<typename Rhs>
        requires store_words_match<Store, Rhs>
    auto operator|(const BitStore<Rhs>& rhs) const {
        auto result = BitVec<word_type>::from(store());
        result |= rhs;
        return result;
    }


    template<typename Rhs>
        requires store_words_match<Store, Rhs>
    void operator+=(const BitStore<Rhs>& rhs) {
        operator^=(rhs);
    }

    template<typename Rhs>
        requires store_words_match<Store, Rhs>
    void operator-=(const BitStore<Rhs>& rhs) {
        operator^=(rhs);
    }

    template<typename Rhs>
        requires store_words_match<Store, Rhs>
    auto operator+(const BitStore<Rhs>& rhs) const {
        auto result = BitVec<word_type>::from(store());
        result += rhs;
        return result;
    }

    template<typename Rhs>
        requires store_words_match<Store, Rhs>
    auto operator-(const BitStore<Rhs>& rhs) const {
        auto result = BitVec<word_type>::from(store());
        result -= rhs;
        return result;
    }


    inline std::string to_string(std::string_view sep = "", std::string_view pre = "",
                                 std::string_view post = "") const {
        return to_binary_string(sep, pre, post);
    }

    inline std::string to_pretty_string() const { return to_binary_string(",", "[", "]"); }

    std::string to_binary_string(std::string_view sep = "", std::string_view pre = "",
                                 std::string_view post = "") const {
        // Edge case ...
        if (is_empty()) { return std::format("{}{}", pre, post); }
 
        // Start with the raw binary string by iterating over the words ...
        // We preallocate enough space to hold the entire string (actually this is usually a bit bigger than we need).
        auto        num_words = words();
        std::string binary_string;
        binary_string.reserve(num_words * bits_per_word);
        for (auto i = 0uz; i < num_words; ++i) {
            auto w = reverse_bits(word(i));
            binary_string.append(gf2::to_binary_string(w));
        }
 
        // The last word may not be fully occupied and we need to remove the spurious zeros
        binary_string.resize(size());
 
        // If there is no separator, return early ...
        if (sep.empty()) { return std::format("{}{}{}", pre, binary_string, post); }
 
        // Build `result` with separators between each character ...
        std::string result;
        result.reserve(pre.size() + size() * (sep.size() + 1) + post.size());
        result.append(pre);
        for (auto i = 0uz; i < size(); ++i) {
            result.push_back(binary_string[i]);
            if (i != size() - 1) { result.append(sep); }
        }
        result.append(post);
        return result;
    }

    std::string to_hex_string() const {
        // Edge case: No bits in the store, return the empty string.
        if (is_empty()) return "";
 
        // The number of digits in the output string. Generally hexadecimal but the last may be to a lower base.
        auto digits = (size() + 3) / 4;
 
        // Preallocate space allowing for a possible lower base on the last digit such as "_2".
        std::string result;
        result.reserve(digits + 2);
 
        //  Reverse the bits in each word to vector-order and get the fully padded hex string representation.
        for (auto i = 0uz; i < words(); ++i) {
            auto w = reverse_bits(word(i));
            result.append(gf2::to_hex_string<word_type>(w));
        }
 
        // Last word may not be fully occupied and padded with spurious zeros so we truncate the output string.
        result.resize(digits);
 
        // Every four elements is encoded by a single hex digit but `size()` may not be a multiple of 4.
        auto k = size() % 4;
        if (k != 0) {
            // That last hex digit should really be encoded to a lower base -- 2, 4 or 8.
            // We compute the number represented by the trailing `k` elements in the bit-vector.
            int num = 0;
            for (auto i = 0uz; i < k; ++i)
                if (get(size() - 1 - i)) num |= 1 << i;
 
            // Convert that number to hex & use it to *replace* the last hex digit in our `result` string.
            // Also append the appropriate base to the output string so that the last digit can be interpreted properly.
            result.resize(result.size() - 1);
            result.append(std::format("{:X}.{}", num, 1 << k));
        }
        return result;
    }

    std::string describe() const {
        std::string result;
        result.reserve(1000);
        result += std::format("binary format:        {}\n", to_binary_string());
        result += std::format("hex format:           {}\n", to_hex_string());
        result += std::format("number of bits:       {}\n", size());
        result += std::format("number of set bits:   {}\n", count_ones());
        result += std::format("number of unset bits: {}\n", count_zeros());
        result += std::format("bits per word:        {}\n", bits_per_word);
        result += std::format("word count:           {}\n", words());
        auto words_vec = to_words();
        result += std::format("words in hex:         {::#x}\n", words_vec);
        ;
        return result;
    }
 
    // The usual output stream operator for a bit-store
    constexpr std::ostream& operator<<(std::ostream& s) const { return s << to_string(); }
};
 
// --------------------------------------------------------------------------------------------------------------------
// Equality operator ...
// -------------------------------------------------------------------------------------------------------------------

// # Example
template<typename Lhs, typename Rhs>
constexpr bool
operator==(const BitStore<Lhs>& lhs, const BitStore<Rhs>& rhs) {
    if constexpr (!store_words_match<Lhs, Rhs>) return false;
    if (&lhs != &rhs) {
        if (lhs.size() != rhs.size()) return false;
        for (auto i = 0uz; i < lhs.words(); ++i)
            if (lhs.word(i) != rhs.word(i)) return false;
    }
    return true;
}
 
// --------------------------------------------------------------------------------------------------------------------
// Other functions ...
// -------------------------------------------------------------------------------------------------------------------

template<typename Lhs, typename Rhs>
auto
join(const BitStore<Lhs>& lhs, const BitStore<Rhs>& rhs) {
    using word_type = typename store_word<Lhs>::word_type;
    auto result = BitVec<word_type>::from(lhs);
    result.append(rhs);
    return result;
}

template<typename Lhs, typename Rhs>
    requires store_words_match<Lhs, Rhs>
constexpr bool
dot(const BitStore<Lhs>& lhs, const BitStore<Rhs>& rhs) {
    gf2_debug_assert_eq(lhs.size(), rhs.size(), "Length mismatch {} != {}", lhs.size(), rhs.size());
    using word_type = typename store_word<Lhs>::word_type;
    auto sum = word_type{0};
    for (auto i = 0uz; i < lhs.words(); ++i) { sum ^= static_cast<word_type>(lhs.word(i) & rhs.word(i)); }
    return count_ones(sum) % 2 == 1;
}

template<typename Lhs, typename Rhs>
    requires store_words_match<Lhs, Rhs>
constexpr bool
operator*(const BitStore<Lhs>& lhs, const BitStore<Rhs>& rhs) {
    return dot(lhs, rhs);
}

template<typename Lhs, typename Rhs>
    requires store_words_match<Lhs, Rhs>
auto
convolve(const BitStore<Lhs>& lhs, const BitStore<Rhs>& rhs) {
    using word_type = typename store_word<Lhs>::word_type;
 
    // Edge case: if either store is empty then the convolution is empty.
    if (lhs.is_empty() || rhs.is_empty()) return BitVec<word_type>{};
 
    // Generally the result will have length `self.size() + other.size() - 1` (could be all zeros).
    auto result = BitVec<word_type>::zeros(lhs.size() + rhs.size() - 1);
 
    // If either vector is all zeros then the convolution is all zeros.
    if (lhs.none() || rhs.none()) return result;
 
    // Only need to consider words in `rhs` up to and including the one holding its final set bit.
    // We have already checked that `rhs` is not all zeros so we know there is a last set bit!
    auto rhs_words_end = word_index<word_type>(rhs.last_set().value()) + 1;
 
    // Initialize `result` by copying the live words from `rhs`
    for (auto i = 0uz; i < rhs_words_end; ++i) result.set_word(i, rhs.word(i));
 
    // Work backwards from the last set bit in `lhs` (which we know exists as we checked `lhs` is not all zeros).
    for (auto i = lhs.last_set().value(); i-- > 0;) {
        word_type prev = 0;
        for (auto j = 0uz; j < result.words(); ++j) {
            auto left = static_cast<word_type>(prev >> (BITS<word_type> - 1));
            prev = result.word(j);
            result.set_word(j, static_cast<word_type>(prev << 1) | left);
        }
 
        if (lhs.get(i)) {
            for (auto j = 0uz; j < rhs_words_end; ++j) result.set_word(j, result.word(j) ^ rhs.word(j));
        }
    }
    return result;
}
 
} // namespace gf2
 
// --------------------------------------------------------------------------------------------------------------------
// Specialise `std::formatter` to handle bit-stores ...
// -------------------------------------------------------------------------------------------------------------------

template<typename Store>
struct std::formatter<gf2::BitStore<Store>> {
 
    // Parse a bit-store format specifier where where we recognize {:p} and {:x}
    constexpr auto parse(const std::format_parse_context& ctx) {
        auto it = ctx.begin();
        while (it != ctx.end() && *it != '}') {
            switch (*it) {
                case 'p': m_pretty = true; break;
                case 'x': m_hex = true; break;
                default: m_error = true;
            }
            ++it;
        }
        return it;
    }
 
    // Push out a formatted bit-vector using the various @c to_string(...) methods in the class.
    template<class FormatContext>
    auto format(const gf2::BitStore<Store>& rhs, FormatContext& ctx) const {
        // Was there a format specification error?
        if (m_error) return std::format_to(ctx.out(), "'UNRECOGNIZED FORMAT SPECIFIER FOR BIT-STORE'");
 
        // Special handling requested?
        if (m_hex) return std::format_to(ctx.out(), "{}", rhs.to_hex_string());
        if (m_pretty) return std::format_to(ctx.out(), "{}", rhs.to_pretty_string());
 
        // Default
        return std::format_to(ctx.out(), "{}", rhs.to_string());
    }
 
    bool m_hex = false;
    bool m_pretty = false;
    bool m_error = false;
};