// Z under two // Because what's some maths code without a bad french pun :) // >>=, <<=, <<, >>, +=, +, *=, *, /, %, &=, &, ==, !=, <, >, = are implemented class zut { unsigned int* number; int size; bool irreducible; public: // Initializers ------------------------------------------------- zut(unsigned int* a, int b){ number=a; size=b; irreducible=false; } zut(initializer_list a){ number=new unsigned int[a.size()]; size=a.size()*32; copy(a.begin(), a.end(), number); irreducible=false; } zut(int b) { size=b; number = new unsigned int[b/32]; irreducible=false; } zut(const zut &b){ size=b.size; number = new unsigned int[b.size/32]; (*this)=b; irreducible=b.irreducible; } void resize(int size){ this->size=size; } // Primitives --------------------------------------------------- // https://stackoverflow.com/questions/2773890/efficient-bitshifting-an-array-of-int void operator<<=(int s){ for (int i = s; i > 0; i -= (i % 32)) { if (!(i % 32)) { for (int y = (size/32)-1; y > 0; y--) number[y] = number[y - 1]; number[0] = 0; i -= 32; } else { for (int y = (size/32)-1; y > 0; y--) number[y] = (number[y] << (i%32)) | ((number[y - 1] >> (32-i%32)) ); number[0] <<= i%32; } } } zut operator<<(int s){ zut res((*this)); res <<= s; return res; } void operator>>=(int s){ for (int i = s; i > 0; i -= (i % 32)) { if (!(i % 32)) { for (int y = (size/32)-1; y >0; y--) number[y] = number[y - 1]; number[0] = 0; i -= 32; } else { for (int y = (size / 32) - 1; y > 0; y--) number[y] = (number[y] >> (i%32)) | ((number[y - 1] << (32 - (i%32)) & 1)); number[0] >>= i%32; } } } zut operator>>(int s){ zut res((*this)); res >>= s; return res; } // i don't think this can be made faster without arch specific code :p // in a z/2z field 1+1=0 and 1+0 = 1, thus it's a simple xor // also i guess you know if you are reading this code but substraction = // addition in Z/2Z so this can be(and is, ie in division) used interchangeably void operator+=(const zut& b){ for (int i = 0; i < size / 32; i++){ if(!(i < b.size / 32)) break; number[i] ^= b.number[i]; } } zut operator+(const zut& b){ zut res((*this)); res+=b; return res; } // https://en.wikipedia.org/wiki/Finite_field_arithmetic#C_programming_example // Ugh, bit by bit comparison, it really looks like a bottleneck to me, but // I can't find a way to make it faster... void operator*=(zut& b){ zut tmp(size); tmp=(*this); (*this)=0; for (int i = 0; i < size; i++) { if(!(i < b.size)) break; // this verifies if the n-th bit inside b is lit if (b.get_bit(i)) (*this)+=tmp; tmp <<=1; } } zut operator*(zut& b){ zut res((*this)); res *= b; return res; } // /!\ primitive might fail if 2 values have different size!! zut operator/(zut& b){ zut r((*this)); zut q(size); q=0; int deg_a=(*this).degree(); int deg_div=deg_a-b.degree(); for (int i = 0; i <= deg_div; i++) { q<<=1; int y=deg_a-i; if (r.get_bit(y)) { q++; r+=(b << (deg_div-i)); } } return q; } zut operator%(zut& b){ zut r((*this)); zut q(size); q=0; int deg_a=(*this).degree(); int deg_div=deg_a-b.degree(); for (int i = 0; i <= deg_div; i++) { q<<=1; int y=deg_a-i; if (r.get_bit(y)) { q++; r+=(b << (deg_div-i)); } } return r; } void operator&=(int b){ for (int i = 0; i < size / 32; i++) number[i] &= b; } zut operator&(int b){ zut res((*this)); res &= b; return res; } void operator&=(const zut& b){ for (int i = 0; i < size / 32; i++){ if(!(i < b.size / 32)) break; number[i] &= b.number[i]; } } zut operator&(const zut& b){ zut res((*this)); res &= b; return res; } void operator++(int a){ number[0]^=1; } // Logical operators -------------------------------------------- bool operator==(const zut& b){ for (int i = 0; i < size / 32; i++){ if(!(i < b.size / 32)) break; if(number[i]!=b.number[i]) return false; } return true; } bool operator==(int b){ if(number[0]!=b) return false; for (int i = 1; i < size / 32; i++) if(number[i]!=0) return false; return true; } bool operator!=(const zut& b){ return !((*this)==b); } bool operator!=(int b){ return !((*this)==b); } bool operator<(int a){ return number[0](int a){ return number[0]>a; } void operator=(const zut& b){ for (int i = 0; i < b.size / 32; i++) number[i] = b.number[i]; for (int i = b.size / 32; i < size / 32; i++) number[i] = 0; } void operator=(int b){ number[0] = b; for (int i = 1; i < size / 32; i++) number[i] = 0; } // Algorithms --------------------------------------------------- // should be able to do this faster... zut get_bits(int start, int end) { zut res(size); res=0; for (int i = start; i < end; i++) { if ((*this).get_bit(i)) res.number[(i-start)/32] |= (1 << (i-start)%32); } return res; } bool get_bit(int bit){ if (((number[(bit / 32)] & (1 << bit % 32)) >> bit % 32) & 1) return true; return false; } /* So, how does this work? * Well to begin with we are dealing with a Z/2Z field and as such have an * equation in the form of sigma(ax^n), for n up to the maximum bit size of * our current number. Knowing a is either 0 or 1 in this field, it can also * be disregarded as d0=0. This means the derivative is in * the form of sigma(nx^(n-1)). This also means that any odd power is * completely disregarded as they would be, by definition, 0. * Thus we only need to lower the power of every power by one and to * drop x^n where n is odd. */ zut derivative() { return ((*this)&0xaaaaaaaa)>>1; // -> 0b101010... -> 1+3 aka mask odd numbers } // https://stackoverflow.com/questions/671815/what-is-the-fastest-most-efficient-way-to-find-the-highest-set-bit-msb-in-an-i int degree() { for (int i = (size / 32) - 1; i >= 0; i--) if (number[i]) return (i * 32) + (31 - __builtin_clz(number[i])); return 0; } // mostly same as derivative: sqrt(x^n) = x^(n/2) (except 1)thus we just shift everything zut sqrt(){ zut sqrt(size); zut bit(size); bit=1; for (int i=0;i>(i>>1)); return sqrt; } // aren't you supposed to learn what this is in middle school? ;) zut gcd(zut b) { zut a((*this)); zut mod=a%b; while(mod>2){ a=b; b=mod; mod=a%b; } if(mod==0) return b; return mod; } // https://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields#Square-free_factorization // ^same as berlekamp, explains it better than i ever could vector sff(vector factors={}) { zut d=(*this).derivative(); // step 2 of sff(wikipedia) if(d==0){ factors=(*this).sqrt().sff(); factors.insert(factors.end(), factors.begin(), factors.end()); } else{ zut gcd=(*this).gcd(d); if(gcd==1){ factors.push_back((*this)); return factors; } factors=gcd.sff(factors); factors=((*this)/gcd).sff(factors); } return factors; } /* * Berlekamp algorithm in all its glory, along with a baked in matrix reduced * echelon form transformation. */ vector bk(){ // first create the matrix Q-I vector subalgebra; vector factors; int deg = (*this).degree(); zut x_2n(size*2); zut fac(size*2); fac=(*this); zut di(size); di=1; vector mat; mat.reserve(deg); for(int i=0; i=0; j--) { int pivrow = -1; for (int i = row; i < mat.size(); i++) { if (mat[i].get_bit(j)) { pivrow = i; break; } } if (pivrow != -1 ) { swap(mat[row], mat[pivrow]); for (int i = 0; i < mat.size(); i++) { if (i != row) { if (mat[i].get_bit(j)) { mat[i]+= mat[row]; } } } row++; } } // then read off berlekamp subalgebra from the null space basis bool empty=true; for(zut &row: mat) { zut row_left = row.get_bits(deg+2, (deg*2)+2); if(row_left==0){ zut row_right = row.get_bits(0, deg); if(row_right!=1){ // bk subalgebra takes half of the matrix thus size/2 row_right.resize(size); subalgebra.push_back(row_right); empty=false; } } } // if we had no correct null basis it is irreductible if(empty){ (*this).irreducible=true; factors.push_back((*this)); return factors; } // otherwise calculates irreducible factor from berlekamp subalgebra factors.push_back((*this)); for(zut &vector: subalgebra) { zut test((*this).gcd(vector)); test=(*this)/test; for(int i=0; i factorize() { vector factors=(*this).sff(); bool done=false; while(!done){ for(int i=0; i< factors.size(); i++) { if(!factors[i].irreducible) { vector factored = factors[i].bk(); factors.erase(factors.begin()+i); factors.insert(factors.end(), factored.begin(), factored.end()); break; } if(i==factors.size()-1){ done=true; } } if(factors.size()==0) done=true; } return factors; } }