#include <NTL/mat_GF2.h>
#include <NTL/vec_long.h>

#include <NTL/new.h>

NTL_START_IMPL


mat_GF2::mat_GF2(const mat_GF2& a)  
{  
   _mat_GF2__numcols = 0;  
   SetDims(a.NumRows(), a.NumCols());  
   _mat_GF2__rep = a._mat_GF2__rep;  
}  
  
mat_GF2& mat_GF2::operator=(const mat_GF2& a)  
{  
   SetDims(a.NumRows(), a.NumCols());  
   _mat_GF2__rep = a._mat_GF2__rep;  
   return *this;
}  
  
  
mat_GF2::mat_GF2(INIT_SIZE_TYPE, long n, long m)  
{  
   _mat_GF2__numcols = 0;  
   SetDims(n, m);  
}  
  
void mat_GF2::kill()  
{  
   _mat_GF2__numcols = 0;  
   _mat_GF2__rep.kill();  
}  
  
void mat_GF2::SetDims(long n, long m)  
{  
   if (n < 0 || m < 0)  
      Error("SetDims: bad args");  
  
   if (m != _mat_GF2__numcols) {  
      _mat_GF2__rep.kill();  
      _mat_GF2__numcols = m;  
   }  
        
   long oldmax = _mat_GF2__rep.MaxLength();  
   long i;  
   _mat_GF2__rep.SetLength(n);  
  
   for (i = oldmax; i < n; i++)  
      _mat_GF2__rep[i].FixLength(m);  
}  
     
        
void conv(mat_GF2& x, const vec_vec_GF2& a)  
{  
   long n = a.length();  
  
   if (n == 0) {  
      x.SetDims(0, 0);  
      return;  
   }  
  
   long m = a[0].length();  
   long i;  
  
   for (i = 1; i < n; i++)  
      if (a[i].length() != m)  
         Error("nonrectangular matrix");  
  
   x.SetDims(n, m);  
   for (i = 0; i < n; i++)  
      x[i] = a[i];  
}  
  
void swap(mat_GF2& X, mat_GF2& Y)  
{  
   swap(X._mat_GF2__numcols, Y._mat_GF2__numcols);  
   swap(X._mat_GF2__rep, Y._mat_GF2__rep);  
}  
  


long operator==(const mat_GF2& a, const mat_GF2& b)  
{  
   if (a.NumCols() != b.NumCols())  
      return 0;  
  
   if (a.NumRows() != b.NumRows())  
      return 0;  
  
   long n = a.NumRows();  
   long i;  
  
   for (i = 0; i < n; i++)  
      if (a[i] != b[i])  
         return 0;  
  
   return 1;  
}  
  
  
long operator!=(const mat_GF2& a, const mat_GF2& b)  
{  
   return !(a == b);  
}  

istream& operator>>(istream& s, mat_GF2& x)  
{  
   vec_vec_GF2 buf;  
   s >> buf;  
   conv(x, buf);  
   return s;  
}  
  
ostream& operator<<(ostream& s, const mat_GF2& a)  
{  
   long n = a.NumRows();  
   long i;  
   s << "[";  
   for (i = 0; i < n; i++) {
      s << a[i];  
      s << "\n";
   }
   s << "]";  
   return s;  
}  


void add(mat_GF2& X, const mat_GF2& A, const mat_GF2& B)  
{  
   long n = A.NumRows();  
   long m = A.NumCols();  
  
   if (B.NumRows() != n || B.NumCols() != m)   
      Error("matrix add: dimension mismatch");  
  
   X.SetDims(n, m);  

   long mw = (m + NTL_BITS_PER_LONG - 1)/NTL_BITS_PER_LONG;
  
   long i;  
   for (i = 0; i < n; i++) {
      _ntl_ulong *xp = X[i].rep.elts();
      const _ntl_ulong *ap = A[i].rep.elts();
      const _ntl_ulong *bp = B[i].rep.elts();
      long j;
      for (j = 0; j < mw; j++)
         xp[j] = ap[j] ^ bp[j];
   }
}  
  
static
void mul_aux(vec_GF2& x, const mat_GF2& A, const vec_GF2& b)  
{  
   long n = A.NumRows();  
   long l = A.NumCols();  
  
   if (l != b.length())  
      Error("matrix mul: dimension mismatch");  
  
   x.SetLength(n);  
  
   long i;  
  
   for (i = 0; i < n; i++) {  
      x.put(i, A[i] * b);
   }  
}  
  
  
void mul(vec_GF2& x, const mat_GF2& A, const vec_GF2& b)  
{  
   if (&b == &x || A.position1(x) != -1) {
      vec_GF2 tmp;
      mul_aux(tmp, A, b);
      x = tmp;
   }
   else
      mul_aux(x, A, b);
}  

static
void mul_aux(vec_GF2& x, const vec_GF2& a, const mat_GF2& B)  
{  
   long n = B.NumRows();  
   long l = B.NumCols();  
  
   if (n != a.length())  
      Error("matrix mul: dimension mismatch");  
  
   x.SetLength(l);  
   clear(x);

   const _ntl_ulong *ap = a.rep.elts();
   _ntl_ulong a_mask = 1;

   _ntl_ulong *xp = x.rep.elts();

   long lw = (l + NTL_BITS_PER_LONG - 1)/NTL_BITS_PER_LONG;
  
   long i;  
   for (i = 0; i < n; i++) {  
      if (*ap & a_mask) {
         const _ntl_ulong *bp = B[i].rep.elts();
         long j;
         for (j = 0; j < lw; j++)
            xp[j] ^= bp[j];
      }

      a_mask <<= 1;
      if (!a_mask) {
         a_mask = 1;
         ap++;
      }
   }  
}  

void mul(vec_GF2& x, const vec_GF2& a, const mat_GF2& B)
{
   if (&a == &x || B.position1(x) != -1) {
      vec_GF2 tmp;
      mul_aux(tmp, a, B);
      x = tmp;
   }
   else
      mul_aux(x, a, B);
}
  
void mul_aux(mat_GF2& X, const mat_GF2& A, const mat_GF2& B)  
{  
   long n = A.NumRows();  
   long l = A.NumCols();  
   long m = B.NumCols();  
  
   if (l != B.NumRows())  
      Error("matrix mul: dimension mismatch");  
  
   X.SetDims(n, m);  
  
   long i;  
  
   for (i = 1; i <= n; i++) {  
      mul_aux(X(i), A(i), B);
   }  
}  
  
  
void mul(mat_GF2& X, const mat_GF2& A, const mat_GF2& B)  
{  
   if (&X == &A || &X == &B) {  
      mat_GF2 tmp;  
      mul_aux(tmp, A, B);  
      X = tmp;  
   }  
   else  
      mul_aux(X, A, B);  
}  
  

     
  
void ident(mat_GF2& X, long n)  
{  
   X.SetDims(n, n);  
   clear(X);
   long i;  
  
   for (i = 0; i < n; i++)  
      X.put(i, i, to_GF2(1));
} 


void determinant(GF2& d, const mat_GF2& M_in)
{
   long k, n;
   long i, j;
   long pos;

   n = M_in.NumRows();

   if (M_in.NumCols() != n)
      Error("determinant: nonsquare matrix");

   if (n == 0) {
      set(d);
      return;
   }

   mat_GF2 M;

   M = M_in;

   long wn = (n + NTL_BITS_PER_LONG - 1)/NTL_BITS_PER_LONG;

   for (k = 0; k < n; k++) {
      long wk = k/NTL_BITS_PER_LONG;
      long bk = k - wk*NTL_BITS_PER_LONG;
      _ntl_ulong k_mask = 1UL << bk;

      pos = -1;
      for (i = k; i < n; i++) {
         if (M[i].rep.elts()[wk] & k_mask) {
            pos = i;
            break;
         }
      }

      if (pos != -1) {
         if (k != pos) {
            swap(M[pos], M[k]);
         }


         _ntl_ulong *y = M[k].rep.elts();

         for (i = k+1; i < n; i++) {
            // M[i] = M[i] + M[k]*M[i,k]

            if (M[i].rep.elts()[wk] & k_mask) {
               _ntl_ulong *x = M[i].rep.elts();

               for (j = wk; j < wn; j++)
                  x[j] ^= y[j];
            }

         }
      }
      else {
         clear(d);
         return;
      }
   }

   set(d);
   return;
}

static
long IsUnitVector(const vec_GF2& a, long i)
{
   long wi = i/NTL_BITS_PER_LONG;
   long bi = i - wi*NTL_BITS_PER_LONG;

   const _ntl_ulong *p = a.rep.elts();
   long wdlen = a.rep.length();

   long j;

   for (j = 0; j < wi; j++)
      if (p[j] != 0) return 0;

   if (p[wi] != (1UL << bi))
      return 0;

   for (j = wi+1; j < wdlen; j++)
      if (p[j] != 0) return 0;

   return 1;
}


long IsIdent(const mat_GF2& A, long n)
{
   if (A.NumRows() != n || A.NumCols() != n)
      return 0;

   if (n == 0) return 1;

   long i;

   for (i = 0; i < n; i++)
      if (!IsUnitVector(A[i], i))
         return 0;

   return 1;
}

void AddToCol(mat_GF2& x, long j, const vec_GF2& a)
// add a to column j of x
// ALIAS RESTRICTION: a should not alias any row of x
{
   long n = x.NumRows();
   long m = x.NumCols();

   if (a.length() != n || j < 0 || j >= m)
      Error("AddToCol: bad args");

   long wj = j/NTL_BITS_PER_LONG;
   long bj = j - wj*NTL_BITS_PER_LONG;
   _ntl_ulong j_mask = 1UL << bj;

   const _ntl_ulong *ap = a.rep.elts();
   _ntl_ulong a_mask = 1;

   long i;
   for (i = 0; i < n; i++) {
      if (*ap & a_mask) 
         x[i].rep.elts()[wj] ^= j_mask;

      a_mask <<= 1;
      if (!a_mask) {
         a_mask = 1;
         ap++;
      }
   }
}


void transpose_aux(mat_GF2& X, const mat_GF2& A)
{
   long n = A.NumRows();
   long m = A.NumCols();

   X.SetDims(m, n);
   clear(X);

   long i;
   for (i = 0; i < n; i++)
      AddToCol(X, i, A[i]);
}
            

void transpose(mat_GF2& X, const mat_GF2& A)
{
   if (&X == &A) {
      mat_GF2 tmp;
      transpose_aux(tmp, A);