# Pseudo Random Number Generation on the GPU

Idea

For a Path Tracing application I need to generate good quality pseudo random numbers in the closed range [0~1]. Because I'm doing this on the GPU (HLSL Shader Model 5) I can only use 32bit variables. My initial approach is the following set-up:

1. Ever Frame the pixel shader receive a good pseudo random number ([0~1]) from the CPU using C++'s std::mt19937 generator and std::uniform_real_distribution.

2. Because for each pixel this number is the same I also use the screen coordinates u and v of each pixel these are also in [0~1].

3. I then call the Multiply With Carry method like below.

Algorithm

// seed is the value given from the cpu
float3 random = Random(seed + u, seed + v)l

// Multiply With Carry, returns 3 floating point values {x, y, z}
// x: the random number
// y, z: new seeds for the next time we need a random number
float3 Random(float seed_a, float seed_b)
{
uint m_z = asuint(seed_a);
uint m_w = asuint(seed_b);

m_z = 36969 * (m_z & 65535) + (m_z >> 16);
m_w = 18000 * (m_w & 65535) + (m_w >> 16);

float r = ((m_z << 16) + m_w) / (float)0xFFFFFFFF;
return float3(r, asfloat(m_z), asfloat(m_w));
}


This produces the following output. (The left part is the random number obtained from the Random method for this pixel, the right part is the visualization of u and v as Red and Green.

Result

As you can see there is clearly a pattern, so the randomness is not good at all. Which hurts the performance of my algorithm tremendously. This is probably due to the fact that the original Multiply with Carry method assumes m_z and m_w are 64 bit integers, not 32 bit ones.

What I want

What I'm looking for is a way to fix my implementation of the Multiply with Carry method so that it produces reasonably good pseudo random numbers and works in the closed [0~1] interval instead of the open [0~1] interval. However since it is quite possible that this method can only work right using 64 bits integers I'm would also be really glad if someone can suggest another pseudo random number generator algorithm that:

• Works with 32 bit numbers
• Produces uniformly distributed results in the closed [0~1] interval
• Does not require too much state information, (this is why I chose MwC since it only needs to store 2 variables) since that is hard on the GPU. 16 32 bit variables would be the maximum I think since I can store that in 4x4 matrix which is easy to pass around.
-
It's badly wrong to re-seed a PRNG each time you need a random number -- that means that the randomness you end up with does not really come from the PRNG itself, but just from how random you can make each new seed (which, in this case, is "not very"). With SIMD parallelism you'll probably want to seed a PRNG for each processing unit once and for all, and then keep carrying that seed forwards through the entire computation. (But I don't know how easy it is to express this if you otherwise let the compiler parallelize your loops for you). – Henning Makholm Mar 22 '13 at 10:35
Unfortunately I cant share information between pixels using a shader so I need to seed each pixel individually. Shaders also do not retain information between frames so that makes everything a lot harder. Still I do believe that the performance I'm getting now is really abysmal. Note that for this example I only take one random number from the PRNG but usually I take hundreds per pixel. – Roy T. Mar 22 '13 at 10:44
So the code you're showing is not what actually runs? Your comment says "new seeds for the next time we need a random number". If you're using "hundreds" of random numbers for each seed, you should have plenty of time to run the initial seed through, say, one of the integer hashes from burtleburtle.net, which ought to get rid of the most visually obvious patterns in the first outputs. – Henning Makholm Mar 22 '13 at 10:51
No wait, that is how MwC works, the two seeds are its state. I mean I only seed them using the number from the CPU and the U and V variables once each frame. – Roy T. Mar 22 '13 at 11:00
You might be interested in this site proposal – Daniel Pendergast Nov 18 '13 at 18:49

There's some Perlin noise code on the following site that produces a number between -1 and +1 that I have used in the past. I am sure it could be changed to produce 0 to 1. http://www.gamedev.net/page/resources/_/technical/game-programming/simple-clouds-part-1-r2085

-
Hey MistyManor, I tried those links, and they indeed give me nice noisy image but it has a really low period which in the end was insufficient for the Monte Carlo simulation which is why I kept searching for something better. I answered my own question now after some searching, but I wish I could accept your answer as well since it is certainly not wrong and will be sufficient for most uses. – Roy T. Mar 24 '13 at 21:02

I queried fellow students and they brought my attention to this paper: http://http.developer.nvidia.com/GPUGems3/gpugems3_ch37.html

Since the paper is quite long and might not be online forever below is the general idea:

Linear Congruent Generators are ideal for use on the GPU because they are simple and do not require much state (only the previously generated number). but they are not 'random enough' for, for example, Monte Carlo based simulation. A generator like a Mersenne Twister would be better but requires way too much state to be stored.

The solution proposed by the paper is combining several LCGs using a combined Tausworthe Generator (as used by the Mersenne Twister) this guarantees a much better randomness without having to store as much state as the Mersenne Twister. The final algorithm looks like this:

struct RandomResult
{
uint4 state;
float value;
};

uint TausStep(uint z, int S1, int S2, int S3, uint M)
{
uint b = (((z << S1) ^ z) >> S2);
return (((z & M) << S3) ^ b);
}

uint LCGStep(uint z, uint A, uint C)
{
return (A * z + C);
}

RandomResult Random(uint4 state)
{
state.x = TausStep(state.x, 13, 19, 12, 4294967294);
state.y = TausStep(state.y, 2, 25, 4, 4294967288);
state.z = TausStep(state.z, 3, 11, 17, 4294967280);
state.w = LCGStep(state.w, 1664525, 1013904223);

RandomResult result;
result.state = state;
result.value = 2.3283064365387e-10 * (state.x ^ state.y ^ state.z ^ state.w);

return result;
}


Note that the initial seed values for state should be larger than 128! (For background information reed the paper) and that you should fill the seed with 4 good random numbers from the CPU + four values unique for that pixel to get a nice result.

-