diff --git a/include/codecrypt.h b/include/codecrypt.h
index 5b56da3..d7e2cbc 100644
--- a/include/codecrypt.h
+++ b/include/codecrypt.h
@@ -4,6 +4,13 @@
 
 #include <vector>
 
+//STL wraparound, because writing (*this)[i] everywhere is clumsy
+#define _ccr_declare_vector_item \
+	inline reference item(size_type n) \
+		{ return (*this)[n]; }; \
+	inline const_reference item(size_type n) const \
+		{ return (*this)[n]; };
+
 namespace ccr
 {
 
@@ -15,10 +22,8 @@ typedef unsigned int uint;
  */
 class bvector : public std::vector<bool>
 {
-	//STL wraparound, because writing (*this)[i] is clumsy
-	inline reference item (size_type n) {
-		return (*this) [n];
-	}
+protected:
+	_ccr_declare_vector_item
 public:
 	uint hamming_weight();
 };
@@ -39,9 +44,8 @@ public:
  */
 class matrix : public std::vector<bvector>
 {
-	inline reference item (size_type n) {
-		return (*this) [n];
-	}
+protected:
+	_ccr_declare_vector_item
 public:
 	matrix operator* (const matrix&);
 
@@ -57,9 +61,8 @@ public:
  */
 class permutation : public std::vector<uint>
 {
-	inline reference item (size_type n) {
-		return (*this) [n];
-	}
+protected:
+	_ccr_declare_vector_item
 public:
 	void compute_inversion (permutation&);
 
@@ -75,13 +78,16 @@ public:
  */
 class polynomial : public bvector
 {
-	inline reference item (size_type n) {
-		return (*this) [n];
-	}
 public:
+	void strip();
+	uint degree() 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&);
+	vector operator<< (uint);
 };
 
 /*
diff --git a/lib/polynomial.cpp b/lib/polynomial.cpp
index 34dc41e..7a69f05 100644
--- a/lib/polynomial.cpp
+++ b/lib/polynomial.cpp
@@ -3,14 +3,100 @@
 
 using namespace ccr;
 
+uint polynomial::degree() const
+{
+	uint r = -1;
+	for (uint i = 0; i < size(); ++i) if (item (i) ) r = i;
+	return r;
+}
+
+void polynomial::strip()
+{
+	resize (degree() + 1);
+}
+
+void polynomial::add (const polynomial&f)
+{
+	uint df = f.degree();
+	if (df > degree() ) resize (df + 1);
+	for (uint i = 0; i <= df; ++i) item (i) = item (i) ^ f[i];
+}
+
+void polynomial::mod (const polynomial&f)
+{
+	uint df = f.degree();
+	uint d;
+	// while there's place to substract, reduce by x^(d-df)-multiply of f
+	while ( (d = degree() ) >= df) {
+		for (uint i = 0; i <= df; ++i)
+			item (i + d - df) = item (i + d - df) ^ f[i];
+	}
+	strip();
+}
+
+void polynomial::mult (const polynomial&b)
+{
+	polynomial a = *this;
+	clear();
+	uint i, j, da, db;
+	da = a.degree();
+	db = b.degree();
+	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];
+}
+
+polynomial polynomial::gcd (polynomial b)
+{
+	polynomial a = *this;
+
+	//eukleides
+	if (a.degree() < 0) return b;
+	for (;;) {
+		if (b.degree() < 0) return a;
+		a.mod (b);
+		if (a.degree() < 0) return b;
+		b.mod (a);
+	}
+	//unreachable
+	return polynomial();
+}
+
 bool polynomial::is_irreducible()
 {
+	//Ben-Or irreducibility test
+	polynomial xi; //x^(2^i) in our case
+	polynomial xmodf, t;
 
-	return false;
+	xmodf.resize (2); //precompute (x mod f)
+	xmodf[0] = 0;
+	xmodf[1] = 1; //x
+	xmodf.mod (*this); //mod f
+
+	uint n = degree();
+	for (uint i = 1; i <= n / 2; ++i) {
+		t = xi;
+		t.mult (xi); //because mult would destroy xi on xi.mult(xi)
+		xi = t;
+		t.add (xmodf);
+
+		t = t.gcd (*this);
+		if (t.degree() != 0) //gcd(f,x^2^i - x mod f) != 1
+			return false;
+	}
+	return true;
 }
 
-void polynomial::generate_random_irreducible (uint size, prng&rng)
+void polynomial::generate_random_irreducible (uint s, 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);
+	}
 }