mce_qd: better parameter selection + fwht speed

src/mce_qd.cpp

@@ -26,25 +26,31 @@ using namespace mce_qd;
 #include <set>
 
 int mce_qd::generate (pubkey&pub, privkey&priv, prng&rng,
-                      uint m, uint T, uint block_discard)
+                      uint m, uint T, uint block_count, uint block_discard)
 {
-	priv.fld.create (m);
-	priv.T = T;
-	uint t = 1 << T;
-
 	//convenience
 	gf2m&fld = priv.fld;
 	std::vector<uint>&essence = priv.essence;
 
+	//initial stuff and sizes
+	fld.create (m);
+	priv.T = T;
+	uint t = 1 << T,
+	     block_size = t,
+	     h_block_count = block_count + block_discard,
+	     n = h_block_count * t;
+
+	if (block_count <= m) return 2; //lower bound on block_count
+	if (n > fld.n / 2) return 2; //n <= q/2
+
 	std::vector<uint> support, Hsig;
 	polynomial g;
 	uint i, j;
 
-	//prepare for data
+	//prepare data arrays
-	Hsig.resize (fld.n / 2);
+	Hsig.resize (n);
-	support.resize (fld.n / 2);
+	support.resize (n);
 	essence.resize (m);
-	//note that q=2^m, algo. n=q/2, log n = m-1
 
 	//retry generating until goppa code is produced.
 	for (;;) {
@@ -58,7 +64,7 @@ int mce_qd::generate (pubkey&pub, privkey&priv, prng&rng,
 		essence[m - 1] = fld.inv (Hsig[0]);
 		//essence[m-1] is now used as precomputed 1/h_0
 
-		for (uint s = 0; s < m - 1; ++s) {
+		for (uint s = 0; (1 << s) < n; ++s) {
 			i = 1 << s; //i = 2^s
 
 			Hsig[i] = choose_random (fld.n, rng, used);
@@ -66,6 +72,7 @@ int mce_qd::generate (pubkey&pub, privkey&priv, prng&rng,
 			used.insert (fld.inv (essence[s]) );
 
 			for (j = 1; j < i; ++j) {
+				if (i + j >= n) break;
 				Hsig[i + j] = fld.inv
 				              (fld.add
 				               (fld.inv (Hsig[i]),
@@ -103,7 +110,7 @@ int mce_qd::generate (pubkey&pub, privkey&priv, prng&rng,
 		if (!consistent) continue; //retry
 
 		//compute the support, retry if it has two equal elements.
-		for (i = 0; i < fld.n / 2; ++i) {
+		for (i = 0; i < n; ++i) {
 			support[i] = fld.add (
 			                 fld.inv (Hsig[i]),
 			                 essence[m - 1]);
@@ -124,12 +131,7 @@ int mce_qd::generate (pubkey&pub, privkey&priv, prng&rng,
 		}
 		if (!consistent) continue; //retry
 
-		//now the blocks.
+		//now the blocks. First assemble blocks to bl
-		uint block_size = 1 << T,
-		     h_block_count = (fld.n / 2) / block_size,
-		     block_count = h_block_count - block_discard;
-
-		//assemble blocks to bl
 		std::vector<polynomial> bl, blp;
 		bl.resize (h_block_count);
 		for (i = 0; i < h_block_count; ++i) {
@@ -211,17 +213,20 @@ int privkey::prepare()
 	uint omega;
 
 	uint block_size = 1 << T,
-	     block_count = hperm.size();
+	     block_count = hperm.size(),
+	     h_block_count = block_perm.size(),
+	     n = h_block_count * block_size;
 
 	//compute H signature from essence
-	Hsig.resize (fld.n / 2);
+	Hsig.resize (n);
 	Hsig[0] = fld.inv (essence[fld.m - 1]);
-	for (s = 0; s < fld.m - 1; ++s) {
+	for (s = 0; (1 << s) < n; ++s) {
 		i = 1 << s; //i = 2^s
 
 		Hsig[i] = fld.inv (fld.add (essence[s], essence[fld.m - 1]) );
 
-		for (j = 1; j < i; ++j)
+		for (j = 1; j < i; ++j) {
+			if (i + j >= n) break;
 			Hsig[i + j] = fld.inv
 			              (fld.add
 			               (fld.inv (Hsig[i]),
@@ -230,6 +235,7 @@ int privkey::prepare()
 			                    essence[fld.m - 1]
 			                ) ) );
 		}
+	}
 
 	//goppa polynomial with omega=0
 	std::set<uint> used;
@@ -249,8 +255,8 @@ int privkey::prepare()
 	}
 
 	//compute the support with omega=0
-	support.resize (fld.n / 2);
+	support.resize (n);
-	for (i = 0; i < fld.n / 2; ++i) {
+	for (i = 0; i < n; ++i) {
 		//don't check discarded support
 		if (block_perm[i / block_size] >= block_count) continue;
 		support[i] = fld.add
@@ -262,7 +268,7 @@ int privkey::prepare()
 		used.insert (support[i]);
 	}
 
-	//choose omega
+	//choose some omega
 	omega = fld.n;
 	for (i = 0; i < fld.n; ++i)
 		if (!used.count (i) ) {
@@ -295,7 +301,7 @@ int privkey::prepare()
 	permuted_support.resize (block_size * block_count);
 
 	//permute support
-	for (i = 0; i < (fld.n / 2) / block_size; ++i) {
+	for (i = 0; i < h_block_count; ++i) {
 		pos = block_perm[i];
 		if (pos >= block_count) continue; //was discarded
 		blk_perm = block_perms[pos];
@@ -376,6 +382,11 @@ int pubkey::encrypt (const bvector & in, bvector & out, const bvector&errors)
 	g.resize (t);
 	r.resize (t);
 
+	std::vector<int> c1, c2, c3;
+	c1.resize (t);
+	c2.resize (t);
+	c3.resize (t);
+
 	for (i = 0; i < qd_sigs.size(); ++i) {
 		//plaintext block
 		in.get_block (i * t, t, p);
@@ -385,7 +396,7 @@ int pubkey::encrypt (const bvector & in, bvector & out, const bvector&errors)
 			qd_sigs[i].get_block (j * t, t, g);
 
 			//block result
-			fwht_dyadic_multiply (p, g, r);
+			fwht_dyadic_multiply (p, g, r, c1, c2, c3);
 			cksum.add_offset (r, t * j);
 		}
 	}

src/mce_qd.h

@@ -98,7 +98,7 @@ public:
 	bool unserialize (sencode*);
 };
 
-int generate (pubkey&, privkey&, prng&, uint m, uint T, uint b);
+int generate (pubkey&, privkey&, prng&, uint m, uint T, uint b, uint bd);
 }
 
 #endif

src/qd_utils.cpp

@@ -31,7 +31,6 @@ static void fwht (vector<int> x, vector<int>&r)
 {
 	int bs, s;
 	s = x.size();
-	r.resize (s);
 	bs = s >> 1;
 	r.swap (x);
 	while (bs) {
@@ -46,18 +45,18 @@ static void fwht (vector<int> x, vector<int>&r)
 	}
 }
 
-//we expect correct parameter size and preallocated out.
+/*
-void fwht_dyadic_multiply (const bvector& a, const bvector& b, bvector& out)
+ * we expect correct parameter size and preallocated out. Last 3 parameters are
+ * used as a cache - just supply the same vectors everytime when you're doing
+ * this multiple times.
+ */
+void fwht_dyadic_multiply (const bvector& a, const bvector& b, bvector& out,
+                           vector<int>&t, vector<int>&A, vector<int>&B)
 {
 
-	//lift everyting to Z.
-	vector<int> t, A, B;
 	uint i;
 
-	t.resize (a.size() );
+	//lift everyting to Z.
-	A.resize (a.size() );
-	B.resize (a.size() );
-
 	for (i = 0; i < a.size(); ++i) t[i] = a[i];
 	fwht (t, A);
 
@@ -109,6 +108,12 @@ bool qd_to_right_echelon_form (std::vector<std::vector<bvector> >&mat)
 
 	bvector tmp;
 	tmp.resize (bs);
+
+	vector<int> c1, c2, c3;
+	c1.resize (bs);
+	c2.resize (bs);
+	c3.resize (bs);
+
 	for (i = 0; i < h; ++i) { //gauss step
 		//first, find a nonsingular matrix in the column
 		for (j = i; j < h; ++j)
@@ -141,13 +146,13 @@ bool qd_to_right_echelon_form (std::vector<std::vector<bvector> >&mat)
 				if (k == w - h + i) continue;
 				fwht_dyadic_multiply
 				(mat[w - h + i][j],
-				 mat[k][j], tmp);
+				 mat[k][j], tmp, c1, c2, c3);
 				mat[k][j] = tmp;
 			}
 			//change the block on the diagonal
 			fwht_dyadic_multiply
 			(mat[w - h + i][j],
-			 mat[w - h + i][j], tmp);
+			 mat[w - h + i][j], tmp, c1, c2, c3);
 			mat[w - h + i][j] = tmp;
 
 			//and zero the column below diagonal
@@ -163,7 +168,7 @@ bool qd_to_right_echelon_form (std::vector<std::vector<bvector> >&mat)
 			//we can safely rewrite the diagonal here (nothing's behind it)
 			fwht_dyadic_multiply
 			(mat[w - i - 1][h - i - 1],
-			 mat[k][h - i - 1], tmp);
+			 mat[k][h - i - 1], tmp, c1, c2, c3);
 			mat[k][h - i - 1] = tmp;
 		}
 
@@ -182,7 +187,7 @@ bool qd_to_right_echelon_form (std::vector<std::vector<bvector> >&mat)
 				fwht_dyadic_multiply
 				(mat[w - i - 1]
 				 [h - j - 1],
-				 mat[k][h - j - 1], tmp);
+				 mat[k][h - j - 1], tmp, c1, c2, c3);
 				mat[k][h - j - 1] = tmp;
 			}
 			//I+I=0

src/qd_utils.h

@@ -25,8 +25,16 @@
 #include "bvector.h"
 #include "prng.h"
 
-//FWHT matrix mult in O(n log n). parameters MUST be of 2^m size.
+/*
-void fwht_dyadic_multiply (const bvector&, const bvector&, bvector&);
+ * FWHT matrix mult in O(n log n). parameters MUST be of 2^m size.
+ *
+ * c1-c3 are caches. Just supply the same vector objects everytime, it's gonna
+ * be a lot faster.
+ */
+void fwht_dyadic_multiply (const bvector&, const bvector&, bvector&,
+                           std::vector<int>& c1,
+                           std::vector<int>& c2,
+                           std::vector<int>& c3);
 
 //create a generator using fwht
 bool qd_to_right_echelon_form (std::vector<std::vector<bvector> >&matrix);