From 53eb902f1b2bf6be584fe2a69c0e9f02813f5b14 Mon Sep 17 00:00:00 2001
From: Mirek Kratochvil <exa.exa@gmail.com>
Date: Mon, 2 Apr 2012 11:14:54 +0200
Subject: [PATCH] matrix operations

---
 include/codecrypt.h |  16 ++++++-
 lib/bvector.cpp     |  16 +++++++
 lib/matrix.cpp      | 106 +++++++++++++++++++++++++++++++++++++++-----
 lib/polynomial.cpp  |   2 +-
 4 files changed, 127 insertions(+), 13 deletions(-)

diff --git a/include/codecrypt.h b/include/codecrypt.h
index e1599e0..c41aec2 100644
--- a/include/codecrypt.h
+++ b/include/codecrypt.h
@@ -26,6 +26,8 @@ protected:
 	_ccr_declare_vector_item
 public:
 	uint hamming_weight();
+	void add (const bvector&);
+	bool operator* (const bvector&); //dot product
 };
 
 /*
@@ -47,12 +49,22 @@ class matrix : public std::vector<bvector>
 protected:
 	_ccr_declare_vector_item
 public:
-	matrix operator* (const matrix&);
+	uint width() const {
+		return size();
+	}
 
+	uint height() const {
+		if (size() ) return item (0).size();
+		return 0;
+	}
+
+	matrix operator* (const matrix&);
+	void mult (const matrix&);
+
+	void compute_transpose (matrix&);
 	bool compute_inversion (matrix&);
 	void generate_random_invertible (uint, prng&);
 	void unit (uint);
-	void compute_transpose (matrix&);
 };
 
 /*
diff --git a/lib/bvector.cpp b/lib/bvector.cpp
index 635c186..8039b6d 100644
--- a/lib/bvector.cpp
+++ b/lib/bvector.cpp
@@ -9,3 +9,19 @@ uint bvector::hamming_weight()
 	return r;
 }
 
+void bvector::add (const bvector&a)
+{
+	if (a.size() > size() ) resize (a.size(), 0);
+	for (uint i = 0; i < size(); ++i)
+		item (i) = item (i) ^ a[i];
+}
+
+bool bvector::operator* (const bvector&a)
+{
+	bool r = 0;
+	uint s = size(), i;
+	if (s > a.size() ) s = a.size();
+	for (i = 0; i < s; ++i) r ^= (item (i) &a[i]);
+	return r;
+}
+
diff --git a/lib/matrix.cpp b/lib/matrix.cpp
index 0e9446a..4492659 100644
--- a/lib/matrix.cpp
+++ b/lib/matrix.cpp
@@ -5,21 +5,107 @@ using namespace ccr;
 
 void matrix::unit (uint size)
 {
-
+	clear();
+	resize (size);
+	for (uint i = 0; i < size; ++i) {
+		item (i).resize (size, 0);
+		item (i) [i] = 1;
+	}
 }
 
-bool matrix::compute_inversion (matrix&r)
+matrix matrix::operator* (const matrix&a)
 {
-
-	return false;
-}
-
-void matrix::generate_random_invertible (uint size, prng&rng)
-{
-
+	matrix r = *this;
+	r.mult (a);
+	return r;
 }
 
 void matrix::compute_transpose (matrix&r)
 {
-
+	uint h = height(), w = width(), i, j;
+	r.resize (h);
+	for (i = 0; i < h; ++i) {
+		r[i].resize (w);
+		for (j = 0; j < w; ++j) r[i][j] = item (j) [i];
+	}
 }
+
+void matrix::mult (const matrix&right)
+{
+	//trivial multiply. TODO strassen algo for larger matrices.
+	matrix leftT;
+	compute_transpose (leftT);
+	uint w = right.width(), h = leftT.width(), i, j;
+	resize (w);
+	for (i = 0; i < w; ++i) {
+		item (i).resize (h);
+		for (j = 0; j < h; ++j) item (i) [j] = leftT[j] * right[i];
+	}
+}
+
+bool matrix::compute_inversion (matrix&res)
+{
+	//gauss-jordan elimination with inversion of the second matrix.
+	//we are computing with transposed matrices for simpler row ops
+
+	uint s = width();
+	if (s != height() ) return false;
+	matrix m, r;
+	r.unit (s);
+	this->compute_transpose (m);
+
+	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].add (m[j]);
+			r[i].add (r[j]);
+		}
+		//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]);
+			}
+
+	r.compute_transpose (res);
+	return true;
+}
+
+void matrix::generate_random_invertible (uint size, prng & rng)
+{
+	matrix lt, ut;
+	uint i, j;
+	// random lower triagonal
+	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 triagonal
+	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);
+	}
+	lt.mult (ut);
+	// permute
+	permutation p;
+	p.generate_random (size, rng);
+	p.permute (lt, *this);
+}
+
diff --git a/lib/polynomial.cpp b/lib/polynomial.cpp
index 8bc926f..7481a1a 100644
--- a/lib/polynomial.cpp
+++ b/lib/polynomial.cpp
@@ -95,7 +95,7 @@ bool polynomial::is_irreducible()
 	for (uint i = 1; i <= n / 2; ++i) {
 		t = xi;
 		t.mult (xi); //because mult would destroy xi on xi.mult(xi)
-		t.mod(*this);
+		t.mod (*this);
 		xi = t;
 		t.add (xmodf);