Function uniformMersenne

src/backend/cuda/kernel/random_engine.hpp:875–929

template<typename T>
__global__ void uniformMersenne(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, size_t 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 elementsPerBlockIteration = (blockDim.x * 4 * sizeof(uint)) / sizeof(T);
    int iter = divup((end - start), elementsPerBlockIteration);

    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;
    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) {
            partialWriteOut128Bytes(out, index + threadIdx.x, o[0], o[1], o[2],
                                    o[3], elements);
        } else {
            writeOut128Bytes(out, index + threadIdx.x, o[0], o[1], o[2], o[3]);
        }
        index += elementsPerBlockIteration;
    }
    state_write(gState, state);
}

template<typename T>
__global__ void normalPhilox(T *out, uint hi, uint lo, uint hic, uint loc,