util_m4ri.c

Go to the documentation of this file.

/************************************************************************ 
 * helper functions for use with m4ri library, 
 * including binary sparse matrices and conversion utilities 
 * author: Leonid Pryadko <leonid.pryadko@ucr.edu> 
 * some code borrowed from various sources 
 ************************************************************************/
#include <m4ri/m4ri.h>
#include <m4ri/mzd.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include "mmio.h"
 
#include "util_m4ri.h"
 
 
#ifdef __M4RI_MAX_MZD_BLOCKSIZE
/* gauss_one_rows() for rows not at equal distances in memory (M4RI before 2024, a matrix spanning several blocks):
 * each row is found with mzd_row().  Kept out of line, so that the inlined gauss_one_rows() stays small. */
int gauss_one_rows_blocks(mzd_t *M, const int idx, const int begrow, const rci_t nrows_eff) {
  const wi_t word_idx = idx >> 6;
  const word bit_mask = (word)1 << (idx & 63);
  const wi_t width = M->width;
  for (rci_t j = begrow; j < nrows_eff; ++j) {
    if (mzd_row(M, j)[word_idx] & bit_mask) {
      if (j != begrow) {
        mzd_row_swap(M, begrow, j);
      }
      const word * const piv_row = mzd_row(M, begrow);
      for (rci_t ii = 0; ii < nrows_eff; ++ii) {
        word * const dst_row = mzd_row(M, ii);
        if (ii != begrow && (dst_row[word_idx] & bit_mask)) {
          for (wi_t w = 0; w < width; ++w) {
            dst_row[w] ^= piv_row[w];
          }
        }
      }
      return 1;
    }
  }
  return 0;
}
#endif /* __M4RI_MAX_MZD_BLOCKSIZE */
 
size_t mzd_weight_naive(const mzd_t *A){
  size_t count = 0;
  for(rci_t i = 0; i < A->nrows; ++i)
    for(rci_t j = 0; j < A->ncols; ++j)
      if(mzd_read_bit(A, i, j))
    ++count;
  return (count);
}
 
size_t mzd_weight_row(const mzd_t * const A, const rci_t i){
  size_t count = 0;
  assert(i >=0 && i < A->nrows);
  if(A->width == 1) {
    for(rci_t j = 0; j < A->ncols; ++j)
      if(mzd_read_bit(A, i, j))
    ++count;
    return (count);
  }
  const word * const truerow = mzd_row_cons(A,i);
  for(wi_t j = 0; j < A->width - 1; j ++){
    count += m4ri_bitcount(truerow[j]);
  }
  for(rci_t j = m4ri_radix * (A-> width - 1) ; j < A->ncols ; ++j)
    if(mzd_read_bit(A, i,  j))
      ++count;
  return count ;
}
 
size_t mzd_weight(const mzd_t *A){
  size_t count = 0;
  if(A->width == 1) {
    for(rci_t i = 0; i < A->nrows; ++i)
      for(rci_t j = 0; j < A->ncols; ++j)
    if(mzd_read_bit(A, i, j))
      ++count;
    return (count);
  }
  for(rci_t i = 0; i < A->nrows; ++i) {
    const word * const truerow = mzd_row_cons(A,i);
    for(wi_t j = 0; j < A->width - 1; j ++){
      count += m4ri_bitcount(truerow[j]);
      //  count += std_bitcount(truerow[j]);
      //      printf(" %d",std_bitcount(truerow[j]));
    }
    for(rci_t j = m4ri_radix * (A-> width - 1) ; j < A->ncols ; ++j)
      if(mzd_read_bit(A, i,  j))
    ++count;
  }
  //  printf(" count=%d\n",count);
  return count ;
}
 

rci_t mzd_gauss_naive(mzd_t *M, mzp_t *q, int full) {
  rci_t startcol=0, startrow = 0;
  rci_t pivots = 0;
#if 0  
  if (q==NULL) 
    q=mzp_init(M->ncols);
  else {
    if ((q->length)!=M->ncols){ 
      mzp_free(q);
      q=mzp_init(M->ncols);
    }
    else 
      mzp_set_ui(q,1);
  }
#endif 
 
  for (rci_t i = startcol; i < M->ncols ; ++i) {
    for(rci_t j = startrow ; j < M->nrows; ++j) {
      if (mzd_read_bit(M, j, i)) {
    mzd_row_swap(M, startrow, j);
    q->values[pivots]=i; /* request cols swap: pivots <--> i */
    ++pivots;
    for(rci_t ii = full ? 0 : startrow + 1;  ii < M->nrows; ++ii) {
      if (ii != startrow) {
        if (mzd_read_bit(M, ii, i)) {
          mzd_row_add_offset(M, ii, startrow, i);
        }
      }
    }
    startrow = startrow + 1;
    break;
      }
    }
  }
  return pivots;
}

int csr_max_row_wght(const csr_t * const p){
  int m=p->rows;
  int wmax=0;
  for(int i=0;i<m;i++){
    int w=(p->p[i+1])-(p->p[i]);
    if (w>wmax)
      wmax=w;
    }
  return wmax;
}
 

csr_t * csr_transpose(csr_t *dst, const csr_t * const p){
  int rows=p->rows, cols=p->cols, nz=p->p[rows];
  if (dst == NULL) 
    dst = csr_init(NULL,cols,rows,nz);
  else if ((dst->cols != rows) || (dst->rows != cols) || (dst->nzmax < MAX(nz, cols+1))) 
    ERROR("Wrong size for return matrix.\n");
  else
    dst->nz=0; /* clear matrix */
  for(int i=0;i<rows;i++)    
    for(int j=p->p[i]; j < p->p[i+1] ; j++){
      dst->p[j]=p->i[j]; // pair format, to be compressed later
      dst->i[j]=i ;
    }
  dst->nz=nz;
  csr_compress(dst);
  return dst;
}
 

mzd_t *mzd_from_csr(mzd_t *dst, const csr_t *p) {
  int m=p->rows, n=p->cols, nz=p->nz;
  if (dst == NULL) 
    dst = mzd_init(m,n);
  else if ((dst->nrows != m) || (dst->ncols != n)) 
    ERROR("Wrong size for return matrix.\n");
  else
    mzd_set_ui(dst,0); /* clear bits */
 
  if(nz==-1){ // binary CSR matrix in Compressed Row Form
    for(int i=0;i<m;i++)
      for(int j=p->p[i]; j < p->p[i+1] ; j++)
    mzd_write_bit(dst, i, p->i[j], 1);
  }
  else{  // binary CSR matrix in Pair Form
    for(int i=0;i<p->nz;i++)
      mzd_write_bit(dst, p->p[i],p->i[i],1);
  }
  return dst;
}

mzd_t *mzd_generator_from_csr(mzd_t *G, const csr_t * const H){  
  rci_t n=H->cols;
  mzd_t *mat=mzd_from_csr(NULL,H); /* convert to dense matrix */
  mzp_t * pivots = mzp_init(n);    /* initialize the permutation */
  rci_t ra = mzd_gauss_naive(mat, pivots, 1);
#if 0  
  if ((prm.debug & 2048) && (n<=80)){ 
    mzd_print(mat);
    mzp_out(pivots);    
    printf("\n");
  }
#endif
  /* now we can check if G is the right size and reallocate if needed */
  if ((G!=NULL) && ((G->ncols!=n) || (G->nrows!=n-ra))){
    mzd_free(G); 
    G=NULL;
  }
  if (G==NULL)
    G=mzd_init(n-ra,n);    
  else
    mzd_set_ui(G,0); /* clear the matrix */
  if (ra > 0 && ra < n) {
    mzd_apply_p_right_trans(mat,pivots); /* permute columns to make std form [ I C ] */
    mzd_t *matC = mzd_submatrix(NULL, mat, 0, ra, ra, n); /* C size ra by n-ra */
    mzd_t * winCT=mzd_init_window(G,0,0,n-ra,ra);
    mzd_transpose(winCT,matC);
    mzd_free(winCT); /* window is no longer needed */
    mzd_free(matC);  /* \todo see why window does not seem to work here */
  }
  for(int i=0; i<n-ra ; i++)
    mzd_write_bit(G,i,i+ra,1); /* identity matrix on the right */
 
  if (ra < n)
    mzd_apply_p_right(G,pivots); /* permute columns back */
  mzp_free(pivots);
  mzd_free(mat);
  return G;
}
 

mzd_t * csr_mzd_mul(mzd_t *C, const csr_t *S, const mzd_t *B, int clear){
  if (S->cols != B->nrows)
    ERROR("column dim of S and row dim of B should match, %d=%d",S->cols,B->nrows);
  if (C == NULL) 
    C = mzd_init(S->rows, B->ncols);
  else {
    if (C->nrows != S->rows || C->ncols != B->ncols) 
      ERROR("Provided return matrix has wrong dimensions.\n");    
  if(clear)
    mzd_set_ui(C,0); 
  }
  rci_t const m = S->rows;
  
  for(rci_t i = 0; i < m; ++i)
    for(int j=S->p[i]; j < S->p[i+1] ; j++)
      mzd_combine(C,i,0, C,i,0, B,S->i[j],0); /* combine the two rows */
  return C;
}
 
void addto(mzd_t *row, const csr_t *spaQ, const int i){
#if 0  
  if (i>=spaQ->rows)
    ERROR("addto: attempt to get row=%d of %d",i,spaQ->rows);
  if (row->ncols != spaQ->cols)
    ERROR("addto: column number mismatch");
  //  mzd_print(row);
#endif
  for (int j=spaQ->p[i]; j<spaQ->p[i+1]; j++)
    mzd_flip_bit(row, 0,spaQ->i[j]); /* flip that bit */
}
 
 
static inline void addto_inline(mzd_t *row, const csr_t *spaQ, const int i){
#ifndef NDEBUG
  if (i>=spaQ->rows)
    ERROR("addto: attempt to get row=%d of %d",i,spaQ->rows);
  if (row->ncols != spaQ->cols)
    ERROR("addto: column number mismatch");
  //  mzd_print(row);
#endif
  for (int j=spaQ->p[i]; j<spaQ->p[i+1]; j++)
    mzd_flip_bit(row, 0,spaQ->i[j]); /* flip that bit */
}
 
 
mzd_t * syndrome_vector(mzd_t *syndrome, mzd_t *row, csr_t *spaQ, int clear){
  if (row->ncols != spaQ->rows)
    ERROR("column dim of row and row dim of spaQ should match, %d=%d",row->ncols,spaQ->rows);
  if (syndrome == NULL) 
    syndrome = mzd_init(row->nrows, spaQ->cols);
  else {
    if (syndrome->nrows != row->nrows || syndrome->ncols != spaQ->cols) 
      ERROR("Provided return matrix has wrong dimensions.\n");    
    if(clear)
      mzd_set_ui(syndrome,0); 
  }
  for (rci_t i=0; i< row->nrows; i++){
    const word * const rawrow = mzd_row_cons(row,i);  
    rci_t j=-1;
    rci_t n=row->ncols;
    do{
      j=nextelement(rawrow,row->width,j+1);
      if((j==-1)||(j>=n)) // no more non-zero elements 
    break; 
      else 
    for (int j1=spaQ->p[j]; j1<spaQ->p[j+1]; j1++)
      mzd_flip_bit(syndrome, i,spaQ->i[j1]); /* flip that bit */
    }
    while( TRUE);
  }
  return syndrome;
}
 
int csr_csr_mul_non_zero(const csr_t * const A, const csr_t * const B){
  if(!A)
    ERROR("matrix A is NULL");
  if(A->nz != -1)
    ERROR("matrix A shold be in compressed form");
  if(!B)
    ERROR("matrix B is NULL");
  if(B->nz != -1)
    ERROR("matrix B should be in compressed form");
  if (A->cols != B->cols) 
    ERROR("col count mismatch: A[%d,%d] and B[%d,%d]",
      A->rows,A->cols, B->rows, B->cols);
 
  /* for each row b of B, the parities of its overlaps with the rows of A, accumulated over the columns j of b from
   * the rows of A with a nonzero in column j (the transpose of A): O(nnz(A) + nnz(B) * column weight of A) operations
   * instead of O(rows(B) * nnz(A)).  A column repeated in a row counts once (the rows of both matrices are sorted), as
   * in mzd_from_csr(), which gives the dense matrices used for the logical operators (Lx_for_CSS_code()). */
  csr_t * const AT = csr_transpose(NULL, A);
  unsigned char * const par = calloc((size_t)A->rows + 1, 1); /* bit 0: the parity, bit 1: touched by the row b */
  int * const touched = malloc(((size_t)A->rows + 1) * sizeof(int));
  if (!par || !touched)
    ERROR("memory allocation");
  int nonzero = 0;
  for (int i = 0; (i < B->rows) && !nonzero; i++) {
    int nt = 0;
    for (int k = B->p[i]; k < B->p[i + 1]; k++) {
      const int j = B->i[k];
      if ((k > B->p[i]) && (j == B->i[k - 1]))
        continue;
      for (int t = AT->p[j]; t < AT->p[j + 1]; t++) {
        const int r = AT->i[t];
        if ((t > AT->p[j]) && (r == AT->i[t - 1]))
          continue;
        if (!(par[r] & 2)) {
          par[r] = 2;
          touched[nt++] = r;
        }
        par[r] ^= 1;
      }
    }
    for (int t = 0; t < nt; t++) {
      nonzero |= par[touched[t]] & 1;
      par[touched[t]] = 0;
    }
  }
  free(touched);
  free(par);
  csr_free(AT);
  return nonzero;
}
 
int syndrome_bit_count(const mzd_t * const row, const csr_t * const spaQ){
  int wei=0;
  int m=spaQ->rows;
  //  mzd_print(row);
  for(int i=0;i<m;i++){
    rci_t bit=0;
    for (int j=spaQ->p[i]; j<spaQ->p[i+1]; j++){
      //      printf("##### here ######## m=%d j=%d %d\n",m,j,spaQ->i[j]);      
      /* printf("checking %d %d: %d\n",i,spaQ->i[j],   */
      /*         mzd_read_bit(row, 0,spaQ->i[j]) ? 1 : 0);   */
      bit ^= mzd_read_bit(row, 0,spaQ->i[j]); 
    }
    wei+=bit;
    //    printf("wei=%d i=%d  bit=%d \n",wei,i,bit);
  }
  return wei;
}

size_t product_weight_csr_mzd(const csr_t *A, const mzd_t *B, int transpose){
  mzd_t *Prod;
  mzd_t *mA=mzd_from_csr(NULL,A);
  if (transpose){
    mzd_t *BT=mzd_transpose(NULL,B);
    Prod=mzd_mul_naive(NULL,mA,BT);
    mzd_free(BT);
  }
  else 
    Prod=mzd_mul_naive(NULL,mA,B);
  mzd_free(mA);
  //  mzd_print(Prod); printf("\n");
  size_t wei=mzd_weight(Prod);
  mzd_free(Prod);
  return wei;
}

int rand_uniform(const int max){
#if 1
  int divisor = RAND_MAX/(max);
  int retval;
  do  
    retval = rand() / divisor;
  while (retval >= max);
  return retval;
#else
  return (int) floor(max * tinymt64_generate_double(&tinymt));
#endif /* 0 */
}
 
mzp_t * mzp_rand_len(mzp_t *q, rci_t length){
  if (q==NULL)
    ERROR("permutation must be initialized!");
  if (length>q->length)
    ERROR("mzp_rand_len: second parameter too large, %d > %d", length, q->length);
  for(int i=0;i<=length-2;i++){
    q->values[i]= i+rand_uniform(length-i);
  }
  for(int i=length-1; i< (q-> length); i++)
    q->values[i]=i; /* no swap for the remaining columns */
  return q;
}

void mzp_out(mzp_t const *p){
  printf("[");
  for(rci_t i=0;i<p->length;i++)
    printf(" %d",p->values[i]);
  printf(" ]\n");
}

mzp_t *perm_p(mzp_t *q, const mzp_t *p,rci_t start){
  assert(p!=NULL);
  if (q==NULL)
    q=mzp_init(p->length);
  else if  (p->length != q->length)
    ERROR("mzp length mismatch %d and %d !\n", p->length , q->length);
  rci_t length = p->length;
  rci_t *perm = q->values; 
  for(rci_t i = start; i < length; ++i) 
    SWAPINT(perm[i], perm[p->values[i]]);
  return q;
}

mzp_t *perm_p_trans(mzp_t *q, const mzp_t *p,const rci_t start){
  assert(p!=NULL);
  if (q==NULL)
    q=mzp_init(p->length);
  else if  (p->length != q->length)
    ERROR("mzp length mismatch %d and %d !\n", p->length , q->length);
  rci_t length = p->length;  
  rci_t *perm = q->values; 
  for(rci_t i = start; i < length; ++i)
    SWAPINT(perm[length-i-1], perm[p->values[length-i-1]]);
  return q;
}
 

csr_t *csr_free(csr_t *p){
  if(p!=NULL){
    free(p->i);
    free(p->p);
    p->nzmax=0;
    p->nz=0;
    p->rows=0;
    p->cols=0;
    free(p);
  }
  return NULL;
}

csr_t *csr_init(csr_t *mat, int rows, int cols, int nzmax){
  if ((mat!=NULL)&&((mat->nzmax < nzmax)||(mat->nzmax < rows+1))){
    // mat=csr_free(mat);  /* allocated size was too small */
    mat->p = realloc(mat->p,MAX(nzmax,(rows+1))*sizeof(int));
    mat->i = realloc(mat->i, nzmax*sizeof(int));
    if ((mat->p==NULL) || (mat->i==NULL))
      ERROR("csr_init: failed to reallocate CSR rows=%d cols=%d nzmax=%d",
            rows,cols,nzmax);
    mat->nzmax=nzmax;
  }
  else if(mat==NULL){
    mat=malloc(sizeof(csr_t));  
    mat->p=calloc(MAX(nzmax,(rows+1)),sizeof(int));
    mat->i=calloc(nzmax,sizeof(int)); 
    if ((mat == NULL) || (mat->p==NULL) || (mat->i==NULL))
      ERROR("csr_init: failed to allocate CSR rows=%d cols=%d nzmax=%d",
            rows,cols,nzmax);
    mat->nzmax=nzmax;
  }
  mat->rows=rows;
  mat->cols=cols;
  mat->nz=0; /* empty */
  return mat;
}
 
#define SORT_NAME int_pair
#define SORT_TYPE int_pair
#define SORT_CMP(x, y) (((x).a != (y).a) ? ((x).a - (y).a) : ((x).b - (y).b))
#define SORT_DEF static inline
#include "sort.h"
 
void csr_compress(csr_t *mat){
  int nz=mat->nz;
  if(nz==-1)
    ERROR("matrix already compressed");
  int_pair *pairs=calloc(nz,sizeof(int_pair));
  for(int i=0;i<nz;i++){
    pairs[i].a=mat->p[i];
    pairs[i].b=mat->i[i];
  }
  int_pair_quick_sort(pairs, nz);
  int i, j=0;
  for(i=0;i<mat->rows;i++){
    mat->p[i]=j;
    while ((j<nz)&&(pairs[j].a == i)){
      mat->i[j]=pairs[j].b;
      j++;
    }
  }
  mat->p[i]=j; /* final value */
  mat->nz=-1; /* indicate compressed form */
  free(pairs);
}
 
csr_t * csr_from_pairs(csr_t *mat, const int nz, int_pair * const prs, const int nrows, const int ncols){
  mat = csr_init(mat, nrows, ncols, nz);
  int_pair_quick_sort(prs, nz);
  int i, j=0;
  for(i=0; i < nrows; i++){
    mat->p[i]=j;
    while ((j<nz) && (prs[j].a == i)){
      mat->i[j]=prs[j].b;
      j++;
    }
  }
  mat->p[i]=j; /* final value */
  mat->nz=-1; /* indicate compressed form */
  return mat;
}
 
void csr_out(const csr_t *mat){
  int nz=mat->nz;
  if(nz==-1){
    printf("# binary CSR matrix (%d x %d) in Compressed Row Form: "
       "nzmax= %d, nz=%d\n",
       mat->rows,mat->cols,mat->nzmax,mat->p[mat->rows]);
    for(int i=0;i<mat->rows;i++){
      printf("%d:",i+1);
      for(int j=mat->p[i]; j < mat->p[i+1] ; j++)
    printf(" %d",mat->i[j]+1);
      printf("\n");      
    }
  }
  else{ 
    printf("# binary CSR matrix (%d x %d) in Pair Form: "
       "nzmax= %d, nz=%d\n",
       mat->rows,mat->cols,mat->nzmax,mat->nz);
    for(int i=0;i<mat->nz;i++)
      printf("%d %d\n",mat->p[i]+1,mat->i[i]+1);
  }
}
 
void csr_print(const csr_t * const smat, const char str[]){
  mzd_t *mH0 = mzd_from_csr(NULL,smat);
  if(str)
    printf("matrix %s:\n",str);
  else
    printf("matrix:\n");
  mzd_print(mH0);
  mzd_free(mH0);
}
 
 

csr_t *csr_mm_read(char *fin, csr_t *mat, int transpose){
  int ret_code;
  MM_typecode matcode;
  FILE *f;
  int M, N, nz;   
 
  if ((f = fopen(fin, "r")) == NULL) 
    ERROR("can't open file %s",fin);
 
  if (mm_read_banner(f, &matcode) != 0)
    ERROR("Could not process Matrix Market banner.");
 
  if (!mm_is_matrix(matcode))
    ERROR("Only Matrix Market 'matrix' object is supported.");
 
  if (!mm_is_integer(matcode) && !mm_is_pattern(matcode)) {
    ERROR("Only integer or pattern matrices are supported. Found: %s", mm_typecode_to_str(matcode));
  }
 
  if (!mm_is_general(matcode) && !mm_is_symmetric(matcode)) {
    ERROR("Only general or symmetric matrices are supported. Found: %s", mm_typecode_to_str(matcode));
  }
 
  /* find out size of matrix .... */
  if (mm_is_coordinate(matcode)) {
    if ((ret_code = mm_read_mtx_crd_size(f, &M, &N, &nz)) != 0)
      ERROR("Cannot read coordinate size in input file %s", fin);
  } else { // array (dense)
    if ((ret_code = mm_read_mtx_array_size(f, &M, &N)) != 0)
      ERROR("Cannot read array size in input file %s", fin);
    if (mm_is_symmetric(matcode)) {
      nz = M * M;
    } else {
      nz = M * N;
    }
  }
 
  int nzmax;
  if (mm_is_coordinate(matcode)) {
    if (mm_is_symmetric(matcode)) {
      nzmax = 2 * nz;
    } else {
      nzmax = nz;
    }
  } else { // array
    if (mm_is_symmetric(matcode)) {
      nzmax = M * M;
    } else {
      nzmax = M * N;
    }
  }
 
  int rows = M;
  int cols = N;
  if (transpose) {
    rows = N;
    cols = M;
  }
 
  mat = csr_init(mat, rows, cols, nzmax);
  int k = 0;
 
  if (mm_is_coordinate(matcode)) {
    for (int i = 0; i < nz; i++) {
      int r, c;
      int v = 1; // default for pattern
      if (mm_is_integer(matcode)) {
        if (fscanf(f, "%d %d %d\n", &r, &c, &v) != 3) {
          ERROR("Failed to read coordinate entry %d in %s", i, fin);
        }
      } else { // pattern
        if (fscanf(f, "%d %d\n", &r, &c) != 2) {
          ERROR("Failed to read pattern entry %d in %s", i, fin);
        }
      }
      r--; c--; // 1-based to 0-based
 
      if (v % 2 != 0) {
        int r_store = transpose ? c : r;
        int c_store = transpose ? r : c;
        mat->p[k] = r_store;
        mat->i[k] = c_store;
        k++;
 
        if (mm_is_symmetric(matcode) && r != c) {
          mat->p[k] = c_store;
          mat->i[k] = r_store;
          k++;
        }
      }
    }
  } else { // array (dense)
    if (mm_is_general(matcode)) {
      for (int c = 0; c < N; c++) {
        for (int r = 0; r < M; r++) {
          int v;
          if (fscanf(f, "%d\n", &v) != 1) {
            ERROR("Failed to read dense general entry at r=%d, c=%d in %s", r, c, fin);
          }
          if (v % 2 != 0) {
            int r_store = transpose ? c : r;
            int c_store = transpose ? r : c;
            mat->p[k] = r_store;
            mat->i[k] = c_store;
            k++;
          }
        }
      }
    } else if (mm_is_symmetric(matcode)) {
      for (int c = 0; c < N; c++) {
        for (int r = c; r < M; r++) {
          int v;
          if (fscanf(f, "%d\n", &v) != 1) {
            ERROR("Failed to read dense symmetric entry at r=%d, c=%d in %s", r, c, fin);
          }
          if (v % 2 != 0) {
            int r_store = transpose ? c : r;
            int c_store = transpose ? r : c;
            mat->p[k] = r_store;
            mat->i[k] = c_store;
            k++;
 
            if (r != c) {
              mat->p[k] = c_store;
              mat->i[k] = r_store;
              k++;
            }
          }
        }
      }
    }
  }
 
  mat->nz = k;
  csr_compress(mat); /* sort entries by row and remove duplicates */
  fclose(f);
  return mat;
}
 

csr_t *csr_apply_perm(csr_t *dst, const csr_t * const src, const mzp_t * const perm){
  int m=src->rows, n=src->cols;
  if (src->nz!=-1)
    ERROR("pair format unsupported, nz=%d; expected \"-1\"",src->nz);
  int nz=src->p[m];
  if(src->cols != perm-> length)
    ERROR("perm length=%d should match cols=%d",perm->length,src->cols);
  if (dst == NULL) 
    dst = csr_init(dst,m,n,nz);
  else if ((dst->rows != m) || (dst->cols != n) || (nz > dst->nzmax)) 
    ERROR("Wrong size for return matrix.\n");
  for(int i=0;i<m;i++)    
    for(int j=src->p[i]; j < src->p[i+1] ; j++){
      dst->p[j]=i; // pair format, to be compressed later
      dst->i[j]= perm->values[src->i[j]];
      //      printf("row %d: mapping %d to %d\n",i,src->i[j],perm->values[src->i[j]]);
    }
  dst->nz=nz;
  //  csr_out(dst);
  csr_compress(dst);
  //csr_out(dst);
  return dst;
}
 
 

 
int do_reduce(mzd_t *row, const mzd_t *matG0, const rci_t rankG0){
  word * const rawrow = mzd_row(row,0);  
  rci_t j=-1;
  rci_t n=row->ncols;
  do{
    j=nextelement(rawrow,row->width,j);
    if(j==-1) // empty line after simplification
      return j; 
    else if (j<rankG0)
      mzd_combine_even_in_place(row,0,0,matG0,j,0);
  } while (j < rankG0);
  if (j<n)
    return j;
  return -1;
}
 
static inline uint64_t rotl(const uint64_t x, int k) {
  return (x << k) | (x >> (64 - k));
}
 
/*** see arXiv:1805.01407 on fast RNGs ****
 
uint64_t s[4]; 
const int A=17, B=45, S=5, R=7, T=9;
const uint64_t result_plus = s[0] + s[3];
const uint64_t result_starstar = rotl(s[1] * S, R) * T;
const uint64_t t = s[1] << A;
s[2] ^= s[0];
s[3] ^= s[1];
s[1] ^= s[2];
s[0] ^= s[3];
s[2] ^= t;
s[3] = rotl(s[3], B);
*/
 
/*
 * generate binary error vector with error probability p 
 */
void make_err(mzd_t *row, double p){
  //  const int max=MIN(10000,RAND_MAX);
  const int cutoff=p*RAND_MAX;
  if (row==NULL)
    ERROR("make_err: row vector must be allocated to correct size!");
  //  mzd_set_ui(row,0);
  //  printf(" rows=%d cols=%d\n",row->nrows,row->ncols);
  for(rci_t i = 0; i < row->nrows; ++i){
    mzd_row_clear_offset(row,i,0);
    printf("i=%d\n",i); // mzd_print(row);
    for(rci_t j = 0; j < row->ncols; ++j){
      //      printf("%d %d\t",j,rand_uniform(max));
      if (rand()<cutoff)
    mzd_write_bit(row,i,j,1);
      //      printf("%d %d -> %d %d\n",i,j,mzd_read_bit(row,i,j),cutoff);
    }
    //    mzd_print(row);
  }
}
 

csr_t * csr_from_mzd(csr_t *mat, const mzd_t * const orig){
  int nz=mzd_weight(orig);
  mat = csr_init(mat, orig->nrows, orig->ncols, nz);
  int i, j=0;
  for(i=0;i < mat->rows; i++){
    mat->p[i]=j;
#if 1 
    int idx=0;
    const word * const rawrow = mzd_row_cons(orig,i);
    while(((idx=nextelement(rawrow,orig->width,idx))!=-1)&&(idx<orig->ncols)){
      mat->i[j++]=idx++;
    }
#else 
    for(int idx=0; idx< orig->ncols; idx++)
      if(mzd_read_bit(orig,i,idx)){
    mat->i[j++]=idx;
    //  printf("i=%d j=%d idx=%d\n",i,j,idx);
      }
#endif /* 0 */    
  }      
  mat->p[i]=j; /* final value */
  assert(j==nz);
  mat->nz=-1; /* indicate compressed form */
  return mat;
}

 
csr_t * Lx_for_CSS_code(const csr_t * const Hx, const csr_t *const Hz){
  rci_t n=Hx->cols;
  if(n!=Hz->cols)
    ERROR("unequal number of columns in matrices Hx[%d,%d] and Hz[%d,%d]",
      Hx->rows, Hx->cols, Hz->rows, Hz->cols);  
  mzd_t *Mx = mzd_from_csr(NULL,Hx); /* convert to dense matrix */
  mzp_t *pivots = mzp_init(n);     /* initialize the permutation */
  rci_t rank = mzd_gauss_naive(Mx, pivots, 1);
  mzd_apply_p_right_trans(Mx,pivots); 
  mzp_t *perm  = perm_p(      NULL, pivots,0);
  mzp_t *permT = perm_p_trans(NULL, pivots,0);
  mzp_free(pivots);
  csr_t* HzPerm=csr_apply_perm(NULL,Hz,permT);
  mzp_free(permT);
#ifndef NDEBUG 
  mzd_t *MxT = mzd_transpose(NULL,Mx); 
  if(product_weight_csr_mzd(HzPerm,MxT,0))
    ERROR("rows of Hx and Hz should be orthogonal ");
  mzd_free(MxT);
#endif
  mzd_t *MzStar = mzd_generator_from_csr(NULL, HzPerm);
  csr_free(HzPerm);
 
  mzd_t *mat = mzd_stack(NULL,Mx,MzStar);
  mzd_free(MzStar);
  //  printf("stacked:\n");  mzd_print(mat);
  int rank1=mzd_echelonize(mat,0);
  //  printf("after gauss %d -> %d:\n", rank, rank1); mzd_print(mat);
 
  mzd_t *window = mzd_init_window(mat,rank,0,rank1,n);
  csr_t *Lx=csr_from_mzd(NULL, window);
  csr_t *ans = csr_apply_perm(NULL,Lx,perm);
  mzd_free(Mx);
  csr_free(Lx);  
  mzd_free(window);
  mzp_free(perm);
  mzd_free(mat);
  return ans;
}
 
mzd_t * mzd_nullspace(const csr_t * const H) {
  if (!H) return NULL;
  mzd_t *mH = mzd_from_csr(NULL, H);
  mzd_t *X = mzd_kernel_left_pluq(mH, 0);
  mzd_free(mH);
  if (!X) return NULL;
  mzd_t *N = mzd_transpose(NULL, X);
  mzd_free(X);
  return N;
}
 
colmask_t * colmask_from_csr(const csr_t * const A) {
  if (!A || A->nz != -1) ERROR("expected a sparse matrix in CSR form");
  colmask_t *cm = malloc(sizeof(colmask_t));
  if (!cm) ERROR("memory allocation");
  cm->rows = A->rows;
  cm->cols = A->cols;
  cm->kw = (A->rows + 63) >> 6;
  const size_t num = (size_t)A->cols * (size_t)cm->kw;
  cm->mask = calloc(num > 0 ? num : 1, sizeof(word));
  if (!cm->mask) ERROR("memory allocation");
  for (int ir = 0; ir < A->rows; ir++) {
    const word bit = (word)1 << (ir & 63);
    for (int iL = A->p[ir]; iL < A->p[ir + 1]; iL++)
      cm->mask[(size_t)A->i[iL] * cm->kw + (ir >> 6)] ^= bit; /* duplicate entries cancel */
  }
  return cm;
}
 
colmask_t * colmask_free(colmask_t *cm) {
  if (cm) {
    free(cm->mask);
    free(cm);
  }
  return NULL;
}
 
cols_ws_t * cols_ws_init(const int m) {
  cols_ws_t *ws = calloc(1, sizeof(cols_ws_t));
  if (!ws) ERROR("memory allocation");
  ws->m = m;
  ws->rmap = malloc((size_t)(m > 0 ? m : 1) * sizeof(int));
  ws->rlist = malloc((size_t)(m > 0 ? m : 1) * sizeof(int));
  if (!ws->rmap || !ws->rlist) ERROR("memory allocation");
  for (int r = 0; r < m; r++) ws->rmap[r] = -1;
  return ws;
}
 
cols_ws_t * cols_ws_free(cols_ws_t *ws) {
  if (ws) {
    free(ws->rmap);
    free(ws->rlist);
    free(ws->bits);
    free(ws);
  }
  return NULL;
}
 
/* The rows of H touched by the columns cols[0..cnt-1], numbered in ws->rmap (ws->rlist: the rows); returns their
 * number.  The caller restores rmap with cols_ws_clear(). */
static int cols_ws_rows(const csr_t * const HT, const int cnt, const int cols[], cols_ws_t * const ws) {
  int nr = 0;
  for (int i = 0; i < cnt; i++) {
    for (int k = HT->p[cols[i]]; k < HT->p[cols[i] + 1]; k++) {
      const int r = HT->i[k];
      if (ws->rmap[r] < 0) {
        ws->rmap[r] = nr;
        ws->rlist[nr++] = r;
      }
    }
  }
  return nr;
}
 
static void cols_ws_clear(cols_ws_t * const ws, const int nr) {
  for (int i = 0; i < nr; i++) ws->rmap[ws->rlist[i]] = -1;
}
 
/* cnt rows of `stride` words in ws->bits, set to zero; the column i in the first words of row i (its bits: the rows
 * touched, see cols_ws_rows()) */
static word * cols_ws_matrix(const csr_t * const HT, const int cnt, const int cols[], cols_ws_t * const ws,
                             const size_t stride) {
  const size_t num = (size_t)cnt * stride;
  if (num > ws->cap) {
    free(ws->bits);
    ws->cap = (num > 2 * ws->cap) ? num : 2 * ws->cap;
    ws->bits = malloc(ws->cap * sizeof(word));
    if (!ws->bits) ERROR("memory allocation");
  }
  word * const M = ws->bits;
  if (num > 0) memset(M, 0, num * sizeof(word));
  for (int i = 0; i < cnt; i++) {
    for (int k = HT->p[cols[i]]; k < HT->p[cols[i] + 1]; k++) {
      const int b = ws->rmap[HT->i[k]];
      M[(size_t)i * stride + (b >> 6)] ^= (word)1 << (b & 63); /* duplicate entries cancel */
    }
  }
  return M;
}
 
/* Gaussian elimination of the cnt rows of M (`stride` words each) on their first nr bits; returns the rank, the
 * independent rows first.  Rows rank..cnt-1 are then zero in the first nr bits (the other words of the rows, e.g.,
 * the tracked combinations, are transformed alike). */
static int cols_ws_eliminate(word * const M, const int cnt, const int nr, const size_t stride) {
  int rank = 0;
  for (int b = 0; b < nr && rank < cnt; b++) {
    const size_t wb = (size_t)(b >> 6);
    const word mb = (word)1 << (b & 63);
    int piv = -1;
    for (int i = rank; i < cnt; i++) {
      if (M[(size_t)i * stride + wb] & mb) {
        piv = i;
        break;
      }
    }
    if (piv < 0) continue;
    word * const pr = M + (size_t)rank * stride;
    if (piv != rank) {
      word * const qr = M + (size_t)piv * stride;
      for (size_t x = 0; x < stride; x++) {
        const word t = pr[x];
        pr[x] = qr[x];
        qr[x] = t;
      }
    }
    /* (the rows rank.. are zero in the bits below b, so that the words below wb need no update) */
    for (int i = rank + 1; i < cnt; i++) {
      word * const ri = M + (size_t)i * stride;
      if (ri[wb] & mb)
        for (size_t x = wb; x < stride; x++) ri[x] ^= pr[x];
    }
    rank++;
  }
  return rank;
}
 
int csr_cols_rank(const csr_t * const HT, const int cnt, const int cols[], cols_ws_t * const ws) {
  if (cnt <= 0) return 0;
  const int nr = cols_ws_rows(HT, cnt, cols, ws);
  const size_t nw = (size_t)((nr + 63) >> 6);
  int rank = 0;
  if (nw > 0) {
    word * const M = cols_ws_matrix(HT, cnt, cols, ws, nw);
    rank = cols_ws_eliminate(M, cnt, nr, nw);
  }
  cols_ws_clear(ws, nr);
  return rank;
}
 
int csr_cols_dependency(const csr_t * const HT, const int cnt, const int cols[], cols_ws_t * const ws, int sub[]) {
  if (cnt <= 1) return 0;
  const int nr = cols_ws_rows(HT, cnt, cols, ws);
  const size_t nw = (size_t)((nr + 63) >> 6), tw = (size_t)((cnt + 63) >> 6);
  const size_t stride = nw + tw; /* the column on the rows touched, then the combination of the columns */
  word * const M = cols_ws_matrix(HT, cnt, cols, ws, stride);
  cols_ws_clear(ws, nr);
  for (int i = 0; i < cnt; i++) M[(size_t)i * stride + nw + (size_t)(i >> 6)] |= (word)1 << (i & 63);
  const int rank = cols_ws_eliminate(M, cnt, nr, stride);
  /* the rows rank..cnt-1: a basis of the kernel (their combinations have zero sums); one of them is neither zero nor
   * the all-ones vector if the kernel has dimension at least 2 */
  for (int i = rank; i < cnt; i++) {
    const word * const t = M + (size_t)i * stride + nw;
    int k = 0;
    for (int j = 0; j < cnt; j++)
      if ((t[j >> 6] >> (j & 63)) & 1) k++;
    if (k == 0 || k == cnt) continue;
    k = 0;
    for (int j = 0; j < cnt; j++)
      if ((t[j >> 6] >> (j & 63)) & 1) sub[k++] = cols[j];
    return k;
  }
  return 0;
}
 
static inline uint64_t util_splitmix64(uint64_t *state) {
  uint64_t z = (*state += 0x9e3779b97f4a7c15ULL);
  z = (z ^ (z >> 30)) * 0xbf58476d1ce4e5b9ULL;
  z = (z ^ (z >> 27)) * 0x94d049bb133111ebULL;
  return z ^ (z >> 31);
}
 
static inline int util_rand_uniform(int max, uint64_t *state) {
  if (max <= 1) return 0;
  return (int)(util_splitmix64(state) % (uint64_t)max);
}
 
void localized_window_perm(mzp_t *perm, int nvar, int seed_col, int kwin, int win_mode,
                           const csr_t *spaH, const csr_t *mHT,
                           int *visited_cols, int *visited_checks, int *col_queue,
                           int *visit_marker, uint64_t *rng_state) {
  if (kwin <= 0 || kwin >= nvar || !visited_cols) {
    for (int i = 0; i < nvar; i++) perm->values[i] = i;
    for (int i = 0; i < nvar - 1; i++) {
      int j = i + util_rand_uniform(nvar - i, rng_state);
      SWAPINT(perm->values[i], perm->values[j]);
    }
    return;
  }
 
  (*visit_marker)++;
  int marker = *visit_marker;
  if (marker <= 0) {
    memset(visited_cols, 0, (size_t)nvar * sizeof(int));
    if (spaH && visited_checks) {
      memset(visited_checks, 0, (size_t)spaH->rows * sizeof(int));
    }
    *visit_marker = marker = 1;
  }
 
  int count = 0;
  if (win_mode == 0 && spaH && mHT && visited_checks && col_queue) {
    int head = 0, tail = 0;
    visited_cols[seed_col] = marker;
    col_queue[tail++] = seed_col;
    perm->values[count++] = seed_col;
 
    while (count < kwin) {
      if (head < tail) {
        int col = col_queue[head++];
        for (int p_idx = mHT->p[col]; p_idx < mHT->p[col + 1] && count < kwin; p_idx++) {
          int chk = mHT->i[p_idx];
          if (visited_checks[chk] != marker) {
            visited_checks[chk] = marker;
            for (int h_idx = spaH->p[chk]; h_idx < spaH->p[chk + 1] && count < kwin; h_idx++) {
              int c_adj = spaH->i[h_idx];
              if (visited_cols[c_adj] != marker) {
                visited_cols[c_adj] = marker;
                col_queue[tail++] = c_adj;
                perm->values[count++] = c_adj;
              }
            }
          }
        }
      } else {
        int r_start = util_rand_uniform(nvar, rng_state);
        int found = -1;
        for (int offset = 0; offset < nvar; offset++) {
          int cand = (r_start + offset) % nvar;
          if (visited_cols[cand] != marker) {
            found = cand;
            break;
          }
        }
        if (found == -1) break;
        visited_cols[found] = marker;
        col_queue[tail++] = found;
        perm->values[count++] = found;
      }
    }
  } else {
    for (int i = 0; i < kwin; i++) {
      int col = (seed_col + i) % nvar;
      visited_cols[col] = marker;
      perm->values[count++] = col;
    }
  }
 
  int win_actual = count;
  for (int i = 0; i < win_actual - 1; i++) {
    int j = i + util_rand_uniform(win_actual - i, rng_state);
    SWAPINT(perm->values[i], perm->values[j]);
  }
 
  int rem_idx = win_actual;
  for (int col = 0; col < nvar; col++) {
    if (visited_cols[col] != marker) {
      perm->values[rem_idx++] = col;
    }
  }
  for (int i = win_actual; i < nvar - 1; i++) {
    int j = i + util_rand_uniform(nvar - i, rng_state);
    SWAPINT(perm->values[i], perm->values[j]);
  }
}
 
void refresh_nullspace_basis(mzd_t *N, const int *cw_arr, int cw_wt, uint64_t *rng_state) {
  if (!N || N->nrows == 0 || N->ncols == 0) return;
  const rci_t nu = N->nrows;
  const rci_t n = N->ncols;
 
  mzp_t *pivs = mzp_init(n);
  mzp_t *perm = mzp_init(n);
  rci_t *piv_cols = malloc((size_t)nu * sizeof(rci_t));
  if (!pivs || !perm || !piv_cols) {
    if (pivs) mzp_free(pivs);
    if (perm) mzp_free(perm);
    free(piv_cols);
    return;
  }
 
  for (int i = 0; i <= (int)n - 2; i++) {
    pivs->values[i] = i + util_rand_uniform(n - i, rng_state);
  }
  pivs->values[n - 1] = n - 1;
 
  mzp_set_ui(perm, 1);
  perm = perm_p_trans(perm, pivs, 0);
 
  rci_t rank = 0;
  for (rci_t i = 0; i < n && rank < nu; i++) {
    int col = perm->values[i];
    if (gauss_one(N, col, rank)) {
      piv_cols[rank++] = col;
    }
  }
 
  int target_row = -1;
  size_t max_row_wt = 0;
  if (cw_arr && cw_wt > 0) {
    for (rci_t r = 0; r < rank; r++) {
      int orig_col = piv_cols[r];
      for (int k = 0; k < cw_wt; k++) {
        if (cw_arr[k] == orig_col) {
          size_t r_wt = mzd_weight_row(N, r);
          if (r_wt > (size_t)cw_wt && r_wt > max_row_wt) {
            max_row_wt = r_wt;
            target_row = r;
          }
          break;
        }
      }
    }
  }
 
  if (target_row >= 0 && cw_arr && cw_wt > 0) {
    mzd_row_clear_offset(N, target_row, 0);
    for (int k = 0; k < cw_wt; k++) {
      if (cw_arr[k] >= 0 && cw_arr[k] < n) {
        mzd_write_bit(N, target_row, cw_arr[k], 1);
      }
    }
  }
 
  free(piv_cols);
  mzp_free(perm);
  mzp_free(pivs);
}
 
int csr_all_ones_in_rowspan(const csr_t * const H) {
  if (!H || H->rows <= 0 || H->cols <= 0 || H->nz != -1) return 0;
  const int m = H->rows, n = H->cols;
  char * const has_nz = calloc((size_t)n, 1);
  if (!has_nz) ERROR("memory allocation");
  int covered = 0;
  for (int t = 0; t < H->p[m]; t++) {
    const int c = H->i[t];
    if (c >= 0 && c < n && !has_nz[c]) {
      has_nz[c] = 1;
      covered++;
    }
  }
  free(has_nz);
  if (covered < n) return 0;
  mzd_t * const mH = mzd_from_csr(NULL, H);
  const rci_t rank_H = mzd_echelonize(mH, 0);
  mzd_t * const m1 = mzd_init(1, n);
  for (int c = 0; c < n; c++) mzd_write_bit(m1, 0, c, 1);
  mzd_t * const mAug = mzd_stack(NULL, mH, m1);
  mzd_free(mH);
  mzd_free(m1);
  const rci_t rank_Aug = mzd_echelonize(mAug, 0);
  mzd_free(mAug);
  return rank_Aug == rank_H;
}