From e42987061f46c12c700ce013b836f9f31c2b542c Mon Sep 17 00:00:00 2001
From: Jiri Pittner <jiri@pittnerovi.com>
Date: Mon, 1 Jan 2024 21:26:35 +0100
Subject: [PATCH] bitvector - some bugfixes and further implementations

---
 bitvector.cc | 182 ++++++++++++++++++++++++++++++++++++++++++++-------
 bitvector.h  |  21 ++++--
 laerror.cc   |  13 ++++
 laerror.h    |   2 +
 t.cc         |  47 ++++++++++++-
 5 files changed, 235 insertions(+), 30 deletions(-)

diff --git a/bitvector.cc b/bitvector.cc
index c6d3dd7..85361fc 100644
--- a/bitvector.cc
+++ b/bitvector.cc
@@ -39,21 +39,48 @@ return s;
 }
 
 
+void bitvector::zero_padding() const
+{
+if(!modulo) return;
+bitvector *p =  const_cast<bitvector *>(this);
+p->v[nn-1] &= (1ULL<<modulo)-1;
+}
+
+
+bitvector& bitvector::operator++()
+{
+copyonwrite();
+zero_padding();
+int i=0;
+while(i<nn) if(++v[i++]) break;
+return *this;
+}
+
+
+bitvector& bitvector::operator--()
+{
+copyonwrite();
+zero_padding();
+int i=0;
+while(i<nn) if(v[i++]--) break;
+return *this;
+}
+
+
+
 //implemented so that vectors of different length are considered different automatically
 bool bitvector::operator!=(const bitvector &rhs) const
 {
-if(nn!=rhs.nn || modulo!=rhs.modulo) return 1;
-if(v==rhs.v) return 0;
-if(!modulo) return memcmp(v,rhs.v,nn*sizeof(bitvector_block));
-if(memcmp(v,rhs.v,(nn-1)*sizeof(bitvector_block))) return 1;
-bitvector_block a=v[nn-1];
-bitvector_block b=rhs.v[nn-1];
-//zero out the irrelevant bits
-bitvector_block mask= ~((bitvector_block)0);
-mask <<=modulo;
-mask = ~mask;
-a&=mask; b&=mask;
-return a!=b;
+if(nn==rhs.nn && modulo==rhs.modulo && v==rhs.v) return false;
+zero_padding();
+rhs.zero_padding();
+int minnn=nn; if(rhs.nn<minnn) minnn=rhs.nn;
+int maxnn=nn; if(rhs.nn>maxnn) maxnn=rhs.nn;
+if(memcmp(v,rhs.v,minnn*sizeof(bitvector_block))) return true;
+if(minnn==maxnn) return false;
+if(nn==minnn) {for(int i=minnn; i<maxnn; ++i) if(rhs.v[i]) return true;}
+if(rhs.nn==minnn){for(int i=minnn; i<maxnn; ++i) if(v[i]) return true;}
+return false;
 }
 
 
@@ -97,6 +124,7 @@ bitvector bitvector::operator~() const
 {
 bitvector r((*this).size());
 for(int i=0; i<nn; ++i) r.v[i] = ~v[i];
+r.zero_padding();
 return r;
 }
 
@@ -197,12 +225,7 @@ void  bitvector::randomize()
 {
 copyonwrite();
 for(int i=0; i<nn; ++i) v[i]=RANDINT64();
-//zero the excess bits
-if(modulo) 
-	{
-	bitvector_block mask = (1ULL<<modulo)-1;
-	v[nn-1] &= mask;
-	}
+zero_padding();
 }
 
 
@@ -322,6 +345,7 @@ bitvector r(size());
 r.clear();
 bitvector tmp(*this);
 tmp.resize(size()+1,true);
+tmp.copyonwrite();
 int rd=rhs.degree();
 for(int i=0; i<=rd; ++i) //avoid making a working copy of rhs and shifting it
         {
@@ -358,7 +382,7 @@ while(!newr.is_zero())
 	t=newt; newt=remainder;
 	}
 
-if(r.degree()>0) laerror("field_inv: polynomial is not irreducible or input is its multiple");
+if(r.degree()>0) laerror("field_inv: polynomial is not irreducible or input is zero modulo the polynomial");
 if(!r[0]) laerror("zero in field_inv");
 
 return t;
@@ -368,12 +392,13 @@ return t;
 
 void bitvector::resize(const unsigned int n, bool preserve)
 {
-int old=size(); 
+if(preserve) zero_padding();
+int oldnn=nn; 
 NRVec<bitvector_block>::resize((n+blockbits-1)/blockbits,preserve); 
 modulo=n%blockbits; 
 if(preserve)  //clear newly allocated memory
 	{
-	for(int i=old; i<nn*blockbits; ++i) this->reset(i);
+	for(int i=oldnn; i<nn; ++i) v[i]=0;
 	} 
 else clear();
 }
@@ -430,6 +455,7 @@ return big;
 
 
 //cf. Brent & Zimmermann ANZMC08t (2008) paper
+//
 bool bitvector::is_irreducible() const
 {
 bitvector tmp(size());
@@ -452,7 +478,8 @@ for(auto p=f.begin(); p!=f.end(); ++p)
 	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;
+	//std::cout << "TEST tmp, ir, gcd, is_one "<<tmp<<" "<<*this<<" "<<g<<" : "<<g.is_one()<<std::endl;
+	if(!g.is_one()) return false;
 	}
 return true;
 }
@@ -508,6 +535,117 @@ if(dimensions)
 NRVec<bitvector_block>::put(fd,dimensions,transp);
 }
 
+static bitvector *tryme;
+static bool irfound;
+static int mynth;
+
+static void irfinder(int nones, int top)
+{
+if(irfound) return;
+if(nones==0) //terminate recursion
+	{
+	bool testit = tryme->is_irreducible();
+	//std::cout <<"candidate = "<< *tryme<< " result = "<<testit<<std::endl;
+	if(testit) 
+		{
+		--mynth;
+		if(!mynth) irfound=true;
+		}
+	return;
+	}
+for(int i=nones; i<=top; ++i) 
+	{
+	tryme->set(i);
+	irfinder(nones-1,i-1);
+	if(irfound) break;
+	else tryme->reset(i);
+	}
+}
+
+
+bitvector find_irreducible(int deg, int pop, int nth)
+{
+if(deg<=0) laerror("illegal degree in find_irreducible");
+if(deg==1) {bitvector r(2); r.set(1); r.reset(0); return r;}
+bitvector r(deg+1);
+if(pop== -1)
+	{
+	do	{
+		r.randomize();
+		r.set(0);
+		r.set(deg);
+		if((r.population()&1)==0) r.flip(1+RANDINT32()%(deg-2));
+		}
+	while(!r.is_irreducible());
+	return r;
+	}
+if(pop<3 || (pop&1)==0) laerror("impossible population of irreducible polynomial requested");
+r.clear();
+r.set(0);
+r.set(deg);
+pop-=2;
+tryme= &r;
+irfound=false;
+mynth=nth;
+irfinder(pop,deg-1);
+if(!irfound) r.clear();
+return r;
+}
+
+
+
+bitvector bitvector::pow(unsigned int n) const
+{
+if(n==0) {bitvector r(size()); r.clear(); r.set(0); return r;}
+if(n==1) return *this;
+bitvector y,z;
+z= *this;
+while(!(n&1))
+        {
+        z = z*z;
+        n >>= 1;
+        }
+y=z;
+while((n >>= 1)/*!=0*/)
+                {
+                z = z*z;
+                if(n&1) y *= z;
+                }
+return y;
+}
+
+bitvector bitvector::field_pow(unsigned int n, const bitvector &ir) const
+{
+if(n==0) {bitvector r(size()); r.clear(); r.set(0); return r;}
+if(n==1) return *this;
+bitvector y,z;
+z= *this;
+while(!(n&1))
+        {
+        z = z.field_mult(z,ir);
+        n >>= 1;
+        }
+y=z;
+while((n >>= 1)/*!=0*/)
+                {
+                z = z.field_mult(z,ir);
+                if(n&1) y = y.field_mult(z,ir);
+                }
+return y;
+}
+
+//sqrt(x) is x^(2^(d-1))
+bitvector bitvector::field_sqrt(const bitvector &ir) const
+{
+int d=ir.degree();
+bitvector r(*this); 
+for(int i=0; i<d-1; ++i) 
+	{
+	r= r.field_mult(r,ir);
+	}
+r.resize(d+1,true);
+return r;
+}
 
 
 
diff --git a/bitvector.h b/bitvector.h
index 8ea03fe..e3abdcf 100644
--- a/bitvector.h
+++ b/bitvector.h
@@ -43,7 +43,7 @@ private:
 	unsigned int modulo;
 public:
 	bitvector() : NRVec<bitvector_block>() {};
-	explicit bitvector (const unsigned int n):NRVec<bitvector_block>((n+blockbits-1)/blockbits) {modulo=n%blockbits;};
+	explicit bitvector (const unsigned int n):NRVec<bitvector_block>((n+blockbits-1)/blockbits) {modulo=n%blockbits; memset(v,0,nn*sizeof(bitvector_block));};
 	bitvector (const bitvector_block a, const unsigned int n):NRVec<bitvector_block>(a,(n+blockbits-1)/blockbits) {modulo=n%blockbits;};
 	bitvector(const bitvector &rhs) : NRVec<bitvector_block>(rhs) {modulo=rhs.modulo;};
 	//operator= seems to be correctly synthetized by the compiler
@@ -63,10 +63,13 @@ public:
 	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));};
-	bool iszero() const {for(int i=0; i<nn; ++i) if(v[i]) return false; return true;};
-	bool is_zero() const {return iszero();};
-	bool is_one() const {if(v[0]!=1) return false; for(int i=1; i<nn; ++i) if(v[i]) return false; return true;};
+	void zero_padding() const; 
+	bool is_zero() const {zero_padding(); for(int i=0; i<nn; ++i) if(v[i]) return false; return true;};
+	bool is_one() const {zero_padding(); if(v[0]!=1) return false; for(int i=1; i<nn; ++i) if(v[i]) return false;return true;};
+	bool iszero() const {return is_zero();};
 	void randomize();
+	bitvector& operator++();
+	bitvector& operator--();
 	bool operator!=(const bitvector &rhs) const;
 	bool operator==(const bitvector &rhs) const {return !(*this != rhs);};
 	bool operator>(const bitvector &rhs) const;
@@ -77,6 +80,9 @@ public:
 	bitvector& operator&=(const bitvector &rhs);
 	bitvector& operator|=(const bitvector &rhs);
 	bitvector& operator^=(const bitvector &rhs);
+	bitvector& operator&=(const bitvector_block rhs) {v[0]&=rhs; return *this;};
+        bitvector& operator|=(const bitvector_block rhs) {v[0]|=rhs; return *this;};
+        bitvector& operator^=(const bitvector_block rhs) {v[0]^=rhs; return *this;};
 	bitvector& operator+=(const bitvector &rhs) {return (*this)^=rhs;}; //addition modulo 2
 	bitvector& operator-=(const bitvector &rhs) {return (*this)^=rhs;}; //subtraction modulo 2
 	bitvector operator&(const bitvector &rhs) const {return bitvector(*this) &= rhs;};
@@ -87,10 +93,13 @@ public:
 	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 pow(unsigned int n) const; 
 	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;
+	bitvector field_pow(unsigned int n, const bitvector &irpolynom) const; 
+	bitvector field_sqrt(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);};
@@ -99,7 +108,7 @@ public:
 	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 population(const unsigned int before=0) const; //number of 1's (Hamming weight)
 	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);};
@@ -117,6 +126,8 @@ public:
 	void write(int fd, bool dimensions=1, bool transp=0);
 	};
 
+extern bitvector find_irreducible(int deg, int pop= -1, int nth=1); //degree and requested Hamming weight or -1 for random trial
+
 //expand to separate bytes or ints
 template <typename T>
 void bitvector_expand(const bitvector &v, NRVec<T> &r)
diff --git a/laerror.cc b/laerror.cc
index 18830c0..7b6186c 100644
--- a/laerror.cc
+++ b/laerror.cc
@@ -205,4 +205,17 @@ extern "C" int cblas_errprn(int ierr, int info, char *form, ...) {
 }
 #endif
 
+
+int endianity()
+{
+static union {int i; char z[sizeof(int)];} u;
+static int jj= -1;
+if(jj<0) /* first call */
+                {
+                u.i=1;
+                jj=0; while(!u.z[jj]) jj++;
+                }
+return jj;
+}
+
 }//namespace
diff --git a/laerror.h b/laerror.h
index 2d24489..63d150c 100644
--- a/laerror.h
+++ b/laerror.h
@@ -53,6 +53,8 @@ s << x.msg;
 return s;
 }
 
+extern int endianity();
+
 }//namespace
 
 #endif
diff --git a/t.cc b/t.cc
index 7e54b08..c53aa97 100644
--- a/t.cc
+++ b/t.cc
@@ -2946,7 +2946,7 @@ cin >>n;
 cout <<factorization(n)<<" phi = "<<eulerphi(n)<<endl;
 }
 
-if(1)
+if(0)
 {
 bitvector ir; cin >>ir;
 if(!ir.is_irreducible()) laerror("input must be an irreducible polynomial");
@@ -2961,8 +2961,49 @@ bitvector cc=a.field_composition(ai,ir);
 cout <<c<<endl;
 cout <<c%ir<<endl;
 cout <<cc<<endl;
+}
+
+if(1)
+{
+int seed;
+int f=open("/dev/random",O_RDONLY);
+if(sizeof(int)!=read(f,&seed,sizeof(int))) laerror("cannot read /dev/random");
+close(f);
+srand(seed);
+
+int d;
+cin >>d;
+bitvector ir=find_irreducible(d,3);
+if(ir.is_zero()) ir=find_irreducible(d,5);
+if(ir.is_zero()) ir=find_irreducible(d,7);
+if(ir.is_zero()) {cout<<"cannot find IR polynomial\n"; exit(1);}
+
+
+cout <<"IR = "<<ir<<endl;
+if(!ir.is_irreducible()) laerror("error in IR finder");
+else cout <<"IS irreducible\n";
+bitvector a(d);
+do { a.randomize(); }while(a.is_zero());
+cout <<"A = "<<a<<endl;
+bitvector aa=a;
+cout <<"A+1 = "<<++aa<<endl;
+cout <<"A+1-1 = "<<--aa<<endl;
+if(a!=aa) laerror("error in inc/dec");
+bitvector g=ir.gcd(a);
+cout<<"GCD A,IR = "<<g<<endl;
+if(!g.is_one() && !a.is_zero()) laerror("ERROR IN is_irreducible\n");
+bitvector ai = a.field_inv(ir);
+cout <<"I = "<<ai<<endl;
+bitvector check = a.field_mult(ai,ir);
+cout<<"check " <<check<<endl;
+if(!check.is_one()) laerror("error in GF(2^n) inversion");
+bitvector s=a.field_sqrt(ir);
+cout <<"test a "<<a<<endl;
+cout <<"tests2 "<<s.field_mult(s,ir)<<endl;
+cout <<"test s "<<s<<endl;
+bitvector dif=a - s.field_mult(s,ir);
+if(!dif.is_zero()) laerror("error in gf 2^n sqrt");
+}
 
 
-}
-
 }