gf2mat.cpp

Go to the documentation of this file.

#include <itpp/base/gf2mat.h>
#include <itpp/base/specmat.h>
#include <itpp/base/matfunc.h>
#include <itpp/base/converters.h>
#include <iostream>

namespace itpp
{

// ====== IMPLEMENTATION OF THE ALIST CLASS ==========

GF2mat_sparse_alist::GF2mat_sparse_alist(const std::string &fname)
    : data_ok(false)
{
  read(fname);
}

void GF2mat_sparse_alist::read(const std::string &fname)
{
  std::ifstream file;
  std::string line;
  std::stringstream ss;

  file.open(fname.c_str());
  it_assert(file.is_open(),
            "GF2mat_sparse_alist::read(): Could not open file \""
            << fname << "\" for reading");

  // parse N and M:
  getline(file, line);
  ss << line;
  ss >> N >> M;
  it_assert(!ss.fail(),
            "GF2mat_sparse_alist::read(): Wrong alist data (N or M)");
  it_assert((N > 0) && (M > 0),
            "GF2mat_sparse_alist::read(): Wrong alist data");
  ss.seekg(0, std::ios::end);
  ss.clear();

  // parse max_num_n and max_num_m
  getline(file, line);
  ss << line;
  ss >> max_num_n >> max_num_m;
  it_assert(!ss.fail(),
            "GF2mat_sparse_alist::read(): Wrong alist data (max_num_{n,m})");
  it_assert((max_num_n >= 0) && (max_num_n <= N) &&
            (max_num_m >= 0) && (max_num_m <= M),
            "GF2mat_sparse_alist::read(): Wrong alist data");
  ss.seekg(0, std::ios::end);
  ss.clear();

  // parse weight of each column n
  num_nlist.set_size(N);
  num_nlist.clear();
  getline(file, line);
  ss << line;
  for (int i = 0; i < N; i++) {
    ss >> num_nlist(i);
    it_assert(!ss.fail(),
              "GF2mat_sparse_alist::read(): Wrong alist data (num_nlist("
              << i << "))");
    it_assert((num_nlist(i) >= 0) && (num_nlist(i) <= M),
              "GF2mat_sparse_alist::read(): Wrong alist data (num_nlist("
              << i << "))");
  }
  ss.seekg(0, std::ios::end);
  ss.clear();

  // parse weight of each row m
  num_mlist.set_size(M);
  num_mlist.clear();
  getline(file, line);
  ss << line;
  for (int i = 0; i < M; i++) {
    ss >> num_mlist(i);
    it_assert(!ss.fail(),
              "GF2mat_sparse_alist::read(): Wrong alist data (num_mlist("
              << i << "))");
    it_assert((num_mlist(i) >= 0) && (num_mlist(i) <= N),
              "GF2mat_sparse_alist::read(): Wrong alist data (num_mlist("
              << i << "))");
  }
  ss.seekg(0, std::ios::end);
  ss.clear();

  // parse coordinates in the n direction with non-zero entries
  nlist.set_size(N, max_num_n);
  nlist.clear();
  for (int i = 0; i < N; i++) {
    getline(file, line);
    ss << line;
    for (int j = 0; j < num_nlist(i); j++) {
      ss >> nlist(i, j);
      it_assert(!ss.fail(),
                "GF2mat_sparse_alist::read(): Wrong alist data (nlist("
                << i << "," << j << "))");
      it_assert((nlist(i, j) > 0) && (nlist(i, j) <= M),
                "GF2mat_sparse_alist::read(): Wrong alist data (nlist("
                << i << "," << j << "))");
    }
    ss.seekg(0, std::ios::end);
    ss.clear();
  }

  // parse coordinates in the m direction with non-zero entries
  mlist.set_size(M, max_num_m);
  mlist.clear();
  for (int i = 0; i < M; i++) {
    getline(file, line);
    ss << line;
    for (int j = 0; j < num_mlist(i); j++) {
      ss >> mlist(i, j);
      it_assert(!ss.fail(),
                "GF2mat_sparse_alist::read(): Wrong alist data (mlist("
                << i << "," << j << "))");
      it_assert((mlist(i, j) > 0) && (mlist(i, j) <= N),
                "GF2mat_sparse_alist::read(): Wrong alist data (mlist("
                << i << "," << j << "))");
    }
    ss.seekg(0, std::ios::end);
    ss.clear();
  }

  file.close();
  data_ok = true;
}

void GF2mat_sparse_alist::write(const std::string &fname) const
{
  it_assert(data_ok,
            "GF2mat_sparse_alist::write(): alist data not ready for writing");

  std::ofstream file(fname.c_str(), std::ofstream::out);
  it_assert(file.is_open(),
            "GF2mat_sparse_alist::write(): Could not open file \""
            << fname << "\" for writing");

  file << N << " " << M << std::endl;
  file << max_num_n << " " << max_num_m << std::endl;

  for (int i = 0; i < num_nlist.length() - 1; i++)
    file << num_nlist(i) << " ";
  file << num_nlist(num_nlist.length() - 1) << std::endl;

  for (int i = 0; i < num_mlist.length() - 1; i++)
    file << num_mlist(i) << " ";
  file << num_mlist(num_mlist.length() - 1) << std::endl;

  for (int i = 0; i < N; i++) {
    for (int j = 0; j < num_nlist(i) - 1; j++)
      file << nlist(i, j) << " ";
    file << nlist(i, num_nlist(i) - 1) << std::endl;
  }

  for (int i = 0; i < M; i++) {
    for (int j = 0; j < num_mlist(i) - 1; j++)
      file << mlist(i, j) << " ";
    file << mlist(i, num_mlist(i) - 1) << std::endl;
  }

  file.close();
}


GF2mat_sparse GF2mat_sparse_alist::to_sparse(bool transpose) const
{
  GF2mat_sparse sbmat(M, N, max_num_m);

  for (int i = 0; i < M; i++) {
    for (int j = 0; j < num_mlist(i); j++) {
      sbmat.set_new(i, mlist(i, j) - 1, bin(1));
    }
  }
  sbmat.compact();

  if (transpose) {
    return sbmat.transpose();
  }
  else {
    return sbmat;
  }
}


// ----------------------------------------------------------------------
// WARNING: This method is very slow. Sparse_Mat has to be extended with
// some extra functions to improve the performance of this.
// ----------------------------------------------------------------------
void GF2mat_sparse_alist::from_sparse(const GF2mat_sparse &sbmat,
                                      bool transpose)
{
  if (transpose) {
    from_sparse(sbmat.transpose(), false);
  }
  else {
    // check matrix dimension
    M = sbmat.rows();
    N = sbmat.cols();

    num_mlist.set_size(M);
    num_nlist.set_size(N);

    // fill mlist matrix, num_mlist vector and max_num_m
    mlist.set_size(0, 0);
    int tmp_cols = 0; // initial number of allocated columns
    for (int i = 0; i < M; i++) {
      ivec temp_row(0);
      for (int j = 0; j < N; j++) {
        if (sbmat(i, j) == bin(1)) {
          temp_row = concat(temp_row, j + 1);
        }
      }
      int trs = temp_row.size();
      if (trs > tmp_cols) {
        tmp_cols = trs;
        mlist.set_size(M, tmp_cols, true);
      }
      else if (trs < tmp_cols) {
        temp_row.set_size(tmp_cols, true);
      }
      mlist.set_row(i, temp_row);
      num_mlist(i) = trs;
    }
    max_num_m = max(num_mlist);

    // fill nlist matrix, num_nlist vector and max_num_n
    nlist.set_size(0, 0);
    tmp_cols = 0; // initial number of allocated columns
    for (int j = 0; j < N; j++) {
      ivec temp_row = sbmat.get_col(j).get_nz_indices() + 1;
      int trs = temp_row.size();
      if (trs > tmp_cols) {
        tmp_cols = trs;
        nlist.set_size(N, tmp_cols, true);
      }
      else if (trs < tmp_cols) {
        temp_row.set_size(tmp_cols, true);
      }
      nlist.set_row(j, temp_row);
      num_nlist(j) = trs;
    }
    max_num_n = max(num_nlist);

    data_ok = true;
  }
}


// ----------------------------------------------------------------------
// Implementation of a dense GF2 matrix class
// ----------------------------------------------------------------------

GF2mat::GF2mat(int i, int j): nrows(i), ncols(j),
    nwords((j >> shift_divisor) + 1)
{
  data.set_size(nrows, nwords);
  data.clear();
}

GF2mat::GF2mat(): nrows(1), ncols(1), nwords(1)
{
  data.set_size(nrows, nwords);
  data.clear();
}

GF2mat::GF2mat(const bvec &x, bool is_column)
{
  if (is_column) {     // create column vector
    nrows = length(x);
    ncols = 1;
    nwords = 1;
    data.set_size(nrows, nwords);
    data.clear();
    for (int i = 0; i < nrows; i++) {
      set(i, 0, x(i));
    }
  }
  else {    // create row vector
    nrows = 1;
    ncols = length(x);
    nwords = (ncols >> shift_divisor) + 1;
    data.set_size(nrows, nwords);
    data.clear();
    for (int i = 0; i < ncols; i++) {
      set(0, i, x(i));
    }
  }
}


GF2mat::GF2mat(const bmat &X): nrows(X.rows()), ncols(X.cols())
{
  nwords = (ncols >> shift_divisor) + 1;
  data.set_size(nrows, nwords);
  data.clear();
  for (int i = 0; i < nrows; i++) {
    for (int j = 0; j < ncols; j++) {
      set(i, j, X(i, j));
    }
  }
}


GF2mat::GF2mat(const GF2mat_sparse &X)
{
  nrows = X.rows();
  ncols = X.cols();
  nwords = (ncols >> shift_divisor) + 1;
  data.set_size(nrows, nwords);
  for (int i = 0; i < nrows; i++) {
    for (int j = 0; j < nwords; j++) {
      data(i, j) = 0;
    }
  }

  for (int j = 0; j < ncols; j++) {
    for (int i = 0; i < X.get_col(j).nnz(); i++)   {
      bin b = X.get_col(j).get_nz_data(i);
      set(X.get_col(j).get_nz_index(i), j, b);
    }
  }
}

GF2mat::GF2mat(const GF2mat_sparse &X, int m1, int n1, int m2, int n2)
{
  it_assert(X.rows() > m2, "GF2mat(): indexes out of range");
  it_assert(X.cols() > n2, "GF2mat(): indexes out of range");
  it_assert(m1 >= 0 && n1 >= 0 && m2 >= m1 && n2 >= n1,
            "GF2mat::GF2mat(): indexes out of range");

  nrows = m2 - m1 + 1;
  ncols = n2 - n1 + 1;
  nwords = (ncols >> shift_divisor) + 1;
  data.set_size(nrows, nwords);

  for (int i = 0; i < nrows; i++) {
    for (int j = 0; j < nwords; j++) {
      data(i, j) = 0;
    }
  }

  for (int i = 0; i < nrows; i++) {
    for (int j = 0; j < ncols; j++)   {
      bin b = X(i + m1, j + n1);
      set(i, j, b);
    }
  }
}


GF2mat::GF2mat(const GF2mat_sparse &X, const ivec &columns)
{
  it_assert(X.cols() > max(columns),
            "GF2mat::GF2mat(): index out of range");
  it_assert(min(columns) >= 0,
            "GF2mat::GF2mat(): column index must be positive");

  nrows = X.rows();
  ncols = length(columns);
  nwords = (ncols >> shift_divisor) + 1;
  data.set_size(nrows, nwords);

  for (int i = 0; i < nrows; i++) {
    for (int j = 0; j < nwords; j++) {
      data(i, j) = 0;
    }
  }

  for (int j = 0; j < ncols; j++) {
    for (int i = 0; i < X.get_col(columns(j)).nnz(); i++)   {
      bin b = X.get_col(columns(j)).get_nz_data(i);
      set(X.get_col(columns(j)).get_nz_index(i), j, b);
    }
  }
}


void GF2mat::set_size(int m, int n, bool copy)
{
  nrows = m;
  ncols = n;
  nwords = (ncols >> shift_divisor) + 1;
  data.set_size(nrows, nwords, copy);
  if (!copy)
    data.clear();
}


GF2mat_sparse GF2mat::sparsify() const
{
  GF2mat_sparse Z(nrows, ncols);
  for (int i = 0; i < nrows; i++) {
    for (int j = 0; j < ncols; j++) {
      if (get(i, j) == 1) {
        Z.set(i, j, 1);
      }
    }
  }

  return Z;
}

bvec GF2mat::bvecify() const
{
  it_assert(nrows == 1 || ncols == 1,
            "GF2mat::bvecify() matrix must be a vector");
  int n = (nrows == 1 ? ncols : nrows);
  bvec result(n);
  if (nrows == 1) {
    for (int i = 0; i < n; i++) {
      result(i) = get(0, i);
    }
  }
  else {
    for (int i = 0; i < n; i++) {
      result(i) = get(i, 0);
    }
  }
  return result;
}


void GF2mat::set_row(int i, bvec x)
{
  it_assert(length(x) == ncols,
            "GF2mat::set_row(): dimension mismatch");
  for (int j = 0; j < ncols; j++) {
    set(i, j, x(j));
  }
}

void GF2mat::set_col(int j, bvec x)
{
  it_assert(length(x) == nrows,
            "GF2mat::set_col(): dimension mismatch");
  for (int i = 0; i < nrows; i++) {
    set(i, j, x(i));
  }
}


GF2mat gf2dense_eye(int m)
{
  GF2mat Z(m, m);
  for (int i = 0; i < m; i++) {
    Z.set(i, i, 1);
  }
  return Z;
}

GF2mat GF2mat::get_submatrix(int m1, int n1, int m2, int n2) const
{
  it_assert(m1 >= 0 && n1 >= 0 && m2 >= m1 && n2 >= n1
            && m2 < nrows && n2 < ncols,
            "GF2mat::get_submatrix() index out of range");
  GF2mat result(m2 - m1 + 1, n2 - n1 + 1);

  for (int i = m1; i <= m2; i++) {
    for (int j = n1; j <= n2; j++) {
      result.set(i - m1, j - n1, get(i, j));
    }
  }

  return result;
}


GF2mat GF2mat::concatenate_vertical(const GF2mat &X) const
{
  it_assert(X.ncols == ncols,
            "GF2mat::concatenate_vertical(): dimension mismatch");
  it_assert(X.nwords == nwords,
            "GF2mat::concatenate_vertical(): dimension mismatch");

  GF2mat result(nrows + X.nrows, ncols);
  for (int i = 0; i < nrows; i++) {
    for (int j = 0; j < nwords; j++) {
      result.data(i, j) = data(i, j);
    }
  }

  for (int i = 0; i < X.nrows; i++) {
    for (int j = 0; j < nwords; j++) {
      result.data(i + nrows, j) = X.data(i, j);
    }
  }

  return result;
}

GF2mat GF2mat::concatenate_horizontal(const GF2mat &X) const
{
  it_assert(X.nrows == nrows,
            "GF2mat::concatenate_horizontal(): dimension mismatch");

  GF2mat result(nrows, X.ncols + ncols);
  for (int i = 0; i < nrows; i++) {
    for (int j = 0; j < ncols; j++) {
      result.set(i, j, get(i, j));
    }
  }

  for (int i = 0; i < nrows; i++) {
    for (int j = 0; j < X.ncols; j++) {
      result.set(i, j + ncols, X.get(i, j));
    }
  }

  return result;
}

bvec GF2mat::get_row(int i) const
{
  bvec result(ncols);
  for (int j = 0; j < ncols; j++) {
    result(j) = get(i, j);
  }

  return result;
}

bvec GF2mat::get_col(int j) const
{
  bvec result(nrows);
  for (int i = 0; i < nrows; i++) {
    result(i) = get(i, j);
  }

  return result;
}


int GF2mat::T_fact(GF2mat &T, GF2mat &U, ivec &perm) const
{
  T = gf2dense_eye(nrows);
  U = *this;

  perm = zeros_i(ncols);
  for (int i = 0; i < ncols; i++) {
    perm(i) = i;
  }

  if (nrows > 250) {   // avoid cluttering output ...
    it_info_debug("Performing T-factorization of GF(2) matrix...  rows: "
                  << nrows << " cols: "  << ncols << " .... " << std::endl);
  }
  int pdone = 0;
  for (int j = 0; j < nrows; j++) {
    // Now working on diagonal element j,j
    // First try find a row with a 1 in column i
    int i1, j1;
    for (i1 = j; i1 < nrows; i1++) {
      for (j1 = j; j1 < ncols; j1++) {
        if (U.get(i1, j1) == 1) { goto found; }
      }
    }

    return j;

  found:
    U.swap_rows(i1, j);
    T.swap_rows(i1, j);
    U.swap_cols(j1, j);

    int temp = perm(j);
    perm(j) = perm(j1);
    perm(j1) = temp;

    // now subtract row i from remaining rows
    for (int i1 = j + 1; i1 < nrows; i1++)  {
      if (U.get(i1, j) == 1) {
        int ptemp = floor_i(100.0 * (i1 + j * nrows) / (nrows * nrows));
        if (nrows > 250 && ptemp > pdone + 10) {
          it_info_debug(ptemp << "% done.");
          pdone = ptemp;
        }
        U.add_rows(i1, j);
        T.add_rows(i1, j);
      }
    }
  }
  return nrows;
}


int GF2mat::T_fact_update_bitflip(GF2mat &T, GF2mat &U,
                                  ivec &perm, int rank,
                                  int r, int c) const
{
  // First, update U (before re-triangulization)
  int c0 = c;
  for (c = 0; c < ncols; c++) {
    if (perm(c) == c0) {
      goto foundidx;
    }
  }
  it_error("GF2mat::T_fact_update_bitflip() - internal error");

foundidx:
  for (int i = 0; i < nrows; i++) {
    if (T.get(i, r) == 1) {
      U.addto_element(i, c, 1);
    }
  }

  // first move column c to the end
  bvec lastcol = U.get_col(c);
  int temp_perm = perm(c);
  for (int j = c; j < ncols - 1; j++) {
    U.set_col(j, U.get_col(j + 1));
    perm(j) = perm(j + 1);
  }
  U.set_col(ncols - 1, lastcol);
  perm(ncols - 1) = temp_perm;

  // then, if the matrix is tall, also move row c to the end
  if (nrows >= ncols) {
    bvec lastrow_U = U.get_row(c);
    bvec lastrow_T = T.get_row(c);
    for (int i = c; i < nrows - 1; i++) {
      U.set_row(i, U.get_row(i + 1));
      T.set_row(i, T.get_row(i + 1));
    }
    U.set_row(nrows - 1, lastrow_U);
    T.set_row(nrows - 1, lastrow_T);

    // Do Gaussian elimination on the last row
    for (int j = c; j < ncols; j++)  {
      if (U.get(nrows - 1, j) == 1) {
        U.add_rows(nrows - 1, j);
        T.add_rows(nrows - 1, j);
      }
    }
  }

  // Now, continue T-factorization from the point (rank-1,rank-1)
  for (int j = rank - 1; j < nrows; j++) {
    int i1, j1;
    for (i1 = j; i1 < nrows; i1++) {
      for (j1 = j; j1 < ncols; j1++) {
        if (U.get(i1, j1) == 1) {
          goto found;
        }
      }
    }

    return j;

  found:
    U.swap_rows(i1, j);
    T.swap_rows(i1, j);
    U.swap_cols(j1, j);

    int temp = perm(j);
    perm(j) = perm(j1);
    perm(j1) = temp;

    for (int i1 = j + 1; i1 < nrows; i1++)  {
      if (U.get(i1, j) == 1) {
        U.add_rows(i1, j);
        T.add_rows(i1, j);
      }
    }
  }

  return nrows;
}

bool GF2mat::T_fact_update_addcol(GF2mat &T, GF2mat &U,
                                  ivec &perm, bvec newcol) const
{
  int i0 = T.rows();
  int j0 = U.cols();
  it_assert(j0 > 0, "GF2mat::T_fact_update_addcol(): dimension mismatch");
  it_assert(i0 == T.cols(),
            "GF2mat::T_fact_update_addcol(): dimension mismatch");
  it_assert(i0 == U.rows(),
            "GF2mat::T_fact_update_addcol(): dimension mismatch");
  it_assert(length(perm) == j0,
            "GF2mat::T_fact_update_addcol(): dimension mismatch");
  it_assert(U.get(j0 - 1, j0 - 1) == 1,
            "GF2mat::T_fact_update_addcol(): dimension mismatch");
  // The following test is VERY TIME-CONSUMING
  it_assert_debug(U.row_rank() == j0,
                  "GF2mat::T_fact_update_addcol(): factorization has incorrect rank");

  bvec z = T * newcol;
  GF2mat Utemp = U.concatenate_horizontal(GF2mat(z, 1));

  // start working on position (j0,j0)
  int i;
  for (i = j0; i < i0; i++) {
    if (Utemp.get(i, j0) == 1) {
      goto found;
    }
  }
  return (false); // adding the new column would not improve the rank

found:
  perm.set_length(j0 + 1, true);
  perm(j0) = j0;

  Utemp.swap_rows(i, j0);
  T.swap_rows(i, j0);

  for (int i1 = j0 + 1; i1 < i0; i1++)  {
    if (Utemp.get(i1, j0) == 1) {
      Utemp.add_rows(i1, j0);
      T.add_rows(i1, j0);
    }
  }

  U = Utemp;
  return (true); // the new column was successfully added
}




GF2mat GF2mat::inverse() const
{
  it_assert(nrows == ncols, "GF2mat::inverse(): Matrix must be square");

  // first compute the T-factorization
  GF2mat T, U;
  ivec perm;
  int rank = T_fact(T, U, perm);
  it_assert(rank == ncols, "GF2mat::inverse(): Matrix is not full rank");

  // backward substitution
  for (int i = ncols - 2; i >= 0; i--) {
    for (int j = ncols - 1; j > i; j--) {
      if (U.get(i, j) == 1) {
        U.add_rows(i, j);
        T.add_rows(i, j);
      }
    }
  }
  T.permute_rows(perm, 1);
  return T;
}

int GF2mat::row_rank() const
{
  GF2mat T, U;
  ivec perm;
  int rank = T_fact(T, U, perm);
  return rank;
}

bool GF2mat::is_zero() const
{
  for (int i = 0; i < nrows; i++) {
    for (int j = 0; j < nwords; j++) {
      if (data(i, j) != 0) {
        return false;
      }
    }
  }
  return true;
}

bool GF2mat::operator==(const GF2mat &X) const
{
  if (X.nrows != nrows) { return false; }
  if (X.ncols != ncols) { return false; }
  it_assert(X.nwords == nwords, "GF2mat::operator==() dimension mismatch");

  for (int i = 0; i < nrows; i++) {
    for (int j = 0; j < nwords; j++) {
      //      if (X.get(i,j)!=get(i,j)) {
      if (X.data(i, j) != data(i, j)) {
        return false;
      }
    }
  }
  return true;
}


void GF2mat::add_rows(int i, int j)
{
  it_assert(i >= 0 && i < nrows, "GF2mat::add_rows(): out of range");
  it_assert(j >= 0 && j < nrows, "GF2mat::add_rows(): out of range");
  for (int k = 0; k < nwords; k++)   {
    data(i, k) ^= data(j, k);
  }
}

void GF2mat::swap_rows(int i, int j)
{
  it_assert(i >= 0 && i < nrows, "GF2mat::swap_rows(): index out of range");
  it_assert(j >= 0 && j < nrows, "GF2mat::swap_rows(): index out of range");
  for (int k = 0; k < nwords; k++) {
    unsigned char temp = data(i, k);
    data(i, k) = data(j, k);
    data(j, k) = temp;
  }
}

void GF2mat::swap_cols(int i, int j)
{
  it_assert(i >= 0 && i < ncols, "GF2mat::swap_cols(): index out of range");
  it_assert(j >= 0 && j < ncols, "GF2mat::swap_cols(): index out of range");
  bvec temp = get_col(i);
  set_col(i, get_col(j));
  set_col(j, temp);
}


void GF2mat::operator=(const GF2mat &X)
{
  nrows = X.nrows;
  ncols = X.ncols;
  nwords = X.nwords;
  data = X.data;
}

GF2mat operator*(const GF2mat &X, const GF2mat &Y)
{
  it_assert(X.ncols == Y.nrows, "GF2mat::operator*(): dimension mismatch");
  it_assert(X.nwords > 0, "Gfmat::operator*(): dimension mismatch");
  it_assert(Y.nwords > 0, "Gfmat::operator*(): dimension mismatch");

  /*
  // this can be done more efficiently?
  GF2mat result(X.nrows,Y.ncols);
  for (int i=0; i<X.nrows; i++) {
  for (int j=0; j<Y.ncols; j++) {
  bin b=0;
  for (int k=0; k<X.ncols; k++) {
  bin x = X.get(i,k);
  bin y = Y.get(k,j);
  b ^= (x&y);
  }
  result.set(i,j,b);
  }
  }
  return result; */

  // is this better?
  return mult_trans(X, Y.transpose());
}

bvec operator*(const GF2mat &X, const bvec &y)
{
  it_assert(length(y) == X.ncols, "GF2mat::operator*(): dimension mismatch");
  it_assert(X.nwords > 0, "Gfmat::operator*(): dimension mismatch");

  /*
  // this can be done more efficiently?
  bvec result = zeros_b(X.nrows);
  for (int i=0; i<X.nrows; i++) {
  // do the nwords-1 data columns first
  for (int j=0; j<X.nwords-1; j++) {
  int ind = j<<shift_divisor;
  unsigned char r=X.data(i,j);
  while (r) {
  result(i) ^= (r & y(ind));
  r >>= 1;
  ind++;
  }
  }
  // do the last column separately
  for (int j=(X.nwords-1)<<shift_divisor; j<X.ncols; j++) {
  result(i) ^= (X.get(i,j) & y(j));
  }
  }
  return result; */

  // is this better?
  return (mult_trans(X, GF2mat(y, 0))).bvecify();
}

GF2mat mult_trans(const GF2mat &X, const GF2mat &Y)
{
  it_assert(X.ncols == Y.ncols, "GF2mat::mult_trans(): dimension mismatch");
  it_assert(X.nwords > 0, "GF2mat::mult_trans(): dimension mismatch");
  it_assert(Y.nwords > 0, "GF2mat::mult_trans(): dimension mismatch");
  it_assert(X.nwords == Y.nwords, "GF2mat::mult_trans(): dimension mismatch");

  GF2mat result(X.nrows, Y.nrows);

  for (int i = 0; i < X.nrows; i++) {
    for (int j = 0; j < Y.nrows; j++) {
      bin b = 0;
      int kindx = i;
      int kindy = j;
      for (int k = 0; k < X.nwords; k++) {
        unsigned char r = X.data(kindx) & Y.data(kindy);
        /* The following can be speeded up by using a small lookup
           table for the number of ones and shift only a few times (or
           not at all if table is large) */
        while (r) {
          b ^= r & 1;
          r >>= 1;
        };
        kindx += X.nrows;
        kindy += Y.nrows;
      }
      result.set(i, j, b);
    }
  }
  return result;
}

GF2mat GF2mat::transpose() const
{
  // CAN BE SPEEDED UP
  GF2mat result(ncols, nrows);

  for (int i = 0; i < nrows; i++) {
    for (int j = 0; j < ncols; j++) {
      result.set(j, i, get(i, j));
    }
  }
  return result;
}

GF2mat operator+(const GF2mat &X, const GF2mat &Y)
{
  it_assert(X.nrows == Y.nrows, "GF2mat::operator+(): dimension mismatch");
  it_assert(X.ncols == Y.ncols, "GF2mat::operator+(): dimension mismatch");
  it_assert(X.nwords == Y.nwords, "GF2mat::operator+(): dimension mismatch");
  GF2mat result(X.nrows, X.ncols);

  for (int i = 0; i < X.nrows; i++) {
    for (int j = 0; j < X.nwords; j++) {
      result.data(i, j) = X.data(i, j) ^ Y.data(i, j);
    }
  }

  return result;
}

void GF2mat::permute_cols(ivec &perm, bool I)
{
  it_assert(length(perm) == ncols,
            "GF2mat::permute_cols(): dimensions do not match");

  GF2mat temp = (*this);
  for (int j = 0; j < ncols; j++) {
    if (I == 0) {
      set_col(j, temp.get_col(perm(j)));
    }
    else {
      set_col(perm(j), temp.get_col(j));
    }
  }
}

void GF2mat::permute_rows(ivec &perm, bool I)
{
  it_assert(length(perm) == nrows,
            "GF2mat::permute_rows(): dimensions do not match");

  GF2mat temp = (*this);
  for (int i = 0; i < nrows; i++) {
    if (I == 0) {
      for (int j = 0; j < nwords; j++) {
        data(i, j) = temp.data(perm(i), j);
      }
    }
    else {
      for (int j = 0; j < nwords; j++) {
        data(perm(i), j) = temp.data(i, j);
      }
    }
  }
}


std::ostream &operator<<(std::ostream &os, const GF2mat &X)
{
  int i, j;
  os << "---- GF(2) matrix of dimension " << X.nrows << "*" << X.ncols
  << " -- Density: " << X.density() << " ----" << std::endl;

  for (i = 0; i < X.nrows; i++) {
    os << "      ";
    for (j = 0; j < X.ncols; j++) {
      os << X.get(i, j) << " ";
    }
    os << std::endl;
  }

  return os;
}

double GF2mat::density() const
{
  int no_of_ones = 0;

  for (int i = 0; i < nrows; i++) {
    for (int j = 0; j < ncols; j++) {
      no_of_ones += (get(i, j) == 1 ? 1 : 0);
    }
  }

  return ((double) no_of_ones) / (nrows*ncols);
}


it_file &operator<<(it_file &f, const GF2mat &X)
{
  // 3 64-bit unsigned words for: nrows, ncols and nwords + rest for char data
  uint64_t bytecount = 3 * sizeof(uint64_t)
                       + X.nrows * X.nwords * sizeof(char);
  f.write_data_header("GF2mat", bytecount);

  f.low_level_write(static_cast<uint64_t>(X.nrows));
  f.low_level_write(static_cast<uint64_t>(X.ncols));
  f.low_level_write(static_cast<uint64_t>(X.nwords));
  for (int i = 0; i < X.nrows; i++) {
    for (int j = 0; j < X.nwords; j++) {
      f.low_level_write(static_cast<char>(X.data(i, j)));
    }
  }
  return f;
}

it_ifile &operator>>(it_ifile &f, GF2mat &X)
{
  it_file::data_header h;

  f.read_data_header(h);
  if (h.type == "GF2mat") {
    uint64_t tmp;
    f.low_level_read(tmp);
    X.nrows = static_cast<int>(tmp);
    f.low_level_read(tmp);
    X.ncols = static_cast<int>(tmp);
    f.low_level_read(tmp);
    X.nwords = static_cast<int>(tmp);
    X.data.set_size(X.nrows, X.nwords);
    for (int i = 0; i < X.nrows; i++) {
      for (int j = 0; j < X.nwords; j++) {
        char r;
        f.low_level_read(r);
        X.data(i, j) = static_cast<unsigned char>(r);
      }
    }
  }
  else {
    it_error("it_ifile &operator>>() - internal error");
  }

  return f;
}

} // namespace itpp