polynomial irreducibility test in GF2

bitvector.cc

@@ -18,6 +18,7 @@
 
 #include "bitvector.h"
 #include <unistd.h>
+#include "numbers.h"
 
 namespace LA {
 
@@ -425,11 +426,68 @@ do      {
         }
 while(! small.is_zero());
 return big;
-
-
-
 }
 
+
+//cf. Brent & Zimmermann ANZMC08t (2008) paper
+bool bitvector::is_irreducible() const
+{
+bitvector tmp(size());
+tmp.clear(); tmp.set(1);
+unsigned int d=degree();
+
+//repeated squaring test
+for(unsigned int j=0; j<d; ++j) tmp = tmp.field_mult(tmp,*this);
+tmp.flip(1);
+if(!tmp.is_zero()) return false;
+
+FACTORIZATION<uint64_t> f = factorization((uint64_t)d);
+if(f.begin()->first==d) return true; //d was prime
+
+//additional tests needed for non-prime degrees
+for(auto p=f.begin(); p!=f.end(); ++p)
+	{
+	unsigned int dm= d / p->first;
+	tmp.clear(); tmp.set(1);
+	for(unsigned int j=0; j<dm; ++j) tmp = tmp.field_mult(tmp,*this);
+	tmp.flip(1);
+	bitvector g=tmp.gcd(*this);
+	if(!g,is_one()) return false;
+	}
+return true;
+}
+
+
+//horner scheme
+bitvector bitvector::composition(const bitvector &x) const
+{
+bitvector r(size());
+r.clear();
+int d=degree();
+for(int i=d; i>0; --i)
+	{
+	if((*this)[i]) r.flip(0);
+	r*=x;
+	}
+if((*this)[0]) r.flip(0);
+return r;
+}
+
+bitvector bitvector::field_composition(const bitvector &x, const bitvector &ir) const
+{
+bitvector r(size());
+r.clear();
+int d=degree();
+for(int i=d; i>0; --i)
+	{
+	if((*this)[i]) r.flip(0);
+	r= r.field_mult(x,ir);
+	}
+if((*this)[0]) r.flip(0);
+return r;
+}
+
+
 void bitvector::read(int fd, bool dimensions, bool transp)
 {
 if(dimensions) 

bitvector.h

@@ -59,6 +59,7 @@ public:
 	int getblocksize() const {return 8*sizeof(bitvector_block);};
 	void set(const unsigned int i) {v[i/blockbits] |= (1UL<<(i%blockbits));};
 	void reset(const unsigned int i) {v[i/blockbits] &= ~(1UL<<(i%blockbits));};
+	void flip(const unsigned int i) {v[i/blockbits] ^= (1UL<<(i%blockbits));};
 	const bool assign(const unsigned int i, const bool r) {if(r) set(i); else reset(i); return r;};
 	void clear() {copyonwrite(true); memset(v,0,nn*sizeof(bitvector_block));};
 	void fill() {memset(v,0xff,nn*sizeof(bitvector_block));};
@@ -85,18 +86,23 @@ public:
 	bitvector operator-(const bitvector &rhs) const {return *this ^ rhs;}; //subtraction modulo 2
 	bitvector multiply(const bitvector &rhs, bool autoresize=true) const; //use autoresize=false only if you know it will not overflow!
 	bitvector operator*(const bitvector &rhs) const {return multiply(rhs,true);}  //multiplication of polynomials over GF(2) NOTE: naive algorithm, does not employ CLMUL nor fft-like approach, only for short vectors!!!
+	bitvector&  operator*=(const bitvector &rhs) {*this = (*this)*rhs; return *this;}
 	bitvector field_mult(const bitvector &rhs, const bitvector &irpolynom) const; //multiplication in GF(2^n)
 	bitvector field_inv(const bitvector &irpolynom) const; //multiplication in GF(2^n)
 	bitvector field_div(const bitvector &rhs, const bitvector &irpolynom) const {return field_mult(rhs.field_inv(irpolynom),irpolynom);};
+	bitvector field_composition(const bitvector &rhs, const bitvector &irpolynom) const;
+	bool is_irreducible() const; //test irreducibility of polynomial over GF2
 	bitvector division(const bitvector &rhs,  bitvector &remainder) const;
 	bitvector operator/(const bitvector &rhs) const {bitvector rem(rhs.size()); return division(rhs,rem);};
 	bitvector operator%(const bitvector &rhs) const {bitvector rem(rhs.size()); division(rhs,rem); return rem;};
 	bitvector gcd(const bitvector &rhs) const; //as a polynomial over GF2
 	bitvector lcm(const bitvector &rhs) const {return (*this)*rhs/this->gcd(rhs);};
+	bitvector composition(const bitvector &rhs) const;
         unsigned int bitdiff(const bitvector &y) const; //number of differing bits (Hamming distance)
 	unsigned int population(const unsigned int before=0) const; //number of 1's
 	unsigned int nlz() const; //number of leading zeroes
 	unsigned int degree() const {if(iszero()) return 0; else return size()-nlz()-1;}; //interprested as a polynomial over GF(2)
+	void truncate(int t=0) {int s=degree()+1; if(t>s) s=t;  resize(s,true);};
 	unsigned int ntz() const;  //number of trailing zeroes
 	//extended, truncated const i.e. not on *this but return new entity, take care of modulo's bits
 	//logical shifts

polynomial.cc

@@ -249,6 +249,12 @@ resize(n,false);
 for(int i=0; i<=n; ++i) (*this)[i] = (T) binom(n,i);
 }
 
+template <typename T>
+Polynomial<T>  Polynomial<T>::composition(const Polynomial &rhs) const
+{
+return value(*this,rhs);
+}
+
 
 /***************************************************************************//**
  * forced instantization in the corresponding object file

polynomial.h

@@ -147,6 +147,7 @@ public:
 			}
 		return r;
 		}
+	Polynomial composition(const Polynomial &rhs) const;
 	Polynomial even_powers() const {int d=degree()/2; Polynomial r(d); for(int i=0; i<=degree(); i+=2) r[i/2] = (*this)[i]; return r;};
 	Polynomial odd_powers() const {int d=(degree()-1)/2; Polynomial r(d); if(degree()==0) {r[0]=0; return r;} for(int i=1; i<=degree(); i+=2) r[(i-1)/2] = (*this)[i]; return r;};
 	void polydiv(const Polynomial &rhs, Polynomial &q, Polynomial &r) const;

t.cc

@@ -2949,12 +2949,19 @@ cout <<factorization(n)<<" phi = "<<eulerphi(n)<<endl;
 if(1)
 {
 bitvector ir; cin >>ir;
+if(!ir.is_irreducible()) laerror("input must be an irreducible polynomial");
 bitvector a; cin >>a;
 bitvector ai = a.field_inv(ir);
 cout<< "inverse = "<<ai<<endl;
 cout<<"check1 " <<(a*ai)%ir<<endl;
 cout<<"check2 " <<a.field_mult(ai,ir)<<endl;
 
+bitvector c=a.composition(ai);
+bitvector cc=a.field_composition(ai,ir);
+cout <<c<<endl;
+cout <<c%ir<<endl;
+cout <<cc<<endl;
+
 
 }