matrix: S and Sinv faster generation step 1

include/codecrypt.h

@@ -41,6 +41,7 @@ protected:
 public:
 	uint hamming_weight();
 	void add (const bvector&);
+	void add_range (const bvector&, uint, uint);
 	void add_offset (const bvector&, uint);
 	bool operator* (const bvector&); //dot product
 	bool zero() const;
@@ -71,6 +72,7 @@ class matrix : public std::vector<bvector>
 protected:
 	_ccr_declare_vector_item
 	_ccr_declare_matrix_item
+
 public:
 	uint width() const {
 		return size();
@@ -81,17 +83,26 @@ public:
 		return 0;
 	}
 
+	void resize2 (uint w, uint h, bool def = 0);
+
 	matrix operator* (const matrix&);
 	void mult (const matrix&); //right multiply - this*param
 
+	void zero ();
 	void unit (uint);
 
 	void compute_transpose (matrix&);
 	bool mult_vecT_left (const bvector&, bvector&);
 	bool mult_vec_right (const bvector&, bvector&);
-	bool compute_inversion (matrix&);
+	bool compute_inversion (matrix&,
+	                        bool upper_tri = false,
+	                        bool lower_tri = false);
 
 	bool set_block (uint, uint, const matrix&);
+	bool add_block (uint, uint, const matrix&);
+	bool set_block_from (uint, uint, const matrix&);
+	bool add_block_from (uint, uint, const matrix&);
+
 	bool get_left_square (matrix&);
 	bool strip_left_square (matrix&);
 	bool get_right_square (matrix&);
@@ -99,9 +110,9 @@ public:
 	void extend_left_compact (matrix&);
 
 	void generate_random_invertible (uint, prng&);
+	void generate_random_with_inversion (uint, matrix&, prng&);
 	bool create_goppa_generator (matrix&, permutation&, prng&);
 	bool create_goppa_generator (matrix&, const permutation&);
-
 };
 
 /*

lib/bvector.cpp

@@ -16,6 +16,13 @@ void bvector::add (const bvector&a)
 		item (i) = item (i) ^ a[i];
 }
 
+void bvector::add_range (const bvector&a, uint b, uint e)
+{
+	if (e > size() ) resize (e, 0);
+	for (uint i = b; i < e; ++i)
+		item (i) = item (i) ^ a[i];
+}
+
 void bvector::add_offset (const bvector&a, uint offset)
 {
 	if (offset + a.size() > size() ) resize (offset + a.size(), 0);

lib/matrix.cpp

@@ -3,6 +3,20 @@
 
 using namespace ccr;
 
+void matrix::resize2 (uint w, uint h, bool def)
+{
+	resize (w);
+	for (uint i = 0; i < w; ++i) item (i).resize (h, def);
+}
+
+void matrix::zero ()
+{
+	uint w = width(), h = height();
+	for (uint i = 0; i < w; ++i)
+		for (uint j = 0; j < h; ++j)
+			item (i, j) = 0;
+}
+
 void matrix::unit (uint size)
 {
 	clear();
@@ -43,7 +57,7 @@ void matrix::mult (const matrix&right)
 	}
 }
 
-bool matrix::compute_inversion (matrix&res)
+bool matrix::compute_inversion (matrix&res, bool upper_tri, bool lower_tri)
 {
 	//gauss-jordan elimination with inversion of the second matrix.
 	//we are computing with transposed matrices for simpler row ops
@@ -57,28 +71,42 @@ bool matrix::compute_inversion (matrix&res)
 	uint i, j;
 
 	//gauss step, create a lower triangular out of m, mirror to r
-	for (i = 0; i < s; ++i) {
-		//we need pivoting 1 at [i][i]. If there's none, get it below.
-		if (m[i][i] != 1) {
-			for (j = i + 1; j < s; ++j) if (m[j][i] == 1) break;
-			if (j == s) return false; //noninvertible
-			m[i].swap (m[j]);
-			r[i].swap (r[j]);
+	if (!upper_tri) for (i = 0; i < s; ++i) {
+			//we need pivoting 1 at [i][i]. If there's none, get it below.
+			if (m[i][i] != 1) {
+				for (j = i + 1; j < s; ++j) if (m[j][i] == 1) break;
+				if (j == s) return false; //noninvertible
+				m[i].swap (m[j]);
+				r[i].swap (r[j]);
+			}
+			//remove 1's below
+			if (lower_tri) {
+				for (j = i + 1; j < s; ++j) if (m[j][i]) {
+						m[j].add_range (m[i], 0, j + 1);
+						r[j].add_range (r[i], 0, j + 1);
+					}
+			} else {
+				for (j = i + 1; j < s; ++j) if (m[j][i]) {
+						m[j].add (m[i]);
+						r[j].add (r[i]);
+					}
+			}
 		}
-		//remove 1's below
-		for (j = i + 1; j < s; ++j) if (m[j][i]) {
-				m[j].add (m[i]);
-				r[j].add (r[i]);
-			}
-	}
 
-	//jordan step (we do it forward because it doesn't matter on GF(2))
-	for (i = 0; i < s; ++i)
-		for (j = 0; j < i; ++j)
-			if (m[j][i]) {
-				m[j].add (m[i]);
-				r[j].add (r[i]);
-			}
+	//jordan step
+	if (!lower_tri) {
+		if (upper_tri) {
+			for (i = s; i > 0; --i)
+				for (j = i - 1; j > 0; --j)
+					if (m[j-1][i-1])
+						r[j-1].add_range (r[i-1], i - 1, s);
+		} else {
+			for (i = s; i > 0; --i)
+				for (j = i - 1; j > 0; --j)
+					if (m[j-1][i-1])
+						r[j-1].add (r[i-1]);
+		}
+	}
 
 	r.compute_transpose (res);
 	return true;
@@ -109,6 +137,32 @@ void matrix::generate_random_invertible (uint size, prng & rng)
 	p.permute (lt, *this);
 }
 
+void matrix::generate_random_with_inversion (uint size, matrix&inversion, prng&rng)
+{
+	matrix lt, ut;
+	uint i, j;
+	// random lower triangular
+	lt.resize (size);
+	for (i = 0; i < size; ++i) {
+		lt[i].resize (size);
+		lt[i][i] = 1;
+		for (j = i + 1; j < size; ++j) lt[i][j] = rng.random (2);
+	}
+	// random upper triangular
+	ut.resize (size);
+	for (i = 0; i < size; ++i) {
+		ut[i].resize (size);
+		ut[i][i] = 1;
+		for (j = 0; j < i; ++j) ut[i][j] = rng.random (2);
+	}
+	*this = lt;
+	this->mult (ut);
+	ut.compute_inversion (inversion,	true, false);
+	lt.compute_inversion (ut, false, true);
+	inversion.mult (ut);
+}
+
+
 bool matrix::get_left_square (matrix&r)
 {
 	uint h = height();
@@ -221,3 +275,46 @@ bool matrix::set_block (uint x, uint y, const matrix&b)
 		for (uint j = 0; j < h; ++j) item (x + i, y + j) = b.item (i, j);
 	return true;
 }
+
+bool matrix::add_block (uint x, uint y, const matrix&b)
+{
+	uint h = b.height(), w = b.width();
+	if (width() < x + w) return false;
+	if (height() < y + h) return false;
+	for (uint i = 0; i < w; ++i)
+		for (uint j = 0; j < h; ++j)
+			item (x + i, y + j) =
+			    item (x + i, y + j)
+			    ^ b.item (i, j);
+	return true;
+}
+
+bool matrix::set_block_from (uint x, uint y, const matrix&b)
+{
+	uint h = b.height(),
+	     w = b.width(),
+	     mh = height(),
+	     mw = width();
+	if (mw > x + w) return false;
+	if (mh > y + h) return false;
+	for (uint i = 0; i < mw; ++i)
+		for (uint j = 0; j < mh; ++j)
+			item (i, j) = b.item (x + i, y + j);
+	return true;
+}
+
+bool matrix::add_block_from (uint x, uint y, const matrix&b)
+{
+	uint h = b.height(),
+	     w = b.width(),
+	     mh = height(),
+	     mw = width();
+	if (mw > x + w) return false;
+	if (mh > y + h) return false;
+	for (uint i = 0; i < mw; ++i)
+		for (uint j = 0; j < mh; ++j)
+			item (i, j) =
+			    item (i, j)
+			    ^ b.item (x + i, y + j);
+	return true;
+}

lib/mce.cpp

@@ -23,8 +23,7 @@ int mce::generate (pubkey&pub, privkey&priv, prng&rng, uint m, uint t)
 
 	//scramble matrix
 	matrix S;
-	S.generate_random_invertible (generator.height(), rng);
-	S.compute_inversion (priv.Sinv);
+	S.generate_random_with_inversion (generator.height(), priv.Sinv, rng);
 
 	//scramble permutation
 	permutation P;

lib/mce_oc.cpp

@@ -35,8 +35,7 @@ int mce_oc::generate (pubkey&pub, privkey&priv,
 
 	//scramble matrix
 	matrix S;
-	S.generate_random_invertible (g.height(), rng);
-	S.compute_inversion (priv.Sinv);
+	S.generate_random_with_inversion (g.height(), priv.Sinv, rng);
 
 	//scramble permutation
 	permutation P;

lib/nd.cpp

@@ -19,8 +19,7 @@ int nd::generate (pubkey&pub, privkey&priv, prng&rng, uint m, uint t)
 
 	//scrambler
 	matrix S;
-	S.generate_random_invertible (h.height(), rng);
-	S.compute_inversion (priv.Sinv);
+	S.generate_random_with_inversion (h.height(), priv.Sinv, rng);
 
 	//permutation
 	priv.Pinv.generate_random (h.width(), rng);

src/main.cpp

@@ -26,7 +26,7 @@ int main()
 
 	ccr::mce_oc::privkey priv;
 	ccr::mce_oc::pubkey pub;
-	ccr::mce_oc::generate (pub, priv, r, 7, 2, 8);
+	ccr::mce_oc::generate (pub, priv, r, 7, 2, 2);
 
 	priv.prepare();