tandem rng tandem-fortran GitHub

Design

Generator layout

tandem_t holds a field-for-field bind(C) mirror of the C struct tandem_rng, so C reads it in place and returns it by value with the C layout. The tests compare the mirror's size and field offsets with the header through tandem_layout_matches.

Bounded integers

below(n) is Lemire's multiply and reject over the 32-bit or 64-bit draw, as Rng::urand of tandem-cuda. A rejected draw is discarded, so below consumes a varying number of draws, and n is an unsigned bit pattern like every integer here. fill_below(x, n) maps draw i of the plain fill to element i, and retries a rejected draw on a fallback generator keyed by the global draw index, the aligned start position over the draw width plus i. It consumes exactly size(x) draws, equals the scalar calls except where a draw is rejected, and a fill cut at any element boundary equals the whole fill.

Weighted choice

tandem_choice_t owns the alias table of Appendix C, which tandem.c builds in exact integers: column j holds mass cut(j) of index j and the rest of the capacity of alias(j). Each call passes the C library a struct that points at the table's arrays, so a copied table stays valid. A draw maps one 64-bit draw, so it never retries. The device fill runs tandem.cuh's fill_choice, the plain 64-bit fill with each draw mapped in the store, over a device copy of the table.

Normals

real64 normals are the 1024-layer ziggurat of Appendix A, one 64-bit draw each. next_normal64 and element i of a real64 fill_normal take that draw. A draw that misses the fast path, 0.43 % of them, continues on split(g) of sub(0x4e524d3634) of the key at position 0, where g is the draw's global index. So a fill equals the scalar calls, and a fill cut at any element equals the whole fill. An empty fill aligns the position to 64 bits. The host and the CUDA C fills are tandem.c's and tandem.cuh's own, bit for bit.

The CUDA Fortran generator tandem_dev_t ports the ziggurat. Its tables come from the spec's tables/normal_f64_zig1024.json through tools/gen_zig_tables.py into cuda/tandem_zig_tables.f90, each Float64 as its bit pattern. The reference logarithm takes libm's fma by bind(C), which nvfortran resolves to the device's in kernels, and keeps products out of sums, so its bits survive nvfortran's default contraction. nvfortran 25.3 has no ieee_fma. transfer in host and device code does not link in the host pass, so the code splits x with fraction and exponent, which give the same m and k exactly.

real32 normals are Box-Muller from two 32-bit float draws. A step gives two normals, the cos half and the sin half. next_normal32 returns the cos half and drops the other, next_normal_pair32 returns both, and fill_normal is the flattened pairs: an odd size keeps the cos half of its last pair and still consumes both uniforms. The host normals are tandem.c's polynomial Box-Muller with explicit fused multiply-adds, the same bits under every compiler.