Twelve posts, twelve codes, one thesis that refused to change. This is the closing summary.

A. The Twelve Codes Together

Every post in this series answered a version of the same question: given a source of positive integers, how do you represent its values compactly as a sequence of bits? The answers differ in shape, in assumptions, and in which distribution each code implicitly expects.

Posts 1 and 2 (Kraft’s Inequality and McMillan’s Converse) established why prefix-free codes are the right unit of analysis. Post 3 (Universal Codes as Priors) named the frame: a code is a hypothesis about the source. Posts 4 through 10 filled in the catalogue. Posts 11 and 12 extended from integer coding to set representation, where the questions shift from “how long is this codeword?” to “how do you store membership?” and “how do you answer rank/select?”

Looking across all twelve, the main lesson is not that one code dominates. It is that the question “which code?” is always empirically answerable given a sample.

B. The Unifying Frame Restated

Post 3 introduced the codes-as-priors thesis with two instances behind it. We now have twelve. The thesis has not changed; it has only become more evidently true.

Every type decomposes structurally. So does its codec.

Codecs as Functors

You have an optional<vector<pair<int, string>>>. The type decomposes structurally: it is an optional of a free monoid of products of an integer and a string. That decomposition is not an observation about memory layout. It is a statement about the algebraic structure of the type.

Now ask: does the codec decompose the same way?

If the answer is yes, you stop writing one-off encoders. You build a codec for optional<T> from a codec for T. You build a codec for vector<T> from a codec for T. The codec for optional<vector<pair<int, string>>> assembles from its parts with no manual layout decisions, no hand-placed length headers, no ad-hoc format negotiation.

This post argues that the answer is always yes, and shows what the machinery looks like. The thesis: codecs are not ad-hoc bit formats. They are constructions on the algebraic structure of types. The algebraic structure of a type determines its codec, the same way it determines its algorithms.

This extends Stepanov’s claim. The peasant algorithm post showed that algorithms arise from algebraic structure. The homomorphism post showed that structure-preserving maps are the natural morphisms. Here, we show the codec itself is a structure-preserving map, and that it lifts from leaf types to compound types by the same algebraic logic.

Bit I/O: The Foundation

Before combinators, we need concrete bit I/O. The approach taken here follows Stepanov’s move in the algorithm posts: state the concept first, then provide a model.

Two concepts govern bit-level I/O:

template<typename T>
concept BitSink = requires(T& s, bool bit) {
    { s.write(bit) } -> std::same_as<void>;
};

template<typename T>
concept BitSource = requires(T& s) {
    { s.read() } -> std::same_as<bool>;
    { s.peek() } -> std::convertible_to<bool>;
};

A BitSink accepts bits. A BitSource supplies them. A codec is an algorithm parameterized over BitSink and BitSource, not a class hierarchy. This is Stepanov’s move at the bit level: require only what the algorithm needs, let anything that satisfies the concept participate.

The standard models are BitWriter and BitReader, which pack bits into byte buffers in LSB-first order:

class BitWriter {
    std::span<std::uint8_t> buf_;
    std::size_t byte_idx_ = 0;
    std::uint8_t byte_ = 0;
    std::uint8_t bit_pos_ = 0;
public:
    explicit BitWriter(std::span<std::uint8_t> buf) noexcept : buf_(buf) {}

    void write(bool bit) noexcept {
        byte_ |= (bit ? std::uint8_t{1} : std::uint8_t{0}) << bit_pos_;
        if (++bit_pos_ == 8) {
            buf_[byte_idx_++] = byte_;
            byte_ = 0;
            bit_pos_ = 0;
        }
    }

    void align() noexcept {
        if (bit_pos_ > 0) {
            buf_[byte_idx_++] = byte_;
            byte_ = 0;
            bit_pos_ = 0;
        }
    }

    [[nodiscard]] std::size_t bytes_written() const noexcept {
        return byte_idx_ + (bit_pos_ > 0 ? 1 : 0);
    }
};

class BitReader {
    std::span<const std::uint8_t> buf_;
    std::size_t byte_idx_ = 0;
    std::uint8_t bit_pos_ = 0;
public:
    explicit BitReader(std::span<const std::uint8_t> buf) noexcept : buf_(buf) {}

    bool read() noexcept {
        bool bit = ((buf_[byte_idx_] >> bit_pos_) & 1) != 0;
        if (++bit_pos_ == 8) {
            ++byte_idx_;
            bit_pos_ = 0;
        }
        return bit;
    }

    [[nodiscard]] bool peek() const noexcept {
        return byte_idx_ < buf_.size();
    }
};

A codec concept rounds out the three-concept core:

The free monoid on a type lifts to bit space. It lifts injectively only when the element codec is prefix-free.

Prefix-Free Codes and the Free Monoid

You have a list of unsigned integers. Encode the list as a single bit string.

Fixed-width encoding wastes space. If you allocate 64 bits per integer, small values like 1 or 7 cost as much as values near \(2^{64}\). Variable-width encoding recovers that space, but immediately raises a harder question: where does one encoded integer end and the next begin?

Two escape routes. First, prefix each encoded item with its length. That works, but the length headers are overhead, and you now need a codec for the lengths as well. Second, choose a code where the structure of the codewords makes boundaries unambiguous without any headers. These are prefix-free codes, and this is the right answer, in a precise categorical sense.

The “precise categorical sense” is what this post develops. Encoding a list as the concatenation of encoded elements is a monoid homomorphism from the free monoid on \(T\) to the monoid of bit strings under concatenation. The universal property of the free monoid guarantees this homomorphism always exists. The question of whether the decoder can invert it comes down to exactly one property of the element codec: whether it is prefix-free.

The Free Monoid, Recalled

A monoid is a set with an associative binary operation and an identity element. The free monoid on a set \(S\) is the set of all finite sequences of elements from \(S\), with concatenation as the operation and the empty sequence as the identity.

“Free” means no equations hold except those forced by the monoid axioms. Nothing is identified with anything else. If you need commutativity or idempotency, you quotient the free monoid by additional equations. But the free monoid itself imposes nothing beyond associativity and identity.

The universal property says: given any monoid \(M\) and any function \(f: S \to M\), there is exactly one monoid homomorphism \(\hat{f}: \text{Free}(S) \to M\) that extends \(f\). That unique extension is fold:

where \(\cdot\) is the operation in \(M\). The free-algebra post develops this in full. For this post, the one fact that matters is that fold is canonical: it is the unique way to extend a per-element map to a list-consuming function that respects the monoid structure.

Huffman codes one symbol at a time. Arithmetic coding encodes the whole sequence as a single number. The difference is a factor of twelve, at least on the right source.

The Last Bit of Redundancy

Huffman coding gets expected codeword length within one bit of entropy. That is the best it can do, because codeword lengths must be integers while entropy is a real number.

The waste is structural. A symbol with probability $p = 0.7$ has optimal (fractional) length $-\log_2(0.7) \approx 0.515$ bits. Huffman rounds that up to 1 bit: 0.485 bits wasted per occurrence. For a nearly-deterministic source with $p_0 = 0.99$ and $p_1 = 0.01$, the entropy is $H \approx 0.081$ bits per symbol. Huffman is stuck at 1 bit per symbol. That is a factor-of-twelve gap, and Huffman cannot close it: a symbol that appears 99% of the time still gets a complete codeword.

Arithmetic coding steps back from per-symbol codewords entirely. It encodes an entire sequence as a single rational number in $[0, 1)$. The bit-length of that number converges to the entropy of the sequence as the sequence grows. No integer rounding, no per-symbol overhead.

This post builds an integer range coder in C++23 and demonstrates the factor-of-twelve improvement on the Bernoulli(0.99) source.

The Continuous View

Start with the unit interval $[0, 1)$. For a two-symbol source, partition it by probability: symbol 0 gets $[0, p_0)$ and symbol 1 gets $[p_0, 1)$.

To encode a sequence, begin with the full interval and narrow it with each symbol. After symbol $s_1$, restrict to the corresponding sub-interval. After $s_2$, apply the same proportional rule inside that sub-interval. After $L$ symbols, the interval has width $\prod_{i=1}^{L} p_{s_i}$.

Any number inside the final interval is a valid encoding. The shortest such number in binary requires approximately $-\log_2(\prod p_{s_i}) = \sum_{i=1}^{L} (-\log_2 p_{s_i})$ bits. As $L \to \infty$, bits per symbol approaches $H(p) = -\sum_k p_k \log_2 p_k$ exactly.

Decoding is the inverse: given the encoded number, determine at each step which sub-interval it falls in, recover the symbol, narrow the interval, and repeat.

The theory is complete. The practice is not. A real interval narrows exponentially fast: after a few dozen symbols you need arbitrary precision. The integer range coder fixes this with 32-bit arithmetic and a renormalization step.

Huffman coding is two things: the optimal length vector for a known distribution, and McMillan’s construction applied to that vector. This post develops both.

From Universal to Optimal

Every code in this series so far has been universal: no prior knowledge of the source distribution required. Elias gamma assigns shorter codewords to smaller integers regardless of which integers actually appear. Fibonacci does the same. VByte packs smaller values into fewer bytes without knowing whether your data clusters at the low end or the high end. Universal codes are defensive: they perform acceptably across a broad class of inputs by committing to none.

Huffman flips that stance. You bring a finite probability distribution. Huffman finds the prefix-free code with minimum expected codeword length for that distribution. The code is tuned to what you provide and will perform poorly on anything else. Call this the move from defensive to distribution-specific coding.

The payoff is real. Shannon’s source coding theorem says no prefix-free code can achieve expected length below $H(p) = -\sum_i p_i \log_2 p_i$. Huffman gets within 1 bit of that bound. For any prefix-free code, expected length satisfies

The upper bound comes from the integer-length constraint: each codeword is a whole number of bits, and $\lceil -\log_2 p_i \rceil \le -\log_2 p_i + 1$. Huffman is optimal subject to this constraint. No prefix-free code does better without abandoning whole-bit codewords.

That last clause points to the limit of this post and the subject of the next. Arithmetic coding breaks the integer-length constraint by assigning fractional bits in effect, reaching entropy exactly in the limit.

The Algorithm

The four steps of Huffman’s construction:

1. Create one leaf node per symbol, weighted by its probability.
2. Push all leaves into a min-priority queue (lowest weight first).
3. While the queue contains more than one node: extract the two lowest-weight nodes, merge them into a new internal node whose weight is their sum, push the internal node back.
4. The remaining node is the root. The path from root to each leaf encodes that leaf’s codeword (“0” for left, “1” for right).

Here is the complete implementation from huffman.hpp:

PFC (Prefix-Free Codecs) is a header-only C++20 library built on a simple observation: data compression and zero-copy access are not contradictory goals, as long as you build on prefix-free codes and generic programming. The library gets 3-10x compression on typical integer distributions while maintaining full STL integration and type safety.

The zero-copy invariant

Traditional compression creates two worlds. Data lives uncompressed in memory (32 bits per integer) and compressed on disk (variable bits). You marshal between them. PFC eliminates that boundary:

// Traditional approach
std::vector<uint32_t> data = {1, 2, 3, 5, 8, 13};
auto compressed = compress(data);     // Marshal to wire format
store_to_disk(compressed);
auto restored = decompress(compressed); // Unmarshal back

// PFC approach
PackedContainer<uint32_t, EliasGamma> data;
data.push_back(1);
data.push_back(2);
data.push_back(3);

// Data is ALREADY compressed in memory
uint32_t value = data[0];  // Decodes from compressed form on access

// Write to disk? Zero copy.
write(fd, data.bytes().data(), data.bytes().size());

The data structure IS the compressed format. There’s no marshaling step.

Prefix-free codes

A code is prefix-free if no codeword is a prefix of another:

Prefix-free:           Not prefix-free:
  1 -> 0                 1 -> 0
  2 -> 10                2 -> 01    ("0" is a prefix of "01")
  3 -> 110               3 -> 011
  4 -> 1110

This matters because prefix-free codes compose naturally. Concatenate encodings without delimiters and decode unambiguously:

encode(1);  // 0
encode(2);  // 10
encode(3);  // 110
// Result: 010110 (self-delimiting)

decode();  // Reads "0" -> 1
decode();  // Reads "10" -> 2
decode();  // Reads "110" -> 3

This enables streaming and composition without any framing overhead.

Universal codes

PFC implements several universal codes, meaning they’re asymptotically optimal for any distribution without knowing the distribution in advance.

Elias Gamma encodes positive integer \(n\) in \(2\lfloor\log_2 n\rfloor + 1\) bits:

n     Binary   Elias Gamma
1     1        1
2     10       010
3     11       011
4     100      00100
5     101      00101
8     1000     0001000

Write \(\lfloor\log_2 n\rfloor\) zeros, then the binary representation of \(n\). Asymptotically optimal for geometric distributions.

Fibonacci encoding uses Zeckendorf representation (sum of non-consecutive Fibonacci numbers) with a terminal “11” marker. Every positive integer has a unique representation.

Rice/Golomb codes are parametric, optimal for geometric distributions with known parameter. Quotient in unary, remainder in binary. Good for values with exponential decay.

Every code in this series so far operates at bit granularity. VByte does not. It gives up bit-level precision for byte-alignment, and in production systems, that trade wins most of the time.

The Practical Question

Every code in this series so far operates at bit granularity. Elias gamma encodes 1 in a single bit. Fibonacci coding uses exactly as many bits as the Zeckendorf representation requires. Bit packing is theoretically attractive because it minimizes the number of bits written, which minimizes the encoded size.

But bit packing is computationally expensive. Reading or writing a single bit requires a shift, a mask, and often a branch to handle byte boundaries. Encoding a sequence of integers this way burns CPU cycles that scale with the number of integers, independent of their values. For high-throughput applications, the overhead of bit manipulation can easily exceed the savings from compact encoding.

VByte (also called Varint in Google’s ecosystem, and LEB128 in the DWARF debug format) trades a small amount of length efficiency for byte-alignment. The idea is simple: encode each integer as a sequence of 7-bit groups, one per byte, with a continuation flag in the high bit of each byte. The result is self-delimiting, compact for small values, and requires no bit-level manipulation to decode.

VByte is the encoding used by Protocol Buffers for all integer fields. It appears in Apache Arrow, Parquet, Snappy’s block format, LevelDB’s metadata, and most production columnar file formats. These are high-throughput systems. Byte-alignment is why VByte is their choice over the more compact universal codes from posts 4 through 7.

The Encoding

VByte splits an integer into 7-bit groups, starting from the least significant bits. Each group occupies one byte where bits 0 through 6 carry 7 data bits and bit 7 is a continuation flag: 1 means more bytes follow, 0 means this is the last byte.

A value in $[0, 127]$ fits in a single byte (continuation flag clear). A value in $[128, 16383]$ requires two bytes (first byte has flag set, second has flag clear). The pattern continues: each additional byte adds 7 bits of capacity.

Every code in this series so far has been fixed. Rice and Golomb are different: they take a parameter, and the parameter is your model of the data.

The First Parametric Code

Every code examined so far in this series has been monolithic. Unary coding is just unary coding. Elias gamma is just Elias gamma. Each one encodes all non-negative integers with a single fixed strategy. You do not get to choose anything about the code beyond whether to use it.

Rice and Golomb codes break this pattern. They are parametric: a single integer parameter, $k$ for Rice or $m$ for Golomb, tunes the code to a specific source distribution. Rice$(k)$ is not one code but a family of codes, one per value of $k$. Each member of the family is optimal for a specific geometric distribution. Choosing $k$ is choosing your prior precisely.

This matters because data sources are rarely uniform. Run-length encodings, inter-frame video differences, and the gap sequences in inverted indexes are all approximately geometrically distributed. If you know the mean of your source, you can pick $k$ so that Rice$(k)$ performs near-optimally, without the overhead of a Huffman table or arithmetic coding.

The key insight: for a geometric source with mean approximately $2^k$, Rice$(k)$ is within a small constant of entropy. No other universal code in this series achieves this. Elias gamma and delta perform well asymptotically but can be far from optimal for a specific geometric distribution with a known mean. Rice exploits that knowledge directly.

Rice Coding

Rice coding splits a non-negative integer $n$ into two parts: a quotient $q = \lfloor n / 2^k \rfloor = n \gg k$ and a remainder $r = n \bmod 2^k = n \mathbin{&} (2^k - 1)$.

The quotient is encoded in unary: $q$ zero bits followed by a stop bit of 1. The remainder is encoded in exactly $k$ bits, MSB first. The total codeword is the concatenation of these two parts.

Codeword examples for $k = 2$ (remainder is always 2 bits):

Codeword length: $(n \gg k) + 1 + k$ bits. The Kraft sum saturates to 1, so Rice is a complete prefix-free code.

Every code in this series so far has optimized expected length under some implied prior. Fibonacci coding does something different: it gives the decoder a way to recover from errors without help from a lower layer.

A Different Design Goal

All the codes in this series have aimed at the same target: assign short codewords to frequent symbols, with length growing roughly as $\log n$ for the $n$-th symbol under some implied prior. Elias gamma minimizes expected length for power-law distributions; delta and omega extend the recursion for heavier tails.

Fibonacci coding has a different goal. It does not optimize for average codeword length under a specific distribution. It optimizes for error resilience. In a stream of gamma-coded integers, a single bit flip in a codeword’s length prefix causes the decoder to misread that codeword’s length, then misread every subsequent codeword. The error propagates without limit until the decoder somehow reacquires sync. On a reliable channel this is a nonissue. On a noisy one, or in stored data that may have silently rotted, it is a serious problem.

Fibonacci coding avoids this. Every Fibonacci codeword ends in two consecutive 1 bits (“11”). This double-one marker appears nowhere else in the codeword. A single bit flip corrupts the codeword it hits, possibly spills into the next codeword, and then the decoder finds the next “11” and resynchronizes. At most two codewords are corrupted per error. The rest of the stream is intact.

The price is length overhead: Fibonacci codewords are approximately $1.44 \times \log_2 n$ bits long, compared to $\log_2 n$ bits for the entropy lower bound. On a reliable channel, that overhead is not worth paying. On a noisy channel, or in a long-running stream where rare bit errors must not lose the entire tail, the self-synchronization property is worth it.

Zeckendorf’s Theorem

Fibonacci numbers starting from $F_2 = 1$: $1, 2, 3, 5, 8, 13, 21, 34, \ldots$

Zeckendorf’s theorem: every positive integer $n$ has a unique representation as a sum of non-consecutive Fibonacci numbers. The greedy algorithm produces it by repeatedly subtracting the largest Fibonacci number that does not exceed $n$.

Elias gamma spends too many bits saying how many bits it will use. Delta fixes that. Omega takes the fix one step further. This post is about what happens when you apply recursion to the length prefix.

Where Gamma Stops Being Good

Elias gamma, from the previous post, encodes a positive integer $n$ in $2\lfloor \log_2 n \rfloor + 1$ bits: a unary count of $\lfloor \log_2 n \rfloor$ zeros, then a stop bit, then the $\lfloor \log_2 n \rfloor$ trailing binary bits of $n$. For small $n$ this is fine. For large $n$, nearly half the bits are spent on the unary prefix alone.

The unary prefix is the bottleneck. It encodes the length $L = \lfloor \log_2 n \rfloor + 1$ in the most wasteful possible way: one bit per unit. For $n = 256$, that is 8 zero bits just to say “the payload is 8 bits long.” The payload itself is also 8 bits, so you are paying a 100% overhead on the length announcement. That is bad, and it gets worse as $n$ grows.

The fix is obvious once you see it: encode $L$ itself in some shorter code instead of unary. Elias delta does exactly this, replacing the unary length prefix with a gamma-coded length. Elias omega takes the idea one step further and applies the recursion to itself, all the way down.

Both codes are universal: they assign finite codewords to every positive integer, and the expected codeword length is within a constant factor of optimal for any source whose probabilities decrease with $n$. The improvement over gamma is real and measurable once $n$ grows past a few dozen.

This post shows both implementations, their implied priors, and the crossover points where each code wins. As in the rest of this series, the code is pedagogical: each header stands alone and the struct-with-encode/decode pattern maps directly onto the PFC library’s EliasDelta and EliasOmega in codecs.hpp.

Elias Delta

Algorithm. Let $L = \lfloor \log_2 n \rfloor + 1$ (the bit-width of $n$, equivalently std::bit_width(n)).

1. Encode $L$ in Elias gamma.
2. Write the $L - 1$ trailing bits of $n$ after its implicit leading 1, MSB first.

Gamma encodes $L$ (a small integer) in $O(\log \log n)$ bits instead of $O(\log n)$ bits for the unary prefix. The payload is identical to gamma’s: the trailing bits of $n$. The total length is $O(\log n + \log \log n)$.

Unary is older than information theory. Elias gamma is its 1975 improvement. Together they span the gap between optimal-but-impractical and practical-but-nearly-optimal. This post derives what each code bets on, and shows numerically what that means.

Unary and Elias Gamma

Unary is the oldest code in this series. It predates information theory by centuries: a shepherd counting sheep on a stick is using unary. Mark one notch per sheep; count the notches to decode. The codeword for \(n\) is \(n\) tally marks. Its information-theoretic justification came later, when Shannon showed it is exactly optimal for a geometric source.

Elias gamma is the 1975 extension by Peter Elias. It brings the codeword length from \(O(n)\) to \(O(\log n)\), making it practical for numbers beyond small single digits, while keeping the prefix-free property that makes self-delimiting streams possible.

Both codes are instances of the claim from Universal Codes as Priors: every prefix-free code is a bet about the source. Unary bets on a geometric distribution with parameter \(1/2\). Gamma bets on a power-law distribution with exponent \(\approx 2\). This post implements both, derives their implied priors, and shows numerically what the bets mean.

Unary: Geometric Prior

The encoding rule for unary is simple: to encode integer \(n \geq 1\), write \((n-1)\) zero bits followed by one 1 bit. The decoder reads bits until it sees the 1; the number of bits read is the decoded value.

Examples: \(1 \to\) 1, \(2 \to\) 01, \(3 \to\) 001, \(4 \to\) 0001.

struct Unary {
    using value_type = std::uint64_t;

    template<BitSink S>
    static void encode(value_type n, S& sink) {
        assert(n >= 1 && "Unary is undefined for n = 0");
        for (value_type i = 1; i < n; ++i) sink.write(false);
        sink.write(true);
    }

    template<BitSource S>
    static value_type decode(S& source) {
        value_type n = 1;
        while (!source.read()) ++n;
        return n;
    }
};

Length analysis. The codeword for \(n\) has length \(n\). The Kraft sum is \(\sum_{n=1}^{\infty} 2^{-n} = 1\): unary saturates Kraft exactly. The implied prior is \(p_n = 2^{-n}\): a geometric distribution with parameter \(1/2\), where each value is half as likely as the previous.

Optimality test. Because the implied prior is dyadic (all probabilities are powers of \(1/2\)) and Kraft saturates, unary achieves entropy exactly on this prior. For a 30-symbol truncation of geometric(1/2), the expected unary length equals the entropy to within the truncation tail (\(\approx 2^{-30}\)):

When you pick a code for integers, you are making a bet about what integers the source will produce. The bet lives in the codeword lengths, not in a separate parameter. This post makes that precise.

Universal Codes as Priors

You want to compress a stream of positive integers. Which code should you use?

The question has more structure than it appears. A code for integers assigns a codeword to each integer. The codeword for 1 is short, for 2 a bit longer, for 100 much longer. The relative lengths encode an implicit bet: what fraction of the stream will be 1s? What fraction will be 100s? If the bet matches the source, the average codeword length will be close to the theoretical minimum, the entropy. If the bet is wrong, you pay an overhead proportional to how wrong you are.

The bet is not a separate parameter. It lives in the codeword lengths themselves. This is the central claim of this post:

Every prefix-free code is a prior over the integers. The codeword lengths determine, up to normalization, a probability distribution. The code is optimal for exactly the sources that match that distribution.

This post makes that claim precise and implements the tools to measure it. The rest of the series (posts 4 through 12) examines ten specific codes and the priors they embody.

The Correspondence: Lengths to Priors

For a prefix-free code with codeword lengths \((l_1, l_2, \ldots, l_n)\), define the unnormalized weight of symbol \(i\) as \(w_i = 2^{-l_i}\). This is the fraction of the Kraft budget consumed by that codeword.

If the code saturates Kraft (meaning \(\sum_i 2^{-l_i} = 1\)), then the weights are already a valid probability distribution: \(p_i = 2^{-l_i}\). If the code does not saturate (meaning \(\sum_i 2^{-l_i} < 1\)), normalize: \(p_i = 2^{-l_i} / \sum_j 2^{-l_j}\).

This is the inverse of Shannon’s prescription. Shannon says: given a distribution \(p_i\), the optimal codeword length is \(\lceil -\log_2 p_i \rceil\) bits. We reverse the direction: given a length \(l_i\), the implied probability is \(2^{-l_i}\).

The function implied_prior computes this map:

inline std::vector<double> implied_prior(const std::vector<std::size_t>& lengths) {
    std::vector<double> probs;
    probs.reserve(lengths.size());
    double total = 0.0;
    for (std::size_t l : lengths) {
        double p = std::ldexp(1.0, -static_cast<int>(l));
        probs.push_back(p);
        total += p;
    }
    // Normalize if Kraft sum is less than 1.
    if (total < 1.0) {
        for (double& p : probs) p /= total;
    }
    return probs;
}

Two examples show the range of priors you get in practice.

Kraft’s inequality is necessary. McMillan’s theorem says it is also sufficient, and the proof is a construction.

McMillan’s Converse

The previous post in this series proved Kraft’s inequality: for any prefix-free binary code with codeword lengths \(l_1, l_2, \ldots, l_n\),

Every prefix-free code satisfies it. No exceptions. But necessity alone is not the useful direction. The question I want answered is the converse: given a length vector that satisfies Kraft, does a prefix-free code with those lengths actually exist?

Yes, and McMillan’s theorem (1956) proves it. Better still, the proof is a construction: given any Kraft-satisfying length vector, you can produce a specific prefix-free code with those exact lengths. No search required. No verification required after the fact. The construction always terminates, always produces a valid code, because Kraft pre-certifies that the budget is sufficient.

This post proves the constructive direction, then goes further. McMillan proved something stronger than just the prefix-free converse. He showed that even uniquely-decodable codes that are not prefix-free must satisfy Kraft. The consequence is worth sitting with: there is no advantage to non-prefix-free designs. If a code can be uniquely decoded, a prefix-free code with the same lengths exists. Prefix-freeness is not a restriction you impose for convenience. It is just the cleanest form of what unique decodability requires.

The Construction

The construction is a left-to-right walk through an imaginary binary trie. Sort the lengths, then assign codewords by taking the next available leaf at each step.

Concretely: fix a counter at zero, and for each length \(l_i\) (in sorted order), emit the binary representation of counter >> (l_max - l_i) left-padded to \(l_i\) bits. Then advance the counter by \(2^{l_{\max} - l_i}\), which skips past the entire subtree rooted at the just-assigned codeword. That advance ensures the next codeword starts at the first unoccupied leaf position in the depth-\(l_{\max}\) trie.

Work through the example from post 1: lengths \(\{1, 2, 3, 3\}\). Sort: \(1, 2, 3, 3\). Take \(l_{\max} = 3\).

• Length 1: counter is 0. Shift right by \(3 - 1 = 2\): emit 0 >> 2 = 0 as a 1-bit string, giving codeword "0". Advance counter by \(2^{3-1} = 4\). Counter is now 4.
• Length 2: counter is 4 (binary 100). Shift right by \(3 - 2 = 1\): emit 4 >> 1 = 2 as a 2-bit string, giving "10". Advance by \(2^{3-2} = 2\). Counter is now 6.
• Length 3: counter is 6 (binary 110). Shift right by \(3 - 3 = 0\): emit 6 as a 3-bit string, giving "110". Advance by \(2^0 = 1\). Counter is 7.
• Length 3: counter is 7 (binary 111). Emit 7 as a 3-bit string: "111". Advance by 1. Counter is 8.

Result: {"0", "10", "110", "111"}. This is exactly the example code from post 1. The construction recovered it directly from the length vector, without any search.

Every prefix-free code satisfies one inequality. That inequality is also sufficient. This post develops the necessary direction.

Kraft’s Inequality

I want a code where each symbol maps to a bit string, and where any concatenation of codewords can be decoded unambiguously. The simplest way to guarantee that is prefix-freeness: no codeword is a prefix of any other. A prefix-free code is self-delimiting. The decoder reads bits left-to-right and knows exactly when each codeword ends, with no lookahead and no length headers.

The question I keep returning to is: which collections of lengths are actually achievable? If I want four codewords of lengths 1, 2, 3, and 3, can I build a prefix-free code with those lengths? What if I want two codewords of length 1? (No: there are only two 1-bit strings, and they are prefixes of everything longer.)

Kraft’s inequality is the answer. A length vector \((l_1, l_2, \ldots, l_n)\) is achievable by a prefix-free binary code only if

This is the constraint you cannot escape. Any prefix-free code satisfies it. Any length vector that violates it cannot be realized as a prefix-free code, full stop.

The converse is also true: any length vector satisfying Kraft is realizable by some prefix-free code. That is McMillan’s theorem, and it is the subject of the next post in this series. This post develops the necessary direction: every prefix-free code satisfies Kraft.

The right tool for understanding why is the binary tree.

The Trie View

Represent each codeword as a path in a binary tree. Start at the root. For each bit, go left (0) or right (1). The codeword ends at a node, which I mark as a terminal. A code is prefix-free if and only if no terminal node has any descendants that are also terminals. Once you reach a terminal on the way down, you stop.

The example code \(\{A \to \texttt{0},\ B \to \texttt{10},\ C \to \texttt{110},\ D \to \texttt{111}\}\) has lengths \((1, 2, 3, 3)\). Its trie looks like this:

         root
        /    \
       0      1
      [A]    / \
            0   1
           [B]  / \
               0   1
              [C] [D]

A is at depth 1, left branch. B is at depth 2, right-then-left. C and D share a parent at depth 2, then split at depth 3. No codeword’s node is an ancestor of another’s: the code is prefix-free.

Twelve posts, twelve codes, one thesis that refused to change. This is the closing summary.

A. The Twelve Codes Together

Every post in this series answered a version of the same question: given a source of positive integers, how do you represent its values compactly as a sequence of bits? The answers differ in shape, in assumptions, and in which distribution each code implicitly expects.

Posts 1 and 2 (Kraft’s Inequality and McMillan’s Converse) established why prefix-free codes are the right unit of analysis. Post 3 (Universal Codes as Priors) named the frame: a code is a hypothesis about the source. Posts 4 through 10 filled in the catalogue. Posts 11 and 12 extended from integer coding to set representation, where the questions shift from “how long is this codeword?” to “how do you store membership?” and “how do you answer rank/select?”

Looking across all twelve, the main lesson is not that one code dominates. It is that the question “which code?” is always empirically answerable given a sample.

B. The Unifying Frame Restated

Post 3 introduced the codes-as-priors thesis with two instances behind it. We now have twelve. The thesis has not changed; it has only become more evidently true.

Huffman codes one symbol at a time. Arithmetic coding encodes the whole sequence as a single number. The difference is a factor of twelve, at least on the right source.

The Last Bit of Redundancy

Huffman coding gets expected codeword length within one bit of entropy. That is the best it can do, because codeword lengths must be integers while entropy is a real number.

The waste is structural. A symbol with probability $p = 0.7$ has optimal (fractional) length $-\log_2(0.7) \approx 0.515$ bits. Huffman rounds that up to 1 bit: 0.485 bits wasted per occurrence. For a nearly-deterministic source with $p_0 = 0.99$ and $p_1 = 0.01$, the entropy is $H \approx 0.081$ bits per symbol. Huffman is stuck at 1 bit per symbol. That is a factor-of-twelve gap, and Huffman cannot close it: a symbol that appears 99% of the time still gets a complete codeword.

Arithmetic coding steps back from per-symbol codewords entirely. It encodes an entire sequence as a single rational number in $[0, 1)$. The bit-length of that number converges to the entropy of the sequence as the sequence grows. No integer rounding, no per-symbol overhead.

This post builds an integer range coder in C++23 and demonstrates the factor-of-twelve improvement on the Bernoulli(0.99) source.

The Continuous View

Start with the unit interval $[0, 1)$. For a two-symbol source, partition it by probability: symbol 0 gets $[0, p_0)$ and symbol 1 gets $[p_0, 1)$.

To encode a sequence, begin with the full interval and narrow it with each symbol. After symbol $s_1$, restrict to the corresponding sub-interval. After $s_2$, apply the same proportional rule inside that sub-interval. After $L$ symbols, the interval has width $\prod_{i=1}^{L} p_{s_i}$.

Any number inside the final interval is a valid encoding. The shortest such number in binary requires approximately $-\log_2(\prod p_{s_i}) = \sum_{i=1}^{L} (-\log_2 p_{s_i})$ bits. As $L \to \infty$, bits per symbol approaches $H(p) = -\sum_k p_k \log_2 p_k$ exactly.

Decoding is the inverse: given the encoded number, determine at each step which sub-interval it falls in, recover the symbol, narrow the interval, and repeat.

The theory is complete. The practice is not. A real interval narrows exponentially fast: after a few dozen symbols you need arbitrary precision. The integer range coder fixes this with 32-bit arithmetic and a renormalization step.

Huffman coding is two things: the optimal length vector for a known distribution, and McMillan’s construction applied to that vector. This post develops both.

From Universal to Optimal

Every code in this series so far has been universal: no prior knowledge of the source distribution required. Elias gamma assigns shorter codewords to smaller integers regardless of which integers actually appear. Fibonacci does the same. VByte packs smaller values into fewer bytes without knowing whether your data clusters at the low end or the high end. Universal codes are defensive: they perform acceptably across a broad class of inputs by committing to none.

Huffman flips that stance. You bring a finite probability distribution. Huffman finds the prefix-free code with minimum expected codeword length for that distribution. The code is tuned to what you provide and will perform poorly on anything else. Call this the move from defensive to distribution-specific coding.

The payoff is real. Shannon’s source coding theorem says no prefix-free code can achieve expected length below $H(p) = -\sum_i p_i \log_2 p_i$. Huffman gets within 1 bit of that bound. For any prefix-free code, expected length satisfies

The upper bound comes from the integer-length constraint: each codeword is a whole number of bits, and $\lceil -\log_2 p_i \rceil \le -\log_2 p_i + 1$. Huffman is optimal subject to this constraint. No prefix-free code does better without abandoning whole-bit codewords.

That last clause points to the limit of this post and the subject of the next. Arithmetic coding breaks the integer-length constraint by assigning fractional bits in effect, reaching entropy exactly in the limit.

The Algorithm

The four steps of Huffman’s construction:

1. Create one leaf node per symbol, weighted by its probability.
2. Push all leaves into a min-priority queue (lowest weight first).
3. While the queue contains more than one node: extract the two lowest-weight nodes, merge them into a new internal node whose weight is their sum, push the internal node back.
4. The remaining node is the root. The path from root to each leaf encodes that leaf’s codeword (“0” for left, “1” for right).

Here is the complete implementation from huffman.hpp:

Every code in this series so far operates at bit granularity. VByte does not. It gives up bit-level precision for byte-alignment, and in production systems, that trade wins most of the time.

The Practical Question

Every code in this series so far operates at bit granularity. Elias gamma encodes 1 in a single bit. Fibonacci coding uses exactly as many bits as the Zeckendorf representation requires. Bit packing is theoretically attractive because it minimizes the number of bits written, which minimizes the encoded size.

But bit packing is computationally expensive. Reading or writing a single bit requires a shift, a mask, and often a branch to handle byte boundaries. Encoding a sequence of integers this way burns CPU cycles that scale with the number of integers, independent of their values. For high-throughput applications, the overhead of bit manipulation can easily exceed the savings from compact encoding.

VByte (also called Varint in Google’s ecosystem, and LEB128 in the DWARF debug format) trades a small amount of length efficiency for byte-alignment. The idea is simple: encode each integer as a sequence of 7-bit groups, one per byte, with a continuation flag in the high bit of each byte. The result is self-delimiting, compact for small values, and requires no bit-level manipulation to decode.

VByte is the encoding used by Protocol Buffers for all integer fields. It appears in Apache Arrow, Parquet, Snappy’s block format, LevelDB’s metadata, and most production columnar file formats. These are high-throughput systems. Byte-alignment is why VByte is their choice over the more compact universal codes from posts 4 through 7.

The Encoding

VByte splits an integer into 7-bit groups, starting from the least significant bits. Each group occupies one byte where bits 0 through 6 carry 7 data bits and bit 7 is a continuation flag: 1 means more bytes follow, 0 means this is the last byte.

A value in $[0, 127]$ fits in a single byte (continuation flag clear). A value in $[128, 16383]$ requires two bytes (first byte has flag set, second has flag clear). The pattern continues: each additional byte adds 7 bits of capacity.

Every code in this series so far has been fixed. Rice and Golomb are different: they take a parameter, and the parameter is your model of the data.

The First Parametric Code

Every code examined so far in this series has been monolithic. Unary coding is just unary coding. Elias gamma is just Elias gamma. Each one encodes all non-negative integers with a single fixed strategy. You do not get to choose anything about the code beyond whether to use it.

Rice and Golomb codes break this pattern. They are parametric: a single integer parameter, $k$ for Rice or $m$ for Golomb, tunes the code to a specific source distribution. Rice$(k)$ is not one code but a family of codes, one per value of $k$. Each member of the family is optimal for a specific geometric distribution. Choosing $k$ is choosing your prior precisely.

This matters because data sources are rarely uniform. Run-length encodings, inter-frame video differences, and the gap sequences in inverted indexes are all approximately geometrically distributed. If you know the mean of your source, you can pick $k$ so that Rice$(k)$ performs near-optimally, without the overhead of a Huffman table or arithmetic coding.

The key insight: for a geometric source with mean approximately $2^k$, Rice$(k)$ is within a small constant of entropy. No other universal code in this series achieves this. Elias gamma and delta perform well asymptotically but can be far from optimal for a specific geometric distribution with a known mean. Rice exploits that knowledge directly.

Rice Coding

Rice coding splits a non-negative integer $n$ into two parts: a quotient $q = \lfloor n / 2^k \rfloor = n \gg k$ and a remainder $r = n \bmod 2^k = n \mathbin{&} (2^k - 1)$.

The quotient is encoded in unary: $q$ zero bits followed by a stop bit of 1. The remainder is encoded in exactly $k$ bits, MSB first. The total codeword is the concatenation of these two parts.

Codeword examples for $k = 2$ (remainder is always 2 bits):

Codeword length: $(n \gg k) + 1 + k$ bits. The Kraft sum saturates to 1, so Rice is a complete prefix-free code.

Every code in this series so far has optimized expected length under some implied prior. Fibonacci coding does something different: it gives the decoder a way to recover from errors without help from a lower layer.

A Different Design Goal

All the codes in this series have aimed at the same target: assign short codewords to frequent symbols, with length growing roughly as $\log n$ for the $n$-th symbol under some implied prior. Elias gamma minimizes expected length for power-law distributions; delta and omega extend the recursion for heavier tails.

Fibonacci coding has a different goal. It does not optimize for average codeword length under a specific distribution. It optimizes for error resilience. In a stream of gamma-coded integers, a single bit flip in a codeword’s length prefix causes the decoder to misread that codeword’s length, then misread every subsequent codeword. The error propagates without limit until the decoder somehow reacquires sync. On a reliable channel this is a nonissue. On a noisy one, or in stored data that may have silently rotted, it is a serious problem.

Fibonacci coding avoids this. Every Fibonacci codeword ends in two consecutive 1 bits (“11”). This double-one marker appears nowhere else in the codeword. A single bit flip corrupts the codeword it hits, possibly spills into the next codeword, and then the decoder finds the next “11” and resynchronizes. At most two codewords are corrupted per error. The rest of the stream is intact.

The price is length overhead: Fibonacci codewords are approximately $1.44 \times \log_2 n$ bits long, compared to $\log_2 n$ bits for the entropy lower bound. On a reliable channel, that overhead is not worth paying. On a noisy channel, or in a long-running stream where rare bit errors must not lose the entire tail, the self-synchronization property is worth it.

Zeckendorf’s Theorem

Fibonacci numbers starting from $F_2 = 1$: $1, 2, 3, 5, 8, 13, 21, 34, \ldots$

Zeckendorf’s theorem: every positive integer $n$ has a unique representation as a sum of non-consecutive Fibonacci numbers. The greedy algorithm produces it by repeatedly subtracting the largest Fibonacci number that does not exceed $n$.

Elias gamma spends too many bits saying how many bits it will use. Delta fixes that. Omega takes the fix one step further. This post is about what happens when you apply recursion to the length prefix.

Where Gamma Stops Being Good

Elias gamma, from the previous post, encodes a positive integer $n$ in $2\lfloor \log_2 n \rfloor + 1$ bits: a unary count of $\lfloor \log_2 n \rfloor$ zeros, then a stop bit, then the $\lfloor \log_2 n \rfloor$ trailing binary bits of $n$. For small $n$ this is fine. For large $n$, nearly half the bits are spent on the unary prefix alone.

The unary prefix is the bottleneck. It encodes the length $L = \lfloor \log_2 n \rfloor + 1$ in the most wasteful possible way: one bit per unit. For $n = 256$, that is 8 zero bits just to say “the payload is 8 bits long.” The payload itself is also 8 bits, so you are paying a 100% overhead on the length announcement. That is bad, and it gets worse as $n$ grows.

The fix is obvious once you see it: encode $L$ itself in some shorter code instead of unary. Elias delta does exactly this, replacing the unary length prefix with a gamma-coded length. Elias omega takes the idea one step further and applies the recursion to itself, all the way down.

Both codes are universal: they assign finite codewords to every positive integer, and the expected codeword length is within a constant factor of optimal for any source whose probabilities decrease with $n$. The improvement over gamma is real and measurable once $n$ grows past a few dozen.

This post shows both implementations, their implied priors, and the crossover points where each code wins. As in the rest of this series, the code is pedagogical: each header stands alone and the struct-with-encode/decode pattern maps directly onto the PFC library’s EliasDelta and EliasOmega in codecs.hpp.

Elias Delta

Algorithm. Let $L = \lfloor \log_2 n \rfloor + 1$ (the bit-width of $n$, equivalently std::bit_width(n)).

1. Encode $L$ in Elias gamma.
2. Write the $L - 1$ trailing bits of $n$ after its implicit leading 1, MSB first.

Gamma encodes $L$ (a small integer) in $O(\log \log n)$ bits instead of $O(\log n)$ bits for the unary prefix. The payload is identical to gamma’s: the trailing bits of $n$. The total length is $O(\log n + \log \log n)$.

Unary is older than information theory. Elias gamma is its 1975 improvement. Together they span the gap between optimal-but-impractical and practical-but-nearly-optimal. This post derives what each code bets on, and shows numerically what that means.

Unary and Elias Gamma

Unary is the oldest code in this series. It predates information theory by centuries: a shepherd counting sheep on a stick is using unary. Mark one notch per sheep; count the notches to decode. The codeword for \(n\) is \(n\) tally marks. Its information-theoretic justification came later, when Shannon showed it is exactly optimal for a geometric source.

Elias gamma is the 1975 extension by Peter Elias. It brings the codeword length from \(O(n)\) to \(O(\log n)\), making it practical for numbers beyond small single digits, while keeping the prefix-free property that makes self-delimiting streams possible.

Both codes are instances of the claim from Universal Codes as Priors: every prefix-free code is a bet about the source. Unary bets on a geometric distribution with parameter \(1/2\). Gamma bets on a power-law distribution with exponent \(\approx 2\). This post implements both, derives their implied priors, and shows numerically what the bets mean.

Unary: Geometric Prior

The encoding rule for unary is simple: to encode integer \(n \geq 1\), write \((n-1)\) zero bits followed by one 1 bit. The decoder reads bits until it sees the 1; the number of bits read is the decoded value.

Examples: \(1 \to\) 1, \(2 \to\) 01, \(3 \to\) 001, \(4 \to\) 0001.

struct Unary {
    using value_type = std::uint64_t;

    template<BitSink S>
    static void encode(value_type n, S& sink) {
        assert(n >= 1 && "Unary is undefined for n = 0");
        for (value_type i = 1; i < n; ++i) sink.write(false);
        sink.write(true);
    }

    template<BitSource S>
    static value_type decode(S& source) {
        value_type n = 1;
        while (!source.read()) ++n;
        return n;
    }
};

Length analysis. The codeword for \(n\) has length \(n\). The Kraft sum is \(\sum_{n=1}^{\infty} 2^{-n} = 1\): unary saturates Kraft exactly. The implied prior is \(p_n = 2^{-n}\): a geometric distribution with parameter \(1/2\), where each value is half as likely as the previous.

Optimality test. Because the implied prior is dyadic (all probabilities are powers of \(1/2\)) and Kraft saturates, unary achieves entropy exactly on this prior. For a 30-symbol truncation of geometric(1/2), the expected unary length equals the entropy to within the truncation tail (\(\approx 2^{-30}\)):

When you pick a code for integers, you are making a bet about what integers the source will produce. The bet lives in the codeword lengths, not in a separate parameter. This post makes that precise.

Universal Codes as Priors

You want to compress a stream of positive integers. Which code should you use?

The question has more structure than it appears. A code for integers assigns a codeword to each integer. The codeword for 1 is short, for 2 a bit longer, for 100 much longer. The relative lengths encode an implicit bet: what fraction of the stream will be 1s? What fraction will be 100s? If the bet matches the source, the average codeword length will be close to the theoretical minimum, the entropy. If the bet is wrong, you pay an overhead proportional to how wrong you are.

The bet is not a separate parameter. It lives in the codeword lengths themselves. This is the central claim of this post:

Every prefix-free code is a prior over the integers. The codeword lengths determine, up to normalization, a probability distribution. The code is optimal for exactly the sources that match that distribution.

This post makes that claim precise and implements the tools to measure it. The rest of the series (posts 4 through 12) examines ten specific codes and the priors they embody.

The Correspondence: Lengths to Priors

For a prefix-free code with codeword lengths \((l_1, l_2, \ldots, l_n)\), define the unnormalized weight of symbol \(i\) as \(w_i = 2^{-l_i}\). This is the fraction of the Kraft budget consumed by that codeword.

If the code saturates Kraft (meaning \(\sum_i 2^{-l_i} = 1\)), then the weights are already a valid probability distribution: \(p_i = 2^{-l_i}\). If the code does not saturate (meaning \(\sum_i 2^{-l_i} < 1\)), normalize: \(p_i = 2^{-l_i} / \sum_j 2^{-l_j}\).

This is the inverse of Shannon’s prescription. Shannon says: given a distribution \(p_i\), the optimal codeword length is \(\lceil -\log_2 p_i \rceil\) bits. We reverse the direction: given a length \(l_i\), the implied probability is \(2^{-l_i}\).

The function implied_prior computes this map:

inline std::vector<double> implied_prior(const std::vector<std::size_t>& lengths) {
    std::vector<double> probs;
    probs.reserve(lengths.size());
    double total = 0.0;
    for (std::size_t l : lengths) {
        double p = std::ldexp(1.0, -static_cast<int>(l));
        probs.push_back(p);
        total += p;
    }
    // Normalize if Kraft sum is less than 1.
    if (total < 1.0) {
        for (double& p : probs) p /= total;
    }
    return probs;
}

Two examples show the range of priors you get in practice.

Kraft’s inequality is necessary. McMillan’s theorem says it is also sufficient, and the proof is a construction.

McMillan’s Converse

The previous post in this series proved Kraft’s inequality: for any prefix-free binary code with codeword lengths \(l_1, l_2, \ldots, l_n\),

Every prefix-free code satisfies it. No exceptions. But necessity alone is not the useful direction. The question I want answered is the converse: given a length vector that satisfies Kraft, does a prefix-free code with those lengths actually exist?

Yes, and McMillan’s theorem (1956) proves it. Better still, the proof is a construction: given any Kraft-satisfying length vector, you can produce a specific prefix-free code with those exact lengths. No search required. No verification required after the fact. The construction always terminates, always produces a valid code, because Kraft pre-certifies that the budget is sufficient.

This post proves the constructive direction, then goes further. McMillan proved something stronger than just the prefix-free converse. He showed that even uniquely-decodable codes that are not prefix-free must satisfy Kraft. The consequence is worth sitting with: there is no advantage to non-prefix-free designs. If a code can be uniquely decoded, a prefix-free code with the same lengths exists. Prefix-freeness is not a restriction you impose for convenience. It is just the cleanest form of what unique decodability requires.

The Construction

The construction is a left-to-right walk through an imaginary binary trie. Sort the lengths, then assign codewords by taking the next available leaf at each step.

Concretely: fix a counter at zero, and for each length \(l_i\) (in sorted order), emit the binary representation of counter >> (l_max - l_i) left-padded to \(l_i\) bits. Then advance the counter by \(2^{l_{\max} - l_i}\), which skips past the entire subtree rooted at the just-assigned codeword. That advance ensures the next codeword starts at the first unoccupied leaf position in the depth-\(l_{\max}\) trie.

Work through the example from post 1: lengths \(\{1, 2, 3, 3\}\). Sort: \(1, 2, 3, 3\). Take \(l_{\max} = 3\).

• Length 1: counter is 0. Shift right by \(3 - 1 = 2\): emit 0 >> 2 = 0 as a 1-bit string, giving codeword "0". Advance counter by \(2^{3-1} = 4\). Counter is now 4.
• Length 2: counter is 4 (binary 100). Shift right by \(3 - 2 = 1\): emit 4 >> 1 = 2 as a 2-bit string, giving "10". Advance by \(2^{3-2} = 2\). Counter is now 6.
• Length 3: counter is 6 (binary 110). Shift right by \(3 - 3 = 0\): emit 6 as a 3-bit string, giving "110". Advance by \(2^0 = 1\). Counter is 7.
• Length 3: counter is 7 (binary 111). Emit 7 as a 3-bit string: "111". Advance by 1. Counter is 8.

Result: {"0", "10", "110", "111"}. This is exactly the example code from post 1. The construction recovered it directly from the length vector, without any search.

Every prefix-free code satisfies one inequality. That inequality is also sufficient. This post develops the necessary direction.

Kraft’s Inequality

I want a code where each symbol maps to a bit string, and where any concatenation of codewords can be decoded unambiguously. The simplest way to guarantee that is prefix-freeness: no codeword is a prefix of any other. A prefix-free code is self-delimiting. The decoder reads bits left-to-right and knows exactly when each codeword ends, with no lookahead and no length headers.

The question I keep returning to is: which collections of lengths are actually achievable? If I want four codewords of lengths 1, 2, 3, and 3, can I build a prefix-free code with those lengths? What if I want two codewords of length 1? (No: there are only two 1-bit strings, and they are prefixes of everything longer.)

Kraft’s inequality is the answer. A length vector \((l_1, l_2, \ldots, l_n)\) is achievable by a prefix-free binary code only if

This is the constraint you cannot escape. Any prefix-free code satisfies it. Any length vector that violates it cannot be realized as a prefix-free code, full stop.

The converse is also true: any length vector satisfying Kraft is realizable by some prefix-free code. That is McMillan’s theorem, and it is the subject of the next post in this series. This post develops the necessary direction: every prefix-free code satisfies Kraft.

The right tool for understanding why is the binary tree.

The Trie View

Represent each codeword as a path in a binary tree. Start at the root. For each bit, go left (0) or right (1). The codeword ends at a node, which I mark as a terminal. A code is prefix-free if and only if no terminal node has any descendants that are also terminals. Once you reach a terminal on the way down, you stop.

The example code \(\{A \to \texttt{0},\ B \to \texttt{10},\ C \to \texttt{110},\ D \to \texttt{111}\}\) has lengths \((1, 2, 3, 3)\). Its trie looks like this:

         root
        /    \
       0      1
      [A]    / \
            0   1
           [B]  / \
               0   1
              [C] [D]

A is at depth 1, left branch. B is at depth 2, right-then-left. C and D share a parent at depth 2, then split at depth 3. No codeword’s node is an ancestor of another’s: the code is prefix-free.

Every type decomposes structurally. So does its codec.

Codecs as Functors

You have an optional<vector<pair<int, string>>>. The type decomposes structurally: it is an optional of a free monoid of products of an integer and a string. That decomposition is not an observation about memory layout. It is a statement about the algebraic structure of the type.

Now ask: does the codec decompose the same way?

If the answer is yes, you stop writing one-off encoders. You build a codec for optional<T> from a codec for T. You build a codec for vector<T> from a codec for T. The codec for optional<vector<pair<int, string>>> assembles from its parts with no manual layout decisions, no hand-placed length headers, no ad-hoc format negotiation.

This post argues that the answer is always yes, and shows what the machinery looks like. The thesis: codecs are not ad-hoc bit formats. They are constructions on the algebraic structure of types. The algebraic structure of a type determines its codec, the same way it determines its algorithms.

This extends Stepanov’s claim. The peasant algorithm post showed that algorithms arise from algebraic structure. The homomorphism post showed that structure-preserving maps are the natural morphisms. Here, we show the codec itself is a structure-preserving map, and that it lifts from leaf types to compound types by the same algebraic logic.

Bit I/O: The Foundation

Before combinators, we need concrete bit I/O. The approach taken here follows Stepanov’s move in the algorithm posts: state the concept first, then provide a model.

Two concepts govern bit-level I/O:

template<typename T>
concept BitSink = requires(T& s, bool bit) {
    { s.write(bit) } -> std::same_as<void>;
};

template<typename T>
concept BitSource = requires(T& s) {
    { s.read() } -> std::same_as<bool>;
    { s.peek() } -> std::convertible_to<bool>;
};

A BitSink accepts bits. A BitSource supplies them. A codec is an algorithm parameterized over BitSink and BitSource, not a class hierarchy. This is Stepanov’s move at the bit level: require only what the algorithm needs, let anything that satisfies the concept participate.

The standard models are BitWriter and BitReader, which pack bits into byte buffers in LSB-first order:

class BitWriter {
    std::span<std::uint8_t> buf_;
    std::size_t byte_idx_ = 0;
    std::uint8_t byte_ = 0;
    std::uint8_t bit_pos_ = 0;
public:
    explicit BitWriter(std::span<std::uint8_t> buf) noexcept : buf_(buf) {}

    void write(bool bit) noexcept {
        byte_ |= (bit ? std::uint8_t{1} : std::uint8_t{0}) << bit_pos_;
        if (++bit_pos_ == 8) {
            buf_[byte_idx_++] = byte_;
            byte_ = 0;
            bit_pos_ = 0;
        }
    }

    void align() noexcept {
        if (bit_pos_ > 0) {
            buf_[byte_idx_++] = byte_;
            byte_ = 0;
            bit_pos_ = 0;
        }
    }

    [[nodiscard]] std::size_t bytes_written() const noexcept {
        return byte_idx_ + (bit_pos_ > 0 ? 1 : 0);
    }
};

class BitReader {
    std::span<const std::uint8_t> buf_;
    std::size_t byte_idx_ = 0;
    std::uint8_t bit_pos_ = 0;
public:
    explicit BitReader(std::span<const std::uint8_t> buf) noexcept : buf_(buf) {}

    bool read() noexcept {
        bool bit = ((buf_[byte_idx_] >> bit_pos_) & 1) != 0;
        if (++bit_pos_ == 8) {
            ++byte_idx_;
            bit_pos_ = 0;
        }
        return bit;
    }

    [[nodiscard]] bool peek() const noexcept {
        return byte_idx_ < buf_.size();
    }
};

A codec concept rounds out the three-concept core:

The free monoid on a type lifts to bit space. It lifts injectively only when the element codec is prefix-free.

Prefix-Free Codes and the Free Monoid

You have a list of unsigned integers. Encode the list as a single bit string.

Fixed-width encoding wastes space. If you allocate 64 bits per integer, small values like 1 or 7 cost as much as values near \(2^{64}\). Variable-width encoding recovers that space, but immediately raises a harder question: where does one encoded integer end and the next begin?

Two escape routes. First, prefix each encoded item with its length. That works, but the length headers are overhead, and you now need a codec for the lengths as well. Second, choose a code where the structure of the codewords makes boundaries unambiguous without any headers. These are prefix-free codes, and this is the right answer, in a precise categorical sense.

The “precise categorical sense” is what this post develops. Encoding a list as the concatenation of encoded elements is a monoid homomorphism from the free monoid on \(T\) to the monoid of bit strings under concatenation. The universal property of the free monoid guarantees this homomorphism always exists. The question of whether the decoder can invert it comes down to exactly one property of the element codec: whether it is prefix-free.

The Free Monoid, Recalled

A monoid is a set with an associative binary operation and an identity element. The free monoid on a set \(S\) is the set of all finite sequences of elements from \(S\), with concatenation as the operation and the empty sequence as the identity.

“Free” means no equations hold except those forced by the monoid axioms. Nothing is identified with anything else. If you need commutativity or idempotency, you quotient the free monoid by additional equations. But the free monoid itself imposes nothing beyond associativity and identity.

The universal property says: given any monoid \(M\) and any function \(f: S \to M\), there is exactly one monoid homomorphism \(\hat{f}: \text{Free}(S) \to M\) that extends \(f\). That unique extension is fold:

where \(\cdot\) is the operation in \(M\). The free-algebra post develops this in full. For this post, the one fact that matters is that fold is canonical: it is the unique way to extend a per-element map to a list-consuming function that respects the monoid structure.

Lists are everywhere in programming. Not because they are convenient. Because they are algebraically universal.

Why Lists?

Arrays are more cache-friendly. Hash maps have better lookup. Yet lists (sequences, vectors, streams) remain the default container in nearly every language. The standard explanation is convention, or ease of construction. The real explanation is algebraic.

A list is the free monoid. It is the most general monoid you can build from a set of generators. And the universal property of free monoids says that fold, the operation that processes a list element by element, is not a design pattern. It is a theorem.

The Free Monoid

Start with a set \(S\) of generators. The free monoid on \(S\) is the set of all finite sequences of elements from \(S\), with concatenation as the operation and the empty sequence as the identity.

“Free” means: no equations hold except those forced by the monoid axioms (associativity and identity). In particular:

• \([a, b] \neq [b, a]\). Commutativity is not imposed.
• \([a, a] \neq [a]\). Idempotency is not imposed.
• \([a, b, c] = [a] \cdot [b] \cdot [c]\). Every sequence is a product of singletons.

In C++:

template<typename T>
class free_monoid {
    std::vector<T> elements_;
public:
    free_monoid() = default;
    explicit free_monoid(T x) : elements_{std::move(x)} {}
    // ...
};

// Monoid operations via ADL
template<typename T>
free_monoid<T> op(const free_monoid<T>& a, const free_monoid<T>& b);  // concatenation

template<typename T>
free_monoid<T> identity(const free_monoid<T>&);  // empty sequence

The free_monoid<int> is the type of finite sequences of integers. Its operation is concatenation. It satisfies the Monoid concept. And it is the most general monoid on int: no structure beyond associativity and identity.

The Universal Property

Here is the key fact. Given any function \(f: S \to M\) where \(M\) is a monoid, there exists a unique monoid homomorphism \(\overline{f}: \text{Free}(S) \to M\) extending \(f\). This homomorphism is defined by:

In code:

template<Monoid M, typename T, typename F>
M extend(F f, const free_monoid<T>& xs) {
    M result = identity(M{});
    for (const auto& x : xs.elements())
        result = op(result, f(x));
    return result;
}

This is fold. The universal property says fold is the only structure-preserving way to interpret a list in a monoid. Any function that respects the monoid structure must agree with fold.

Fold Is a Theorem

When you write std::accumulate or std::reduce, you are invoking the universal property. The homomorphism condition:

A homomorphism is a function that preserves algebraic structure. This post shows that fold, sum, length, and even the logarithm are all the same idea.

Structures and Maps

The series so far has built up algebraic structures: monoids in the peasant post, rings in the modular post, Euclidean domains in the polynomial post, product monoids in the accumulator post. But structures alone are half the story. The other half is the maps between them.

A homomorphism is a function \(f: A \to B\) between two structures of the same kind that preserves the operation:

The operation on the left is in \(A\). The operation on the right is in \(B\). The function \(f\) “commutes” with the operation. That is the entire definition.

The Concept

We need a monoid concept for this post. A monoid has an associative binary operation and an identity element. As always, the operations are ADL free functions:

template<typename M>
concept Monoid = std::semiregular<M> &&
    requires(M a, M b) {
        { op(a, b) } -> std::convertible_to<M>;
        { identity(a) } -> std::convertible_to<M>;
    };

And a runtime check that a function preserves the monoid structure:

template<Monoid A, Monoid B, typename F>
bool is_homomorphism(F f, const A& a1, const A& a2) {
    return f(op(a1, a2)) == op(f(a1), f(a2));
}

This tests one pair of inputs. It is not a proof, but it catches violations.

Examples Everywhere

Length. The string monoid (concatenation, empty string) maps to the integers under addition. The map is length. And length is a homomorphism:

The length of a concatenation is the sum of the lengths. This is not a coincidence. It is the homomorphism property.

Sum. Lists of integers under concatenation form a monoid. The integers under addition form a monoid. The map sum is a homomorphism:

Product. Same source monoid, different target. Now the integers under multiplication:

Logarithm. The positive reals under multiplication form a monoid. The reals under addition form a monoid. The logarithm maps one to the other:

This is the defining property of the logarithm, stated as algebra. The logarithm is a homomorphism from \((\mathbb{R}^+, \times, 1)\) to \((\mathbb{R}, +, 0)\).

Count. For any type \(T\), the function that sends a list to its length is a homomorphism from the list monoid to \((\mathbb{Z}, +, 0)\). This is the same as the string length example, generalized.

Lattices have two operations of a different kind than rings. The structure determines a fixed-point algorithm.

Two Operations, Different Rules

Monoids have one binary operation. Rings have two (addition and multiplication) linked by distributivity. Lattices also have two operations, but with different laws entirely.

A lattice is a set with two operations:

• meet (\(\wedge\)): greatest lower bound
• join (\(\vee\)): least upper bound

Both are idempotent, commutative, and associative. And they satisfy the absorption laws:

Absorption is what distinguishes lattices from a pair of unrelated monoids. It ties meet and join together: knowing one constrains the other.

A bounded lattice adds a least element (bottom, \(\bot\)) and a greatest element (top, \(\top\)). Bottom is the identity for join, top is the identity for meet.

In C++20 concepts, with ADL free functions:

template<typename L>
concept Lattice = std::semiregular<L> &&
    requires(L a, L b) {
        { meet(a, b) } -> std::convertible_to<L>;
        { join(a, b) } -> std::convertible_to<L>;
    };

template<typename L>
concept BoundedLattice = Lattice<L> &&
    requires(L a) {
        { bottom(a) } -> std::convertible_to<L>;
        { top(a) } -> std::convertible_to<L>;
    };

Four Examples

Sign lattice. Abstract signs of integers: bottom (unreachable), negative, zero, positive, top (unknown). Meet is greatest lower bound, join is least upper bound in the Hasse diagram. This is the classic abstract interpretation domain. You can define abstract arithmetic on it: pos * neg = neg, neg + neg = neg, pos + neg = top.

Intervals. Closed intervals \([a, b]\) ordered by inclusion. Meet is intersection. Join is the smallest enclosing interval. Bottom is the empty interval. Top is the full range. This is the foundation of interval arithmetic.

Divisors. Positive integers ordered by divisibility. Meet is gcd, join is lcm. Bottom is 1 (divides everything), top is 0 (everything divides 0). Lattice structure appearing in number theory.

Power sets. Subsets of \(\{0, \ldots, N-1\}\). Meet is intersection (bitwise AND), join is union (bitwise OR). Bottom is the empty set, top is the full set.

All four satisfy BoundedLattice. All four satisfy the same laws. The concept constrains the interface; the laws constrain the semantics.

The Algorithm: Tarski’s Fixed-Point Theorem

Here is the payoff. Tarski’s theorem: any monotone function on a complete lattice has a least fixed point, computable by iterating from bottom.

The peasant post showed that power() works on any monoid. But what happens when you have two operations instead of one?

From Monoids to Semirings

A monoid is one operation with an identity element. The peasant algorithm exploits this: give it any monoid and it computes powers by repeated squaring. The accumulator post used the same structure for streaming statistics.

A semiring is two monoids on the same set, linked by a compatibility condition. Formally, a semiring \((S, +, \times, 0, 1)\) satisfies:

1. \((S, +, 0)\) is a commutative monoid
2. \((S, \times, 1)\) is a monoid
3. \(\times\) distributes over \(+\): \(a \times (b + c) = a \times b + a \times c\)
4. \(0\) annihilates: \(0 \times a = a \times 0 = 0\)

In C++20, using the same ADL free functions as the peasant post:

template<typename S>
concept Semiring = std::semiregular<S> &&
    requires(S a, S b) {
        { a + b } -> std::convertible_to<S>;
        { a * b } -> std::convertible_to<S>;
        { zero(a) } -> std::convertible_to<S>;
        { one(a) } -> std::convertible_to<S>;
    };

The concept captures syntax. The axioms (associativity, distributivity, annihilation) are semantic requirements that the programmer must ensure.

Five Semirings

The ordinary integers are a semiring. But there are others, and each one corresponds to a different graph problem.

The naming of the tropical semirings is counterintuitive but standard. In the tropical min semiring, the “addition” is min and the “multiplication” is ordinary addition. This matters because matrix multiplication uses both operations, and we need the algebraic structure to be correct: the inner product of a row and column computes the best path through an intermediate node.

Each semiring is a small struct with operator+, operator*, and ADL functions zero() and one():

struct boolean_semiring {
    bool val;
    constexpr boolean_semiring operator+(boolean_semiring rhs) const {
        return boolean_semiring(val || rhs.val);
    }
    constexpr boolean_semiring operator*(boolean_semiring rhs) const {
        return boolean_semiring(val && rhs.val);
    }
};
constexpr boolean_semiring zero(boolean_semiring) { return boolean_semiring(false); }
constexpr boolean_semiring one(boolean_semiring)  { return boolean_semiring(true); }

Matrices Over a Semiring

Matrix multiplication requires addition and multiplication. That is exactly what a semiring provides. If \(S\) is a semiring, then \(n \times n\) matrices over \(S\) form a semiring too, with entry-wise addition and the usual row-times-column product (using \(S\)’s operations).

Accumulators are monoids. The same algebraic structure from the peasant post, in a different domain.

Accumulators as Monoids

An accumulator processes a stream of values, maintaining state that can be queried at any point. Write a class with operator+= for each statistic you need. Sum, mean, variance, min, max. Five statistics, five classes.

The problem is combinations. Sum and min? Write a sixth class. Sum, min, and max? A seventh. Every new combination requires new code.

But every accumulator has the same structure:

1. Process a value incrementally: operator+=(value)
2. Combine with another accumulator of the same type: operator+=(accumulator)
3. Extract a result: .eval()

Default construction gives you an empty accumulator: the identity element. Combination via += is associative. Together, a monoid. The peasant post used the same structure for exponentiation. Here we use it for streaming computation.

In C++20 concepts:

template<typename A>
concept Accumulator = std::semiregular<A> &&
    requires(A a, A b, typename A::value_type v) {
        typename A::value_type;
        { a += v } -> std::same_as<A&>;   // process one value
        { a += b } -> std::same_as<A&>;   // combine two accumulators
        { a.eval() };                       // extract result
    };

KBN: Compensated Summation

The simplest accumulator is a sum. But naive floating-point summation accumulates O(n) rounding error:

double sum = 0.0;
sum += 1.0;
for (int i = 0; i < 1'000'000; ++i)
    sum += 1e-10;
// Expected: 1.0001    Actual: ~1.00009999999...8

When you add a tiny number to a large one, the tiny number’s low-order bits get dropped. After a million additions, these losses add up.

Kahan-Babuska-Neumaier (KBN) summation tracks what gets lost:

template<std::floating_point T>
class kbn_sum {
    T sum_ = T(0);
    T comp_ = T(0);   // compensation for lost bits

public:
    using value_type = T;

    constexpr kbn_sum& operator+=(const T& v) {
        T t = sum_ + v;
        comp_ += abs_(sum_) >= abs_(v) ? (sum_ - t) + v
                                       : (v - t) + sum_;
        sum_ = t;
        return *this;
    }

    constexpr T eval() const { return sum_ + comp_; }
};

The correction term comp_ recovers the bits that floating-point addition drops. O(1) error instead of O(n), regardless of sequence length.

kbn_sum is a monoid:

• Identity: kbn_sum{} (sum=0, compensation=0)
• Operation: a += b (combine two compensated sums)

Welford: Online Mean and Variance

Computing the mean is just sum/count. Variance is harder. The textbook formula \(\sigma^2 = \frac{1}{n}\sum(x_i - \bar{x})^2\) requires two passes: one for the mean, one for the deviations.

Welford’s algorithm computes both in a single pass:

welford& operator+=(const T& v) {
    ++n_;
    T delta = v - mean_;
    mean_ += delta / static_cast<T>(n_);
    T delta2 = v - mean_;   // uses *updated* mean
    m2_ += delta * delta2;
    return *this;
}

delta uses the old mean, delta2 uses the new mean. Their product accumulates into m2_, the sum of squared deviations. At any point, variance = m2_ / n.

What makes this a monoid is the combination formula. Given two independent Welford accumulators with means \(\bar{x}_A, \bar{x}_B\) and counts \(n_A, n_B\), Chan et al. showed how to merge them:

Many structures come in pairs. Recognizing duality lets you transfer insights between domains.

The Motivating Example

This collection includes two approaches to automatic differentiation:

• Forward mode (in dual): Propagate derivatives alongside values, from inputs toward outputs
• Reverse mode (in autodiff): Build a graph during forward evaluation, then propagate gradients backward from outputs toward inputs

These aren’t just two implementations of the same idea. They’re duals, mirror images with complementary strengths.

Forward mode computes one column of the Jacobian per pass. If \(f: \mathbb{R}^n \to \mathbb{R}^m\), computing the full Jacobian takes \(n\) passes. Reverse mode computes one row per pass, \(m\) passes for the full Jacobian.

For neural network training, we have many inputs (millions of parameters) and one output (the loss). Reverse mode wins overwhelmingly: one backward pass gives all gradients. This is why backpropagation dominates deep learning.

For sensitivity analysis with few parameters and many outputs, forward mode wins. Same algorithm structure, opposite traversal direction, complementary use cases.

The mathematical explanation: forward mode computes Jacobian-vector products (\(Jv\)); reverse mode computes vector-Jacobian products (\(v^T J\)). These are transposes of each other. Duality is transposition.

Push vs Pull

Consider two ways to traverse a sequence:

Pull (iterator/consumer controls):

for (auto it = seq.begin(); it != seq.end(); ++it) {
    process(*it);  // Consumer pulls each element
}

Push (producer controls):

seq.for_each([](auto x) {
    process(x);  // Producer pushes each element
});

Same traversal. Same elements processed. But control flow is reversed:

C++ ranges are pull-based: view | filter | transform creates an iterator that pulls through the pipeline. Reactive streams (Rx) are push-based: events flow through a pipeline of observers.

These are duals. Given a pull-based algorithm, you can mechanically derive its push-based counterpart by reversing who initiates each step. The transformation preserves correctness because it’s just changing direction, not content.

Encode vs Decode

Compression algorithms come in pairs:

// Encoder: structure -> bits
auto encode(const Document& doc) -> Bitstream;

// Decoder: bits -> structure
auto decode(const Bitstream& bits) -> Document;

These must be inverses: decode(encode(x)) == x. But their implementations are often strikingly different:

A reflection on eleven explorations in generic programming

The Question Behind the Code

What do these computations have in common?

• Computing the millionth Fibonacci number
• Finding the shortest path between cities in a weighted graph
• Calculating compound interest over thirty years
• Composing ten 3D rotations into one
• Repeating a string n times

The answer: they’re all computed by the same twenty lines of code.

template<typename T>
constexpr T power(T const& base, T exp) {
    if (exp == zero(exp)) return one(exp);
    if (exp == one(exp))  return base;

    return even(exp)
        ? square(power(base, half(exp)))
        : product(base, power(base, decrement(exp)));
}

This shouldn’t work. Fibonacci numbers involve integer sequences. Shortest paths involve graphs. Rotations involve 3D geometry. Different domains, different mathematics.

Yet they share structure. Once you see it, a single algorithm serves them all.

This collection of eleven blog posts is an extended meditation on one idea: algorithms arise from algebraic structure. The posts cover different domains (number theory, calculus, linear algebra, polymorphism) but they circle the same insight. Recognize the structure; the algorithm follows.

The Principle

Alex Stepanov articulated this most clearly in Elements of Programming: “Generic programming is about abstracting and classifying algorithms and data structures.” But the deeper point is how to abstract. Not by common syntax or superficial similarity, but by the algebraic laws a type obeys.

Why does structure appear everywhere? Because reality has structure. The algebraic structures we discover in programming (groups, rings, monoids) are the same structures physicists discover in nature. Rotations form a group. Spacetime transformations form a group. This isn’t coincidence. We’re uncovering patterns that exist.

Noether’s theorem makes this precise: every continuous symmetry corresponds to a conservation law. Time-translation symmetry gives conservation of energy. Space-translation symmetry gives conservation of momentum. Rotational symmetry gives conservation of angular momentum. The symmetry groups of physics are algebraic structures.

When we recognize “this is a monoid” in our code, we’re tapping into the same mathematical substrate that governs physical law. The algorithms follow because the structure constrains what’s possible, both in computation and in nature.

Consider the power() function above. What does it require?

• An associative binary operation (so we can regroup: \((a \cdot b) \cdot c = a \cdot (b \cdot c)\))
• An identity element (so \(1 \cdot x = x \cdot 1 = x\))
• Halving and parity testing on the exponent

That’s it. Any type providing these operations, with these laws, can use this algorithm. The requirements are algebraic, not syntactic.

The rest of this series asks: given a structure, what algorithms does it support? This post inverts the question.

The Flip Side

The peasant post showed that power() works on any monoid. The semirings post showed that matrix multiplication over different semirings solves six graph problems. The thread running through the whole series is: algorithms arise from algebraic structure.

But there’s a flip side that we haven’t addressed directly.

Sometimes you’re stuck with an expensive algorithm not because the problem is hard, but because you’re working in the wrong algebra. Change the algebra, and the algorithm becomes trivial. The cost shows up somewhere else, always. But if the cheap operation is the one you actually need, you win.

This is a very old idea. Napier invented logarithms in 1614 to turn multiplication into addition. What’s worth noticing is that logarithms, odds ratios, tropical semirings, and quaternions are all doing the same thing.

The Pattern

A computational basis transform takes values from one domain and represents them in another, where different operations are cheap:

Each row follows the same structure. A transform \(\varphi: D \to D’\) makes some operations cheaper and others more expensive. There is no free lunch.

This is not a deep theorem. It’s almost tautological: if you could make everything cheaper by relabeling, the labels would already be the standard ones. But making the pattern explicit helps you recognize when you’re paying for an operation you don’t need.

Three Examples

Log Space

The most familiar example. You have a million small probabilities and need their product.

// Standard: underflows to 0 after ~30 terms
double product = 1.0;
for (double p : probs) product *= p;  // 0.0

// Log domain: addition instead of multiplication
mutatio::lgd product(1.0);
for (double p : probs) product = product * mutatio::lgd(p);
// product.log() is finite. product.value() would overflow,
// but you stay in log space.

The algebra changed from \((\mathbb{R}^+, \times)\) to \((\mathbb{R}, +)\). Multiplication became addition. The isomorphism is \(\log\).

Tropical Semirings

The semirings post showed this already. Replace \((+, \times)\) with \((\min, +)\) and matrix multiplication becomes shortest-path computation. That’s the same move: you changed the semiring to make the algorithm you wanted (all-pairs shortest paths) fall out of a generic operation (matrix power) that you already had.

A synthesis of three earlier posts, comparing forward-mode AD, reverse-mode AD, and numerical differentiation.

Computing derivatives shows up everywhere: optimization, machine learning, physics simulation, numerical analysis. This series has explored three distinct approaches:

1. Forward-mode AD via dual numbers
2. Reverse-mode AD via computational graphs
3. Numerical differentiation via finite differences

Each has different strengths. The right choice depends on the shape of your problem.

The Landscape

The key point: problem structure determines the best method.

Forward-Mode AD: Dual Numbers

Forward-mode AD extends numbers with an infinitesimal \(\varepsilon\) where \(\varepsilon^2 = 0\). The derivative falls out of the arithmetic for free:

// f(x) = x^3 - 3x + 1
// f'(x) = 3x^2 - 3

auto x = dual<double>::variable(2.0);  // x = 2, dx = 1
auto f = x*x*x - 3.0*x + 1.0;

std::cout << f.value() << "\n";       // 3.0
std::cout << f.derivative() << "\n";  // 9.0

Strengths:

• Simple implementation (operator overloading)
• No memory overhead
• Naturally composable for higher derivatives
• Works with any function of overloaded operators

When to use:

• Single input variable (or few inputs)
• Computing Jacobian-vector products
• Higher-order derivatives via nesting
• Sensitivity analysis along one direction

Complexity: One forward pass per input variable. For f: R^n -> R^m, computing the full Jacobian requires n passes.

Reverse-Mode AD: Computational Graphs

Reverse-mode AD builds a computational graph during the forward pass, then propagates gradients backward via the chain rule:

auto f = [](const auto& x) {
    return sum(pow(x, 2.0));  // f(x) = sum(x^2)
};

auto df = grad(f);  // Returns gradient function
auto gradient = df(x);  // One backward pass for all partials

Strengths:

• O(1) backward passes regardless of input dimension
• Powers modern deep learning (backpropagation)
• Efficient for loss functions: f: R^n -> R

When to use:

• Many inputs, scalar output (neural networks)
• Computing vector-Jacobian products
• Optimization where you need the full gradient

Complexity: One forward pass to build the graph, one backward pass to compute all gradients. Memory scales with the number of operations because you have to store intermediate values.

Numerical Differentiation: Finite Differences

Approximate the derivative using the limit definition:

// Central difference: f'(x) ~ (f(x+h) - f(x-h)) / 2h
double df = central_difference(f, x);

Strengths:

Numerical integration (quadrature) for C++20.

Overview

The definite integral is the signed area under a curve:

Most functions do not have closed-form antiderivatives, so we approximate integrals numerically using quadrature rules: weighted sums of function evaluations.

Different rules choose different nodes x_i and weights w_i. The tradeoff is always accuracy vs. computational cost.

Quick Start

#include <integration/integrate.hpp>
#include <cmath>
#include <iostream>

int main() {
    using namespace integration;

    // integral from 0 to pi of sin(x) dx = 2
    double result = integrate([](double x) { return std::sin(x); }, 0.0, 3.14159265);

    std::cout << "integral of sin(x) dx = " << result << "\n";  // ~2.0
}

The Quadrature Zoo

Basic Rules

double m = midpoint_rule(f, a, b);
double t = trapezoidal_rule(f, a, b);
double s = simpsons_rule(f, a, b);

Composite Rules

Divide [a,b] into n subintervals and apply the basic rule to each:

// Error: O(h^2) where h = (b-a)/n
double m = composite_midpoint(f, a, b, 100);
double t = composite_trapezoidal(f, a, b, 100);

// Error: O(h^4) - much more accurate!
double s = composite_simpsons(f, a, b, 100);  // n must be even

Gauss-Legendre Quadrature

The optimal choice: n points exactly integrate polynomials of degree 2n-1.

// 5-point Gauss-Legendre: exact for degree <= 9
double g = gauss_legendre<5>(f, a, b);

Adaptive Integration

Automatically refines where the function is difficult:

// Recommended for general use
double result = integrate(f, a, b);              // Default tolerance 1e-10
double result = integrate(f, a, b, 1e-12);       // Custom tolerance

// With error estimate
auto [value, error, evals] = integrate_with_error(f, a, b);

Deriving Simpson’s Rule

Taylor expand \(f(x)\) around the midpoint \(m = (a+b)/2\). Odd powers vanish by symmetry when integrating from \(a\) to \(b\):

Simpson’s rule is the unique combination of endpoint and midpoint values that cancels the \(h^2\) error:

It also cancels the h^3 term. Simpson gets a “bonus degree” of accuracy for free. This is one of those happy accidents in numerical analysis.

Why Gauss-Legendre is Optimal

With \(n\) evaluation points, we have \(2n\) free parameters (\(n\) nodes + \(n\) weights). We can match \(2n\) conditions: exact integration of \(1, x, x^2, \ldots, x^{2n-1}\).

The nodes turn out to be roots of the \(n\)-th Legendre polynomial \(P_n(x)\). Orthogonal polynomials arise naturally from the optimization. This is not a coincidence. It is the same reason Legendre polynomials show up in approximation theory: they are the optimal basis for polynomial approximation on [-1,1] with the right inner product.

Forward-mode automatic differentiation via dual numbers for C++20.

Overview

Dual numbers are a simple yet powerful technique for computing exact derivatives. The key insight: if we extend our number system with an element epsilon where epsilon^2 = 0, then evaluating f(x + epsilon) yields f(x) + epsilon * f'(x). The derivative emerges automatically from the algebra.

Quick Start

#include <dual/dual.hpp>
#include <iostream>

int main() {
    using namespace dual;

    // Create a dual variable at x = 2
    auto x = dual<double>::variable(2.0);

    // Compute f(x) = x^3 - 3x + 1
    auto f = x*x*x - 3.0*x + 1.0;

    std::cout << "f(2) = " << f.value() << "\n";       // 3.0
    std::cout << "f'(2) = " << f.derivative() << "\n"; // 9.0
}

The Mathematics

A dual number has the form \(a + b\varepsilon\) where \(\varepsilon^2 = 0\). Arithmetic follows naturally:

Notice how the \(bd\varepsilon^2\) term vanishes because \(\varepsilon^2 = 0\).

For a function \(f\), Taylor expansion gives:

If we set \(b = 1\) (marking \(x\) as “the variable we’re differentiating with respect to”), then:

The derivative appears as the coefficient of epsilon!

API Reference

dual

The core dual number type.

// Create a variable for differentiation
auto x = dual<double>::variable(3.0);  // x = 3, dx = 1

// Create a constant
auto c = dual<double>::constant(2.0);  // c = 2, dc = 0

// Access values
double val = x.value();       // 3.0
double deriv = x.derivative(); // 1.0

// Arithmetic operators: +, -, *, /
auto y = sin(x*x) + exp(-x);

// Convenience function
auto [value, deriv] = differentiate([](auto x) { return x*x; }, 3.0);

Mathematical Functions

All standard math functions are supported with correct derivative propagation:

• Basic: sqrt, cbrt, abs
• Exponential: exp, exp2, expm1, log, log2, log10, log1p
• Trigonometric: sin, cos, tan, asin, acos, atan, atan2
• Hyperbolic: sinh, cosh, tanh, asinh, acosh, atanh
• Power: pow, hypot
• Special: erf, erfc

Higher-Order Derivatives

Second derivatives with dual2:

auto result = differentiate2([](auto x) { return sin(x); }, 1.0);
// result.value  = sin(1)
// result.first  = cos(1)
// result.second = -sin(1)

Arbitrary order with jets:

// Compute f, f', f'', f''', f'''' at x = 1
auto derivs = derivatives<4>([](auto x) { return exp(x); }, 1.0);
// All derivatives of e^x at x=1 equal e

Forward vs Reverse Mode

This library implements forward mode AD:

The same GCD algorithm works for integers and polynomials. That’s not a coincidence. It’s because both are Euclidean domains.

The Observation

// For integers: gcd(48, 18) = 6
// For polynomials: gcd(x^3 - 1, x^2 - 1) = x - 1

// Same algorithm, different types
template<euclidean_domain E>
E gcd(E a, E b) {
    while (b != E(0)) {
        a = std::exchange(b, a % b);
    }
    return a;
}

That template compiles and works correctly for both integers and polynomials. The reason it works is algebraic: both types support division with remainder, and the remainder is always “smaller” than the divisor in a well-defined sense.

Quick Start

#include <polynomials/polynomial.hpp>
#include <iostream>

using namespace poly;

int main() {
    // Create polynomial x^2 - 1 = (x-1)(x+1)
    auto p = polynomial<double>{-1, 0, 1};

    // Create polynomial x^3 - 1 = (x-1)(x^2+x+1)
    auto q = polynomial<double>{-1, 0, 0, 1};

    // GCD should be (x - 1)
    auto g = gcd(p, q);

    std::cout << "gcd(x^2-1, x^3-1) has degree " << g.degree() << "\n";  // 1

    // Find roots of x^2 - 1
    auto roots = find_roots(p, -10.0, 10.0);
    for (double r : roots) {
        std::cout << "Root: " << r << "\n";  // -1 and 1
    }
}

API Reference

Creating Polynomials

// From dense coefficients (a[i] = coefficient of x^i)
polynomial<double> p{1, -2, 1};  // 1 - 2x + x^2

// Monomial: coefficient * x^degree
auto m = polynomial<double>::monomial(3.0, 4);  // 3x^4

// The variable x
auto x = polynomial<double>::x();  // x

// Constant
polynomial<double> c{5.0};  // 5

Arithmetic

auto sum = p + q;
auto diff = p - q;
auto prod = p * q;
auto [quot, rem] = divmod(p, q);  // Division with remainder
auto quot_only = p / q;
auto rem_only = p % q;

GCD and Related

auto g = gcd(p, q);                    // Greatest common divisor
auto [g, s, t] = extended_gcd(p, q);   // Bezout: g = p*s + q*t
auto l = lcm(p, q);                    // Least common multiple
bool d = divides(p, q);                // Does p divide q?

Evaluation and Calculus

double val = evaluate(p, x);           // p(x)
auto dp = derivative(p);               // p'(x)
auto integral = antiderivative(p);     // integral of p
auto roots = find_roots(p, -10, 10);   // All real roots in interval
auto crit = stationary_points(p, -10, 10);  // Where p'(x) = 0

The Euclidean Domain Structure

What makes this work is shared algebraic structure. A Euclidean domain has a norm function and a division algorithm where the remainder is always smaller than the divisor:

The GCD algorithm doesn’t care which type it’s operating on. It only needs the division-with-remainder property. Stepanov’s whole point is exactly this: algorithms arise from algebraic structure. When you recognize that polynomials and integers share the same abstract structure, you immediately get:

Floating-point lies to you.

double x = 0.1 + 0.2;
std::cout << (x == 0.3);  // Prints 0 (false!)

The number 0.1 has no exact binary representation, for the same reason 1/3 has no exact decimal representation. Floating-point represents numbers as m x 2^e, and most decimal fractions don’t land on a power of two.

Rational arithmetic fixes this. 1/3 stays exactly 1/3.

The Representation

A rational number is a pair (numerator, denominator) kept in lowest terms:

template<std::integral T>
class rat {
    T num_;  // numerator (carries sign)
    T den_;  // denominator (always positive)

    void reduce() {
        T g = std::gcd(abs(num_), den_);
        num_ /= g;
        den_ /= g;
    }
};

Three invariants, always maintained:

1. The denominator is positive (sign lives in the numerator)
2. GCD(|num|, den) = 1 (always reduced)
3. Zero is uniquely 0/1

Arithmetic

Addition needs a common denominator:

Then reduce. In code:

rat operator+(rat const& rhs) const {
    return rat(num_ * rhs.den_ + rhs.num_ * den_,
               den_ * rhs.den_);
}

The constructor calls reduce() automatically.

Multiplication is simpler:

Division multiplies by the reciprocal:

Exact Comparison

No floating-point fuzziness. Two reduced rationals are equal iff their numerators and denominators match:

bool operator==(rat const& rhs) const {
    return num_ == rhs.num_ && den_ == rhs.den_;
}

For ordering, cross-multiply (valid because denominators are positive):

The Mediant

The mediant of a/b and c/d is (a+c)/(b+d). It’s not the average. It has different, more interesting properties:

rat mediant(rat const& a, rat const& b) {
    return rat(a.numerator() + b.numerator(),
               a.denominator() + b.denominator());
}

If a/b < c/d, then a/b < mediant < c/d. The mediant is always in lowest terms when a/b and c/d are neighbors in the Stern-Brocot tree. And mediants generate all positive rationals exactly once.

The Stern-Brocot Tree

Start with 0/1 and 1/0 (representing infinity). Repeatedly take mediants:

Level 0:     0/1                     1/0
Level 1:     0/1       1/1           1/0
Level 2:     0/1   1/2   1/1   2/1   1/0
Level 3:  0/1 1/3 1/2 2/3 1/1 3/2 2/1 3/1 1/0

Every positive rational appears exactly once. The path from root to any node encodes its continued fraction. This connects to best rational approximations and Farey sequences.

GCD Ties Everything Together

Reducing fractions requires GCD. The algorithm is Euclid’s, from around 300 BCE:

T gcd(T a, T b) {
    while (b != 0) {
        a = a % b;
        std::swap(a, b);
    }
    return a;
}

The same algorithm works for any Euclidean domain. That’s not a coincidence. It’s a consequence of the algebraic structure.

Rational numbers form a field: every non-zero element has a multiplicative inverse (the reciprocal). The requirement that denominators be non-zero and fractions reduced comes from this algebraic structure, not from arbitrary convention.

The problem is combinatorial. You have N algorithms (sort, search, find, copy) and M containers (array, list, tree, hash table). The naive approach: implement each algorithm for each container. That is NxM implementations.

The insight is to interpose an abstraction layer.

The Iterator Abstraction

Instead of algorithms knowing about containers directly, we define iterator categories, capabilities that algorithms require and containers provide:

Input: Single-pass read. You can advance (++) and dereference (*), but once you move forward, you cannot go back. Stream-like.

Forward: Multi-pass. You can iterate multiple times; begin() always gives the same starting point.

Bidirectional: Can go backward (--). Enables algorithms like reverse iteration.

Random-access: Can jump anywhere (+n, []). Enables binary search, sorting.

This is a hierarchy of requirements. Each level adds capabilities and enables more algorithms. An algorithm declares the weakest category it needs, and any container providing at least that category works.

A True Input Iterator

The input iterator category exists for a reason. Here is a working example that reads entropy from /dev/urandom:

#include <fstream>
#include <iterator>
#include <cstdint>

struct entropy_iterator {
    using iterator_category = std::input_iterator_tag;
    using value_type        = uint8_t;
    using difference_type   = std::ptrdiff_t;
    using pointer           = const uint8_t*;
    using reference         = uint8_t;  // returns by value, not reference

    std::ifstream* source = nullptr;
    uint8_t byte = 0;

    entropy_iterator() = default;  // sentinel (end iterator)

    explicit entropy_iterator(std::ifstream& s) : source(&s) {
        ++(*this);  // prime the first byte
    }

    uint8_t operator*() const { return byte; }

    entropy_iterator& operator++() {
        if (source && source->good()) {
            source->read(reinterpret_cast<char*>(&byte), 1);
            if (!source->good()) source = nullptr;
        }
        return *this;
    }

    entropy_iterator operator++(int) {
        auto tmp = *this;
        ++(*this);
        return tmp;
    }

    bool operator==(const entropy_iterator& other) const {
        return source == other.source;
    }
};

Use it like any input iterator:

int main() {
    std::ifstream urandom("/dev/urandom", std::ios::binary);
    entropy_iterator it(urandom);

    // generate 16 random bytes
    std::vector<uint8_t> key(16);
    std::copy_n(it, 16, key.begin());

    // or use in algorithms
    int sum = 0;
    for (int i = 0; i < 1000; ++i, ++it) {
        sum += *it;
    }
    // sum ≈ 127500 (mean of uniform [0,255] × 1000)
}

Each ++ consumes a fresh entropy byte from the kernel. You literally cannot iterate twice over the same sequence. This is why the input iterator category exists: some sources are inherently single-pass. Claiming forward iterator capabilities would be a lie.

The same pattern applies to network streams, sensor readings, and any source where data is consumed by reading it.

The Payoff

Now binary_search does not need to know about vectors, deques, or sorted arrays. It only needs random-access iterators. The algorithm expresses its requirements; the container provides capabilities. They compose through the iterator abstraction.

The Miller-Rabin primality test and the mathematics of certainty

The Problem

Given a large number n, is it prime? Trial division up to sqrt(n) is too slow for cryptographic-sized numbers. We need something faster, and we are willing to accept “probably prime” with quantifiable certainty.

Fermat’s Little Theorem

For prime \(p\) and any \(a\) not divisible by \(p\):

This suggests a test: pick random a, compute a^(n-1) mod n. If the result is not 1, n is definitely composite. But if it is 1, n might be prime, or might be a Carmichael number that fools this test.

The Miller-Rabin Improvement

Miller and Rabin observed something stronger. For odd prime \(p\), write \(p-1 = 2^r \cdot d\) (factor out all 2s). Then the sequence:

must either:

1. Start with 1, or
2. Contain -1 (i.e., p-1) somewhere before reaching 1

Why? Because the only square roots of 1 mod p are plus or minus 1. If we ever see 1 without first seeing -1, we have found a non-trivial square root of 1, proving n is composite.

The Witness Test

bool witness_test(int64_t n, int64_t a) {
    // Write n-1 = 2^r x d
    int64_t d = n - 1;
    int r = 0;
    while ((d & 1) == 0) { d >>= 1; r++; }

    // Compute x = a^d mod n
    int64_t x = mod_pow(a, d, n);

    if (x == 1 || x == n - 1) return true;  // Probably prime

    // Square r-1 times, looking for n-1
    for (int i = 1; i < r; i++) {
        x = (x * x) % n;
        if (x == n - 1) return true;
    }

    return false;  // Definitely composite
}

If witness_test(n, a) returns false, n is definitely composite. The value a is a “witness” to compositeness.

Error Bounds

Here is the part I find most satisfying. For any composite n, at least 3/4 of all possible witnesses a in [2, n-2] will detect it. Each random witness has at most 1/4 chance of failing to detect a composite.

With \(k\) independent witnesses:

The error drops exponentially. 40 witnesses gives you a false positive probability smaller than the chance of a cosmic ray flipping a bit in your RAM during the computation.

Parameterizing by Error

Rather than asking “how many iterations?”, ask “what error rate is acceptable?”:

Finite algebraic structures and what they teach us about algorithms

The Stepanov Perspective

Stepanov’s central insight: algorithms arise from algebraic structure. The same algorithm that works on integers works on matrices, polynomials, and modular integers, not by accident, but because they share algebraic properties.

Integers modulo N form a ring: a set with addition and multiplication satisfying familiar laws. When N is prime, it is a field, meaning every non-zero element has a multiplicative inverse. Understanding these structures tells you which algorithms apply where.

The Ring Z/NZ

Integers modulo N, written Z/NZ, are equivalence classes:

• 0 = {…, -N, 0, N, 2N, …}
• 1 = {…, -N+1, 1, N+1, 2N+1, …}
• …

Operations are inherited from integers:

[a] + [b] = [a + b]
[a] x [b] = [a x b]

The implementation keeps one representative per class, in [0, N):

template<int64_t N>
struct mod_int {
    int64_t v;  // Always in [0, N)

    static constexpr int64_t normalize(int64_t x) {
        x %= N;
        return x < 0 ? x + N : x;
    }

    constexpr mod_int(int64_t x) : v(normalize(x)) {}
};

Ring Axioms

A ring (R, +, x) satisfies:

Addition forms an abelian group:

• Associative: (a + b) + c = a + (b + c)
• Identity: a + 0 = a
• Inverses: a + (-a) = 0
• Commutative: a + b = b + a

Multiplication is a monoid:

• Associative: (a x b) x c = a x (b x c)
• Identity: a x 1 = 1 x a = a

Distributive:

• a x (b + c) = a x b + a x c
• (b + c) x a = b x a + c x a

These axioms enable algorithms. Power-by-squaring works because multiplication is associative. The extended GCD works because of the ring structure.

Fermat’s Little Theorem

When N is prime, something special happens: every non-zero element has a multiplicative inverse. The set of non-zero elements forms a multiplicative group of order N-1.

Fermat’s Little Theorem: for prime \(p\) and \(a\) not congruent to \(0 \pmod{p}\):

This gives us the inverse:

So \(a^{p-2}\) is the multiplicative inverse of \(a\):

constexpr mod_int inverse() const {
    return pow(N - 2);  // Using repeated squaring
}

The Connection to Peasant

The power function uses the same peasant algorithm from the previous post:

How the Russian peasant algorithm reveals the universal structure of exponentiation

The Algorithm

Russian peasants had a clever method for multiplication that does not require memorizing times tables. To compute 23 x 17:

23    17
11    34     (halve, double)
 5    68
 2   136
 1   272

Add the right column wherever the left is odd: 17 + 34 + 68 + 272 = 391. That is 23 x 17.

Why does this work? Because we are really computing:

23 x 17 = (16 + 4 + 2 + 1) x 17 = 16x17 + 4x17 + 2x17 + 17

The algorithm only needs three operations on the multiplier:

• half(n), integer division by 2
• even(n), test if divisible by 2
• Addition on the result

From Multiplication to Exponentiation

Here is the insight that makes this interesting: the same algorithm computes powers.

Replace “add to accumulator” with “multiply into accumulator” and “double the multiplicand” with “square the base”:

T power(T base, int exp) {
    T result = 1;
    while (exp > 0) {
        if (!even(exp)) result = result * base;
        base = base * base;
        exp = half(exp);
    }
    return result;
}

This is O(log n) multiplications instead of O(n). Computing 2^1000 takes about 10 multiplications, not 1000.

The Monoid Connection

The peasant algorithm works whenever you have:

1. An associative binary operation *
2. An identity element 1 where 1 * x = x * 1 = x

This structure is called a monoid. The algorithm computes x * x * ... * x (n times) using O(log n) operations.

What makes this powerful is that many things form monoids:

Why Associativity Unlocks Efficiency

Why does the peasant algorithm achieve O(log n) instead of O(n)? The answer lies in a single algebraic law: associativity.

Associativity says \((a \cdot b) \cdot c = a \cdot (b \cdot c)\). This looks innocuous, but it means we can restructure computation without changing results. Consider computing \(a^8\):

Naive:     a x a x a x a x a x a x a x a     (7 multiplications)
Peasant:   ((a^2)^2)^2                        (3 multiplications)

Both produce the same answer because we can freely regroup. The peasant algorithm exploits this freedom systematically: instead of accumulating one factor at a time, it squares intermediate results and combines them.

I keep building the same thing.

Not the same product — the products are different. One indexes a Hugo blog. One indexes AI conversations. One consolidates medical records from three hospitals. One catalogs a hundred git repositories. But underneath, they all have the same skeleton. After the fifth time, I think the skeleton deserves a name.

The pattern

Domain files (ground truth)
    ↓ index
SQLite database (read-only cache, FTS5)
    ↓ expose
MCP server (tools + resources → AI assistant)

That’s it. Three layers. The domain files are always canonical — the database is a disposable cache you can rebuild from them at any time. SQLite gives you structured queries, full-text search, and JSON extraction over data that was previously trapped in flat files. MCP exposes it to an AI assistant that can write SQL, retrieve content, and (in some cases) create new content.

Here’s the inventory:

Five projects. Five completely different domains. One architecture.

Why SQLite

SQLite is the most deployed database in history. It’s on every phone, every browser, every Python installation. But that’s not why I use it.

I use it because it solves three problems at once:

Structured queries over unstructured data. Hugo front matter is YAML trapped inside markdown files. Medical records are scattered across three incompatible EHR export formats. AI conversations are JSON trees with branching paths. SQLite turns all of these into tables you can JOIN, GROUP BY, and aggregate. json_extract() handles the long tail of fields that don’t fit a fixed schema.

Full-text search. FTS5 with porter stemming and unicode61 tokenization gives you relevance-ranked search across any text corpus. No Elasticsearch, no external service, no running daemon. Just a virtual table that lives in the same database file.

Time is finite in ways I can’t ignore. That changes which questions about code feel important.

I read a post arguing that the most valuable programming skill in 2026 is deleting code. The thesis: AI generates code faster than anyone can review it, so the real value is in curation and subtraction. Code is a liability, not an asset.

I agree with the observation. I disagree with the prescription.

I keep building the same thing.

Not the same product — the products are different. One indexes a Hugo blog. One indexes AI conversations. One consolidates medical records from three hospitals. One catalogs a hundred git repositories. But underneath, they all have the same skeleton. After the fifth time, I think the skeleton deserves a name.

The pattern

Domain files (ground truth)
    ↓ index
SQLite database (read-only cache, FTS5)
    ↓ expose
MCP server (tools + resources → AI assistant)

That’s it. Three layers. The domain files are always canonical — the database is a disposable cache you can rebuild from them at any time. SQLite gives you structured queries, full-text search, and JSON extraction over data that was previously trapped in flat files. MCP exposes it to an AI assistant that can write SQL, retrieve content, and (in some cases) create new content.

Here’s the inventory:

Five projects. Five completely different domains. One architecture.

Why SQLite

SQLite is the most deployed database in history. It’s on every phone, every browser, every Python installation. But that’s not why I use it.

I use it because it solves three problems at once:

Structured queries over unstructured data. Hugo front matter is YAML trapped inside markdown files. Medical records are scattered across three incompatible EHR export formats. AI conversations are JSON trees with branching paths. SQLite turns all of these into tables you can JOIN, GROUP BY, and aggregate. json_extract() handles the long tail of fields that don’t fit a fixed schema.

Full-text search. FTS5 with porter stemming and unicode61 tokenization gives you relevance-ranked search across any text corpus. No Elasticsearch, no external service, no running daemon. Just a virtual table that lives in the same database file.

I have cancer. My oncologist is at one hospital system (Siteman/BJC), my primary care doctor at another, and my earlier treatment history lives at a third (Anderson, where my first oncologist practiced). Patient portals are fine for browsing, but they don’t answer questions. They show you your data one lab result at a time, one note at a time, one visit at a time.

I wanted to run queries against my medical records. Correlate lab trends with treatment changes. Generate structured question lists before oncology visits. Ask “what changed since my last appointment” and get a real answer. That means getting the data out of the portal and into something programmable.

Chartfold loads EHR exports into SQLite and exposes them to Claude via MCP.

Time is finite in ways I can’t ignore. That changes which questions about code feel important.

I read a post arguing that the most valuable programming skill in 2026 is deleting code. The thesis: AI generates code faster than anyone can review it, so the real value is in curation and subtraction. Code is a liability, not an asset.

I agree with the observation. I disagree with the prescription.

I keep building the same thing.

Not the same product — the products are different. One indexes a Hugo blog. One indexes AI conversations. One consolidates medical records from three hospitals. One catalogs a hundred git repositories. But underneath, they all have the same skeleton. After the fifth time, I think the skeleton deserves a name.

The pattern

Domain files (ground truth)
    ↓ index
SQLite database (read-only cache, FTS5)
    ↓ expose
MCP server (tools + resources → AI assistant)

That’s it. Three layers. The domain files are always canonical — the database is a disposable cache you can rebuild from them at any time. SQLite gives you structured queries, full-text search, and JSON extraction over data that was previously trapped in flat files. MCP exposes it to an AI assistant that can write SQL, retrieve content, and (in some cases) create new content.

Here’s the inventory:

Five projects. Five completely different domains. One architecture.

Why SQLite

SQLite is the most deployed database in history. It’s on every phone, every browser, every Python installation. But that’s not why I use it.

I use it because it solves three problems at once:

Structured queries over unstructured data. Hugo front matter is YAML trapped inside markdown files. Medical records are scattered across three incompatible EHR export formats. AI conversations are JSON trees with branching paths. SQLite turns all of these into tables you can JOIN, GROUP BY, and aggregate. json_extract() handles the long tail of fields that don’t fit a fixed schema.

Full-text search. FTS5 with porter stemming and unicode61 tokenization gives you relevance-ranked search across any text corpus. No Elasticsearch, no external service, no running daemon. Just a virtual table that lives in the same database file.

I keep building the same thing.

Not the same product — the products are different. One indexes a Hugo blog. One indexes AI conversations. One consolidates medical records from three hospitals. One catalogs a hundred git repositories. But underneath, they all have the same skeleton. After the fifth time, I think the skeleton deserves a name.

The pattern

Domain files (ground truth)
    ↓ index
SQLite database (read-only cache, FTS5)
    ↓ expose
MCP server (tools + resources → AI assistant)

That’s it. Three layers. The domain files are always canonical — the database is a disposable cache you can rebuild from them at any time. SQLite gives you structured queries, full-text search, and JSON extraction over data that was previously trapped in flat files. MCP exposes it to an AI assistant that can write SQL, retrieve content, and (in some cases) create new content.

Here’s the inventory:

Five projects. Five completely different domains. One architecture.

Why SQLite

SQLite is the most deployed database in history. It’s on every phone, every browser, every Python installation. But that’s not why I use it.

I use it because it solves three problems at once:

Structured queries over unstructured data. Hugo front matter is YAML trapped inside markdown files. Medical records are scattered across three incompatible EHR export formats. AI conversations are JSON trees with branching paths. SQLite turns all of these into tables you can JOIN, GROUP BY, and aggregate. json_extract() handles the long tail of fields that don’t fit a fixed schema.

Full-text search. FTS5 with porter stemming and unicode61 tokenization gives you relevance-ranked search across any text corpus. No Elasticsearch, no external service, no running daemon. Just a virtual table that lives in the same database file.

I keep building the same thing.

Not the same product — the products are different. One indexes a Hugo blog. One indexes AI conversations. One consolidates medical records from three hospitals. One catalogs a hundred git repositories. But underneath, they all have the same skeleton. After the fifth time, I think the skeleton deserves a name.

The pattern

Domain files (ground truth)
    ↓ index
SQLite database (read-only cache, FTS5)
    ↓ expose
MCP server (tools + resources → AI assistant)

That’s it. Three layers. The domain files are always canonical — the database is a disposable cache you can rebuild from them at any time. SQLite gives you structured queries, full-text search, and JSON extraction over data that was previously trapped in flat files. MCP exposes it to an AI assistant that can write SQL, retrieve content, and (in some cases) create new content.

Here’s the inventory:

Five projects. Five completely different domains. One architecture.

Why SQLite

SQLite is the most deployed database in history. It’s on every phone, every browser, every Python installation. But that’s not why I use it.

I use it because it solves three problems at once:

Structured queries over unstructured data. Hugo front matter is YAML trapped inside markdown files. Medical records are scattered across three incompatible EHR export formats. AI conversations are JSON trees with branching paths. SQLite turns all of these into tables you can JOIN, GROUP BY, and aggregate. json_extract() handles the long tail of fields that don’t fit a fixed schema.

Full-text search. FTS5 with porter stemming and unicode61 tokenization gives you relevance-ranked search across any text corpus. No Elasticsearch, no external service, no running daemon. Just a virtual table that lives in the same database file.

I have around 100 git repos. When I’m working with Claude Code on one of them, the AI has deep knowledge of that repo but zero awareness of the rest. Questions like “which of my repos already has a fuzzy search implementation?” or “what other projects use this pattern?” require me to go dig around manually.

repoindex fixes that.

The Idea

Separation of concerns:

Claude Code (deep work on ONE repo)
         |
         |  "What else do I have?"
         |  "Which repos need X?"
         v
    repoindex (collection awareness)
         |
         +-- repo://...     -> what exists
         +-- tags://...     -> organization
         +-- stats://...    -> aggregations
         +-- events://...   -> what happened

Claude Code works inside repositories. repoindex knows about repositories: metadata, tags, status, relationships. Together they give you full portfolio awareness.

MCP Server Integration

The most useful part is the MCP (Model Context Protocol) server. Add it to your Claude Code configuration and the AI can query your collection directly:

• “Which of my Python repos don’t have a LICENSE file?”
• “What repos have I updated in the last week?”
• “Show me all projects tagged with ml”

The server exposes resources like repo://, tags://, stats://, and events:// that Claude Code reads to understand your portfolio.

Core Features

Tag-Based Organization. Hierarchical tags for categorizing repos. Tags can be explicit (repoindex tag add myproject topic:ml) or implicit (derived automatically from language, directory, features).

Query Language. Filter repos with expressions:

repoindex query "language == 'Python' and 'ml' in tags"
repoindex query "stars > 10 and has:docs"

Event Tracking. What happened across your collection:

repoindex events --since 7d --pretty

New releases, tags, PyPI publishes, all in one view.

JSONL Output. Every command outputs newline-delimited JSON by default, so it plays well with Unix pipelines:

repoindex status | jq 'select(.status.clean == false)'

Installation

Available on PyPI:

pip install repoindex

Configure your repository directories and start indexing:

repoindex config generate
repoindex list --pretty

Why the Rename?

This was previously called ghops. The new name is more honest about what it does: it indexes repositories. The old name implied GitHub-specific operations, but the tool works with any git repo.

I had a bookmark manager. Then an ebook library manager. Then a chat history manager. Each started with the standard CRUD CLI:

btk add https://example.com --tags python,tutorial
btk list --tag python
btk search "async"
btk delete 1234

ebk import book.pdf --author "Knuth"
ebk list --author Knuth
ebk search "algorithms"

This works fine until you have 10,000+ bookmarks organized with hierarchical tags like programming/python/async, research/ml/transformers, work/clients/acme. Your ebook library has similar structure. Your exported chat conversations from Claude, ChatGPT, and Copilot are piling up.

Traditional CRUD commands become unwieldy:

btk list --tag programming/python/async/io --format json | jq '.[].title'
ebk list --category "Computer Science/Algorithms/Graph Theory" --limit 50
ctk search "machine learning" --source ChatGPT --date-from 2024-01-01

Each command requires precise arguments. Each tool has different flag conventions. You can’t navigate your data. You can only query it. And queries require knowing exactly what you’re looking for.

The insight: everything is a file

When I have thousands of source files organized in directories, I don’t run:

list-files --path /src/components/auth --extension .tsx

I run:

cd src/components/auth
ls *.tsx

The difference matters. With a filesystem, I can navigate incrementally (cd from general to specific), explore (ls to see what’s there), compose (cat file | grep pattern | wc -l), and use familiar tools (find, grep, xargs, pipes, redirection).

What if my bookmarks, ebooks, and chat histories were filesystems?

The pattern

Over the past year, I built six Python tools that all follow the same architecture:

Each tool provides:

1. A stateless CLI for scripting: btk bookmark add URL, ebk import book.pdf
2. An interactive shell with a virtual filesystem: btk shell, ebk shell, ctk chat
3. POSIX-like commands: cd, ls, pwd, cat, mv, cp, rm, find, grep
4. Unix pipeline support: most commands output JSONL by default for piping

The interesting part is the shell.

Navigating 10,000 bookmarks

Live recording captured with asciinema. You can pause, copy text, and replay. The entire recording is 78KB of text.

I should be upfront about what happened here. I did not compute coprime Ramsey numbers. I did not write 92 Python modules or 5,922 tests. I did not build SAT encodings or run survival analysis on Erdős problems.

Claude Code (Opus 4.6) did all of that. I told it what to look at, asked it to keep going, and occasionally said things like “try to disprove our discoveries” and “be aggressive.” The AI did the rest. 131 subagents, 78,000 lines of code, three minted DOIs. In one session.

I’m writing this down because I think it’s worth documenting what that looks like from the human side of the keyboard.

The Setup

Terence Tao maintains a database of 1,183 Erdős problems on GitHub. Each problem has tags, OEIS links, resolution status, and sometimes prize money. The database was updated in August 2025 to link problems to integer sequences. Since then, 213 problems have been solved, many with AI assistance.

I had been poking at this database on and off for a few months. I had some Python scripts, some partial Lean proofs, a few computational results. Nothing organized. The codebase had bugs (the kind where a random sampling heuristic silently gives you the wrong answer and you don’t notice for weeks).

I started a Claude Code session intending to fix those bugs. Then I said “iterate.” Then I kept saying “iterate.”

What Claude Found

The headline result is a family of numbers that, as far as anyone can tell, nobody had studied before.

Take the integers 1 through n. Connect every coprime pair with an edge. This is the coprime graph. Now 2-color every edge. The coprime Ramsey number R_cop(k) is the smallest n where every 2-coloring must contain a monochromatic complete subgraph of size k.

Classical Ramsey: R(3,3) = 6. Coprime Ramsey: R_cop(3) = 11.

The value R_cop(4) = 59 required SAT solving (Glucose4 via pysat). A random sampling heuristic had said 20. It was off by a factor of three. The SAT solver finds avoiding colorings instantly at every n up to 58. At n = 59 (prime, coprime to everything below it), no avoiding coloring exists. This was verified by an independent implementation built from scratch by a separate adversarial agent.

Facts and Degrees

In classical logic, something is true or false. The cat is on the mat, or it is not. A patient has a fever, or they do not. There is no middle ground.

Fuzzy logic adds a dial.

Instead of true/false, every statement carries a degree of belief – a number between 0 and 1. A degree of 1.0 means certainty. A degree of 0.0 means we have no belief at all. And everything in between is fair game.

Here is the simplest possible fuzzy fact:

# A fuzzy fact: "Rex has hair" with 85% confidence
engine.add_fact("has-hair", ["rex"], 0.85)

The predicate is has-hair. The argument is rex. The degree is 0.85. Maybe we observed Rex from a distance, or the photo was blurry. We are fairly sure Rex has hair, but not certain.

This is the building block of everything that follows. A fuzzy knowledge base is just a collection of these facts, each with its own degree. Some facts we are sure about (deg=1.0). Others are tentative guesses (deg=0.3). The engine treats them all the same way – it just pays attention to the number.

One important detail: when two sources assert the same fact with different degrees, the engine keeps the higher one. This is called fuzzy-OR. If one sensor says has-hair(rex) at 0.85 and another says it at 0.92, the engine stores 0.92. Optimistic, but reasonable – the stronger evidence wins.

engine.add_fact("has-hair", ["rex"], 0.85)
engine.add_fact("has-hair", ["rex"], 0.92)  # fuzzy-OR: keeps 0.92

In the widget below, you can create fuzzy facts and drag the degree slider to see how the degree changes the visual representation. A fact at 1.0 is solid and bright. A fact at 0.1 is faded, barely there. This is not just decoration – it is the engine’s uncertainty, made visible.

Rules

Facts alone are inert. To reason, we need rules – if-then statements that produce new facts from existing ones.

A fuzzy rule looks like this: “If X has hair, then X is a mammal.” In code:

engine.add_rule(
    name="mammal-rule",
    conditions=[{"pred": "has-hair", "args": ["?x"], "degVar": "?d"}],
    actions=[{
        "type": "add",
        "fact": {"pred": "is-mammal", "args": ["?x"], "deg": ["*", 0.95, "?d"]}
    }],
    priority=60,
)