Hailstones on a Mac Air

3n+1 through GMP, hand-rolled NEON, and a Metal GPU kernel on a MacBook Air. Two results ran opposite to what I predicted.

created 2026-05-03·updated 2026-07-08·explainer, explainer, design, writing

#rust #simd #gpu #collatz #performance #metal

The Collatz map:

T(n) = n/2      if n even
       3n + 1   if n odd

Iterate it on any starting n and, conjecturally, you land on 1. I wanted to know how fast a MacBook Air could run that iteration on numbers too big for a u64, and whether the GPU was worth reaching for. Two of the results ran opposite to what I predicted going in.

Attempt 1: the loop, rewritten three times

Started with rug (Rust bindings to GMP) and the most direct loop: bit-test, halve-or-triple, repeat. Rewrote it twice more, benchmarking each version with Criterion before touching the next line:

original: 1.87ms

optimized: 840µs

batch_shift: 60.6µs

Same GMP Integer, same algorithm — only the shape of the loop changed.

optimized stopped reallocating a fresh Integer every iteration: 2.2×. batch_shift added another 13.9× on top of that, 30.9× total, from one change: instead of looping n >>= 1 one bit at a time until n turns odd again, count the trailing zeros once and shift by that count:

n ↦ (3n + 1) ≫ v₂(3n + 1)

where v₂ is the 2-adic valuation — the trailing zero bit count. GMP exposes it as find_one(0). This one change outweighed everything I did to the multiply itself.

Attempt 2: does a narrower type beat GMP?

GMP is a general bignum — it multiplies two arbitrary-size numbers against each other. This loop never does that. It only does n × 3 + 1 (a known constant) and a shift by a count it just computed. Narrower problem, so I wrote CollatzInt: a Vec<u64> of little-endian limbs with exactly two operations, mul3_add1 and shr_assign.

Multiplying a limb by a fixed 3 makes the carry-out fully determined by that limb’s own value, no dependency on its neighbors:

carry(ℓ) = ⌊3ℓ / 2⁶⁴⌋ ∈ 2, for ℓ ∈ [0, 2⁶⁴)

which two threshold comparisons pin down exactly (ℓ ≥ ⌈2⁶⁴/3⌉, and ℓ ≥ ⌈2·2⁶⁴/3⌉). Every limb computes its own carry in parallel — a NEON vcgeq_u64 op across all limbs at once, phase one of mul3_add1. Phase two is a short scalar carry-propagate for the +1, bounded to {0, 1} since phase one already absorbed the ×3 overflow:

limbs[]

×3 (NEON, all limbs parallel)

carry detect (2 threshold cmps)

+1 scalar carry propagate

trailing_zeros()

shr_assign(tz)

One fused Collatz step. Phase 1 has no cross-limb dependency, so it runs in parallel; phase 2 is a short, branch-free carry chain.

Steps-only comparison (power computation excluded from both):

	time	vs. GMP
`gmp`	66.99µs	—
`custom_seq` (1 thread)	57.58µs	1.16×
`custom_par` (rayon, 8-way)	36.51µs	1.83×

Single-threaded and restricted to two operations, CollatzInt already beats GMP. Parallelism adds more, but there isn’t much work to spread across 8 threads at this scale.

Attempt 3: throwing it at the GPU

shaders/collatz.metal runs the same fused step as a fixed 4-limb (256-bit) struct, batched over 65,536 starting values at once. One-thread-per-value stalls here: Collatz sequences vary wildly in length, so a SIMD group idles waiting on its slowest lane. MetalCollatz::run_batch instead dispatches a fixed pool of 4,096 persistent threads that atomically claim indices from a shared counter as they finish — a thread that draws a short sequence immediately steals more work rather than idling.

metal_worker.rs — a small, persistent thread pool claiming work from a shared atomic counter, instead of one GPU thread per input value.

GPU speedup isn’t constant. On a uniform batch of nearby starting values (short, similar-length sequences), the GPU edges the CPU-parallel path by 1.05×. On a batch with real variance in descent length, at 3-limb width, the GPU reaches 19.4× over sequential CPU, 3.7× over CPU-parallel. Same kernel, same core loop. The GPU pays for itself once there’s enough divergent work for the persistent pool’s work-stealing to earn its keep — at small or uniform scale, dispatch overhead erases most of the advantage over a plain rayon loop.

Pre-sorting by magnitude did nothing. The intuitive fix for SIMD-group divergence is to sort the batch by expected sequence length before dispatch, so similar-difficulty values land in the same claimed block and no group waits on one straggler. I predicted this would help. Ran it both ways on the same divergent 3-limb batch:

unsorted: 3.327ms
sorted_by_magnitude: 3.319ms

0.2% — noise. The persistent-thread work-stealing pool was already the fix for divergence: threads claim work continuously as they finish rather than being locked to a static block for the whole dispatch, so sorting solves a problem the scheduler had already solved.

Open question

Does the 19.4× hold as the batch grows past 65,536, or does the persistent pool’s fixed size start working against it at some larger scale? Haven’t run that yet.