tandem rng tandem-jax GitHub

Design

Fills

XLA path

Without jax_enable_x64 the 32-bit multiplies are built from 16-bit halves, which also serves backends without 64-bit integers.

CPU fills

The package in cpu/ builds the tandem-c fills as an XLA FFI handler for the CPU, under the same target name and operands as the CUDA one. tandem_jax loads the extension on import when it is installed. On the CPU, stream for every dtype but bool, uniform, normal, randint with scalar bounds, their stream_* forms and jax.random.bits then run on the tandem-c fills, at any chunk length. Fills of 256 KiB or more split into parts at stream row boundaries, one per thread of XLA's intra-op pool, and the call completes when the last part does. jit and vmap work: a batch writes its rows back to back. Stream draws, uniforms, bounded integers and float64 normals equal the XLA path bit for bit. float32 normals equal tandem-c bit for bit and the XLA path within the tolerance of Appendix A, since XLA's log, sin and cos are not tandem-c's polynomials.

CUDA kernels

tandem_jax loads the extension on import when it is installed. On a CUDA device, stream, uniform, normal, randint, their stream_* forms and the jax.random functions on a Tandem key then run on the kernels. lax.platform_dependent makes the choice when XLA lowers the code, so a CPU device in the same program keeps the XLA path. The kernels write 32- and 64-bit draws directly. Narrower dtypes come from the 32-bit word fill and one XLA pass. They need a chunk length of 8 or more, and randint with per-element bounds keeps the XLA path. jit and vmap work: a batch of keys or positions runs as one kernel launch. The values are the XLA path's: bit for bit for the stream, uniforms, bounded integers and float64 normals, and within the tolerance of Appendix A for float32 normals, which use the fast __sincosf of tandem-cuda, within 16 ulps + 1e-6. The float64 normals take the ziggurat in one kernel that continues each miss where it finds it, as tandem-cuda's short fills do.

Bounded integers

Bounded integers use Lemire's method on draw i for element i. The draw width follows the range, 32 bits up to 2^32 and 64 bits above, so the dtype does not change the values, and width names it as the u32 and u64 fills of the C and CUDA ports do. A rejected draw retries on the fallback stream split(g) of sub(0x424c573332) (0x424c573634 for 64 bits) from position 0, with g the index of the draw in the key's stream, so a fill cut at any element boundary equals the whole fill. The retry loop runs only when a draw was rejected.

Normals

float64 normals are the 1024-layer ziggurat of Appendix A: element i from 64-bit draw i and the spec's tables, which src/tandem_jax holds as the spec's file normal_f64_zig1024.json. A draw that misses the inner rectangles, 0.43 % of them, continues on the fallback stream split(g) of sub(0x4e524d3634), with g the index of the draw in the key's stream. XLA fuses a product into the sum that reads it, so the XLA path emulates the reference logarithm's fused multiply-adds exactly and keeps every other product out of a sum.

The XLA path writes each miss as NaN and packs the misses into 32-bit words. A binary search over the running count of the words then finds miss k, where nonzero over the fill took longer than the fill. The slow path takes a batch of the expected misses plus six standard deviations, and more misses take further batches. Two of its rounds settle all but a few dozen misses, and the rest finish on a gathered set of 256. XLA fuses an unrolled F into each consumer of its eight words, so the slow path runs F as a loop.

float32 normals are Box-Muller pairs: elements 2j and 2j + 1 are r cos 2 pi b and r sin 2 pi b from uniform draws 2j and 2j + 1, computed in float32 in one fused jit.