ilyamironov: Fix sampling bug and add unit test for sampled distribution

A unit test for sampling routines.

- Constant- and variable-time sampling fixed; - Chi_squared test added

- make test command added - chi squared test is run 10,000 times

Restored and fixed comparator routine from BCNS. Unit test added to rlwe_test.c.

Fixing merger with master.

Incorporating Joppe's comments.

Makefile tweaks

Hide compiler warnings.

Pretty print.

.gitignore

@@ -29,3 +29,4 @@
 *.hex
 rlwe_benchmark
 rlwe_main
+rlwe_test

Makefile

@@ -18,19 +18,23 @@ PRNG_LDFLAGS=-lcrypto
 PRNG_CCFLAGS=-I/usr/local/opt/openssl/include
 PRNG_LDFLAGS=-L/usr/local/opt/openssl/lib -lcrypto
 
-CCFLAGS=-O3 -Wall -Wextra -std=c99 -Wno-deprecated-declarations -D$(PRNG) $(PRNG_CCFLAGS)
+CCFLAGS=-O3 -Wall -Wextra -std=c99 -Wno-deprecated-declarations -D$(PRNG) $(PRNG_CCFLAGS) #-DCONSTANT_TIME
 LDFLAGS=$(PRNG_LDFLAGS)
 
 all:
 	$(CC) $(CCFLAGS) -Wno-unused-function -c fft.c
 	$(CC) $(CCFLAGS) -Wno-unused-function -c rlwe.c
 	$(CC) $(CCFLAGS) -Wno-unused-function -c rlwe_kex.c
 	$(CC) $(CCFLAGS) -Wno-unused-function -Wno-unused-parameter -c rlwe_rand.c
-	$(CC) $(CCFLAGS) -o rlwe_main -lcrypto rlwe_main.c fft.o rlwe.o rlwe_kex.o rlwe_rand.o $(LDFLAGS)
-	$(CC) $(CCFLAGS) -o rlwe_benchmark -lcrypto rlwe_benchmark.c fft.o rlwe.o rlwe_kex.o rlwe_rand.o $(LDFLAGS)
+	$(CC) $(CCFLAGS) -o rlwe_main rlwe_main.c fft.o rlwe.o rlwe_kex.o rlwe_rand.o $(LDFLAGS)
+	$(CC) $(CCFLAGS) -o rlwe_benchmark rlwe_benchmark.c fft.o rlwe.o rlwe_kex.o rlwe_rand.o $(LDFLAGS)
+	$(CC) $(CCFLAGS) -Wno-unused-function -o rlwe_test rlwe_test.c fft.o rlwe_kex.o rlwe_rand.o $(LDFLAGS)
 
 clean:
-	rm fft.o rlwe.o rlwe_kex.o rlwe_rand.o rlwe_main rlwe_benchmark
+	rm fft.o rlwe.o rlwe_kex.o rlwe_rand.o rlwe_main rlwe_benchmark rlwe_test
+
+test:
+	./rlwe_test
 
 prettyprint:
 	astyle --style=java --indent=tab --pad-header --pad-oper --align-pointer=name --align-reference=name --suffix=none *.c *.h

rlwe.c

@@ -112,25 +112,21 @@ static int cmplt_ct(uint64_t *a, uint64_t *b) {
 	uint64_t r = 0; /* result */
 	uint64_t m = 0; /* mask   */
 	int i;
-	for (i = 0; i < 3; ++i) {
+	for (i = 2; i >= 0; --i) {
 		r |= ct_lt_u64(a[i], b[i]) & ~m;
 		m |= ct_mask_u64(ct_ne_u64(a[i], b[i])); /* stop when a[i] != b[i] */
 	}
 	return r & 1;
 }
 
 static uint32_t single_sample(uint64_t *in) {
-	uint32_t lower_index = 0, this_index = 32, upper_index = 64;
-	int i;
-	for (i = 0; i < 6; i++) {
-		if (cmplt_ct(in, rlwe_table[this_index])) {
-			upper_index = this_index;
-		} else {
-			lower_index = this_index;
-		}
-		this_index = (lower_index + upper_index) / 2;
+	size_t i = 0;
+
+	while (cmplt_ct(rlwe_table[i], in)) { // ~3.5 comparisons in expectation
+		i++;
 	}
-	return lower_index;
+
+	return i;
 }
 
 /* We assume that e contains two random bits in the two
@@ -148,7 +144,7 @@ static uint64_t dbl(const uint32_t in, int32_t e) {
 static uint32_t single_sample_ct(uint64_t *in) {
 	uint32_t index = 0, i;
 	for (i = 0; i < 52; i++) {
-		index = ct_select_u64(index, i, cmplt_ct(in, rlwe_table[i]));
+		index = ct_select_u64(index, i + 1, cmplt_ct(in, rlwe_table[i]));
 	}
 	return index;
 }
@@ -166,7 +162,7 @@ void rlwe_sample_ct(uint32_t *s, RAND_CTX *rand_ctx) {
 			r >>= 1;
 			// use the constant time version single_sample
 			s[i * 64 + j] = single_sample_ct(rnd);
-			t = 0xFFFFFFFF - s[i * 64 + j];
+			t = (uint32_t) - s[i * 64 + j];
 			s[i * 64 + j] = ct_select_u64(t, s[i * 64 + j], ct_eq_u64(m, 0));
 		}
 	}
@@ -226,10 +222,9 @@ void rlwe_sample(uint32_t *s, RAND_CTX *rand_ctx) {
 			RANDOM192(rnd, rand_ctx);
 			m = (r & 1);
 			r >>= 1;
-			m = 2 * m - 1;
 			s[i * 64 + j] = single_sample(rnd);
-			if (m == -1) {
-				s[i * 64 + j] = 0xFFFFFFFF - s[i * 64 + j];
+			if (m) {
+				s[i * 64 + j] = (uint32_t) - s[i * 64 + j];
 			}
 		}
 	}

rlwe_test.c

@@ -0,0 +1,140 @@
+/* This is free and unencumbered software released into the public domain.
+ *
+ * Anyone is free to copy, modify, publish, use, compile, sell, or
+ * distribute this software, either in source code form or as a compiled
+ * binary, for any purpose, commercial or non-commercial, and by any
+ * means.
+ *
+ * See LICENSE for complete information.
+ */
+
+#include <stdio.h>
+#include <stdint.h>
+#include <string.h>
+#include <stdlib.h>
+
+#include "rlwe_rand.h"
+
+#include "rlwe.c"
+
+const int ROUNDS = 10000;
+const int MAX_NOISE = 52;
+
+static uint64_t test_values[][3] = {
+	{0ULL, 0ULL, 0ULL}, // 0
+	{1ULL, 0ULL, 0ULL}, // 1
+	{2ULL, 0ULL, 0ULL}, // 2
+	{0xFFFFFFFFFFFFFFFFULL, 0ULL, 0ULL}, // 2^64-1
+	{0ULL, 1ULL, 0ULL}, // 2^64
+	{1ULL, 1ULL, 0ULL}, // 2^64+1
+	{0xFFFFFFFFFFFFFFFF, 0xFFFFFFFFFFFFFFFF, 0ULL}, // 2^128-1
+	{0ULL, 0ULL, 1ULL}, // 2^128
+	{1ULL, 0ULL, 1ULL}, // 2^128 + 1
+	{0x6a09e667f3bcc907ULL, 0xb2fb1366ea957d3eULL,  0x3adec17512775099ULL}, // -1
+	{0x6a09e667f3bcc908ULL, 0xb2fb1366ea957d3eULL,  0x3adec17512775099ULL}, // digits of sqrt{2}
+	{0x6a09e667f3bcc909ULL, 0xb2fb1366ea957d3eULL,  0x3adec17512775099ULL}, // +1
+	{0ULL, 0ULL, 0xFFFFFFFFFFFFFFFFULL},
+	{0ULL, 0xFFFFFFFFFFFFFFFFULL, 0xFFFFFFFFFFFFFFFFULL},
+	{0xFFFFFFFFFFFFFFFFULL, 0xFFFFFFFFFFFFFFFFULL, 0xFFFFFFFFFFFFFFFFULL}  // = 2^192-1
+};
+
+int test_comparer() {
+	int i, j;
+	int n = sizeof(test_values) / (sizeof(uint64_t) * 3);
+
+	printf("Testing the comparison function... ");
+
+	for (i = 0; i < n; i++) {
+		for (j = 0; j < n; j++) {
+			if ((i < j) != cmplt_ct(test_values[i], test_values[j])) {
+				printf("Test values (%d, %d) are compared incorrectly.\n", i, j);
+				return -1;
+			}
+		}
+	}
+	printf("OK\n");
+	return 1;
+}
+
+
+typedef void FuncSampler(uint32_t *s, RAND_CTX *rand_ctx);
+
+int chi_squared(FuncSampler sampler, RAND_CTX *rand_ctx) {
+	int ret = -1;
+	uint32_t s[1024];
+	uint32_t count[2 * MAX_NOISE];
+	memset(count, 0, sizeof(uint32_t) * 2 * MAX_NOISE);
+
+	int i, j;
+	for (i = 0; i < ROUNDS; i++) {
+		sampler(s, rand_ctx);
+		for (j = 0; j < 1024; j++) {
+			if (abs(s[j]) >= MAX_NOISE) {
+				printf("Noise is out of bounds\n");
+				goto err;
+			}
+			count[s[j] + MAX_NOISE]++;
+		}
+	}
+	double sum = 0;
+	int df = 0; // degrees of freedom
+	double chi_sqr = 0;
+	for (j = 1; j < 2 * MAX_NOISE; j++) { // skip count[0], which must be zero anyways
+		int v = abs(j - MAX_NOISE);
+		double point_mass = v > 0 ? (rlwe_table[v][2] - rlwe_table[v - 1][2]) * .5 : (double) rlwe_table[0][2];
+		double expect = ROUNDS * 1024 * (point_mass /  0xFFFFFFFFFFFFFFFFULL);
+
+		if (expect == 0) {
+			if (count[j] == 0) {
+				continue; // nothing to do
+			} else {
+				printf("Unexpected value (%d) has been output.\n", j - MAX_NOISE);
+				goto err;
+			}
+		}
+		printf("%2d: observed %7d, expected %g\n", j - MAX_NOISE, count[j], expect);
+		sum += expect;
+
+		chi_sqr += (count[j] - expect) * (count[j] - expect) / expect;
+		df++;
+	}
+	printf("The chi-squared statistic = %f (df = %d)\n", chi_sqr, df);
+
+	if (chi_sqr > 2 * df) {
+		printf("Chi-squared test failed.\n");
+		goto err;
+	}
+
+	ret = 1;
+
+err:
+	return ret;
+}
+
+int main() {
+	int ret = -1;
+
+	RAND_CTX rand_ctx;
+
+	if (!RAND_CHOICE_init(&rand_ctx)) {
+		printf("Cannot initialize PRNG.\n");
+		return -1;
+	}
+
+	if (test_comparer() <= 0) {
+		goto cleanup;
+	}
+
+#ifdef CONSTANT_TIME
+	printf("Testing constant time sampling...\n");
+	chi_squared(rlwe_sample_ct, &rand_ctx);
+#else
+	printf("Testing variable time sampling...\n");
+	chi_squared(rlwe_sample, &rand_ctx);
+#endif
+
+	ret = 0; // success!
+cleanup:
+	RAND_CHOICE_cleanup(&rand_ctx);
+	return ret;
+}