feat(gpu): add circulant matrix for one vs many poly product

[andrei-stoian-zama]

To support encrypted GLWE x clear matrix product, all polys of the GLWE are multiplied with the clear matrix. For each poly of the GLWE, this PR builds a circulant matrix that is multiplied with the clear matrix to obtain the polynomial product. Sample N of this product contains the matrix product of the original clear vector that was encrypted and the clear matrix.

Now, with make test_core_crypto_gpu there are two new tests:

test core_crypto::gpu::algorithms::test::glwe_dot_product_with_clear::test_gpu_glwe_dot_product_with_clear_test_params_4_bits_native_u64 ... ok

test core_crypto::gpu::algorithms::test::glwe_dot_product_with_clear::test_gpu_poly_product_with_clear_test_params_4_bits_native_u64 ... ok

[andrei-stoian-zama]

this test checks the poly product without noise

[andrei-stoian-zama]

I added a stride_C parameter, since I use this function to write the output matrix in a bigger buffer (a GLWE list) that has a bigger stride than the width of C which is N.

[andrei-stoian-zama]

here is the python equivalent of the algorithm:
with bi=blockIdx.x, bj=blockIdx.y, ti=threadIdx.x,tj=threadIdx.y

def make_circulant_transpose_cuda(v1):
    N = v1.shape[0]
    result = np.zeros((N * N,), dtype=np.uint64)

    BS = 4
    for bi in range(0, N // BS):
        for bj in range(0, N // BS):
            buf = np.zeros((2 * BS - 1), dtype=np.uint64)

            block_start = bi * BS * N + bj * BS

            for ti in range(BS):
                for tj in range(BS):
                    tid = ti * BS + tj
                    if tid < 2 * BS - 1:
                        read_idx_start = (bj - bi) * BS + tid - BS + 1
                        if read_idx_start < 0:
                            read_idx_start = N + read_idx_start
                        buf[tid] = v1[read_idx_start]

           # Sync threads

            for ti in range(BS):
                for tj in range(BS):
                    fact = 1
                    if bi * BS + ti > bj * BS + tj:
                        fact = -1
                    result[block_start + ti * N + tj] = buf[tj - ti + BS - 1] * fact

    return result.reshape((N, N))

[andrei-stoian-zama]

only the 2 rows of threads in the block read data since we only need to read 2*block_tile-1 values for a block of block_tile x block_tile threads

[andrei-stoian-zama]

sometimes check with n_polys=poly_size, sometimes check with arbitrary number of polys