title: Generating random numbers from a gaussian distribution in CUDA
body: <p>I've searched a lot over the internet to find a way to generate random numbers on my CUDA device, within a kernel. The numbers must come from a gaussian distribution.</p>
  
  <p><a href="http://http.developer.nvidia.com/GPUGems3/gpugems3_ch37.html" rel="nofollow">The best thing I found</a> was from NVIDIA itself. It is the Wallace algorithm, that uses a uniform distribution to build a gaussian one. But the code samples they give lack explanation and I really need to understand how the algorithm goes, especially on the device. For example, they give:</p>
  
  <pre><code> __device__ void generateRandomNumbers_wallace(  
  unsigned seed,  // Initialization seed  
   float *chi2Corrections,  // Set of correction values  
   float *globalPool,  // Input random number pool  
   float *output  // Output random numbers  
  
  
      unsigned tid=threadIdx.x;  
      // Load global pool into shared memory.  
       unsigned offset = __mul24(POOL_SIZE, blockIdx.x);  
      for( int i = 0; i &lt; 4; i++ )  
        pool[tid+THREADS*i] = globalPool[offset+TOTAL_THREADS*i+tid];  
      __syncthreads();  
        const unsigned lcg_a=241;  
        const unsigned lcg_c=59;  
        const unsigned lcg_m=256;  
        const unsigned mod_mask = lcg_m-1;  
        seed=(seed+tid)&amp;mod_mask ;  
        // Loop generating outputs repeatedly  
       for( int loop = 0; loop &lt; OUTPUTS_PER_RUN; loop++ )  
        {  
          Transform();  
          unsigned intermediate_address;  
          i_a = __mul24(loop,8*TOTAL_THREADS)+8*THREADS *  
            blockIdx.x + threadIdx.x;  
          float chi2CorrAndScale=chi2Corrections[  
            blockIdx.x * OUTPUTS_PER_RUN + loop];  
          for( i = 0; i &lt; 4; i++ )  
            output[i_a + i*THREADS]=chi2CorrAndScale*pool[tid+THREADS*i];  
      }  
  </code></pre>
  
  <p>First of all, many of the variables declared aren't even used in the function! And I really don't get what the "8" is for in the second loop. I understand the "4" in the other loops have something to do with the 4x4 orthogonal matrix block, am I right? Could anyone give me a better idea of what is going on here?</p>
  
  <p>Anyway, does anyone have any good code samples I could use? Or does anyone have another way of generating random gaussian numbers in a CUDA kernel? Code samples will be much appreciated.</p>
  
  <p>Thanks!</p>