I’m writing this to tell you about an elegant and simple idea that might make your code go faster. Depending on how you use GPUs, your code might already be using this idea: it has gone mainstream in ML frameworks since I first heard about it and wrote a library ten years ago in 2016. The idea is not mine! I just find it aesthetically pleasing, and so might you.
(Don’t use my old library! It used unsupported internal cuRAND functions. There are better options now, like cuRANDDx.)
Here’s the idea: if you’re generating random numbers on a GPU, by default you should avoid RNGs with internal state that’s updated for every number produced. Instead, prefer a stateless RNG that works by encrypting or hashing a counter. There are some excellent options so you don’t need to roll your own, called Philox and Threefry. PyTorch, JAX and cuRANDDx all use these counter-based RNGs.
You can tell you’re using a stateful RNG if you have to allocate memory for per-thread RNG state, and then read and write that state.
Why counters are nice: both VRAM and memory bandwidth are scarce resources on the GPU. It’s likely that your GPU algorithm is memory bandwidth-bound, or perhaps VRAM-bound. By using a counter-based RNG, each GPU thread can generate its own sequence of random numbers using pure computation, without having to read or write VRAM. You also avoid any stalls caused by waiting for those memory operations. The amount of state can be large! Mersenne Twister needs 2.5 kB of state per thread, which is why no-one uses it unmodified on GPUs.
When this advice doesn’t apply: counter-based RNGs do more computation, so if your benchmarks reveal that your code is compute-bound, you might want to spend memory to save compute. In WGSL, Philox is particularly slow because it needs 64-bit multiplies that need to be emulated. Or if you initialize an RNG and use it for a single kernel launch and never need random numbers after that, then you don’t need to persist state to memory.
If you’ve read this far, I think I’ve communicated the main idea. If you like, you can keep reading for more explanation and some tips, or you can just play with these demos!
Demo time!
All these demos use WebGPU. They compare three methods: a fast stateful RNG (XORWOW), a counter-based RNG but with counter state stored in VRAM, and a counter-based RNG with no stored state. The latter two give a clean measurement of the cost of reading and writing the state. You’ll see that the stateless RNG doesn’t always win! Philox needs the top 32 bits of a 32x32-bit multiply. cuda has an instruction for that (__umulhi), but WGSL doesn’t. So these browser demos somewhat exacerbate Philox’s higher compute intensity because one of the key operations needs to be emulated using several 16-bit multiplies.
The first demo is really simple: it estimates pi by throwing lots of random darts in a square and counting the ones that land inside the circle. With Chrome on my 2021 M1 Pro, stateless Philox is about 70% faster than XORWOW:
This next one shows ink diffusion. That sounds fancy! But this demo is super simple, it has no fluid dynamics or anything like that. It models each particle of ink as a jiggling wiggling random walker, being jostled by tiny molecules. (For fancier and more accurate physics simulations, see some of my other projects). So at each timestep, each ink particle gets a small random gaussian added to its position in each dimension. On my machine, stateless Philox is about 75% faster than XORWOW on this demo:
The two demos above use a lot of random numbers and do little computation otherwise. Here’s something with more computation per random number sampled. It’s a WebGPU port of a simple path tracer I wrote in cuda. It renders an image by sending multiple rays through each pixel, each one slightly jittered, and letting them bounce around several times. By doing this we get physical effects like soft shadows and bounced colored light. On my laptop, the stateful XORWOW RNG performs the best on this workload:
Your results will be different depending on your hardware and your browser’s WebGPU implementation; Philox performance is worse on Firefox on my laptop. But the stateless version will always use less memory (zero memory for RNG state).
Why do people need random numbers on the GPU?
These days, the most common reason is probably machine learning. Most famously, LLM inference uses random numbers to sample from the next-token distribution. LLM training often uses dropout, which needs a lot of random numbers: every forward pass needs a fresh random mask over activations. For image and audio models, the training set is sometimes augmented on-device with random crops, flips and color jitter, so that models learn to be invariant to those transformations. In diffusion models, training adds random noise at a random timestep to each example, and sampling starts from pure Gaussian noise and might inject more during denoising steps. Reinforcement learning strategies often use huge amounts of random numbers. And on and on. Modern ML is all on GPU and uses random numbers everywhere.
The path tracer demo above is a simple example of physically based rendering, which needs multiple random numbers for every pixel.
GPU random numbers are pervasive in physics simulations. I’ve used them for diffusion monte carlo for high-accuracy quantum chemistry calculations, and also for setting the temperature in molecular dynamics simulations, where you need to sample velocities from the Maxwell-Boltzmann distribution.
Another place I’ve used random numbers is MCMC sampling of a Bayesian model. For example, I designed this website by pitting pairs of designs against each other and picking the winner, many times over. To pick the right pair to compare I used MCMC to sample from the posterior distribution, after including my votes. I’ll write about this elsewhere but yes, that’s right! Somehow even frontend design can use large amounts of random numbers! (My implementation used numpy on the CPU, so it wasn’t, strictly speaking, GPU random numbers. But Bayesian inference on the GPU is common.)
What’s wrong with using normal RNGs on GPUs?
Most RNGs have two main parts: state transition and output. There’s some internal state S, a state transition function t(S), and an output function o(S). The way you use the RNG is, you call o(S) to get the random numbers you actually want, then you call t(S) to move the state along:
S = initialize_rng(seed)
while (...):
r = o(S)
S = t(S)
# Use the random number r for whatever
On a GPU with many threads, each thread needs its own independent sequence of random numbers. (If each thread generated the same sequence, the whatever you’re doing with the random numbers would break, because they’re totally not random any more.) This looks something like:
# threadIdx is some integer that identifies this thread.
def concat_bits(a,b):
return (a << 32) + b
def initialize(S):
# Each thread gets initialized with a different value, so it'll
# generate a different sequence.
S[threadIdx] = initialize_rng(concat_bits(threadIdx, seed))
def execute_step(S):
S_local = S[threadIdx]
r = o(S_local)
S[threadIdx] = t(S_local)
# Use the random number r for whatever
def main():
S = allocate_on_gpu(num_threads)
launch_kernel(initialize, S)
while (...):
launch_kernel(execute_step, S)
The code looks so innocuous! And for a serial machine like a CPU, it really is. But there’s a big trap if you’re using this kind of RNG on the GPU. Do you see it? Our programming languages hide what’s often the most expensive part of a GPU program, the memory operations. Here’s the trap:
S_local = S[threadIdx] # needs to read from VRAM
S[threadIdx] = t(S_local) # needs to write to VRAM
On a GPU, memory bandwidth is almost everything. These memory operations are one of the nicer kinds because all the indices are right next to each other, analogous to sequential ops with DRAM. But it would be better to have no memory ops at all.
Additionally, S needs to be allocated on VRAM, so it takes space that you could use to hold more model parameters or particles or whatever.
Solution: get rid of the state
Why do we need to keep state at all? It’s because the sequence of random numbers needs to pick up from where it left off, and the way we’ve set things up (a t and an o function), there’s no cheap way to skip ahead to the nth random number in a sequence.
What if the transition function could skip ahead? That is, t now takes two arguments, the state S and an integer n, and returns the result of applying the original t function n times.
Then we wouldn’t need to store state at all, and instead could do something like this:
# threadIdx is some integer that identifies this thread.
def concat_bits(a,b):
return (a << 32) + b
def execute_step(stepIdx):
S_local = initialize_rng(concat_bits(threadIdx, seed))
S_local = t(S_local, stepIdx)
r = o(S_local)
# Use the random number r for whatever
def main():
# No need for initialization or allocation.
for stepIdx in range(num_steps):
launch_kernel(execute_step, stepIdx)
The simplest way to do this is for S to just be a counter, and t to just be “add 1”.
The elegant part I really like
But if the state is just a counter, how on earth are we going to get random numbers out? The key insight is to use fruits of research in cryptography. Any encryption or hash function worth the name can take orderly input and produce essentially random bits as output. A precise way to state this is that for a good encryption or hashing algorithm, if any bit in the input changes, then every bit in the output flips with probability one half. If this weren’t so, then an attacker could try to discover the key by perturbing inputs to the algorithm and seeing how the output changes.
(What I’ve said above is a simplification. Real cryptographic algorithms are subjected to lots of tests and might not achieve exactly the property I’ve described above. But they get close.)
Well, that property is very useful for making a random number generator as well! Simply encrypt a counter with some key (or hash it), and every time the counter changes, a new set of random bits has been generated. Usually, the seed becomes the key that each counter is encrypted with.
As far as I can tell, this idea was first made practical by Salmon et al. in 2011 in the paper Parallel Random Numbers: As Easy as 1, 2, 3. (I used to work with John Salmon at D. E. Shaw Research, where he did this work to generate large amounts of random numbers for molecular dynamics simulations on the Anton family of supercomputers). People had previously proposed using ciphers like this but they were too slow. Salmon et al. designed cheaper bijections that still pass the BigCrush test suite for randomness. This approach has since gone mainstream: PyTorch, JAX and the new cuRANDDx library all use counter-based random number generators. But the old cuRAND library still exposes a stateful API, probably for backwards compatibility.
Note that this kind of RNG implies a different API. There’s no need to talk of a state transition function, and no need to read or write state. There simply is no state. You simply ask for the Nth random number. If you need a million random numbers in parallel, assign each thread its own range for it to run its counters. Since every thread already has some index variables that distinguish it from every other thread, this is a very natural setup for GPU programs.
How to use these RNGs correctly
The caveat is that if you ask for the Nth random number twice, you’ll get the same result twice. Stateful RNG libraries hide this from you by automatically doing the state transition out of your sight. With counter-based RNGs, instead of explicitly managing state, you must implicitly ensure that the RNG is always called with different input bits in a given run of your program. Passing in the same bits will give you exactly the same output, and incrementing them in sequence will recapitulate a sequence you have seen before. This can be helpful because deterministic programs are easier to test and debug, but also a potential pitfall.
So here are some rules of thumb to avoid this problem and the sticky statistical issues that come up with it:
-
Use application-level indices as counters. For example, if you have a 32-bit index for each thread and 32-bit outer index for each iteration, then combine these to get 64 unique bits per thread per iteration.
-
Differentiate between calls on different source text lines using constants. So when you call the RNG, pass in something like
(iterIdx << 64) + (threadIdx << 32) + M, whereMis different for every place in the code where you call it. -
Ask for random numbers in sets of four. One Philox call produces four 32-bit numbers, so try to use all those bits of randomness before asking for a fresh set. The computation isn’t cheap, so make it go as far as you can.
A very easy bug to write is to pass in iterIdx + threadIdx + M. In a given thread this is correct and would produce a sequence of good random numbers. But notice what’s happening in the thread next door. When iterIdx = 1 and threadIdx = 0, we get the same random number as when iterIdx = 0 and threadIdx = 1, so thread t+1’s numbers are thread t’s numbers shifted by one step. The whole grid of threads is using the same sequence with different offsets. Guess why I know about this bug :)
That’s why you see all that bit-shifting I’m doing everywhere, for example in concat_bits. It’s to construct a single counter of 96 or 64 bits by combining two smaller counters. Real APIs will often just let you pass in several different counters instead of asking for a single big type that’s 160 bits or something. Either way, the right way to think about this is that every time the RNG is called the input must be a globally unique sequence of bits.
What do libraries actually use? Is it really crypto-strength encryption and hash algorithms?
No, that would be wasteful. As of November 2026, nvidia’s cuRANDDx and libcu++ use Philox, PyTorch uses Philox on the GPU, JAX uses Threefry, Triton uses something based on Philox and C++26 just added a std::philox_engine.
Philox is a fast bijection built in the style of a cipher. Threefry is the Threefish cipher (the one inside the Skein hash) with fewer rounds. WebGPU doesn’t have a built-in random number generator in WGSL, so you’ll need to port one of the existing ones. In practice, a lot of WGSL code uses a cheap hash like PCG, which is the same idea with weaker mixing, because WGSL doesn’t have 32x32->64-bit integer multiply.
None of these RNG schemes produces cryptographically secure random numbers! They could be designed to do that but for most applications that would be overkill. So it’s faster to steal the ideas from cryptography but change the specifics to suit the problem at hand.
Bonus: reasons to reuse random numbers
I’d like to leave you with some interesting reasons to reuse random numbers on purpose. You can do this easily with counter-based RNGs by simply arranging to pass in the same seed and the same counter values.
Suppose you’re computing something that uses random numbers as well as some fixed parameters, for example pricing options. You can numerically compute a derivative with respect to some parameter by changing the value of the parameter and rerunning the computation with the same sequence of random numbers. If you’d used fresh random numbers, the noise due to randomness might swamp the difference that’s due to the change in the parameter. Another way to think of this is that you might need a lot more random numbers to converge the calculation enough to where the noise is smaller than the actual difference you’re trying to measure. By reusing random numbers, you can cancel out most of the noise.
When training a model with dropout, you want to drop out the same units in both the forward and backward passes. You might think this means you need to store which units are dropped. That takes memory. Instead, you can replay the random numbers exactly so you end up choosing the same units to drop in the backward pass. This is what FlashAttention does: it saves the RNG state instead of the dropout mask, and rebuilds the mask in the backward pass.
By using counter-based RNGs, simulations can be made precisely deterministic even if the underlying work was split differently on the hardware. For example, if you run the same simulation on a bigger GPU, you might use more threads for each kernel launch. If the random number generation is done carefully, for example by tying the counter to particle ids rather than thread ids, you can have both simulations produce exactly the same trajectory. This is important so that scientific findings can be reproduced by different labs.