     egblas_shuffle_seed(n, x, incx, rd());
 }
 
+template <typename T>
+void egblas_par_shuffle_one(size_t n, T* x, T * y, size_t seed) {
+    // Allocate room for the states
+    curandState* states;
+    cuda_check(cudaMalloc((void**)&states, 64 * 64 * sizeof(curandState)));
+
+    if (n < 1000) {
+        setup_kernel<<<8, 1>>>(states, seed);
+
+        par_shuffle_one_kernel_states<8, 1, T><<<8, 1>>>(states, n, x, y);
+        par_shuffle_one_kernel_states<8, 2, T><<<8, 1>>>(states, n, x, y);
+        par_shuffle_one_kernel_states<8, 1, T><<<8, 1>>>(states, n, x, y);
+        par_shuffle_one_kernel_states<8, 2, T><<<8, 1>>>(states, n, x, y);
+    } else if (n < 50000) {
+        setup_kernel<<<64, 1>>>(states, seed);
+
+        par_shuffle_one_kernel_states<64, 1, T><<<64, 1>>>(states, n, x, y);
+        par_shuffle_one_kernel_states<64, 2, T><<<64, 1>>>(states, n, x, y);
+        par_shuffle_one_kernel_states<64, 1, T><<<64, 1>>>(states, n, x, y);
+        par_shuffle_one_kernel_states<64, 2, T><<<64, 1>>>(states, n, x, y);
+    } else {
+        setup_kernel<<<64, 64>>>(states, seed);
+
+        par_shuffle_one_kernel_states<64 * 64, 1, T><<<64, 64>>>(states, n, x, y);
+        par_shuffle_one_kernel_states<64 * 64, 2, T><<<64, 64>>>(states, n, x, y);
+        par_shuffle_one_kernel_states<64 * 64, 1, T><<<64, 64>>>(states, n, x, y);
+        par_shuffle_one_kernel_states<64 * 64, 2, T><<<64, 64>>>(states, n, x, y);
+    }
+
+    cuda_check(cudaFree(states));
+}
+
 void egblas_par_shuffle_seed(size_t n, void* x, size_t incx, void* y, size_t incy, size_t seed){
     std::default_random_engine g(seed);
 


    distribution_t dist;

    // Optimized version for arrays of 4B
    if (incx == 4) {
        uint32_t* x_flat = reinterpret_cast<uint32_t*>(x);

        // Allocate room for the states
        curandState* states;
        cuda_check(cudaMalloc((void**)&states, 64 * 64 * sizeof(curandState)));

        if (n < 1000) {
            setup_kernel<<<8, 1>>>(states, seed);

            shuffle_4_kernel_states<8, 1, uint32_t><<<8, 1>>>(states, n, x_flat);
            shuffle_4_kernel_states<8, 2, uint32_t><<<8, 1>>>(states, n, x_flat);
            shuffle_4_kernel_states<8, 1, uint32_t><<<8, 1>>>(states, n, x_flat);
            shuffle_4_kernel_states<8, 2, uint32_t><<<8, 1>>>(states, n, x_flat);
        } else if (n < 50000) {
            setup_kernel<<<64, 1>>>(states, seed);

            shuffle_4_kernel_states<64, 1, uint32_t><<<64, 1>>>(states, n, x_flat);
            shuffle_4_kernel_states<64, 2, uint32_t><<<64, 1>>>(states, n, x_flat);
            shuffle_4_kernel_states<64, 1, uint32_t><<<64, 1>>>(states, n, x_flat);
            shuffle_4_kernel_states<64, 2, uint32_t><<<64, 1>>>(states, n, x_flat);
        } else {
            setup_kernel<<<64, 64>>>(states, seed);

            shuffle_4_kernel_states<64 * 64, 1, uint32_t><<<64, 64>>>(states, n, x_flat);
            shuffle_4_kernel_states<64 * 64, 2, uint32_t><<<64, 64>>>(states, n, x_flat);
            shuffle_4_kernel_states<64 * 64, 1, uint32_t><<<64, 64>>>(states, n, x_flat);
            shuffle_4_kernel_states<64 * 64, 2, uint32_t><<<64, 64>>>(states, n, x_flat);
        }

        cuda_check(cudaFree(states));

#ifdef EGBLAS_SYNCHRONIZE
    cudaDeviceSynchronize();


} // end of anonymous namespace

template <typename T>
void egblas_shuffle_one(size_t n, T * x, size_t seed) {
    // Allocate room for the states
    curandState* states;
    cuda_check(cudaMalloc((void**)&states, 64 * 64 * sizeof(curandState)));

    if (n < 1000) {
        setup_kernel<<<8, 1>>>(states, seed);

        shuffle_one_kernel_states<8, 1, T><<<8, 1>>>(states, n, x);
        shuffle_one_kernel_states<8, 2, T><<<8, 1>>>(states, n, x);
        shuffle_one_kernel_states<8, 1, T><<<8, 1>>>(states, n, x);
        shuffle_one_kernel_states<8, 2, T><<<8, 1>>>(states, n, x);
    } else if (n < 50000) {
        setup_kernel<<<64, 1>>>(states, seed);

        shuffle_one_kernel_states<64, 1, T><<<64, 1>>>(states, n, x);
        shuffle_one_kernel_states<64, 2, T><<<64, 1>>>(states, n, x);
        shuffle_one_kernel_states<64, 1, T><<<64, 1>>>(states, n, x);
        shuffle_one_kernel_states<64, 2, T><<<64, 1>>>(states, n, x);
    } else {
        setup_kernel<<<64, 64>>>(states, seed);

        shuffle_one_kernel_states<64 * 64, 1, T><<<64, 64>>>(states, n, x);
        shuffle_one_kernel_states<64 * 64, 2, T><<<64, 64>>>(states, n, x);
        shuffle_one_kernel_states<64 * 64, 1, T><<<64, 64>>>(states, n, x);
        shuffle_one_kernel_states<64 * 64, 2, T><<<64, 64>>>(states, n, x);
    }

    cuda_check(cudaFree(states));
}

void egblas_shuffle_seed(size_t n, void* x, size_t incx, size_t seed){
    std::default_random_engine g(seed);

