     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 && incy == 4){
        uint32_t* x_flat = reinterpret_cast<uint32_t*>(x);
        uint32_t* y_flat = reinterpret_cast<uint32_t*>(y);

        // 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_4_kernel_states<8, 1, uint32_t><<<8, 1>>>(states, n, x_flat, y_flat);
            par_shuffle_4_kernel_states<8, 2, uint32_t><<<8, 1>>>(states, n, x_flat, y_flat);
            par_shuffle_4_kernel_states<8, 1, uint32_t><<<8, 1>>>(states, n, x_flat, y_flat);
            par_shuffle_4_kernel_states<8, 2, uint32_t><<<8, 1>>>(states, n, x_flat, y_flat);
        } else if (n < 50000) {
            setup_kernel<<<64, 1>>>(states, seed);

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

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

        cuda_check(cudaFree(states));

#ifdef EGBLAS_SYNCHRONIZE
    cudaDeviceSynchronize();
#endif

        return;
    }


}

template <size_t Threads, size_t Distance, typename T>
__global__ void shuffle_one_kernel_states(curandState* states, size_t n, T* x) {
    const auto base_index  = threadIdx.x + blockIdx.x * blockDim.x;

    const size_t split_size = n / Threads;
