Collatz search

This program has no render pass, so nothing is ever drawn here. Press Play and it reports numbers instead.

PlayDownload PNG

Every other page in this series draws something. This one draws nothing, and that is the whole point of it. Press Play: the frame above stays empty, because the program contains no render pass and never touches the canvas. Four numbers appear over it instead, and those numbers were computed on your GPU.

The rule being searched is Collatz’s. Take a whole number; if it is even, halve it; if it is odd, triple it and add one; repeat. Every number anyone has ever tried eventually reaches 1, and nobody has proved that every number must. This program takes the 131072 starting numbers below 2^17, counts how many steps each one needs to reach 1, and reports the number that took the longest, how long it took, and the total across all of them.

It is here because it reframes what a pngine file is. SJON, the S-expression format pngine compiles, is usually read as a way to describe a picture: each form is one WebGPU resource, and the shaders inside the forms are WGSL, WebGPU’s shading language. But nothing in the format requires a picture. What a compiled file actually contains is a description of GPU work, stored as PNGB bytecode in the PNG’s pNGb chunk alongside the executor, the small WebAssembly interpreter that turns that bytecode into GPU commands. Here the work happens to produce a number. The file you can download under the player is a picture of nothing that computes something: an image that is also a kernel.

examples/samples/30_collatz_search.sjon

; Collatz search. A program with nothing to draw: two compute passes and a
; buffer copy, no render pass and no canvas. The GPU runs the Collatz sequence
; for every starting number below 131072, a second pass reduces the 131072
; results to a winner and a sum inside one workgroup, and a `(queue …)` copies
; that result buffer into a `map-read` buffer, which is host memory the GPU may
; only copy into. The browser runtime maps it after each submit and hands the
; contents to the page through the `onQueryResult` option:
;
;   const p = await pngine(url, { canvas, onQueryResult: (lanes) => {
;     const [bound, winner, steps, total] = lanes;
;     hud.textContent = `${winner} takes ${steps} steps; mean ${total / bound}`;
;   }});
;
; It prints 106239 / 353 / 110.30, which is checkable against the literature.
; A payload cannot consume its own readback, so what to do with the answer is
; the host page's job. The native `--frame` renderer runs both passes and the
; copy but maps nothing, so it writes a fully transparent PNG: the picture is
; empty because the program never had a picture.
;
; SEARCH_MAX stops at 131072 because WGSL has no 64-bit integer type: 159487's
; sequence climbs past 2^32 and would wrap silently in `u32`. Every start below
; 131072 peaks at 2482111348, which still fits.
;
; The reduction packs each candidate into one `u32` as `(steps << 20) | n`, so
; a plain `max` finds the longest sequence and carries the number that produced
; it. n needs 17 bits and steps needs 9, so both fit either side of the split.
;
; SEARCH_MAX sizes the buffer only; `(define …)` values never reach WGSL, so
; each shader restates the bound as a local `const`. Keep them in step.

(define :name SEARCH_MAX :value 131072)

; steps[n] is the number of Collatz steps n takes to reach 1.
(buffer :name steps :size (* SEARCH_MAX 4) :usage [storage])

; Four 64-bit lanes: search bound, winning number, its step count, total steps.
(buffer :name result :size 32 :usage [storage copy-src])

; The readback: `map-read` pairs with `copy-dst` and nothing else.
(buffer :name readback :size 32 :usage [map-read copy-dst])

(shader-module :name climbModule :code """
const SEARCH_MAX: u32 = 131072u;

@group(0) @binding(0) var<storage, read_write> steps: array<u32>;

// One thread per starting number. 0 and 1 are defined to take no steps.
@compute @workgroup_size(64)
fn main(@builtin(global_invocation_id) gid: vec3u) {
  let n = gid.x;
  if (n >= SEARCH_MAX) { return; }

  var v = n;
  var count = 0u;
  // The longest sequence under SEARCH_MAX is 353 steps, so 1024 is a ceiling
  // no start reaches. It is here so a thread cannot loop forever on a GPU
  // with no way to interrupt it.
  for (var i = 0u; i < 1024u; i += 1u) {
    if (v <= 1u) { break; }
    if ((v & 1u) == 0u) { v = v >> 1u; } else { v = 3u * v + 1u; }
    count += 1u;
  }
  steps[n] = count;
}
""")

(shader-module :name foldModule :code """
const SEARCH_MAX: u32 = 131072u;
const THREADS: u32 = 256u;

@group(0) @binding(0) var<storage, read> steps: array<u32>;
@group(0) @binding(1) var<storage, read_write> result: array<vec2u>;

var<workgroup> best: array<u32, 256>;
var<workgroup> total: array<u32, 256>;

// A single workgroup reduces the whole array: each thread strides across it,
// then the 256 partial answers collapse pairwise in log2(256) = 8 rounds.
@compute @workgroup_size(256)
fn main(@builtin(local_invocation_id) lid: vec3u) {
  var b = 0u;
  var t = 0u;
  for (var n = lid.x; n < SEARCH_MAX; n += THREADS) {
    let s = steps[n];
    t += s;
    b = max(b, (s << 20u) | n);
  }
  best[lid.x] = b;
  total[lid.x] = t;
  workgroupBarrier();

  for (var stride = THREADS / 2u; stride > 0u; stride = stride >> 1u) {
    if (lid.x < stride) {
      best[lid.x] = max(best[lid.x], best[lid.x + stride]);
      total[lid.x] = total[lid.x] + total[lid.x + stride];
    }
    workgroupBarrier();
  }

  // The readback channel is 64 bits wide per lane, so each answer goes in the
  // low half of a vec2u and the high half stays zero.
  if (lid.x == 0u) {
    result[0] = vec2u(SEARCH_MAX, 0u);
    result[1] = vec2u(best[0] & 0xFFFFFu, 0u);
    result[2] = vec2u(best[0] >> 20u, 0u);
    result[3] = vec2u(total[0], 0u);
  }
}
""")

(compute-pipeline :name climbPipeline :layout auto
  (compute :module climbModule :entry main))

(compute-pipeline :name foldPipeline :layout auto
  (compute :module foldModule :entry main))

(bind-group :name climbBindGroup :layout climbPipeline :group 0
  (entry :binding 0 :buffer steps))

(bind-group :name foldBindGroup :layout foldPipeline :group 0
  (entry :binding 0 :buffer steps)
  (entry :binding 1 :buffer result))

(compute-pass :name climb :pipeline climbPipeline :bind-groups [climbBindGroup]
  (dispatch :workgroups [(ceil (/ SEARCH_MAX 64))]))

(compute-pass :name fold :pipeline foldPipeline :bind-groups [foldBindGroup]
  (dispatch :workgroups [1]))

(queue :name readOut
  (copy-buffer-to-buffer :source result :source-offset 0
    :destination readback :destination-offset 0 :size 32))

(frame :name main :perform [climb fold readOut])

Show all 134 lines

examples/samples/30_collatz_search.sjon in the pngine repository.

How it works

A frame with no render pass

(frame :name main :perform [climb fold readOut]) is the whole program’s schedule, and it names two compute passes and one queue:

(frame :name main :perform [climb fold readOut])

There is no (render-pass …) form anywhere in the document, no (color-attachment …), and no mention of context-current-texture, the canvas texture the drawing samples render into. Nothing in the compiler requires one. The canvas the player hands the runtime is simply never written.

Pass one: one thread per number

climbModule is a compute entry point declared @compute @workgroup_size(64). Its @builtin(global_invocation_id) gives each thread a unique index across the whole dispatch, and this shader reads that index as the starting number itself: thread 106239 is the one that runs the sequence for 106239. It walks the sequence in a for loop and writes the step count into steps[n].

I cap the loop at 1024 iterations rather than looping until the value reaches 1. The longest sequence in this range is 353 steps, so the cap is never reached; it is there because a GPU has no way to interrupt a thread that never finishes, and an unbounded loop in a shader is a hang, not an error message.

The pass dispatches (dispatch :workgroups [(ceil (/ SEARCH_MAX 64))]), which is 2048 workgroups of 64 threads: exactly 131072 threads, one per number. (ceil (/ SEARCH_MAX 64)) is an SJON expression over the (define :name SEARCH_MAX :value 131072) at the top of the file, evaluated when the document is compiled. WebGPU guarantees at least 65535 workgroups per dimension, so 2048 is comfortable.

Why the search stops at 131072

WGSL has no 64-bit integer type, so the sequence runs in u32, and 3 * v + 1 can carry a value far above the number it started from. The bound is not arbitrary: 159487 is the first starting number whose sequence climbs past 2^32. Every start below 131072 stays under 2482111348, which is the highest value reached in this whole range (from 113383), and that still fits.

One number past the bound and the arithmetic wraps. It does not fault, it does not warn, and pngine validate cannot see it: the shader is well formed and the answer is simply wrong. If you raise SEARCH_MAX, this is the reason to stop at 159486.

Pass two: one workgroup reduces everything

foldModule collapses 131072 results into four. It runs as a single workgroup of 256 threads and uses the workgroup address space, memory shared by the threads of one workgroup and gone when it finishes:

var<workgroup> best: array<u32, 256>;
var<workgroup> total: array<u32, 256>;

That is 2 KiB, against a limit WebGPU guarantees to be at least 16 KiB. Each thread first strides through the array 256 elements apart, accumulating its own running maximum and sum, so all 131072 values are read exactly once. Then the 256 partial answers collapse pairwise: thread i folds in thread i + stride, the stride halves, and after eight rounds thread 0 holds the answer for the whole array.

Between every round sits a workgroupBarrier(), which no thread passes until all of them arrive, so nobody reads a slot the round before was still writing. The barrier is placed after the if (lid.x < stride) block rather than inside it, because a barrier has to be reached by every thread in the workgroup: WGSL’s uniformity analysis rejects a barrier that only some threads would run.

Packing two answers into one max

A reduction like this needs an operator that can be applied in any grouping, which max is. “Which number produced the maximum” is not: two threads comparing step counts have no way to carry the number along. The shader gets both out of one max by shifting them into a single u32:

b = max(b, (s << 20u) | n);

one packed u32

bit    31                   20 19                            0
       +----------------------+------------------------------+
       |     step count s     |      starting number n       |
       +----------------------+------------------------------+
        12 bits (353 needs 9)   20 bits (131071 needs 17)

Comparing the packed values compares step counts first and carries the winning number for free; the split works because the two quantities are small in the right way, each fitting its side of bit 20 with room to spare. Thread 0 unpacks the winner with best[0] & 0xFFFFFu and best[0] >> 20u.

Getting the answer out: copy, then map

Three buffers carry the work, and their usage flags are what make the readback possible:

(define :name SEARCH_MAX :value 131072)
(buffer :name steps :size (* SEARCH_MAX 4) :usage [storage])
(buffer :name result :size 32 :usage [storage copy-src])
(buffer :name readback :size 32 :usage [map-read copy-dst])

steps is half a megabyte of intermediate state that never leaves the GPU. result is the four-slot answer, STORAGE so a shader can write it and COPY_SRC so it can be copied out of. readback is the one that crosses back to the CPU.

MAP_READ is a usage with almost no company: WebGPU allows it to pair with COPY_DST and nothing else, and pngine rejects any other combination at compile time. That restriction is the shape of the feature. A mappable buffer is host memory the GPU may only copy into or out of, never memory a shader binds, so getting a result to the CPU is always two steps: compute into a storage buffer, then copy that into a mappable one.

The copy is an ordinary queue operation:

(queue :name readOut
  (copy-buffer-to-buffer :source result :source-offset 0
    :destination readback :destination-offset 0 :size 32))

which compiles to copyBufferToBuffer. After each submit the runtime notices that a MAP_READ buffer was written, maps it, and hands the contents to the page through the onQueryResult option:

const p = await pngine(url, { canvas, onQueryResult: (lanes) => {
  const [bound, winner, steps, total] = lanes;
  hud.textContent = `${winner} takes ${steps} steps; mean ${total / bound}`;
}});

That channel was built for timestamp and occlusion queries, whose results are 64-bit, so it reads the mapped range as 64-bit lanes. This program has no queries, but it inherits the shape: each answer is written as vec2u(value, 0u), the low half of a lane, with the high half left zero. Four lanes, 32 bytes, four numbers.

A payload cannot consume its own readback. Deciding what to do with the answer is the host page’s job, which is why the numbers above are drawn by this page in HTML and not by the program.

What the four numbers say

131072, 106239, 353, 14457822. The search covered 131072 starting numbers; the longest sequence in that range begins at 106239 and takes 353 steps to reach 1; and the 131072 sequences together take 14457822 steps, an average of about 110. Both of the middle numbers are checkable against the literature: 106239 is a record holder, meaning no smaller starting number takes longer.

The poster is the output

The still image behind the player, and the empty card in the gallery, is not a placeholder. It is what pngine <file> --frame produced: a 640 by 400 frame of a program with nothing to render, every byte of it zero, which is why the area reads as blank in either theme. The same file still carries the compiled bytecode and the executor in its pNGb chunk, so it runs. Download it and it is a fully transparent PNG that searches Collatz sequences.

Things to try

• Raise SEARCH_MAX. It appears three times on purpose: once as the (define …) that sizes the buffer and the dispatch, and once as a const in each shader, because a (define …) value feeds SJON expressions and never reaches WGSL text. Stop at 159486, for the reason above.
• Add a fifth lane. Widen both 32-byte buffers to 40, write result[4] = vec2u(something, 0u), and add a label to the page. The readback is just bytes.
• Delete the (queue …) from :perform. Both compute passes still run, the GPU still does all the work, and nothing ever comes back.

Related

• Game of Life and Particle fountain are the other direction: simulations that live on the GPU and never send anything back.
• The readback channel is onQueryResult in the JavaScript API; the copy is a (queue …) form and the buffers are (buffer …) forms.