Bit-parallel sequence alignment kernels that Zig derives at compile time from the scoring scheme you give it.
Fast bit-vector aligners exist, but each one was worked out by hand for one cost model. Myers did it for unit edit costs in 1999, and BitPAl later designed a construction for linear-gap match, mismatch and gap weights that its generator instantiates per weight set. However, if your scheme falls outside what someone already derived, you are back to the scalar dynamic program. bitdp takes the recurrence itself and produces the bit-parallel program with comptime, so a new scheme costs a recompile instead of a paper.
Read the DP matrix one text column at a time, with the pattern rows packed into the bits of a 64-bit word. Down a column, each cell hands the next one its horizontal score difference. When differences are bounded, that makes the column a finite-state machine whose state is the difference, and for min-plus recurrences every step of that machine is monotone in its state.
Monotonicity is the key fact. Write the state in thermometer code (one bit per threshold, [d >= t]). Each threshold bit of the next state is then either a constant or a copy of one threshold bit of the previous state. A threshold that copies itself behaves like a carry chain with generate, propagate and kill positions, and a carry chain across a machine word is exactly what integer addition computes. A threshold that copies a different threshold is a shift and some bitwise logic. When the copies between different thresholds form no cycle, which holds for every linear-gap scheme, the whole column becomes a short cascade of additions and boolean operations.
Myers' algorithm comes out of this as the special case with one carry chain. Nothing about it is hard-coded here, and the derived kernel uses 14 word operations per column where Myers' uses 15.
Wide score ranges need one more idea. Most planes ask whether a difference of two thermometer-coded numbers clears a threshold, and on thermometer codes, addition is merging, since the merged sequence of two sorted bit strings is the thermometer code of their sum. A Batcher odd-even merging network (one OR and one AND per comparator) therefore delivers every threshold at once. The algebra also fixes a single cut-off, the highest substitution cost minus the gap cost, above which no level needs a carry chain. Levels below it keep their chains and the rest of the column comes out of two merges.
For small score ranges the compiler also tries exhaustive truth-table synthesis and keeps whichever program is shorter.
In addition, any integer substitution matrix works, not only match and mismatch, because the cost itself becomes one more thermometer-coded input and the same rules apply.
Research code, not yet stable. Word operations per column, as derived:
| Scheme (match, mismatch, gap) | Distinct differences | Additions | Word ops | Check |
|---|---|---|---|---|
| 0, 1, 1 (edit distance) | 3 | 1 | 14 | 10^6 random pairs, bit-exact |
| 0, 2, 1 (indel distance) | 2 | 1 | 6 | 2 x 10^5 pairs, bit-exact |
| 0, 1, 2 | 5 | 1 | 28 | 2 x 10^5 pairs, bit-exact |
| 0, 3, 2 | 5 | 3 | 57 | 2 x 10^5 pairs, bit-exact |
| -2, 3, 5 | 13 | 5 | 181 | 10^5 pairs, bit-exact |
| -3, 4, 6 | 16 | 7 | 253 | 10^5 pairs, bit-exact |
| -4, 5, 9 | 23 | 9 | 399 | 10^5 pairs, bit-exact |
| -4, 7, 11 | 27 | 11 | 509 | 10^5 pairs, bit-exact |
| DNA: transition 1, transversion 2, gap 2 | 5 | 2 | 51 | 10^5 pairs, bit-exact |
| BLOSUM62 costs, gap 4 | 20 | 15 | 799 | 10^5 protein pairs, bit-exact |
Rows five to eight are the weight sets that BitPAl (Loving, Hernandez and Benson, 2014) benchmarks, written as costs. A score scheme (M, I, G) becomes costs (-M, -I, -G), which gives the same optimal alignments. For (2, -3, -5), the one set whose operation counts the BitPAl paper reports, BitPAl needs 265 operations per 64-bit word and its packed variant 166. The derived kernel needs 181, fewer than the first and still more than the second. Operation counts include every AND, OR, XOR, NOT, shift and addition, and a + b + 1 counts as two.
All numbers are from one core of a Ryzen 7 8840U, with every tool's total cost checked against the others; bench/RESULTS.md has the full tables and bench/README.md the way to reproduce them.
The closest tool is BGSA, which runs Myers' and BitPAl's hand-derived kernels with one alignment per SIMD lane. In its own setting (all queries against all subjects of equal length) and at the same vector width, the derived kernels land close to it. For edit distance, bitdp reaches 0.73 to 0.93 of BGSA's speed with AVX2 and 0.90 to 1.20 with AVX-512. On BitPAl's (2, -3, -5) weights with AVX2, bitdp ties or beats BGSA's standard BitPAl (8.6 against 8.6 and 11.7 against 9.2 GCUPS) and trails its packed variant (8.6 against 13.0 and 11.7 against 14.0). Schemes BGSA cannot express still run at the same order of speed, for example 45.8 GCUPS with AVX2 for transition/transversion costs on 1 kbp sequences.
On independent pairs, where each lane has its own text, batched bitdp is faster than edlib, ksw2 and parasail's correct kernels on every DNA workload measured (for instance 35.2 against 21.7, 2.2 and 5.2 GCUPS for edit distance on 1 kbp pairs with AVX2). The same holds on real Illumina reads checked against the reference span a mapper placed them on. There the margin over edlib shrinks to 15 %, since those reads carry only 0.2 edits each. For BLOSUM62 it is several times slower than parasail, since fifteen cost classes make the column program too long to pay off.
On a real genome, 32 CRISPR-style 20 nt guides scanned along E. coli K-12 (4.6 Mbp) with transition/transversion costs return the same 78 sites as the plain DP in 0.30 s against 13.3 s. A 20 nt guide leaves most of a 64-bit word empty, so Packed puts three guides in each word and a spacer row after each guide to stop carries and shifts. That raises the scan from 12 to 27 GCUPS with AVX-512 and from 8 to 20 with AVX2. Hyyrö, Fredriksson and Navarro packed patterns this way for unit costs in 2005; here the same spacer rule works for every derived kernel.
parasail's striped global kernels, which would otherwise be its fastest, scored some pairs below their optimum (426 of 100 000 on one workload), so the comparison leaves them out.
- Global and search (semi-global) modes, linear gap costs only. With affine gaps the value carried down a column can shift both up and down between thresholds, the copy graph acquires cycles, and the cascade above no longer applies.
- At most 32 distinct score differences and 24 distinct substitution costs.
- The carry-chain part still grows quadratically with the spread between substitution costs. This is where BitPAl's packed variant stays ahead, and why protein matrices are slow.
Packedhandles search mode only, with patterns of at most 63 characters.- Bytes outside the scheme's alphabet (N, for instance) count as the costliest substitution.
- Deriving a scheme happens inside the Zig compiler. Small schemes take seconds, while BLOSUM62 takes about 15 seconds and 375 MB of compiler memory.
Requires Zig 0.16.0.
const bitdp = @import("bitdp");
const Edit = bitdp.Kernel(.{ .match = 0, .mismatch = 1, .gap = 1 });
// Patterns up to 62 characters, no allocation.
const d = Edit.distance("ACGTTGCA", "ACGTGCA"); // 1
// Any pattern length: prepare it once, align it against many texts.
var a = try Edit.Aligner.init(gpa, long_pattern);
defer a.deinit(gpa);
const d2 = a.distance(text);
// Many pairs, one per SIMD lane.
try Edit.distances(gpa, patterns, texts, out);
// Up to Edit.lanes patterns against the same text, one per lane.
var g = try Edit.Group.init(gpa, subjects);
defer g.deinit(gpa);
const costs = g.distances(query);
// Search mode: the pattern against its best-matching stretch of the text,
// or every end position with cost at most k.
const Find = bitdp.Kernel(.{ .sub = &bitdp.schemes.tsTvCost, .gap = 2, .mode = .search });
var guides = try Find.Group.init(gpa, guide_list);
defer guides.deinit(gpa);
try guides.scan(gpa, genome, 4, &hits);
// Short patterns, several per lane word. `count` says how many it took.
const pk = Find.Packed.init(guide_list);
try pk.scan(gpa, genome, 4, &hits); // Hit.lane is the pattern index
// Any substitution cost function over a declared alphabet.
const TsTv = bitdp.Kernel(bitdp.schemes.ts_tv);
const Blosum = bitdp.Kernel(bitdp.schemes.blosum62Linear(4));zig build test
zig build mve -Doptimize=ReleaseFast # prints derived programs, verifies, timesSee REFERENCES.md.
MIT