ResumeWritingProjects
Writing

Modernizing SFMT-19937 for Zig

Moving SFMT-19937 to Zig vector types made it slower. The vector types were not the problem, and the Debug benchmark reporting it could not have said so.

View on GitHub

The first version of SFMT-19937 in zig-prng was correct, and it still thought like C: a 128-bit word carried as a struct of four u32 values, loops walked with manual indices, and scalar helper functions imitating operations that the algorithm expects the CPU to do as SIMD. Replacing that state with @Vector(4, u32) made it slower. That is the wrong result in a useful way, because the vector types were not the problem and the Debug benchmark reporting it could not have said so.

SFMT is not a scalar generator with a SIMD paint job. The name is literal: SIMD-oriented Fast Mersenne Twister. If a modern implementation makes the vector path slower, something specific has gone wrong, and here it was one operation out of the recurrence still routing through scalar arithmetic.

ReleaseFast throughput

The takeoff only appears when the optimizer can see the SIMD shape.

5.0x

Struct baselineOriginal C-shaped state representation332 M/s
@Vector rewriteByte-lane shifts expressed with @shuffle1,618 M/s
Documented finalCorrectness held against SFMT 1.5.1 vectors1,675 M/s

Debug mode told the wrong story: the vector rewrite looked slower because the conversions and helper calls were still visible to the compiler.

Baseline82 M/sNaive vector48 M/sPointer rollback57 M/s@shuffle70 M/s
One billion generated u32 values per run. The last value stayed fixed at 0x7c7d388d.

The invariant came first

Before touching the implementation I added a benchmark that generates one billion u32 values and reports throughput. The first measurement was boring, which is what a baseline should be: 82 M/s in Debug, reference vectors passing.

Correctness was the number that mattered, and two tests pin it: they read SFMT 1.5.1’s own SFMT.19937.out.txt, load the expected stream for init_gen_rand and init_by_array, and compare it word for word, so a refactor was free to move code around as long as the sequence came out the same. The benchmark keeps a cruder guardrail of its own: it prints the last of the billion values, 0x7c7d388d.

Two passes, and the number went the wrong way

The first pass swapped C-style while loops for Zig range for loops wherever the iteration count was fixed. The benchmark did not move. That was expected, since the bottleneck lay somewhere else entirely, and the change bought readability instead of speed: the recurrence now reads like Zig rather than a translation unit with different punctuation.

I left initByArray alone, because its loops carry interdependent counters and % size modular arithmetic, which the explicit while form states more clearly than a range would.

The second pass replaced the 128-bit word struct with @Vector(4, u32), which on paper was the right move, since the state finally had a type that matched the algorithm. Debug throughput dropped to 48 M/s. My first guess was extra copying from passing 128-bit vectors by value, and moving the recursion helper back to pointer parameters recovered some of it at 57 M/s, still under the struct baseline. The signature was not it either.

The byte-lane shift was the leak

The recurrence uses two kinds of shift, and they sit one line apart:

inline fn doRecursion(r: *W128, a: *const W128, b: *const W128, c: *const W128, d: *const W128) void {
    const x = lshift128(a.*);
    const y = rshift128(c.*);
    const sr1_splat: @Vector(4, u5) = @splat(SR1);
    const sl1_splat: @Vector(4, u5) = @splat(SL1);
    r.* = a.* ^ x ^ ((b.* >> sr1_splat) & MSK) ^ y ^ (d.* << sl1_splat);
}

SR1 and SL1 are element shifts, 11 and 18 bits, and they are plain >> and << on the vector. SR2 and SL2 are byte shifts of the whole 128-bit register. Those need lshift128 and rshift128, and they are not the same operation, however alike they look in the recurrence.

The Zig code was simulating the byte shift with scalar u64 composition: pack halves of the register, shift them, unpack the result. Every bit survived that. What the scalar composition destroyed was the compiler’s view of a single vector operation, which is the whole reason the vector rewrite measured slower than the struct it replaced.

Reinterpret the vector as sixteen bytes, hand @shuffle a zero vector as its second input, and let the negative mask lanes pull zeros into the vacated bytes:

inline fn rshift128(in: W128) W128 {
    const bytes: @Vector(16, u8) = @bitCast(in);
    const mask = @Vector(16, i8){ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, -1 };
    const shifted = @shuffle(u8, bytes, @as(@Vector(16, u8), @splat(0)), mask);
    return @bitCast(shifted);
}

That is _mm_srli_si128(c, 1) written in Zig, with no scalar arithmetic in the path. lshift128 is the same shape with the mask running the other way. Once the code matched the algorithm, ReleaseFast reported what Debug had been hiding: 332 M/s for the struct version, 1,618 M/s for the vector implementation, and 1,675 M/s after the final cleanup.

What Debug mode could and could not tell me

Debug caught regressions quickly and was the wrong place to judge a vector rewrite. The conversions and helper calls had not been inlined away, so the benchmark was measuring scaffolding that ReleaseFast removes. The false starts were still worth having. 48 M/s and 57 M/s told me something was wrong; they could not tell me what.

Correctness is what hid it: the byte-lane shift was bit-exact from the first version through the last, and an operation that produces the right answer is not an operation you go looking at.

The finished version holds the SFMT 1.5.1 reference vectors for both init paths at 1,675 M/s.