jschmied · Oct 12, 2026 · Oct 13, 2026 · Oct 13, 2026 · Oct 13, 2026 · Oct 13, 2026
diff --git a/docs/recipes/README.md b/docs/recipes/README.md
@@ -17,6 +17,7 @@ For agent-assisted contributions, read the [TensorFold contributor skill](SKILL.
 | Gemma 4 26B-A4B | [MLX, fused one-row decode](gemma-4.md) |
 | DeepSeek-V4-Flash | [MLX on a 256 GB Mac, DSpark and MTP drafts](deepseek-v4-flash.md) |
 | Qwen3.6-35B-A3B | [MLX with MTP drafts, one-GPU CUDA](qwen3.6-moe.md) |
+| Kolibri-1 | [Native CUDA, block FP8, research checkpoint](kolibri-1.md) |
 
 ## Capability floor
 

diff --git a/docs/recipes/kolibri-1.md b/docs/recipes/kolibri-1.md
@@ -0,0 +1,86 @@
+# Kolibri-1
+
+This page describes the native engine. `kolibri1` serves Aleph Alpha's Kolibri-1, a research checkpoint, on one
+NVIDIA GPU through the lane host. It is a draft family: plain decoding and copy drafts from the prompt's own text are
+in; the EAGLE-3 drafter is the next step and is not loaded here.
+
+## Checkpoint
+
+Tested revision: `Aleph-Alpha/Kolibri-1` at `e52eb4627d11516b0c01de49210ab5a4e4061444` (32 shards, 74 GB). The
+weights are block FP8: e4m3 codes in 128 x 128 blocks with fp32 `weight_scale_inv`, the routers, norms, embeddings and
+head in bf16. Publisher sha256 of `model.safetensors.index.json`:
+`c6769a34820c9057b404b2a5cae1ac24fd7b4a9ca08b84eb0a06e7b5bad3005e`.
+
+```bash
+tensorfold-native serve /path/to/Kolibri-1 --name kolibri-1 --context 32768 --parallel 8
+```
+
+The model has 50 layers, 2,560 wide, 48 query heads over 4 key heads of 128. Every fifth layer (4, 9, ..., 49) attends
+to every key without positions; the other 40 attend to a 513-key window with RoPE. Each layer routes by sigmoid with an
+expert bias to 6 of 384 experts, weighted by the unnormalised sigmoid, plus one ungated shared expert.
+
+## What the family adds
+
+It brings no kernels of its own. Every launch is a shared CUDA operator:
+
+- `cuda/fp8.zig` and `qlinear` for the attention projections (block FP8, `fp8g` lane tiles);
+- `cuda/experts.zig` for the routed and shared experts (`fp8g`, the shared expert stacked last at weight 1);
+- `qk_rows` (QK norm, RoPE, the cache write), `norms`, `embed_rows`, `moe_glue` (sigmoid top-k with bias),
+  `dense_rows` (router and head);
+- `window_attention` for prompt rows of every layer (the full layers with no window) and for sliding decode rows;
+- `tree_attention` for the full layers' decode rows over committed keys and each window's chain.
+
+The family supplies its config check, the tensor map (116,303 tensors, none left over), the packing into those
+layouts, the layer graph and the lane backend. One-row steps replay a CUDA graph a stream; wider windows and shared
+rounds run eager. Both give a row the same bits.
+
+## Serving
+
+- Prompts run in passes of up to 8,192 rows, as the Python engine ran them, one a step, so other streams' rounds run
+  between them (a full pass holds them for about 3.7 s). A row's bits do not depend on its chunk: 1,000-, 1,536- and
+  2,048-row cuts of a 5,000- and a 20,000-token prompt give the 8,192-row passes' logits byte for byte.
+- Concurrent streams share rounds of up to 64 rows. Each stream's caches are allocated at load, one per `--parallel`
+  stream: at a 32k window 1.29 GiB a stream, at 64k 1.92 GiB. Each sliding layer's ring holds a pass and its window
+  (8,768 keys).
+- The prompt cache keeps a stream's state at any row: full layers' keys up to it and each sliding layer's last 512.
+- Drafts: the lane core's copy drafts verify in windows of up to 16 rows. There is no draft head yet, so a request
+  for head drafts and `--drafter` are refused by name.
+- Refused at load: other quantizations, `norm_topk_prob`, other layer kinds, a config with a vision tower, shapes the
+  shared kernels do not take, `--vision`, `--precision` and a captured kernel set (`TENSORFOLD_CUDA_KERNELS`).
+
+## Checks
+
+```bash
+zig build -Doptimize=ReleaseFast -Dnvcc=/usr/local/cuda/bin/nvcc
+zig-out/bin/tf-cuda-test kolibri-reference MODEL IDS   # layer 0 and the head against float64 from the checkpoint
+zig-out/bin/tf-cuda-test kolibri-layers MODEL IDS      # the 513-key window edge, layer 4's full reach
+zig-out/bin/tf-cuda-test kolibri-full MODEL IDS 128    # steps, drafted, sampled, shared rounds, cuts, kept states
+zig-out/bin/tf-cuda-test kolibri-greedy MODEL PROMPTS N OUT   # greedy ids and top-2 margins a step
+zig-out/bin/tf-cuda-test kolibri-logits MODEL IDS OUT  # teacher-forced logits of 2,048 rows, for quality
+```
+
+`IDS` is a file of real-text token ids (20,000 or more for `kolibri-full`, 3,000 for `kolibri-layers`). Random ids give
+flat logits and say nothing about quality.
+
+The native rows are not byte-equal to the Python engine's: its glue was Triton, these are the shared CUDA operators,
+and their fp32 roundings differ by an ulp here and there. Kolibri-1's top-6-of-384 routing turns such differences into
+different experts on near-tied rows, so greedy replies part where two tokens are nearly tied. On 2,048 rows of real
+text, teacher-forced: next-token top-1 agrees on 92.3% of rows, KL 0.06; NLL 3.574 here against 3.578 for the
+Python engine, top-1 35.96% against 35.47%, mean entropy 2.354 against 2.358.
+
+## Speed
+
+GB10, greedy, served by `tensorfold-native` with copy drafts on:
+
+| | tok/s |
+| --- | ---: |
+| One stream, code / chat (`bench_openai.py`, 256 tokens) | 51.9 / 49.8 |
+| 4 streams, aggregate | 124.0 / 126.0 |
+| 8 streams, aggregate | 175.5 / 183.1 |
+
+Cold prompts (`bench_cache.py --require-cold --repeats 3`, `--prompt-cache-gib 0`, median of the cold requests):
+2,048 tokens 0.97 s, 8,192 tokens 3.72 s, 32,768 tokens 16.8 s, 64,000 tokens 37.8 s to the first token.
+
+## Next
+
+- The EAGLE-3 drafter through the lane core's drafter interface.
diff --git a/zig/build/cuda.zig b/zig/build/cuda.zig
@@ -45,8 +45,16 @@ const kernels = [_]Kernel{
     .{ .name = "embed_rows", .flags = &.{"-O3"} }, // ours: 4-bit affine and bf16 token rows (cuda/embed_rows.zig)
     .{ .name = "lfm2", .flags = glue }, // ours: LFM2's short convolution, QK norm and rotary, gates and label rows
     .{ .name = "vit", .flags = glue }, // ours: vision towers' LayerNorm, Linear epilogues, activations, unshuffle
+    .{ .name = "act_rows", .flags = &.{"-O3"} }, // ours: SwiGLU, output gates, group sums (cuda/act_rows.zig)
+    .{ .name = "gdn_prep", .flags = &.{"-O3"} }, // ours: Gated DeltaNet conv, norms, g and beta (cuda/gdn_prep.zig)
+    .{ .name = "qk_rows", .flags = &.{"-O3"} }, // ours: QK norm, RoPE and the cache write (cuda/qk_rows.zig)
     .{ .name = "qmmf", .flags = &.{"-O3"} }, // nvfp4/qmmf.cu's device code (FP8G, FP4), tensorfold_nvfp4_v3
     .{ .name = "prompt16", .flags = &.{"-O3"} }, // nvfp4/prompt.cu's FP4 and FP8 tile 4, tensorfold_nvfp4_prompt_v1
+    .{ .name = "tree_attention", .flags = &.{"-O3"} }, // ours: committed keys and paths (cuda/tree_attention.zig)
+    .{ .name = "block_attention", .flags = &.{"-O3"} }, // ours: a drafter's block over its context (cuda/block_attention.zig)
+    .{ .name = "window_attention", .flags = &.{"-O3"} }, // ours: sliding-window attention (cuda/window_attention.zig)
+    .{ .name = "moe_glue", .flags = &.{"-O3"} }, // ours: sigmoid top-k and the weighted combine (cuda/moe_glue.zig)
+    .{ .name = "dense_rows", .flags = &.{"-O3"} }, // ours: bf16 rows times weight rows into fp32 (cuda/dense_rows.zig)
     .{ .name = "modelopt_pack", .flags = &.{"-O3"} }, // ours: ModelOpt tensors into the lane, prompt and expert layouts
     .{ .name = "glm_kda", .flags = &.{ "-O3", "--fmad=false" } }, // GLM: the Python family's glm5_next/cuda/kda.cu
     .{ .name = "glm_glue", .flags = &.{"-O3"} }, // ours: the Python family's Triton glue by hand, its PTX op for op
@@ -154,6 +162,7 @@ pub fn targets(b: *std.Build, target: std.Build.ResolvedTarget, optimize: std.bu
     for (exlr8_build.denseModules(b, target, optimize, exlr8, exlr8_tests[0])) |m| runner.addImport(m.name, m.module);
     const glm_cuda = glmCuda(b, target, optimize, cuda, mods.core, mods.lanes, mods.heat);
     runner.addImport("glm_cuda", glm_cuda);
+    runner.addImport("kolibri1", kolibri1Module(b, target, optimize, cuda, mods.core, mods.lanes));
     const runner_exe = b.addExecutable(.{ .name = "tf-cuda-test", .root_module = runner });
     b.installArtifact(runner_exe);
     const controls = if (nvcc) |tool| fatbin(b, tool, version.?, .{ .name = "contract_controls", .flags = &.{"-O3"} }, sms) else null;
@@ -173,6 +182,11 @@ fn glmCuda(b: *std.Build, target: std.Build.ResolvedTarget, optimize: std.builti
     return m;
 }
 
+/// The Kolibri-1 CUDA family over the runtime, the core and the lane core.
+fn kolibri1Module(b: *std.Build, target: std.Build.ResolvedTarget, optimize: std.builtin.OptimizeMode, cuda: *std.Build.Module, core: *std.Build.Module, lanes: *std.Build.Module) *std.Build.Module {
+    return b.createModule(.{ .root_source_file = b.path("zig/src/families/kolibri1/cuda.zig"), .target = target, .optimize = optimize, .link_libc = true, .imports = &.{ .{ .name = "cuda", .module = cuda }, .{ .name = "core", .module = core }, .{ .name = "lanes", .module = lanes } } });
+}
+
 /// LFM2 (d1's decision model) on CUDA: config, weights, kernels and the forward that scores labels.
 fn lfm2Module(b: *std.Build, target: std.Build.ResolvedTarget, optimize: std.builtin.OptimizeMode, cuda: *std.Build.Module, core: *std.Build.Module) *std.Build.Module {
     return b.createModule(.{ .root_source_file = b.path("zig/src/families/lfm2/lfm2.zig"), .target = target, .optimize = optimize, .link_libc = true, .imports = &.{ .{ .name = "cuda", .module = cuda }, .{ .name = "core", .module = core } } });
@@ -199,6 +213,7 @@ fn engines(b: *std.Build, target: std.Build.ResolvedTarget, optimize: std.builti
             .{ .name = "lfm2", .module = lfm2Module(b, target, optimize, cuda, core) },
             .{ .name = "core", .module = core },
             .{ .name = "glm_cuda", .module = glmCuda(b, target, optimize, cuda, core, lanes, heat) },
+            .{ .name = "kolibri1", .module = kolibri1Module(b, target, optimize, cuda, core, lanes) },
         },
     });
     return .{ .api = api, .engines = mod };
@@ -269,9 +284,14 @@ pub fn hostTests(b: *std.Build, draft_ids: *std.Build.Module, build_options: *st
     }
     // the tokenizer's own tests (memo, steps, regex, unicode) run here: no other test binary reaches them
     for ([_]*std.Build.Module{ cuda, mods.core, mods.lanes, mods.nemotron, lfm2Module(b, host, .debug, cuda, mods.core), mods.radio, mods.heat, mods.tokenizer, native, stagger(b, host, .debug) }) |m| step.dependOn(&b.addRunArtifact(b.addTest(.{ .root_module = m })).step);
+    // the contract's attention reference (tests/cuda/contract_attn.zig): its allocations under failure
+    const attn_ref = b.createModule(.{ .root_source_file = b.path("zig/tests/cuda/contract_attn.zig"), .target = host, .optimize = .debug, .imports = &.{.{ .name = "core", .module = mods.core }} });
+    step.dependOn(&b.addRunArtifact(b.addTest(.{ .root_module = attn_ref })).step);
     // the GLM-5.3-Flash CUDA family's host side
     const glm = glmCuda(b, host, .debug, cuda, mods.core, mods.lanes, mods.heat);
     step.dependOn(&b.addRunArtifact(b.addTest(.{ .root_module = glm })).step);
+    // the Kolibri-1 CUDA family's host side
+    step.dependOn(&b.addRunArtifact(b.addTest(.{ .root_module = kolibri1Module(b, host, .debug, cuda, mods.core, mods.lanes) })).step);
     // the shares a GLM-5.3-Flash rank takes (cluster/canon.zig against the Python family's)
     const shares = b.createModule(.{ .root_source_file = b.path("zig/src/cluster/canon_glm_test.zig"), .target = host, .optimize = .debug });
     step.dependOn(&b.addRunArtifact(b.addTest(.{ .root_module = shares })).step);

diff --git a/zig/kernels/cuda/act_rows.cu b/zig/kernels/cuda/act_rows.cu
@@ -0,0 +1,143 @@
+// Activation rows for every CUDA family (float64 references: core/act_ref.zig): SwiGLU, output gates, group sums.
+
+// A thread 8 values with 16-byte loads; a row's bits depend on that row alone.
+
+#include <cuda_bf16.h>
+#include <stdint.h>
+
+namespace tf_act {
+
+constexpr int NT = 256;
+
+__device__ __forceinline__ void unpack8(const uint4& u, float* f) {
+    const __nv_bfloat162* h = reinterpret_cast<const __nv_bfloat162*>(&u);
+#pragma unroll
+    for (int i = 0; i < 4; ++i) {
+        const float2 t = __bfloat1622float2(h[i]);
+        f[2 * i] = t.x;
+        f[2 * i + 1] = t.y;
+    }
+}
+
+__device__ __forceinline__ uint4 pack8(const float* f) {
+    uint4 u;
+    __nv_bfloat162* h = reinterpret_cast<__nv_bfloat162*>(&u);
+#pragma unroll
+    for (int i = 0; i < 4; ++i) h[i] = __floats2bfloat162_rn(f[2 * i], f[2 * i + 1]);
+    return u;
+}
+
+// Exactly rounded division on purpose: the same bits on every GPU, where div.full and ex2.approx are not.
+__device__ __forceinline__ float sigmoid(float x) { return __fdiv_rn(1.f, __fadd_rn(1.f, expf(-x))); }
+
+// A thread's 8 stored values in order, then the group's 8 lanes in a fixed tree (every lane of the warp calls).
+__device__ __forceinline__ float group_sum8(const uint4& packed, bool live) {
+    float f[8];
+    unpack8(packed, f);
+    float s = 0.f;
+    if (live)
+#pragma unroll
+        for (int i = 0; i < 8; ++i) s = __fadd_rn(s, f[i]);
+    s = __fadd_rn(s, __shfl_xor_sync(0xFFFFFFFFu, s, 4));
+    s = __fadd_rn(s, __shfl_xor_sync(0xFFFFFFFFu, s, 2));
+    return __fadd_rn(s, __shfl_xor_sync(0xFFFFFFFFu, s, 1));
+}
+
+}  // namespace tf_act
+
+using namespace tf_act;
+
+// out[r, i] = bf16(silu(gate[r * g_ld + i]) * up[r * u_ld + i]) for i < n (n % 64 == 0); xs (or null) out's sums.
+extern "C" __global__ void __launch_bounds__(NT) tf_swiglu_rows(const __nv_bfloat16* __restrict__ gate, long long g_ld,
+                                                                const __nv_bfloat16* __restrict__ up, long long u_ld,
+                                                                __nv_bfloat16* __restrict__ out, float* __restrict__ xs,
+                                                                int n) {
+    const int64_t r = blockIdx.x;
+    const int c = blockIdx.y * blockDim.x + threadIdx.x; // 8 values
+    const bool live = 8 * c < n;
+    float g[8], u[8], o[8];
+    if (live) {
+        unpack8(*reinterpret_cast<const uint4*>(gate + r * g_ld + 8 * c), g);
+        unpack8(*reinterpret_cast<const uint4*>(up + r * u_ld + 8 * c), u);
+    }
+#pragma unroll
+    for (int i = 0; i < 8; ++i) o[i] = live ? __fmul_rn(__fmul_rn(g[i], sigmoid(g[i])), u[i]) : 0.f;
+    const uint4 op = pack8(o);
+    if (live) *reinterpret_cast<uint4*>(out + r * n + 8 * c) = op;
+    if (xs != nullptr) {
+        const float s = group_sum8(op, live);
+        if (live && (c & 7) == 0) xs[r * (n / 64) + c / 8] = s;
+    }
+}
+
+// out[r, h * d + e] = bf16(o[r, h * d + e] * sigmoid(g[r * g_ld + h * g_hs + e])) (d % 64 == 0); xs (or null) its sums.
+extern "C" __global__ void __launch_bounds__(NT) tf_gate_rows(const __nv_bfloat16* __restrict__ o,
+                                                              const __nv_bfloat16* __restrict__ g, long long g_ld,
+                                                              long long g_hs, __nv_bfloat16* __restrict__ out,
+                                                              float* __restrict__ xs, int heads, int d) {
+    const int64_t r = blockIdx.x;
+    const int c = blockIdx.y * blockDim.x + threadIdx.x; // 8 values
+    const int n = heads * d;
+    const bool live = 8 * c < n;
+    const int h = 8 * c / d, e = 8 * c % d;
+    float ov[8], gv[8], y[8];
+    if (live) {
+        unpack8(*reinterpret_cast<const uint4*>(o + r * n + 8 * c), ov);
+        unpack8(*reinterpret_cast<const uint4*>(g + r * g_ld + h * g_hs + e), gv);
+    }
+#pragma unroll
+    for (int i = 0; i < 8; ++i) y[i] = live ? __fmul_rn(ov[i], sigmoid(gv[i])) : 0.f;
+    const uint4 yp = pack8(y);
+    if (live) *reinterpret_cast<uint4*>(out + r * n + 8 * c) = yp;
+    if (xs != nullptr) {
+        const float s = group_sum8(yp, live);
+        if (live && (c & 7) == 0) xs[r * (n / 64) + c / 8] = s;
+    }
+}
+
+// xs[r, j] = the sum of x[r * ld + 64 j ..][0 .. 64] (k % 64 == 0), in group_sum8's order.
+extern "C" __global__ void __launch_bounds__(NT) tf_group_sums_rows(const __nv_bfloat16* __restrict__ x, long long ld,
+                                                                    float* __restrict__ xs, int k) {
+    const int64_t r = blockIdx.x;
+    const int c = blockIdx.y * blockDim.x + threadIdx.x;
+    const bool live = 8 * c < k;
+    const uint4 v = live ? *reinterpret_cast<const uint4*>(x + r * ld + 8 * c) : make_uint4(0, 0, 0, 0);
+    const float s = group_sum8(v, live);
+    if (live && (c & 7) == 0) xs[r * (k / 64) + c / 8] = s;
+}
+
+// DFlash2's two-tap conv: out = bf16(x k0 + x[r - 1] k1) [+ res], no r - 1 at a seg block's start, k = bf16(b + dyn).
+extern "C" __global__ void __launch_bounds__(NT) tf_dconv_rows(const __nv_bfloat16* __restrict__ x,
+                                                               const __nv_bfloat16* __restrict__ dyn,
+                                                               const __nv_bfloat16* __restrict__ base,
+                                                               const __nv_bfloat16* __restrict__ res,
+                                                               __nv_bfloat16* __restrict__ out, int d, int gs, int branch,
+                                                               int seg) {
+    const int64_t r = blockIdx.x;
+    const int c = blockIdx.y * blockDim.x + threadIdx.x;
+    if (8 * c >= d) return;
+    const int g = d / gs, grp = 8 * c / gs;
+    float xv[8], pv[8], b0[8], b1[8], y[8];
+    unpack8(*reinterpret_cast<const uint4*>(x + r * d + 8 * c), xv);
+    if (r % seg > 0) unpack8(*reinterpret_cast<const uint4*>(x + (r - 1) * d + 8 * c), pv);
+    else
+#pragma unroll
+        for (int i = 0; i < 8; ++i) pv[i] = 0.f;
+    unpack8(*reinterpret_cast<const uint4*>(base + (int64_t)(branch * 2) * d + 8 * c), b0);
+    unpack8(*reinterpret_cast<const uint4*>(base + (int64_t)(branch * 2 + 1) * d + 8 * c), b1);
+    const float d0 = __bfloat162float(dyn[((r * 2 + branch) * 2) * g + grp]);
+    const float d1 = __bfloat162float(dyn[((r * 2 + branch) * 2 + 1) * g + grp]);
+#pragma unroll
+    for (int i = 0; i < 8; ++i) {
+        const float k0 = __bfloat162float(__float2bfloat16_rn(__fadd_rn(b0[i], d0)));
+        const float k1 = __bfloat162float(__float2bfloat16_rn(__fadd_rn(b1[i], d1)));
+        y[i] = __fmaf_rn(pv[i], k1, __fmul_rn(xv[i], k0));
+    }
+    if (res != nullptr) {
+        float rv[8];
+        unpack8(*reinterpret_cast<const uint4*>(res + r * d + 8 * c), rv);
+#pragma unroll
+        for (int i = 0; i < 8; ++i) y[i] = __fadd_rn(rv[i], __bfloat162float(__float2bfloat16_rn(y[i])));
+    }
+    *reinterpret_cast<uint4*>(out + r * d + 8 * c) = pack8(y);
+}