one day i'll finish to write this article

README.md

@@ -0,0 +1,4 @@
+zut
+=====
+zut(Z under two) is a small and efficient C++ class for dealing with the
+factorisation of large Z/2Z polynomials.

zut.cpp

@@ -0,0 +1,443 @@
+// 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<unsigned int> 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 :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]<a;
+}
+
+bool operator>(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<size;i++,bit<<=1)
+        sqrt+=((bit&(*this))>>(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<zut> sff(vector<zut> 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;
+}
+
+/*
+ * Honestly I was about to write a whole paragraph on how this works but,
+ * honestly, wikipedia is a much better resource than me for that. I separated
+ * and commented briefly each subsection of the algorithm. It's not _that_
+ * complex plus, honestly, I have no idea who I'd write that for, as I'm sure
+ * you already knew about that before making this problem.
+ */
+vector<zut> bk(){
+    // first create the matrix Q-I
+    vector<zut> subalgebra;
+    vector<zut> factors;
+    int deg = (*this).degree();
+    zut x_2n(size*2);
+    zut fac(size*2); fac=(*this);
+    zut di(size); di=1;
+    vector<zut> mat; mat.reserve(deg);
+    for(int i=0; i<deg;i++){
+        x_2n=1; x_2n <<=i*2;
+        zut row(size*2); 
+        // Q = (x_2n%fac)+di), I=di
+        row=(((x_2n%fac)+di)<<(deg+1))+di; mat.push_back(row); 
+        di<<=1;
+    }
+
+    // then compute the row reduced echelon form of the matrix     
+    int row=0;
+    for (int j = deg*2; j>=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<factors.size(); i++){
+           if((factors[i]%test)==0){
+            if((factors[i]/test)!=1){
+               zut help((factors[i]/test));
+               factors[i]=test;
+               factors.push_back(help);
+               empty=false;
+               break;
+            }
+           }
+       }
+    }
+    if(empty)
+        factors[0].irreducible=true;
+    return factors;
+}
+
+// we get square free factors, then recurse with berlekamp to get irreducible
+// forms of all possible factors
+vector<zut> factorize() {
+    vector<zut> factors=(*this).sff();
+
+    bool done=false;
+    while(!done){
+
+        for(int i=0; i< factors.size(); i++) {
+            if(!factors[i].irreducible) {
+                vector<zut> 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;
+}
+
+}