bitvector - some bugfixes and further implementations

This commit is contained in:
Jiri Pittner 2024-01-01 21:26:35 +01:00
parent 9bceebdd29
commit e42987061f
5 changed files with 235 additions and 30 deletions

View File

@ -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 //implemented so that vectors of different length are considered different automatically
bool bitvector::operator!=(const bitvector &rhs) const bool bitvector::operator!=(const bitvector &rhs) const
{ {
if(nn!=rhs.nn || modulo!=rhs.modulo) return 1; if(nn==rhs.nn && modulo==rhs.modulo && v==rhs.v) return false;
if(v==rhs.v) return 0; zero_padding();
if(!modulo) return memcmp(v,rhs.v,nn*sizeof(bitvector_block)); rhs.zero_padding();
if(memcmp(v,rhs.v,(nn-1)*sizeof(bitvector_block))) return 1; int minnn=nn; if(rhs.nn<minnn) minnn=rhs.nn;
bitvector_block a=v[nn-1]; int maxnn=nn; if(rhs.nn>maxnn) maxnn=rhs.nn;
bitvector_block b=rhs.v[nn-1]; if(memcmp(v,rhs.v,minnn*sizeof(bitvector_block))) return true;
//zero out the irrelevant bits if(minnn==maxnn) return false;
bitvector_block mask= ~((bitvector_block)0); if(nn==minnn) {for(int i=minnn; i<maxnn; ++i) if(rhs.v[i]) return true;}
mask <<=modulo; if(rhs.nn==minnn){for(int i=minnn; i<maxnn; ++i) if(v[i]) return true;}
mask = ~mask; return false;
a&=mask; b&=mask;
return a!=b;
} }
@ -97,6 +124,7 @@ bitvector bitvector::operator~() const
{ {
bitvector r((*this).size()); bitvector r((*this).size());
for(int i=0; i<nn; ++i) r.v[i] = ~v[i]; for(int i=0; i<nn; ++i) r.v[i] = ~v[i];
r.zero_padding();
return r; return r;
} }
@ -197,12 +225,7 @@ void bitvector::randomize()
{ {
copyonwrite(); copyonwrite();
for(int i=0; i<nn; ++i) v[i]=RANDINT64(); for(int i=0; i<nn; ++i) v[i]=RANDINT64();
//zero the excess bits zero_padding();
if(modulo)
{
bitvector_block mask = (1ULL<<modulo)-1;
v[nn-1] &= mask;
}
} }
@ -322,6 +345,7 @@ bitvector r(size());
r.clear(); r.clear();
bitvector tmp(*this); bitvector tmp(*this);
tmp.resize(size()+1,true); tmp.resize(size()+1,true);
tmp.copyonwrite();
int rd=rhs.degree(); int rd=rhs.degree();
for(int i=0; i<=rd; ++i) //avoid making a working copy of rhs and shifting it 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; 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"); if(!r[0]) laerror("zero in field_inv");
return t; return t;
@ -368,12 +392,13 @@ return t;
void bitvector::resize(const unsigned int n, bool preserve) 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); NRVec<bitvector_block>::resize((n+blockbits-1)/blockbits,preserve);
modulo=n%blockbits; modulo=n%blockbits;
if(preserve) //clear newly allocated memory 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(); else clear();
} }
@ -430,6 +455,7 @@ return big;
//cf. Brent & Zimmermann ANZMC08t (2008) paper //cf. Brent & Zimmermann ANZMC08t (2008) paper
//
bool bitvector::is_irreducible() const bool bitvector::is_irreducible() const
{ {
bitvector tmp(size()); 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); for(unsigned int j=0; j<dm; ++j) tmp = tmp.field_mult(tmp,*this);
tmp.flip(1); tmp.flip(1);
bitvector g=tmp.gcd(*this); 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; return true;
} }
@ -508,6 +535,117 @@ if(dimensions)
NRVec<bitvector_block>::put(fd,dimensions,transp); 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;
}

View File

@ -43,7 +43,7 @@ private:
unsigned int modulo; unsigned int modulo;
public: public:
bitvector() : NRVec<bitvector_block>() {}; 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_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;}; bitvector(const bitvector &rhs) : NRVec<bitvector_block>(rhs) {modulo=rhs.modulo;};
//operator= seems to be correctly synthetized by the compiler //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;}; 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 clear() {copyonwrite(true); memset(v,0,nn*sizeof(bitvector_block));};
void fill() {memset(v,0xff,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;}; void zero_padding() const;
bool is_zero() const {return iszero();}; bool is_zero() const {zero_padding(); for(int i=0; i<nn; ++i) if(v[i]) return false; return true;};
bool is_one() const {if(v[0]!=1) return false; for(int i=1; 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(); void randomize();
bitvector& operator++();
bitvector& operator--();
bool operator!=(const bitvector &rhs) const; bool operator!=(const bitvector &rhs) const;
bool operator==(const bitvector &rhs) const {return !(*this != rhs);}; bool operator==(const bitvector &rhs) const {return !(*this != rhs);};
bool operator>(const bitvector &rhs) const; 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 &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;}; //addition modulo 2
bitvector& operator-=(const bitvector &rhs) {return (*this)^=rhs;}; //subtraction modulo 2 bitvector& operator-=(const bitvector &rhs) {return (*this)^=rhs;}; //subtraction modulo 2
bitvector operator&(const bitvector &rhs) const {return bitvector(*this) &= rhs;}; 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 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) 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& 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_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_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_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_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 bool is_irreducible() const; //test irreducibility of polynomial over GF2
bitvector division(const bitvector &rhs, bitvector &remainder) const; 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()); return division(rhs,rem);};
@ -99,7 +108,7 @@ public:
bitvector lcm(const bitvector &rhs) const {return (*this)*rhs/this->gcd(rhs);}; bitvector lcm(const bitvector &rhs) const {return (*this)*rhs/this->gcd(rhs);};
bitvector composition(const bitvector &rhs) const; bitvector composition(const bitvector &rhs) const;
unsigned int bitdiff(const bitvector &y) const; //number of differing bits (Hamming distance) 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 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) 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);}; 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); 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 //expand to separate bytes or ints
template <typename T> template <typename T>
void bitvector_expand(const bitvector &v, NRVec<T> &r) void bitvector_expand(const bitvector &v, NRVec<T> &r)

View File

@ -205,4 +205,17 @@ extern "C" int cblas_errprn(int ierr, int info, char *form, ...) {
} }
#endif #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 }//namespace

View File

@ -53,6 +53,8 @@ s << x.msg;
return s; return s;
} }
extern int endianity();
}//namespace }//namespace
#endif #endif

47
t.cc
View File

@ -2946,7 +2946,7 @@ cin >>n;
cout <<factorization(n)<<" phi = "<<eulerphi(n)<<endl; cout <<factorization(n)<<" phi = "<<eulerphi(n)<<endl;
} }
if(1) if(0)
{ {
bitvector ir; cin >>ir; bitvector ir; cin >>ir;
if(!ir.is_irreducible()) laerror("input must be an irreducible polynomial"); 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<<endl;
cout <<c%ir<<endl; cout <<c%ir<<endl;
cout <<cc<<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");
}
}
} }