Design

Fills

On x86-64 with GCC 12+ or clang, tandem.c compiles the row loop, the ziggurat and the normal and exponential loops a second time for AVX2 and FMA and picks that copy at run time. A plain -O2 build with no -m flags so gets 256-bit rows and hardware fused multiply-adds on Haswell, Zen and newer, and still runs on older CPUs. With -mavx2 -mfma only the AVX2 copy remains. Define TANDEM_NO_AVX2 to leave it out. Both copies give the same bits.

The eight chunks of a row step together in registers. With GCC 12+ or clang the step is written with vector extensions and compiles to NEON, SSE2 or, in the AVX2 copy, one 256-bit vector per word. Define TANDEM_NO_SIMD for the scalar version. After alignment every integer fill is one byte stream, so one routine serves all widths and floats convert in place.

Bounded integers

Bounded integers tandem_u32_below and tandem_u64_below (Lemire) go beyond the specification and return the same values as tandem-cuda. The scalar bounded draws loop on rejection. The bounded fills use the parallel rule of tandem-cuda instead: element i maps draw i of the plain fill, which keeps the SIMD speed, and a rejected draw retries on a fallback generator keyed by the global draw index, so a fill consumes exactly len draws and a fill cut anywhere equals the whole fill.

Normals

The f64 normals are the 1024-layer ziggurat of Appendix A of the specification. Element i of a fill comes from draw i of the u64 fill: the low 10 bits pick the layer, bit 10 the sign and the top 53 bits the abscissa, which one table product scales. 99.57 % of the draws land in the inner rectangles and are done. A miss continues on its own fallback stream, split(g) of sub(0x4e524d3634) of the key at position 0, where g is the global index of the draw, so a fill consumes exactly n draws, a scalar draw is element 0 of a fill, and a fill cut anywhere equals the whole fill. The slow path takes Marsaglia's tail beyond R and compares ln y with -x^2 / 2 in the wedges, so the polynomial logarithm below is its only function. The tables in tandem_normal_tables.h are generated by make tables from the spec's tables/normal_f64_zig1024.json, whose SHA-256 the script checks.

The fill runs a branch-free table pass over 512 draws at a time, in L1, which writes every element and lists the misses. The misses queue across blocks, and their fallback keys and first blocks are computed eight at a time in the lanes of the row loop, at a sixth to a fifth of the cost of separate seeds. The fill runs at a little over 40 % of the speed of the u64 fill, almost all of it in the u64 fill and the table pass.

The f32 normals are Box-Muller. tandem_normal2_f32 returns both halves of a pair (cos half first), tandem_normal_f32 returns the cos half, and tandem_fill_normal_f32 writes the flattened pairs, so an odd count still consumes 2 * ceil(n / 2) uniforms. They draw two f32 uniforms and compute in float, so they agree across ports to a few ulps because float libm functions differ. One loop without libm calls computes them: polynomials for the logarithm of the exponent-split argument and for the sine and cosine of the angle after an exact quarter-turn reduction, which the compiler vectorizes. make accuracy compares the f32 fills with libm on 5e7 pairs, within 3.3 ulps. make bench reports the normal fills too.

Every multiply-add of the logarithm and of the f32 loop is an explicit fma, and the code is built with floating point contraction off, so every compiler and target produces the same bits: tests/test_normal_bits.c checks hashes of f64 and f32 normals against the values from the M4. Where the build has no FMA instruction, as the base copy on x86-64, fma and fmaf would be library calls that keep the loops scalar. There tandem.c rounds each multiply-add once by an exact emulation instead, so the base copy vectorizes and gives the same bits. The float form takes the product exactly in double and rounds the sum to odd before the rounding to float. The double form adds Dekker's exact product by two error-free sums, the last rounded to odd (Boldo and Melquiond, IEEE Trans. Computers 57, 2008).

Exponentials

Exponentials tandem_exponential_f64 and _f32 and the fills tandem_fill_exponential_f64 and _f32 return -ln(1 - u) from one uniform u each, as Appendix A of the specification describes: element i of a fill comes from uniform i, so a fill equals the scalar draws and a fill cut anywhere equals the whole fill. They use no libm call, and tandem-cuda runs the same arithmetic on the device, so exponentials are bit exact across compilers, targets and devices. f64 uses the logarithm of the normals. f32 carries u = (2 - 2m) / (m + 1) as two floats and adds k ln 2 by an exact two-sum, because an error near 1 ulp moves 1 - exp(-x) to a neighbouring 2^-24 grid point. The maximum error is 1.1e-15 relative for f64 and 0.571 ulp for f32, checked over all 2^24 f32 inputs, and every f32 draw maps back to its own grid point.