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/assert.h>
#include <gf2/Unsigned.h>
#include <gf2/RNG.h>
 
#include <algorithm>
#include <bitset>
#include <concepts>
#include <format>
#include <iterator>
#include <optional>
#include <ranges>
#include <string>
#include <string_view>
#include <utility>
#include <vector>
#include <type_traits>
 
// --------------------------------------------------------------------------------------------------------------------
// The BitStore concept ...
// --------------------------------------------------------------------------------------------------------------------
 
namespace gf2 {

template<typename Store>
concept BitStore = requires(Store store, const Store const_store) {
    typename Store::word_type;
    requires Unsigned<typename Store::word_type>;

    { const_store.size() } -> std::same_as<usize>;

    { const_store.store() } -> std::same_as<const typename Store::word_type*>;

    { store.store() } -> std::convertible_to<const typename Store::word_type*>;

    { const_store.offset() } -> std::same_as<u8>;

    { const_store.words() } -> std::same_as<usize>;

    { const_store.word(usize{}) } -> std::same_as<typename Store::word_type>;

    { store.set_word(usize{}, typename Store::word_type{}) } -> std::same_as<void>;
};
 
} // namespace gf2
 
// --------------------------------------------------------------------------------------------------------------------
// Forward Declarations of BitStore Related Types.
// --------------------------------------------------------------------------------------------------------------------
 
namespace gf2 {
 
// clang-format off
template<usize N, Unsigned Word> class BitArray;
template<Unsigned Word> class BitVector;
template<Unsigned Word> class BitSpan;
template<BitStore Store> class BitRef;
template<BitStore Store, bool is_const> class Bits;
template<BitStore Store> class SetBits;
template<BitStore Store> class UnsetBits;
template<BitStore Store> class Words;
// clang-format on
 
} // namespace gf2
 
// --------------------------------------------------------------------------------------------------------------------
// Free Functions for BitStore Types.
// --------------------------------------------------------------------------------------------------------------------
 
namespace gf2 {

template<BitStore Store>
constexpr bool
get(Store const& store, usize i) {
    gf2_debug_assert(i < store.size(), "Index {} is out of bounds for store of length {}", i, store.size());
    auto [word_index, word_mask] = index_and_mask<typename Store::word_type>(i);
    return store.word(word_index) & word_mask;
}

// # Panics
template<BitStore Store>
constexpr auto
ref(Store& store, usize i) {
    gf2_debug_assert(i < store.size(), "Index {} is out of bounds for store of length {}", i, store.size());
    return BitRef{&store, i};
}

template<BitStore Store>
constexpr bool
front(Store const& store) {
    gf2_debug_assert(!is_empty(store), "The store is empty!");
    return get(store, 0);
}

template<BitStore Store>
constexpr bool
back(Store const& store) {
    gf2_debug_assert(!is_empty(store), "The store is empty!");
    return get(store, store.size() - 1);
}

template<BitStore Store>
constexpr void
set(Store& store, usize i, bool value = true) {
    gf2_debug_assert(i < store.size(), "Index {} is out of bounds for store of length {}", i, store.size());
    auto [word_index, word_mask] = index_and_mask<typename Store::word_type>(i);
    auto word_value = store.word(word_index);
    auto bit_value = (word_value & word_mask) != 0;
    if (bit_value != value) store.set_word(word_index, word_value ^ word_mask);
}

template<BitStore Store>
constexpr void
flip(Store& store, usize i) {
    gf2_debug_assert(i < store.size(), "Index {} is out of bounds for store of length {}", i, store.size());
    auto [word_index, word_mask] = index_and_mask<typename Store::word_type>(i);
    auto word_value = store.word(word_index);
    store.set_word(word_index, word_value ^ word_mask);
}

template<BitStore Store>
constexpr void
swap(Store& store, usize i0, usize i1) {
    gf2_debug_assert(i0 < store.size(), "index {} is out of bounds for a store of length {}", i0, store.size());
    gf2_debug_assert(i1 < store.size(), "index {} is out of bounds for a store of length {}", i1, store.size());
    if (i0 != i1) {
        auto [w0, mask0] = index_and_mask<typename Store::word_type>(i0);
        auto [w1, mask1] = index_and_mask<typename Store::word_type>(i1);
        auto word0 = store.word(w0);
        auto word1 = store.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
                store.set_word(w0, word0 ^ mask0 ^ mask1);
            } else {
                // The two bits are in different words
                store.set_word(w0, word0 ^ mask0);
                store.set_word(w1, word1 ^ mask1);
            }
        }
    }
}


template<BitStore Store>
constexpr bool
is_empty(Store const& store) {
    return store.size() == 0;
}

template<BitStore Store>
constexpr bool
any(Store const& store) {
    for (auto i = 0uz; i < store.words(); ++i)
        if (store.word(i) != 0) return true;
    return false;
}

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

template<BitStore Store>
constexpr bool
none(Store const& store) {
    return !any(store);
}


template<BitStore Store>
constexpr void
set_all(Store& store, bool value = true) {
    using word_type = typename Store::word_type;
    auto word_value = value ? MAX<word_type> : word_type{0};
    for (auto i = 0uz; i < store.words(); ++i) store.set_word(i, word_value);
}

template<BitStore Store>
constexpr void
flip_all(Store& store) {
    using word_type = typename Store::word_type;
    for (auto i = 0uz; i < store.words(); ++i) store.set_word(i, static_cast<word_type>(~store.word(i)));
}


template<Unsigned Src, BitStore Store>
constexpr void
copy(Src src, Store& store) {
    constexpr auto src_bits = BITS<Src>;
    gf2_assert(store.size() == src_bits, "Lengths do not match: {} != {}.", store.size(), src_bits);
 
    using dst_type = typename Store::word_type;
    constexpr auto dst_bpw = BITS<dst_type>;
 
    if constexpr (src_bits <= dst_bpw) {
        // The source fits into a single store word. We can pad the source to get it to the correct size.
        // It happens that casting zero-pads --- but `set_word` only ever looks at the low `src_bits` piece anyway.
        store.set_word(0, static_cast<dst_type>(src));
    } else {
        // The source spans some integer number of store words (we assume it is an integer number).
        // Proceed by nibbling off destination-sized words from the larger source until its exhausted.
        constexpr auto dst_words = src_bits / dst_bpw;
        constexpr auto dst_mask = MAX<dst_type>;
        for (auto i = 0uz; i < dst_words; ++i, src >>= dst_bpw) {
            auto src_word = static_cast<dst_type>(src & dst_mask);
            store.set_word(i, src_word);
        }
    }
}

template<typename Iter, BitStore Store>
    requires std::is_unsigned_v<typename std::iterator_traits<Iter>::value_type>
constexpr void
copy(Iter src_begin, Iter src_end, Store& store) {
 
    // What are the two word types (without any const qualifiers)?
    using dst_type = std::remove_const_t<typename Store::word_type>;
    using src_type = typename std::iterator_traits<Iter>::value_type;
 
    // Numbers of bits per respective word.
    constexpr auto dst_bpw = BITS<dst_type>;
    constexpr auto src_bpw = BITS<src_type>;
 
    // Size of the source in words
    auto src_words = static_cast<std::size_t>(std::distance(src_begin, src_end));
 
    // The sizes in bits must match ...
    gf2_assert(store.size() == src_words * src_bpw, "Lengths do not match: {} != {}.", store.size(),
               src_words * src_bpw);
 
    // What we do next depends on the size of the respective words.
    if constexpr (src_bpw == dst_bpw) {
        // Same sized words -- we can copy whole words directly ...
        auto src_iter = src_begin;
        for (auto i = 0uz; i < store.words(); ++i, ++src_iter) store.set_word(i, static_cast<dst_type>(*src_iter));
    } else {
        // Otherwise proceed by copying one source word at a time into the appropriate span in the destination store.
        // The actual copying is done by the earlier method that copies any unsigned into the store,
        auto src_iter = src_begin;
        auto dst_begin = 0uz;
        for (auto i = 0uz; i < src_words; ++i, ++src_iter, dst_begin += src_bpw) {
            auto dst_span = span(store, dst_begin, dst_begin + src_bpw);
            copy(*src_iter, dst_span);
        }
    }
}

template<BitStore Src, BitStore Dst>
constexpr void
copy(Src const& src, Dst& dst) {
    // The sizes must match ...
    gf2_assert(dst.size() == src.size(), "Lengths do not match: {} != {}.", dst.size(), src.size());
 
    // What are the word types (without any const qualifiers)?
    using dst_type = std::remove_const_t<typename Dst::word_type>;
    using src_type = std::remove_const_t<typename Src::word_type>;
 
    // Numbers of bits per respective words
    constexpr auto dst_bpw = BITS<dst_type>;
    constexpr auto src_bpw = BITS<src_type>;
 
    // What we do next depends on the size of the respective words.
    if constexpr (dst_bpw == src_bpw) {
        // Same sized words holding the same number of bits => same number of words in each ...
        for (auto i = 0uz; i < dst.words(); ++i) {
            auto value = static_cast<dst_type>(src.word(i));
            dst.set_word(i, value);
        }
    } else if constexpr (dst_bpw > src_bpw) {
        // The destination word type is larger -- some number of source words fit into each destination word.
        // We assume that the destination word size is an integer multiple of the source word size.
        constexpr auto ratio = dst_bpw / src_bpw;
 
        auto src_words = src.words();
        for (auto i = 0uz; i < dst.words(); ++i) {
            // Combine `ratio` source words into a single destination word being careful not to run off the end of `src`
            auto value = ZERO<dst_type>;
            for (auto j = 0uz, src_index = i * ratio; j < ratio && src_index < src_words; ++j, ++src_index) {
                auto src_word = static_cast<dst_type>(src.word(src_index));
                value |= (src_word << (j * src_bpw));
            }
            dst.set_word(i, value);
        }
    } else {
        // The destination word type is smaller -- each source word becomes multiple destination words.
        // We assume that the source word size is an integer multiple of the destination word size.
        constexpr auto ratio = src_bpw / dst_bpw;
 
        auto dst_words = dst.words();
        auto dst_mask = MAX<dst_type>;
        for (auto i = 0uz; i < src.words(); ++i) {
            // Split each source word into `ratio` destination words being careful not to run off the end of `dst`
            auto src_word = src.word(i);
            for (auto j = 0uz, dst_index = i * ratio; j < ratio && dst_index < dst_words; ++j, ++dst_index) {
                dst.set_word(dst_index, static_cast<dst_type>(src_word & dst_mask));
                src_word >>= dst_bpw;
            }
        }
    }
}

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

template<BitStore Store>
constexpr void
copy(Store& store, std::invocable<usize> auto f) {
    set_all(store, false);
    for (auto i = 0uz; i < store.size(); ++i)
        if (f(i) == true) set(store, i);
}

template<BitStore Store>
constexpr void
fill_random(Store& store, 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) {
        set_all(store, false);
        return;
    }
 
    // 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) {
        set_all(store);
        return;
    }
 
    // 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(store, false);
    for (auto i = 0uz; i < store.size(); ++i)
        if (rng.u64() < scaled_p) set(store, i);
 
    // Restore the old seed if necessary.
    if (seed != 0) rng.set_seed(old_seed);
}


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

template<BitStore Store>
constexpr usize
count_zeros(Store const& store) {
    return store.size() - count_ones(store);
}

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

template<BitStore Store>
constexpr usize
trailing_zeros(Store const& store) {
    if (is_empty(store)) return 0;
 
    // The last word may not be fully occupied so we need to subtract the number of unused bits.
    auto bits_per_word = BITS<typename Store::word_type>;
    auto tail_bits = store.size() % bits_per_word;
    auto unused_bits = tail_bits == 0 ? 0 : bits_per_word - tail_bits;
    auto num_words = store.words();
    auto i = num_words;
    while (i--) {
        auto w = store.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 store.size();
}


template<BitStore Store>
constexpr std::optional<usize>
first_set(Store const& store) {
    // 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.
    auto bits_per_word = BITS<typename Store::word_type>;
    for (auto i = 0uz; i < store.words(); ++i) {
        if (auto loc = lowest_set_bit(store.word(i))) return i * bits_per_word + loc.value();
    }
    return {};
}

template<BitStore Store>
constexpr std::optional<usize>
last_set(Store const& store) {
    // 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.
    auto bits_per_word = BITS<typename Store::word_type>;
    for (auto i = store.words(); i--;) {
        if (auto loc = highest_set_bit(store.word(i))) return i * bits_per_word + loc.value();
    }
    return {};
}

template<BitStore Store>
constexpr std::optional<usize>
next_set(Store const& store, usize index) {
    // 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 >= store.size()) return {};
 
    // Where is that starting index located in the word store?
    using word_type = typename Store::word_type;
    auto bits_per_word = BITS<word_type>;
    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 < store.words(); ++i) {
        auto w = store.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 {};
}

template<BitStore Store>
constexpr std::optional<usize>
previous_set(Store const& store, usize index) {
    // Edge case: If the store is empty or we are at the start, there are no previous set bits.
    if (is_empty(store) || index == 0) return {};
 
    // Silently fix large indices and also adjust the index down a slot.
    if (index > store.size()) index = store.size();
    index--;
 
    // Where is that starting index located in the word store?
    using word_type = typename Store::word_type;
    auto bits_per_word = BITS<word_type>;
    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 = store.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 {};
}


template<BitStore Store>
constexpr std::optional<usize>
first_unset(Store const& store) {
    // Iterate forward looking for a word with a unset bit and use the lowest of those ...
    using word_type = typename Store::word_type;
    auto bits_per_word = BITS<word_type>;
    for (auto i = 0uz; i < store.words(); ++i) {
        auto w = store.word(i);
        if (i == store.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>(store.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 {};
}

template<BitStore Store>
constexpr std::optional<usize>
last_unset(Store const& store) {
    // Iterate backwards looking for a word with a unset bit and use the highest of those ...
    using word_type = typename Store::word_type;
    auto bits_per_word = BITS<word_type>;
    for (auto i = store.words(); i--;) {
        auto w = store.word(i);
        if (i == store.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>(store.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 {};
}

template<BitStore Store>
constexpr std::optional<usize>
next_unset(Store const& store, usize index) {
    // 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 >= store.size()) return {};
 
    // Where is that starting index located in the word store?
    using word_type = typename Store::word_type;
    auto bits_per_word = BITS<word_type>;
    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 < store.words(); ++i) {
        auto w = store.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 == store.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>(store.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 {};
}

template<BitStore Store>
constexpr std::optional<usize>
previous_unset(Store const& store, usize index) {
    // Edge case: If the store is empty or we are at the start, there are no previous set bits.
    if (is_empty(store) || index == 0) return {};
 
    // Silently fix large indices and also adjust the index down a slot.
    if (index > store.size()) index = store.size();
    index--;
 
    // Where is that starting index located in the word store?
    using word_type = typename Store::word_type;
    auto bits_per_word = BITS<word_type>;
    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 = store.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 == store.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>(store.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 {};
}


template<BitStore Store>
constexpr auto
bits(Store const& store) {
    return Bits<const Store, true>(&store, 0);
}

template<BitStore Store>
constexpr auto
bits(Store& store) {
    return Bits<Store, false>(&store, 0);
}

template<BitStore Store>
constexpr auto
set_bits(Store const& store) {
    return SetBits<Store>(&store);
}

template<BitStore Store>
constexpr auto
unset_bits(Store const& store) {
    return UnsetBits<Store>(&store);
}

template<BitStore Store>
constexpr auto
store_words(Store const& store) {
    return Words<Store>(&store);
}

template<BitStore Store>
constexpr void
to_words(Store const& store, std::output_iterator<typename Store::word_type> auto out) {
    for (auto i = 0uz; i < store.words(); ++i, ++out) *out = store.word(i);
}

template<BitStore Store>
constexpr auto
to_words(Store const& store) {
    return std::ranges::to<std::vector>(store_words(store));
}


template<BitStore Store>
static constexpr auto
span(Store const& store, usize begin, usize end) {
    // Check that the range is correctly formed & does not extend beyond the end of the bit-store.
    gf2_assert(begin <= end, "Span range [{}, {}) is mis-ordered.", begin, end);
    gf2_assert(end <= store.size(), "Span end {} extends beyond the store end {}.", end, store.size());
 
    // Get the zero-based word index and bit offset in that word for the begin location.
    using word_type = typename Store::word_type;
    auto bits_per_word = BITS<word_type>;
    auto [data_index, bit_offset] = index_and_offset<word_type>(begin);
 
    // If the store is already a span, we need to adjust those values.
    bit_offset += store.offset();
    if (bit_offset >= bits_per_word) {
        data_index += 1;
        bit_offset %= bits_per_word;
    }
 
    // Return the read-only bit-span over the underlying word data.
    return BitSpan<const word_type>(store.store() + data_index, bit_offset, end - begin);
}

template<BitStore Store>
static constexpr auto
span(Store& store, usize begin, usize end) {
    // Check that the range is correctly formed & does not extend beyond the end of the bit-store.
    gf2_assert(begin <= end, "Span range [{}, {}) is mis-ordered.", begin, end);
    gf2_assert(end <= store.size(), "Span end {} extends beyond the store end {}.", end, store.size());
 
    // Get the zero-based word index and bit offset in that word for the begin location.
    using word_type = typename Store::word_type;
    auto bits_per_word = BITS<word_type>;
    auto [data_index, bit_offset] = index_and_offset<word_type>(begin);
 
    // If the store is already a span, we need to adjust those values.
    bit_offset += store.offset();
    if (bit_offset >= bits_per_word) {
        data_index += 1;
        bit_offset %= bits_per_word;
    }
 
    // Return the read-only bit-span over the underlying word data.
    return BitSpan<word_type>(store.store() + data_index, bit_offset, end - begin);
}


template<BitStore Store>
constexpr auto
sub(Store const& store, usize begin, usize end) {
    return BitVector<typename Store::word_type>::from(gf2::span(store, begin, end));
}


template<BitStore Store>
constexpr void
split(Store const& store, usize at, BitVector<typename Store::word_type>& left,
      BitVector<typename Store::word_type>& right) {
    auto sz = store.size();
    gf2_assert(at <= sz, "Oops, split point {} is beyond the end of the bit-store {}", at, sz);
    left.resize(at);
    right.resize(sz - at);
    copy(span(store, 0, at), left);
    copy(span(store, at, sz), right);
}

template<BitStore Store>
constexpr auto
split(Store const& store, usize at) {
    BitVector<typename Store::word_type> left, right;
    split(store, at, left, right);
    return std::pair{left, right};
}

template<BitStore Lhs, BitStore Rhs>
auto
join(Lhs const& lhs, Rhs const& rhs) {
    auto lhs_size = lhs.size();
    auto rhs_size = rhs.size();
 
    using word_type = typename Lhs::word_type;
    BitVector<word_type> result(lhs_size + rhs_size);
 
    result.span(0, lhs_size).copy(lhs);
    result.span(lhs_size, lhs_size + rhs_size).copy(rhs);
    return result;
}


template<BitStore Store>
constexpr void
riffle(Store const& store, BitVector<typename Store::word_type>& dst) {
    // Nothing to do if the store is empty or has just one element. With two bits `ab` we fill `dst` with `a0b`.
    auto sz = store.size();
    if (sz < 2) {
        dst.resize(sz);
        dst.copy(store);
        return;
    }
 
    // Make sure `dst` is large enough to hold the riffled bits (a bit too big but we fix that below).
    dst.resize(2 * sz);
    auto dst_words = dst.words();
 
    // Riffle each word in the store into two adjacent words in `dst`.
    for (auto i = 0uz; i < store.words(); ++i) {
        auto [lo, hi] = riffle(store.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();
}

template<BitStore Store>
constexpr auto
riffle(Store const& store) {
    BitVector<typename Store::word_type> result;
    riffle(store, result);
    return result;
}


template<BitStore Store>
static std::string
to_binary_string(Store const& store, std::string_view sep = "", std::string_view pre = "", std::string_view post = "") {
    // Edge case ...
    if (is_empty(store)) { 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 = store.words();
    std::string binary_string;
    binary_string.reserve(num_words * BITS<typename Store::word_type>);
    for (auto i = 0uz; i < num_words; ++i) {
        auto w = reverse_bits(store.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(store.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() + store.size() * (sep.size() + 1) + post.size());
    result.append(pre);
    for (auto i = 0uz; i < store.size(); ++i) {
        result.push_back(binary_string[i]);
        if (i != store.size() - 1) { result.append(sep); }
    }
    result.append(post);
    return result;
}

template<BitStore Store>
static std::string
to_string(Store const& store, std::string_view sep = "", std::string_view pre = "", std::string_view post = "") {
    return to_binary_string(store, sep, pre, post);
}

template<BitStore Store>
static std::string
to_pretty_string(Store const& store) {
    return to_binary_string(store, ",", "[", "]");
}

template<BitStore Store>
static std::string
to_hex_string(Store const& store) {
    // Edge case: No bits in the store, return the empty string.
    if (is_empty(store)) return "";
 
    // The number of digits in the output string. Generally hexadecimal but the last may be to a lower base.
    auto digits = (store.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.
    using word_type = typename Store::word_type;
    for (auto i = 0uz; i < store.words(); ++i) {
        auto w = reverse_bits(store.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 = store.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(store, store.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;
}

template<BitStore Store>
static std::string
describe(Store const& store) {
    using word_type = typename Store::word_type;
    auto        bits_per_word = BITS<word_type>;
    std::string result;
    result.reserve(1000);
    result += std::format("binary format:        {}\n", to_binary_string(store));
    result += std::format("hex format:           {}\n", to_hex_string(store));
    result += std::format("number of bits:       {}\n", store.size());
    result += std::format("number of set bits:   {}\n", count_ones(store));
    result += std::format("number of unset bits: {}\n", count_zeros(store));
    result += std::format("bits per word:        {}\n", bits_per_word);
    result += std::format("word count:           {}\n", store.words());
    result += "words in hex:         [";
    for (auto i = 0uz; i < store.words(); ++i) {
        auto w = store.word(i);
        result += std::format("{:0{}X}", w, (bits_per_word + 3) / 4);
        result += std::format("{:x}", store.word(i));
        if (i + 1 < store.words()) result += ", ";
    }
    result += "]\n";
    return result;
}

template<BitStore Store>
inline std::ostream&
operator<<(Store const& store, std::ostream& s) {
    return s << to_string(store);
}


// # Example
template<BitStore Lhs, BitStore Rhs>
constexpr bool
operator==(Lhs const& lhs, Rhs const& rhs) {
    if constexpr (!std::same_as<typename Lhs::word_type, typename Rhs::word_type>) 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;
}


template<BitStore Store>
constexpr void
operator<<=(Store& store, usize shift) {
    // Edge cases where there is nothing to do.
    if (shift == 0 || store.size() == 0) return;
 
    // Perhaps we are shifting all the bits out and are left with all zeros.
    if (shift >= store.size()) {
        set_all(store, false);
        return;
    }
 
    // For larger shifts, we can efficiently shift by whole words first.
    using word_type = typename Store::word_type;
    auto  bits_per_word = BITS<word_type>;
    usize word_shift = shift / bits_per_word;
    usize end_word = store.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) store.set_word(i, store.word(i + word_shift));
 
        // Fill in the high order words with zeros.
        for (auto i = end_word; i < store.words(); ++i) store.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>(store.word(i) >> shift);
                auto hi = static_cast<word_type>(store.word(i + 1) << shift_complement);
                store.set_word(i, lo | hi);
            }
        }
 
        // Do the last word.
        auto value = store.word(end_word - 1);
        store.set_word(end_word - 1, static_cast<word_type>(value >> shift));
    }
}

template<BitStore Store>
constexpr void
operator>>=(Store& store, usize shift) {
    // Edge cases where there is nothing to do.
    if (shift == 0 || store.size() == 0) return;
 
    // Perhaps we are shifting all the bits out and are left with all zeros.
    if (shift >= store.size()) {
        set_all(store, false);
        return;
    }
 
    // For larger shifts, we can efficiently shift by whole words first.
    using word_type = typename Store::word_type;
    auto  bits_per_word = BITS<word_type>;
    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 = store.words(); i-- > word_shift;) store.set_word(i, store.word(i - word_shift));
 
        // Fill in the low order words with zeros.
        for (auto i = 0uz; i < word_shift; ++i) store.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 = store.words(); i-- > word_shift + 1;) {
            auto lo = static_cast<word_type>(store.word(i - 1) >> shift_complement);
            auto hi = static_cast<word_type>(store.word(i) << shift);
            store.set_word(i, lo | hi);
        }
    }
 
    // Do the "first" word.
    auto value = store.word(word_shift);
    store.set_word(word_shift, static_cast<word_type>(value << shift));
}


template<BitStore Store>
constexpr auto
operator<<(Store const& store, usize shift) {
    auto result = BitVector<typename Store::word_type>::from(store);
    result <<= shift;
    return result;
}

template<BitStore Store>
constexpr auto
operator>>(Store const& store, usize shift) {
    auto result = BitVector<typename Store::word_type>::from(store);
    result >>= shift;
    return result;
}


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

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

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


template<BitStore Store>
constexpr auto
operator~(Store const& store) {
    auto result = BitVector<typename Store::word_type>::from(store);
    result.flip_all();
    return result;
}

template<BitStore Lhs, BitStore Rhs>
    requires std::same_as<typename Lhs::word_type, typename Rhs::word_type>
constexpr auto
operator^(Lhs const& lhs, Rhs const& rhs) {
    auto result = BitVector<typename Lhs::word_type>::from(lhs);
    result ^= rhs;
    return result;
}

template<BitStore Lhs, BitStore Rhs>
    requires std::same_as<typename Lhs::word_type, typename Rhs::word_type>
constexpr auto
operator&(Lhs const& lhs, Rhs const& rhs) {
    auto result = BitVector<typename Lhs::word_type>::from(lhs);
    result &= rhs;
    return result;
}

template<BitStore Lhs, BitStore Rhs>
    requires std::same_as<typename Lhs::word_type, typename Rhs::word_type>
constexpr auto
operator|(Lhs const& lhs, Rhs const& rhs) {
    auto result = BitVector<typename Lhs::word_type>::from(lhs);
    result |= rhs;
    return result;
}


template<BitStore Lhs, BitStore Rhs>
    requires std::same_as<typename Lhs::word_type, typename Rhs::word_type>
constexpr void
operator+=(Lhs& lhs, Rhs const& rhs) {
    lhs ^= rhs;
}

template<BitStore Lhs, BitStore Rhs>
    requires std::same_as<typename Lhs::word_type, typename Rhs::word_type>
constexpr void
operator-=(Lhs& lhs, Rhs const& rhs) {
    lhs ^= rhs;
}


template<BitStore Lhs, BitStore Rhs>
    requires std::same_as<typename Lhs::word_type, typename Rhs::word_type>
constexpr auto
operator+(Lhs const& lhs, Rhs const& rhs) {
    auto result = BitVector<typename Lhs::word_type>::from(lhs);
    result += rhs;
    return result;
}

template<BitStore Lhs, BitStore Rhs>
    requires std::same_as<typename Lhs::word_type, typename Rhs::word_type>
constexpr auto
operator-(Lhs const& lhs, Rhs const& rhs) {
    auto result = BitVector<typename Lhs::word_type>::from(lhs);
    result -= rhs;
    return result;
}


template<BitStore Lhs, BitStore Rhs>
    requires std::same_as<typename Lhs::word_type, typename Rhs::word_type>
constexpr bool
dot(Lhs const& lhs, Rhs const& rhs) {
    gf2_debug_assert_eq(lhs.size(), rhs.size(), "Length mismatch {} != {}", lhs.size(), rhs.size());
    using word_type = typename 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<BitStore Lhs, BitStore Rhs>
    requires std::same_as<typename Lhs::word_type, typename Rhs::word_type>
constexpr bool
operator*(Lhs const& lhs, Rhs const& rhs) {
    return dot(lhs, rhs);
}


template<BitStore Lhs, BitStore Rhs>
    requires std::same_as<typename Lhs::word_type, typename Rhs::word_type>
auto
convolve(Lhs const& lhs, Rhs const& rhs) {
    using word_type = typename Lhs::word_type;
    auto bits_per_word = BITS<word_type>;
 
    // Edge case: if either store is empty then the convolution is empty.
    if (lhs.is_empty() || rhs.is_empty()) return BitVector<word_type>{};
 
    // Set up the result bit-vector.
    auto result = BitVector<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_per_word - 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

template<gf2::BitStore Store>
struct std::formatter<Store> {
    // Parse a bit-store format specifier where where we recognize {:p} and {:x}
    constexpr auto parse(std::format_parse_context const& 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;
    }
 
    // Format the bit-store according using the appropriate string conversion function.
    template<class FormatContext>
    auto format(Store const& 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(), "{}", to_hex_string(rhs));
        if (m_pretty) return std::format_to(ctx.out(), "{}", to_pretty_string(rhs));
 
        // Default
        return std::format_to(ctx.out(), "{}", to_string(rhs));
    }
 
    bool m_hex = false;
    bool m_pretty = false;
    bool m_error = false;
};