gf2m arithmetic phase

include/codecrypt.h

@@ -14,7 +14,13 @@
 namespace ccr
 {
 
+/*
+ * typedefs. uint and sint should be able to comfortably hold the field
+ * elements of underlying calculations (esp. with polynomials. Switching to
+ * 64bits is adviseable when computing with n=64K and larger.
+ */
 typedef unsigned int uint;
+typedef int sint;
 
 /*
  * vector over GF(2). We rely on STL's vector<bool> == bit_vector
@@ -85,22 +91,42 @@ public:
 };
 
 /*
- * polynomial over GF(2) is effectively a vector with a_n binary values
+ * galois field of 2^m elements. Stored in an integer, for convenience.
+ */
+
+class gf2m
+{
+public:
+	uint poly;
+	uint n, m;
+
+	bool create (uint m);
+
+	uint add (uint, uint);
+	uint mult (uint, uint);
+	uint exp (uint, int);
+	uint inv (uint);
+};
+
+/*
+ * polynomial over GF(2^m) is effectively a vector with a_n binary values
  * with some added operations.
  */
-class polynomial : public bvector
+class polynomial : public std::vector<uint>
 {
+protected:
+	_ccr_declare_vector_item
 public:
 	void strip();
 	int degree() const;
 	bool zero() const;
-	void add (const polynomial&);
-	void mod (const polynomial&);
-	void mult (const polynomial&);
-	polynomial gcd (polynomial);
-	bool is_irreducible();
-	void generate_random_irreducible (uint n, prng&);
-	void compute_mod_squaring_matrix (matrix&);
+	void add (const polynomial&, gf2m&);
+	void mod (const polynomial&, gf2m&);
+	void mult (const polynomial&, gf2m&);
+	polynomial gcd (polynomial, gf2m&);
+	bool is_irreducible (gf2m&);
+	void generate_random_irreducible (uint s, gf2m&, prng&);
+	void compute_square_root_matrix (std::vector<polynomial>&, gf2m&);
 };
 
 /*

lib/gf2m.cpp

@@ -0,0 +1,118 @@
+
+#include "codecrypt.h"
+
+using namespace ccr;
+
+#include <iostream>
+using namespace std;
+
+/*
+ * helpful stuff for arithmetic in GF(2^m) - polynomials over GF(2).
+ */
+
+int gf2p_degree (uint p)
+{
+	int r = -1;
+	for (int i = 0; p; p >>= 1, ++i) r = i;
+	return r;
+}
+
+uint gf2p_mod (uint a, uint p)
+{
+	if (!p) return 0;
+	int t, degp = gf2p_degree (p);
+	while ( (t = gf2p_degree (a) ) >= degp)
+		a ^= p << (t - degp);
+	return a;
+}
+
+uint gf2p_gcd (uint a, uint b)
+{
+	uint c;
+	if (!a) return b;
+	while (b) {
+		c = gf2p_mod (a, b);
+		a = b;
+		b = c;
+	}
+	return a;
+}
+
+uint gf2p_modmult (uint a, uint b, uint p)
+{
+	a = gf2p_mod (a, p);
+	b = gf2p_mod (b, p);
+	uint r = 0;
+	uint d = 1 << gf2p_degree (p);
+	while (a) {
+		if (a & 1) r ^= b;
+		a >>= 1;
+		b <<= 1;
+		if (b <= d) b ^= p;
+	}
+	return r;
+}
+
+bool is_irreducible_gf2_poly (uint p)
+{
+	if (!p) return false;
+	int d = gf2p_degree (p) / 2;
+	uint test = 2; //x^1+0
+	for (int i = 0; i < d; ++i) {
+		test = gf2p_modmult (test, test, p);
+
+		if (gf2p_gcd (test ^ 2 /* test - x^1 */, p) != 1)
+			return false;
+	}
+	return true;
+}
+
+bool gf2m::create (uint M)
+{
+	if (M < 1) return false; //too small.
+	m = M;
+	n = 1 << m;
+	if (!n) return false; //too big.
+	for (uint t = 1 + (1 << m), e = 1 << (1 + m); t < e; t += 2)
+		if (is_irreducible_gf2_poly (t) ) {
+			poly = t;
+			return true;
+		}
+	return false;
+}
+
+/*
+uint gfn_mult(uint a, uint b, uint n)
+{
+	uint irp=0;
+	while(n) { irp=(irp<<1)|1; n>>=1;}
+	uint r=a*b;
+	//TODO probably move this to own file
+}
+
+uint gfn_inv (uint a, uint n);
+
+uint gfn_exp (uint a, sint k, uint n)
+{
+	if (!a) return 0;
+	if (a == 1) return 1;
+	if (k < 0) {
+		a = gfn_inv (a, n);
+		k = -k;
+	}
+	uint r = 1;
+	while (k) {
+		if (k & 1) r=gfn_mult(r,a,n);
+		a=gfn_mult(a,a,n);
+		k >>= 2;
+	}
+	return r;
+}
+
+uint gfn_inv (uint a, uint n)
+{
+	if (n == 2) return a;
+	return gfn_exp (a, ( (sint) n) - 2, n);
+}
+
+*/

lib/polynomial.cpp

@@ -1,17 +1,16 @@
 
+#if 0
 #include "codecrypt.h"
 
 using namespace ccr;
 
-#if 0
 #include <iostream>
 using namespace std;
 void dump (const polynomial&t)
 {
-	for (uint i = 0; i < t.size(); ++i) cout << t[i];
+	for (uint i = 0; i < t.size(); ++i) cout << t[i] << ' ';
 	cout << endl;
 }
-#endif
 
 int polynomial::degree() const
 {
@@ -31,26 +30,37 @@ bool polynomial::zero() const
 	return true;
 }
 
-void polynomial::add (const polynomial&f)
+void polynomial::add (const polynomial&f, gf2m&fld)
 {
 	int df = f.degree();
 	if (df > degree() ) resize (df + 1);
 	for (int i = 0; i <= df; ++i) item (i) = item (i) ^ f[i];
 }
 
-void polynomial::mod (const polynomial&f)
+void polynomial::mod (const polynomial&f, gf2m&fld)
 {
 	int df = f.degree();
 	int d;
+	uint hi = fld.inv (f[df]);
+	cout << "mod by inv " << hi << endl;
+	dump (*this);
+	dump (f);
 	// while there's place to substract, reduce by x^(d-df)-multiply of f
-	while ( (d = degree() ) >= df) {
-		for (int i = 0; i <= df; ++i)
-			item (i + d - df) = item (i + d - df) ^ f[i];
-	}
+	for (d = degree(); d >= df; --d)
+		if (item (d) ) {
+			uint t = fld.mult (item (d), hi);
+			cout << "mult " << t << endl;
+			for (int i = 0; i <= df; ++i)
+				item (i + d - df) = fld.add (item (i + d - df)
+				                             , fld.mult (t, f[i]) );
+			cout << "now ";
+			dump (*this);
+		}
+	cout << "end mod" << endl;
 	strip();
 }
 
-void polynomial::mult (const polynomial&b)
+void polynomial::mult (const polynomial&b, gf2m&fld)
 {
 	polynomial a = *this;
 	clear();
@@ -60,10 +70,11 @@ void polynomial::mult (const polynomial&b)
 	resize (da + db + 1, 0);
 	for (i = 0; i <= da; ++i)
 		if (a[i]) for (j = 0; j <= db; ++j)
-				item (i + j) = item (i + j) ^ b[j];
+				item (i + j) = fld.add (item (i + j),
+				                        fld.mult (a[i], b[j]) );
 }
 
-polynomial polynomial::gcd (polynomial b)
+polynomial polynomial::gcd (polynomial b, gf2m&fld)
 {
 	polynomial a = *this;
 
@@ -71,15 +82,17 @@ polynomial polynomial::gcd (polynomial b)
 	if (a.degree() < 0) return b;
 	for (;;) {
 		if (b.zero() ) return a;
-		a.mod (b);
+		dump (a);
+		a.mod (b, fld);
 		if (a.zero() ) return b;
-		b.mod (a);
+		dump (b);
+		b.mod (a, fld);
 	}
 	//unreachable
 	return polynomial();
 }
 
-bool polynomial::is_irreducible()
+bool polynomial::is_irreducible (gf2m&fld)
 {
 	//Ben-Or irreducibility test
 	polynomial xi; //x^(2^i) in our case
@@ -89,36 +102,41 @@ bool polynomial::is_irreducible()
 	xmodf[0] = 0;
 	xmodf[1] = 1; //x
 	xi = xmodf;
-	xmodf.mod (*this); //mod f
+	xmodf.mod (*this, fld); //mod f
 
-	uint n = degree();
-	for (uint i = 1; i <= n / 2; ++i) {
+	uint d = degree();
+	for (uint i = 1; i <= d / 2; ++i) {
 		t = xi;
-		t.mult (xi); //because mult would destroy xi on xi.mult(xi)
-		t.mod (*this);
+		t.mult (xi, fld); //because mult would destroy xi on xi.mult(xi)
+		t.mod (*this, fld);
 		xi = t;
-		t.add (xmodf);
+		t.add (xmodf, fld);
 
-		t = t.gcd (*this);
+		t = t.gcd (*this, fld);
 		if (t.degree() != 0) //gcd(f,x^2^i - x mod f) != 1
 			return false;
 	}
 	return true;
 }
 
-void polynomial::generate_random_irreducible (uint s, prng & rng)
+void polynomial::generate_random_irreducible (uint s, gf2m&fld, prng & rng)
 {
 	resize (s + 1);
 	item (s) = 1; //degree s
-	item (0) = 1; //not divisible by x^1
-	for (uint i = 1; i < s; ++i) item (i) = rng.random (2);
-	while (!is_irreducible() ) {
-		uint pos = 1 + rng.random (s - 1);
-		item (pos) = !item (pos);
+	item (0) = 1 + rng.random (fld.n - 1); //not divisible by x^1
+	for (uint i = 1; i < s; ++i) item (i) = rng.random (fld.n);
+	cout << "start ";
+	dump (*this);
+	while (!is_irreducible (fld) ) {
+		cout << "retry ";
+		dump (*this);
+		uint pos = rng.random (s);
+		item (pos) = pos == 0 ?
+		             (1 + rng.random (fld.n - 1) ) : rng.random (fld.n);
 	}
 }
 
-void polynomial::compute_mod_squaring_matrix (matrix&r)
+void polynomial::compute_square_root_matrix (vector<polynomial>&r, gf2m&fld)
 {
 	int d = degree();
 	if (d < 0) return;
@@ -129,9 +147,12 @@ void polynomial::compute_mod_squaring_matrix (matrix&r)
 		col.resize (i + 1, 0);
 		col[i] = 1;
 		t = col;
-		col.mult (t);
-		col.mod (*this);
+		col.mult (t, fld);
+		col.mod (*this, fld);
 		col.resize (d, 0);
 		r[i] = col;
 	}
+
+	//TODO gauss
 }
+#endif