Function normalMersenne

src/backend/cuda/kernel/random_engine.hpp:976–1031

template<typename T>
__global__ void normalMersenne(T *const out, uint *const gState,
                               const uint *const pos_tbl,
                               const uint *const sh1_tbl,
                               const uint *const sh2_tbl, uint mask,
                               const uint *const g_recursion_table,
                               const uint *const g_temper_table,
                               uint elementsPerBlock, uint elements) {
    __shared__ uint state[STATE_SIZE];
    __shared__ uint recursion_table[TABLE_SIZE];
    __shared__ uint temper_table[TABLE_SIZE];
    uint start = blockIdx.x * elementsPerBlock;
    uint end   = start + elementsPerBlock;
    end        = (end > elements) ? elements : end;
    int iter = divup((end - start) * sizeof(T), blockDim.x * 4 * sizeof(uint));

    uint pos = pos_tbl[blockIdx.x];
    uint sh1 = sh1_tbl[blockIdx.x];
    uint sh2 = sh2_tbl[blockIdx.x];
    state_read(state, gState);
    read_table(recursion_table, g_recursion_table);
    read_table(temper_table, g_temper_table);
    __syncthreads();

    uint index                    = start;
    int elementsPerBlockIteration = blockDim.x * 4 * sizeof(uint) / sizeof(T);
    uint o[4];
    int offsetX1 = (STATE_SIZE - N + threadIdx.x) % STATE_SIZE;
    int offsetX2 = (STATE_SIZE - N + threadIdx.x + 1) % STATE_SIZE;
    int offsetY  = (STATE_SIZE - N + threadIdx.x + pos) % STATE_SIZE;
    int offsetT  = (STATE_SIZE - N + threadIdx.x + pos - 1) % STATE_SIZE;
    int offsetO  = threadIdx.x;

    for (int i = 0; i < iter; ++i) {
        for (int ii = 0; ii < 4; ++ii) {
            uint r = recursion(recursion_table, mask, sh1, sh2, state[offsetX1],
                               state[offsetX2], state[offsetY]);
            state[offsetO] = r;
            o[ii]          = temper(temper_table, r, state[offsetT]);
            offsetX1       = (offsetX1 + blockDim.x) % STATE_SIZE;
            offsetX2       = (offsetX2 + blockDim.x) % STATE_SIZE;
            offsetY        = (offsetY + blockDim.x) % STATE_SIZE;
            offsetT        = (offsetT + blockDim.x) % STATE_SIZE;
            offsetO        = (offsetO + blockDim.x) % STATE_SIZE;
            __syncthreads();
        }
        if (i == iter - 1) {
            partialBoxMullerWriteOut128Bytes(out, index + threadIdx.x, o[0],
                                             o[1], o[2], o[3], elements);
        } else {
            boxMullerWriteOut128Bytes(out, index + threadIdx.x, o[0], o[1],
                                      o[2], o[3]);
        }
        index += elementsPerBlockIteration;
    }
    state_write(gState, state);
}

template<typename T>