C++ version

```
#include
#include
#include
#include
#include

#include
#include

#define RESTRICT __restrict
#define ITER 10

__global__
void dice(float * RESTRICT out,
const unsigned * RESTRICT A,
const unsigned * RESTRICT B,
int asize, int bsize)
{
__shared__ unsigned shared_A[16 * 32];
__shared__ unsigned shared_B[16 * 32];

int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;
//if (i > asize || j > bsize || i > j)
if (i > asize || j > bsize || blockIdx.x > blockIdx.y)
return;

// only the first row of threads are in charge of loading to shared memory
//if (threadIdx.y == 0)
// for (int index = threadIdx.x; index < 16 * 32; index += 16)
// {
// shared_A[index] = A[blockIdx.x * blockDim.x * 32 + index];
// shared_B[index] = B[blockIdx.y * blockDim.y * 32 + index];
// }
for (int index = threadIdx.y * 16 + threadIdx.x; index < 16 * 32; index += 256)
{
shared_A[index] = A[blockIdx.x * blockDim.x * 32 + index];
shared_B[index] = B[blockIdx.y * blockDim.y * 32 + index];
}
__syncthreads();

unsigned a, b, a_and_b;
unsigned pop_a = 0, pop_b = 0, pop_ab = 0;

for (int k = 0; k < 32; ++k) {
a = shared_A[threadIdx.x * 32 + k];
b = shared_B[threadIdx.y * 32 + k];
a_and_b = a & b;

pop_a += __popc(a);
pop_b += __popc(b);
pop_ab += __popc(a_and_b);
}

out[i * bsize + j] = 2.0 * pop_ab / (pop_a + pop_b);
}

using namespace std;

int main()
{
srand(time(NULL));

int encoding_size_in_bits = 1024;
int bflen = encoding_size_in_bits / 32;
int total_n = 1 << 20;
int total_out_size = total_n / bflen;

size_t input_size = sizeof(unsigned) * total_n;
size_t out_size = sizeof(float) * total_out_size * total_out_size;
unsigned * A = (unsigned *)malloc(input_size);
unsigned * B = (unsigned *)malloc(input_size);
float * out = (float *)malloc(out_size);
unsigned * A_dev, * B_dev;
float * out_dev;
cudaMalloc(&A_dev, input_size);
cudaMalloc(&B_dev, input_size);
cudaMalloc(&out_dev, out_size);

for (int i = 0; i < total_n; ++i)
{
A[i] = rand();
B[i] = rand();
}
printf("%d %d %d %d %d\n", A[0], A[1], A[2], A[3], A[4]);
printf("%d %d %d %d %d\n", B[0], B[1], B[2], B[3], B[4]);
cudaMemcpy(A_dev, A, total_n, cudaMemcpyHostToDevice);
cudaMemcpy(B_dev, B, total_n, cudaMemcpyHostToDevice);

cudaError_t err;
dim3 grid(total_out_size / 16, total_out_size / 16);
dim3 block(16, 16);

for (int i = 0; i < ITER; ++i)
{
auto start = chrono::high_resolution_clock::now();
dice<<>>(out_dev, A_dev, B_dev, total_out_size, total_out_size);
cudaStreamSynchronize(0);
auto end = chrono::high_resolution_clock::now();
err = cudaGetLastError();
if (err != cudaSuccess)
printf("Error in kernel\n");
cout << "Duration " << chrono::duration_cast(end - start).count() << endl;
}

cudaMemcpy(out, out_dev, out_size, cudaMemcpyDeviceToHost);
printf("%f %f %f %f %f\n", out[0], out[1], out[2], out[3], out[4]);

cudaDeviceSynchronize();
cudaProfilerStop();

return 0;
}
```