X = tmp;
   }
   else
      transpose_aux(X, A);
}

   

void solve(GF2& d, vec_GF2& X, const mat_GF2& A, const vec_GF2& b)

{
   long n = A.NumRows();
   if (A.NumCols() != n)
      Error("solve: nonsquare matrix");

   if (b.length() != n)
      Error("solve: dimension mismatch");

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

   long i, j, k, pos;

   mat_GF2 M;
   M.SetDims(n, n+1);

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

   AddToCol(M, n, b);

   long wn = ((n+1) + 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;
      }
   }

   vec_GF2 XX;
   XX.SetLength(n+1);
   XX.put(n, 1);

   for (i = n-1; i >= 0; i--) {
      XX.put(i, XX*M[i]);
   }

   XX.SetLength(n);
   X = XX;

   set(d);
   return;
}



void inv(GF2& d, mat_GF2& X, const mat_GF2& A)
{
   long n = A.NumRows();
   if (A.NumCols() != n)
      Error("solve: nonsquare matrix");

   if (n == 0) {
      X.SetDims(0, 0);
      set(d);
   }

   long i, j, k, pos;

   mat_GF2 M;
   M.SetDims(n, 2*n);

   vec_GF2 aa;
   aa.SetLength(2*n);


   for (i = 0; i < n; i++) {
      aa = A[i];
      aa.SetLength(2*n);
      aa.put(n+i, 1);
      M[i] = aa;
   }

   long wn = ((2*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;
      }
   }

   vec_GF2 XX;
   XX.SetLength(2*n);

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

   for (j = 0; j < n; j++) {
      XX.SetLength(n+j+1);
      clear(XX);
      XX.put(n+j, to_GF2(1));
      
      for (i = n-1; i >= 0; i--) {
         XX.put(i, XX*M[i]);
      }
   
      XX.SetLength(n);
      AddToCol(X, j, XX);
   }

   set(d);
   return;
}





long gauss(mat_GF2& M, long w)
{
   long k, l;
   long i, j;
   long pos;

   long n = M.NumRows();
   long m = M.NumCols();

   if (w < 0 || w > m)
      Error("gauss: bad args");

   long wm = (m + NTL_BITS_PER_LONG - 1)/NTL_BITS_PER_LONG;

   l = 0;
   for (k = 0; k < w && l < 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 = l; i < n; i++) {
         if (M[i].rep.elts()[wk] & k_mask) {
            pos = i;
            break;
         }
      }

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

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

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

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

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

         l++;
      }
   }
   
   return l;
}

long gauss(mat_GF2& M)
{
   return gauss(M, M.NumCols());
}


void image(mat_GF2& X, const mat_GF2& A)
{
   mat_GF2 M;
   M = A;
   long r = gauss(M);
   M.SetDims(r, M.NumCols());
   X = M;
}

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

   mat_GF2 M;
   long r;

   transpose(M, A);
   r = gauss(M);

   X.SetDims(m-r, m);
   clear(X);

   long i, j, k;

   vec_long D;
   D.SetLength(m);
   for (j = 0; j < m; j++) D[j] = -1;

   j = -1;
   for (i = 0; i < r; i++) {
      do {
         j++;
      } while (M.get(i, j) == 0); 

      D[j] = i;
   }

   for (k = 0; k < m-r; k++) {
      vec_GF2& v = X[k];
      long pos = 0;
      for (j = m-1; j >= 0; j--) {
         if (D[j] == -1) {
            if (pos == k) {
               v[j] = 1;
               // v.put(j, to_GF2(1));
            }
            pos++;
         }
         else {
            v[j] = v*M[D[j]];
            // v.put(j, v*M[D[j]]);
         }
      }
   }
}

   
void mul(mat_GF2& X, const mat_GF2& A, GF2 b)
{
   X = A;
   if (b == 0)
      clear(X);
}

void diag(mat_GF2& X, long n, GF2 d)  
{  
   if (d == 1)
      ident(X, n);
   else {
      X.SetDims(n, n);
      clear(X);
   }
} 

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

   if (d == 1)
      return IsIdent(A, n);
   else
      return IsZero(A);
}


long IsZero(const mat_GF2& a)
{
   long n = a.NumRows();
   long i;

   for (i = 0; i < n; i++)
      if (!IsZero(a[i]))
         return 0;

   return 1;
}

void clear(mat_GF2& x)
{
   long n = x.NumRows();
   long i;
   for (i = 0; i < n; i++)
      clear(x[i]);
}


mat_GF2 operator+(const mat_GF2& a, const mat_GF2& b)
{
   mat_GF2 res;
   add(res, a, b);
   NTL_OPT_RETURN(mat_GF2, res);
}

mat_GF2 operator*(const mat_GF2& a, const mat_GF2& b)
{
   mat_GF2 res;
   mul_aux(res, a, b);
   NTL_OPT_RETURN(mat_GF2, res);
}

mat_GF2 operator-(const mat_GF2& a, const mat_GF2& b)
{
   mat_GF2 res;
   add(res, a, b);
   NTL_OPT_RETURN(mat_GF2, res);
}


vec_GF2 operator*(const mat_GF2& a, const vec_GF2& b)
{
   vec_GF2 res;
   mul_aux(res, a, b);
   NTL_OPT_RETURN(vec_GF2, res);
}

vec_GF2 operator*(const vec_GF2& a, const mat_GF2& b)
{
   vec_GF2 res;
   mul_aux(res, a, b);
   NTL_OPT_RETURN(vec_GF2, res);
}


void inv(mat_GF2& X, const mat_GF2& A)
{
   GF2 d;
   inv(d, X, A);
   if (d == 0) Error("inv: non-invertible matrix");
}

void power(mat_GF2& X, const mat_GF2& A, const ZZ& e)
{
   if (A.NumRows() != A.NumCols()) Error("power: non-square matrix");

   if (e == 0) {
      ident(X, A.NumRows());
      return;
   }

   mat_GF2 T1, T2;
   long i, k;

   k = NumBits(e);
   T1 = A;

   for (i = k-2; i >= 0; i--) {
      sqr(T2, T1);
      if (bit(e, i))
         mul(T1, T2, A);
      else
         T1 = T2;
   }

   if (e < 0)
      inv(X, T1);
   else
      X = T1;
}

NTL_END_IMPL