     return c;
 }
 
+template < class RealType, class IntType >
+__device__ void fill_random_pool( curandStateMtgp32_t * state, volatile RealType * rands, RealType bin_start, RealType bin_width, IntType offset, IntType N ) {
+    N += (!!(N & (THREAD_NUM - 1))) * (THREAD_NUM - (N & (THREAD_NUM - 1))); // pad N st N = 256 n; assumes THREAD_NUM = 2^{p}
+
+    // over generate random numbers per bin
+    // avoids branch divergence
+    while(offset < N ) {
+        RealType event = curand_uniform( state );
+
+        rands[ offset ] = bin_start + event * bin_width;
+
+        offset += THREAD_NUM;
+    }
+}
 
 /**
  *  This version assumes that allele_list is organized into units of WARP_SIZE (32) alleles

        unsigned int * seq = sequences + (seq_idx * sequence_width);

        // use all threads to generate random numbers
        event_bins[ tid ] = curand_poisson( &states[blockIdx.x], rho );
        __syncthreads();

#ifdef LOG_RANDOM_EVENTS

        evt_list += BLOCK_NUM_MAX * THREAD_NUM;
#endif  // LOG_RANDOM_EVENTS

        int_type evt_lo, evt_hi;
        int_type max_rounds = warp_event_range( event_bins[lane_id], lane_id, evt_lo, evt_hi);
        __syncthreads();

        // within a warp each thread (bin) has [evt_lo, evt_hi) events

        // over generate random numbers per bin
        // avoids branch divergence
        fill_random_pool( states, rand_pool, bin_start, bin_range, lane_id, (max_rounds << 5) );
        __syncthreads();

#ifdef LOG_RANDOM_EVENTS
        i = tid;
        while( i < 1024 ) {
            pool[ i ] = rand_pool[ i ];
            i += THREAD_NUM;
        }
        pool += BLOCK_NUM_MAX * 1024;
        __syncthreads();
#endif  // LOG_RANDOM_EVENTS

        i = tid;
