From 9bceebdd2937b86875e7c9118c783c32127f54fa Mon Sep 17 00:00:00 2001 From: Jiri Pittner Date: Mon, 1 Jan 2024 10:58:30 +0100 Subject: [PATCH] polynomial irreducibility test in GF2 --- bitvector.cc | 64 ++++++++++++++++++++++++++++++++++++++++++++++++--- bitvector.h | 6 +++++ polynomial.cc | 6 +++++ polynomial.h | 1 + t.cc | 7 ++++++ 5 files changed, 81 insertions(+), 3 deletions(-) diff --git a/bitvector.cc b/bitvector.cc index ca91541..c6d3dd7 100644 --- a/bitvector.cc +++ b/bitvector.cc @@ -18,6 +18,7 @@ #include "bitvector.h" #include +#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 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; j0; --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) diff --git a/bitvector.h b/bitvector.h index c23d212..8ea03fe 100644 --- a/bitvector.h +++ b/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 diff --git a/polynomial.cc b/polynomial.cc index 8ff084a..220a424 100644 --- a/polynomial.cc +++ b/polynomial.cc @@ -249,6 +249,12 @@ resize(n,false); for(int i=0; i<=n; ++i) (*this)[i] = (T) binom(n,i); } +template +Polynomial Polynomial::composition(const Polynomial &rhs) const +{ +return value(*this,rhs); +} + /***************************************************************************//** * forced instantization in the corresponding object file diff --git a/polynomial.h b/polynomial.h index 81bdff2..49ff8d5 100644 --- a/polynomial.h +++ b/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; diff --git a/t.cc b/t.cc index 7581412..7e54b08 100644 --- a/t.cc +++ b/t.cc @@ -2949,12 +2949,19 @@ cout <>a; bitvector ai = a.field_inv(ir); cout<< "inverse = "<