PrimeStore, or how to store all primes ≤ 10¹² in at most 19.05 GiByte.

Storing primality for every integer below a bound $N$ is a solved problem in the boring sense: a bitmap, one bit per number, $N/8$ bytes. It works, and it wastes most of its space on numbers whose compositeness arithmetic hands you for free.

Today I’m open sourcing PrimeStore: a frozen container format and a security-hardened C23 decoder library that answers “is $k$ prime?” for any $k < N$, from a store between 2.0× and 6.1× smaller than that bitmap — 2.0× is the worst wheel on the format’s menu (the bare odd-only store), 6.1× the best, both asymptotic rates. Measured at $N = 10^{10}$ the largest wheel comes in at 6.25×, 7.3× with zstd layered over the payload. The decode semantics are written down as a Lean 4 specification and machine-checked against the shipped code.

What’s in the box

libprimestore.so — the decoder. A tiny API: ps_store_create decodes and validates a container, ps_lookup answers primality queries, ps_store_destroy releases the store. Stores are immutable after creation; ps_lookup is a pure function of the store and $k$, so it is thread-safe with no global state. The decoder’s only dynamic dependency is libc.
psencode — the standalone encoder that builds containers, with prime enumeration delegated to the excellent primesieve library. The decoder deliberately never links primesieve — a rule that is checked at link time, so an untrusted-container decoder has no sieve attack surface at all.
proof/ — the format’s frozen v0 specification in Lean 4 (with Mathlib): the wheel arithmetic, the bit layout, bitCount, rank by inclusion–exclusion, and the exact lookup semantics the C decoder implements.
saw/ and vst/ — machine-checked connections between that spec and the shipped code: Galois SAW proofs over the clang-emitted LLVM bitcode, and VST/CompCert source-level proofs of the C itself. More on this below, including the trust tables that say exactly what is proved and what is not.
The usual security apparatus — threat model, security requirements, fuzzing harnesses (libFuzzer and AFL++ builds), golden byte-identity fixtures, scan-build/annocheck gates, and a SECURITY.md with severity classes and response SLAs.

The mathematics behind it

The format is less designed than derived: take the first fact everyone learns about primes, formalize it, and keep pulling until the returns run out. Every headline number — the ratios, the container sizes, the entropy behavior — falls out of that pull, and each step of it is a machine-checked theorem in proof/, not a comment in a header.

Step zero: you only ever check the odds

Apart from 2 itself, no even number is prime. That single fact retires half of all integers from consideration: past 2, primality lives entirely on the odd numbers. A bitmap storing one bit per integer spends half its space on numbers whose parity already condemns them.

That observation is a wheel — the smallest one. Keep every second residue, store one bit per two integers, lose nothing. In the format’s menu this is $n = 1$: wheel $P = 2$, one surviving class, ratio 2.00×. The whole format is the answer to: what happens if you keep pulling this thread?

Residue classes: the honest way to say “skip”

To pull the thread you need the right language for “skipping”. Fix a modulus $m$. Every integer has a remainder $r \in \{0, 1, \dots, m-1\}$ mod $m$, and the $m$ residue classes partition the integers. A class can be dropped without loss exactly when some fixed prime divides all of its members — which happens iff $\gcd(r, m) > 1$: any prime dividing both $r$ and $m$ divides every $k \equiv r \pmod m$. The classes worth keeping — those coprime to $m$ — are the units (they are exactly the invertible elements of the ring $\mathbb{Z}/m\mathbb{Z}$, hence the name; Dirichlet’s theorem guarantees every unit class produces primes forever, so nothing else may be dropped).

The sieve lemma, stated in this language:

Every prime is either one of the primes dividing $m$ — a seed — or a unit. Every $k \ge m$ in a non-unit class has a seed as a proper divisor and is therefore composite.

Mod 2, this says “primes are odd”: class 0 dies, class 1 lives. Mod 6 ($m = 2 \cdot 3$): classes 1 and 5 live, and two-thirds of all integers are retired before a single bit is stored. Mod 12 ($2^2 \cdot 3$): classes 1, 5, 7, 11 live — the same survivors as mod 6, with the pattern pointlessly stretched to period 12. Hold that thought.

$\varphi$: counting what survives

How many classes survive a modulus $m$? The count is Euler’s totient: $\varphi(m) = \#\{\, r < m : \gcd(r, m) = 1 \,\}$. A wheel mod $m$ keeps a $\varphi(m)/m$ fraction of the integers, spends that many bits per integer, and therefore beats the plain bitmap by $m/\varphi(m)$. The decisive property is

$$ \frac{\varphi(m)}{m} \;=\; \prod_{p \,\mid\, m} \left(1 - \frac{1}{p}\right) $$

— only which primes divide $m$ appear; exponents never do. A wheel on 12 keeps exactly the same integers as a wheel on 6 and repeats its pattern twice as often: the squared factor bought nothing and doubled the period. So prime powers are dead weight, every candidate wheel is squarefree, and the only real choice left is which primes.

Which wheel? The primorial, Pareto-optimally

Adding a prime $p$ to a squarefree wheel multiplies the ratio by exactly $p/(p-1)$ and the period by $p$. Small primes are the efficient ones: 2 buys a factor of 2 outright, 3 buys another 1.5×, by the time you reach 23 a new prime buys 1.046×. Now the exchange argument that forces the choice: suppose your wheel contains the prime $q$ but skips a smaller $p$. Swap them. The ratio strictly improves ($p/(p-1) > q/(q-1)$, since $x/(x-1)$ shrinks as $x$ grows) and the period strictly shrinks ($\times\, p/q < 1$) — the swap dominates on both axes. A wheel that is not downward-closed in the primes is therefore dominated, full stop. The undominated wheels are exactly the primorials $P = p_n\# = \prod_{i \le n} p_i$.

Concretely: the wheel mod 10 (primes 2, 5) retires multiples of 2 and 5 — ratio 2.5×, period 10. The wheel mod 6 (primes 2, 3) retires multiples of 2 and 3 — ratio 3.0×, period 6. Smaller period, more removal; there is no axis on which mod 10 wins. Only the primorial menu is free of such dominated choices, and its ratios follow Mertens:

$$ \frac{P}{\varphi(P)} \;=\; \prod_{i=1}^{n} \frac{p_i}{p_i - 1} \;\sim\; e^{\gamma} \ln p_n \quad \text{(Mertens)} $$

— 2.00, 3.00, 3.75, 4.38, 4.81, 5.21, 5.54, 5.85, 6.11 for $n = 1 \dots 9$ (the asymptote $e^\gamma \ln 23 \approx 5.6$ still undershoots the finite product; Mertens converges slowly). Lean pins the coarse guarantee phi_le_half: every menu wheel at least halves the store.

Why not pull the thread forever? Because the costs compound faster than the gains: the ratio grows only like $\ln p_n$ while the period grows like $e^{p_n}$ — the logarithm of a primorial is the sum of $\ln p$ over its seeds, and the prime number theorem (stated below) evaluates that sum as $\sim p_n$. In practice: each rung doubles the inclusion–exclusion work per lookup below ($2^n$ terms), the free first block dominates the store until $N \gtrsim 100\,P$, and the format’s u32 arithmetic caps the wheel at $p_9\# = 223{,}092{,}870 < 2^{32}$ (the next rung, $p_{10}\# \approx 6.5 \times 10^9$, overflows). That is the $n = 9$ ceiling and the 6.11× asymptote.

The container, assembled

The derivation lands on the concrete format. The store keeps one bit per unit residue in $[P, N)$, packed MSB-first into a byte-aligned bitmap, under a 12-byte header — "PRM" magic, a version/n nibble pair, $N$ as uint64 little-endian. The stored region starts at $P$: block 0 and everything below it, seeds included, is answered at decode time by the seed table and sqrt-bounded trial division, at zero stored cost. The wheel parameter $n \in [1, 9]$ picks the rung; $N$ states the query contract.

Finding a bit: rank by inclusion–exclusion

Because only units are stored, $k$’s bit is not at position $k - P$. The frozen spec (Basic.lean, bitIndexOf) puts it at

$$ \operatorname{bitIndex}(k) \;=\; \big(\lfloor k / P \rfloor - 1\big) \cdot \varphi(P) \;+\; \operatorname{rank}_P(k \bmod P) $$

where $\operatorname{rank}_P(r) = \#\{\, u < r : \gcd(u, P) = 1 \,\}$ counts units strictly below the residue. While $P \le 510\,510$ ($n \le 7$) a flat per-residue table answers it in one hop; for $n = 8, 9$ the decoder instead evaluates the Möbius inversion over the seeds directly (RankIE.lean, rankOf_eq_inclusion_exclusion, proved for every $n \ge 1$):

$$ \operatorname{rank}_P(r) \;=\; \sum_{S \subseteq \{p_1, \dots, p_n\}} (-1)^{|S|} \left\lfloor \frac{r - 1}{\prod_{s \in S} s} \right\rfloor $$

— $2^n = 512$ terms at $n = 9$ (the empty subset contributes $r - 1$), pure integer floors, no table. This is also the identity whose $r = 0$ edge the SAW campaign caught: in $\mathbb{N}$, $r - 1$ truncates to 0; in u64 it wraps.

The size of a store, exactly

The same rank function counts the whole payload (Compression.lean, bitCount_eq): the stored region $[P, N)$ holds

$$ \big(\lfloor N / P \rfloor - 1\big) \cdot \varphi(P) \;+\; \operatorname{rank}_P(N \bmod P) $$

bits — full blocks plus the truncated tail — and the container is $12 + \lceil \text{bits} / 8 \rceil$ bytes with at most 7 padding bits (encode_size). Read backwards, that is the decoder’s validation contract: any buffer whose length is not exactly this is rejected before a payload bit is touched.

Evaluated at $n = 9$, $N = 10^{10}$, the closed form gives a 199,923,336-byte container — computable before a single bit is encoded, and enforced by the decoder on every input. Against the 1.25 GB plain bitmap that is 6.25×, sitting above the 6.11× asymptote exactly as the tail factor $N/(N-P) \approx 1.023$ (at $N = 44.8\,P$) predicts; at $n = 5$ the same $N$ measures 4.81×, already at its asymptote. Generating the $10^{10}$ container is honest work: psencode needs about 31 minutes, single-threaded.

Why the payload is noise: Dirichlet, the PNT, and the entropy floor

Two classical facts about how primes are distributed. The prime number theorem: the prime-counting function $\pi(x)$ — the number of primes below $x$ — satisfies

$$ \pi(x) \;\sim\; \frac{x}{\ln x} \qquad (x \to \infty) $$

so primes near $x$ occur with density $\approx 1/\ln x$. Dirichlet’s theorem on arithmetic progressions supplies the second half: those primes spread evenly across the $\varphi(P)$ unit residue classes. Combine the two and a unit slot near $x$ is set with probability $\approx (P/\varphi(P))/\ln x$. At $n = 9$, $N = 10^9$ that predicts a density of $\approx 0.32$ at the start of the stored range, falling to $\approx 0.29$ at the end — the measured block densities run 0.3118, 0.3036, 0.2986, 0.2958.

A Bernoulli coin of bias $p$ carries $h(p) = -p \log_2 p - (1-p) \log_2 (1-p)$ bits per toss; the payload’s mean bias gives $h(0.3034) = 0.886$ bits per bit, against a measured order-0 entropy of 0.885. The payload is, to measurement error, exactly what an iid coin at local density produces. That one sentence is the whole story of the entropy floor: generic compressors saturate because there is nothing behind the coin; block-position awareness can only harvest the slow $1/\ln x$ drift ($\le 0.18\%$); and the ceiling works out to $1/0.886 = 1.13\times$ — the study’s measured order-0 bound for that container is 1.129×, and the best generic tool reaches 1.122×. The math, the bound, and the measurement agree to the third digit.

The compressibility study measured this properly before I believed it: every wheel setting swept across $N$ up to $10^{10}$, five generic compressors against entropy bounds. The payload sits at 98.3–99.4% of its order-0 bound on every configuration measured, and the tools converge on each other — the signature of data that is already noise. Don’t design a cleverer coder; there is nothing left for it to find.

What a generic compressor can still collect is the global bit density: zstd on the payload lifts the $n = 9$, $N = 10^{10}$ store from 6.25× to 7.30× end-to-end — free at decode time, zero decoder state. That layer stays outside the container for now: every v0 container is a literal bitmap, though the spec reserves bit 63 of the bound field as a compression flag and an order-0 rANS codec sits in-tree as a non-v0 spike; the study sizes a zstd alternative at ~266 KB of decoder if a compressed variant ever ships.

The proof stack

I built PrimeStore because a compact primality store is genuinely useful — but a store is only as useful as your reason to trust its decoder, and a decoder that eats untrusted buffers is exactly where bugs hurt. So the library ships with the verification stack I want every binary-format decoder to have:

A frozen spec in Lean 4. proof/PrimeStore/Basic.lean defines the container byte-for-byte: what the header must contain, how many bits the payload must have, and what lookup answers for every query. The format is frozen — changing it means a version bump, not an edit.
Bitcode proofs with SAW. Each decoder function is symbolically executed with SAW over the clang-emitted LLVM IR and proved to implement its Lean spec term: get_bit, primorial, phi, is_unit, the rank table path, and more. Reproduce with make verify-saw.
Source-level proofs with VST. The same specs, hand-mirrored into Rocq, are proved against the real C translation units via CompCert’s Clight AST and VST semax: ps_primorial, ps_phi, and ps_gcd are PROVED over the full input domains, zero admits. Notably ps_gcd — the Euclid loop with a mathematical loop invariant — exceeds the SAW route, which had to assume the very contract VST went on to prove.

Two things make this more than proof-theater. First, the trust tables: both saw/README.md and vst/README.md list every function with an explicit verdict — PROVED, ASSUMED, or RESIDUAL — plus every residual assumption. What isn’t proved says so, in those words, next to what covers it empirically (golden fixtures, sanitizers, 1.5×10¹⁰ fuzzed bytes without a crash).

Second, the tooling found a real bug. SAW divergence-hunting turned up a wraparound in the inclusion–exclusion rank path: at residue $r = 0$ the C computed below − 1 in u64 and got 2⁶⁴−1, while the Lean identity truncates to 0. A wrong-number bug, not a memory-safety bug — exactly the class that ships silently when your spec lives in a comment instead of a proof assistant. It’s fixed, and the fix is pinned by proof and oracle tests.

Hardening

The build treats the compiler as part of the verification stack: C23 with -Wall -Wextra -Werror, full CET/CF protection, stack protector, _FORTIFY_SOURCE=3, thin LTO, full RELRO, noexecstack, annobin annotations, and mimalloc for allocations; debug builds add ASan/UBSan/LSan with use-after-scope. The threat model treats the container buffer as fully attacker-controlled — the decoder does no I/O, reads no environment, and validates the exact expected length before touching a payload bit.

How it was built

In the spirit of my earlier posts: PrimeStore was developed AI-assisted — and with that help, it took a single weekend to go from specification to fully tested implementation. The arc is git log, so I’ll just quote it: the v0 format was specified, proved, and frozen on Friday — the main correctness theorem, lookup (encode k) = decide (Nat.Prime) for every k < N, landed Friday afternoon, and the first decoder sources that evening. The golden byte-identity gates against the Lean fixtures, the libFuzzer campaign, and the first SAW bitcode proofs were green by Saturday night. Sunday went to the VST source-level rungs and the compression-flag amendment. Fifty commits, three days.

My role was operator: I wrote the specifications and the rulings; coding agents wrote the C, the Lean, the SAW scripts, and the Rocq proofs — and the honesty of the trust tables is part of what’s under test. The claims above are the checkable ones: every proof reruns from make verify-saw / make verify-vst, every measurement regenerates from the committed scripts, the dates above are the commit history, and the fuzz corpus policy is written down.

Try it

console

$ git clone https://gitlab.com/the_coh/prime_store
$ cd prime_store && make DEBUG=0 all && make test   # unit + golden byte-identity tests
$ build/linux-release/bin/psencode 9 1000000000 primes.prm

Linking the decoder is three calls and a pkg-config primestore: create a store from the container bytes, loop ps_lookup over your queries, destroy the store. The header primestore.h documents the exact contracts — including the one borrowed-buffer rule the API hangs on.

Status: v0.1.0, format v0 frozen, released under the Apache-2.0 license. Use it where you’d otherwise ship a plain prime bitmap: it answers the same queries from a store several times smaller, decodes untrusted input without third-party dependencies, and tells you — row by row — exactly how much of it is proven. Issues, review of the residual-assumption rows, and benchmark numbers from other machines are all welcome.

The next pull of the thread

The derivation above exposes v0’s one real design compromise: the format commits to a single $n$ for the whole store, at encode time. That commit is a guess about scale — the study’s crossover tables exist precisely because the best $n$ depends on $N$, and the tail factor punishes any store where $N \lesssim 100\,P$.

The sketch I like most for a v1 removes the guess: an adaptive wheel that grows as the store climbs. Rung boundaries are the primorials themselves, and the wheel steps up the moment the values being stored reach the current one — as soon as the encoder is laying down primes at or beyond $p_n\#$, the region it is writing already uses the $p_{n+1}\#$ wheel, monotonically, forever upward. Three things fall out:

No encoder guess, ever. The menu choice becomes positional instead of global: every region of the store rides a wheel sized for its own scale, encode streams (append forever, $N$ never needs to be known in advance), and the v0 empty-container pathology — $N \le P$ — cannot occur, because a rung that would store nothing simply never opens.
Every number $\ge 2$ is stored exactly once, seeds included. The seed primes live in lower rungs; trial division shrinks to answering 0 and
1. A lookup binary-searches the ~15 primorial boundaries, runs the unit test, ranks, and reads one bit — the same machinery, one rung-dispatch in front of it.
The store rides the top of the menu. Each rung spans at least twice the volume of everything below it combined, so the overall ratio tracks the top rung’s — and with u64 wheel arithmetic the menu extends from $n = 9$ to $n = 15$ ($p_{15}\# \approx 6.1 \times 10^{17} < 2^{63}$; $p_{16}\#$ overflows), lifting the asymptote from 6.11× to ≈ 7.2×.

The growth law stays what Mertens says it must be — logarithmic twice removed:

$$ r(N) \;\approx\; e^{\gamma} \ln p_m \;\sim\; e^{\gamma} \ln \ln N, \qquad p_m\# \;\le\; N \;<\; p_{m+1}\# $$

($\ln P \sim p_n$ by the PNT, so the largest primorial below $N$ has $p_m \approx \ln N$.) And the race from the derivation runs to its end: each rung buys $\approx 1/p$ in ratio and costs a doubling of the inclusion–exclusion work — $2^{15} = 32{,}768$ terms per rank at the ceiling — so the returns remain structural, not a matter of effort. The study’s other lever, the optional compressed payload riding the already-reserved header flag (zstd at ~266 KB of decoder, or the in-tree rANS spike), composes multiplicatively with all of it; its golden gate would degrade gracefully from byte-identity to decompress-then-compare.

None of this touches v0 — the format is frozen, and the house rule is that the Lean spec moves first, the C follows, and the trust tables grow a row per rung. That discipline is what makes the store worth shipping, not a demo of it. Step zero was “you only ever check the odds”. The adaptive wheel is what happens when you refuse to stop pulling.