util_m4ri.h

Go to the documentation of this file.

#ifndef UTIL_M4RI_H
#define UTIL_M4RI_H
 
/************************************************************************ 
 * 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 
 ************************************************************************/
 
#define SWAPINT(a,b) do{ int t=a; a=b; b=t; } while(0)
 
#define ERROR(fmt,...) \
  do{ \
    fprintf (stderr, "%s:%d: *** ERROR in function '%s()' ***\n", __FILE__, __LINE__, __FUNCTION__); \
    fprintf(stderr, " ␛[31;1m " fmt " ␛[0m\n",##__VA_ARGS__); \
    exit(-1); \
  } \
  while(0)
 
static inline word const * mzd_row_cons(const mzd_t * mat, const int row){
  //  return mat->rows[row] ;
#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wdiscarded-qualifiers"
  return mzd_row(mat,row);
#pragma GCC diagnostic pop
}
 
 
#define SETWD(pos) ((pos)>>6)
#define SETBT(pos) ((pos)&0x3F)
#define TIMESWORDSIZE(w) ((w)<<6)    /* w*WORDSIZE */
 
#define FIRSTBIT(x) __builtin_ctzll(x) // number of trailing zeros 
 
// #ifdef __POPCNT__
#if 1
 
/* 
 * optimized code copied verbatim from https://danluu.com/assembly-intrinsics/ 
 */
/* uint32_t builtin_popcnt_unrolled_errata_manual(const uint64_t* buf, int len) { */
/*   assert(len % 4 == 0); */
/*   uint64_t cnt[4]; */
/*   for (int i = 0; i < 4; ++i) { */
/*     cnt[i] = 0; */
/*   } */
 
/*   for (int i = 0; i < len; i+=4) { */
/*     __asm__( */
/*         "popcnt %4, %4  \n\t" */
/*         "add %4, %0     \n\t" */
/*         "popcnt %5, %5  \n\t" */
/*         "add %5, %1     \n\t" */
/*         "popcnt %6, %6  \n\t" */
/*         "add %6, %2     \n\t" */
/*         "popcnt %7, %7  \n\t" */
/*         "add %7, %3     \n\t" // +r means input/output, r means intput */
/*         : "+r" (cnt[0]), "+r" (cnt[1]), "+r" (cnt[2]), "+r" (cnt[3]) */
/*         : "r"  (buf[i]), "r"  (buf[i+1]), "r"  (buf[i+2]), "r"  (buf[i+3])); */
/*   } */
/*   return cnt[0] + cnt[1] + cnt[2] + cnt[3]; */
/* } */
 
static inline int m4ri_bitcount(word w){
  return __builtin_popcountll(w);  
}
 
#else /* no __POPCNT__ */
 
#define MASK(c)    (((uint64_t)(-1)) / (__M4RI_TWOPOW(__M4RI_TWOPOW(c)) + 1))
#define COUNT(x,c) ((x) & MASK(c)) + (((x) >> (__M4RI_TWOPOW(c))) & MASK(c))
 
static inline int m4ri_bitcount(word w)  {
   uint64_t n = __M4RI_CONVERT_TO_UINT64_T(w);
   n = COUNT(n, 0);
   n = COUNT(n, 1);
   n = COUNT(n, 2);
   n = COUNT(n, 3);
   n = COUNT(n, 4);
   n = COUNT(n, 5);
   return (int)n;
}
 
static inline int std_bitcount ( uint64_t x) {
  uint64_t c1 = UINT64_C (0x5555555555555555 );
  uint64_t c2 = UINT64_C (0x3333333333333333 );
  uint64_t c4 = UINT64_C (0x0F0F0F0F0F0F0F0F );
  x -= (x >> 1) & c1;
  x = (( x >> 2) & c2) + (x & c2);
  x = ( x + (x >> 4) ) & c4;
  x *= UINT64_C (0x0101010101010101 );
  return (int) (x >> 56);
}
 
#endif /* __POPCNT__ */
 
 
typedef struct{    /*  */
  int rows ;        /* number of rows */
  int cols ;        /* number of columns */
  int nz ;      /* # of entries in triplet matrix */
  int nzmax ;       /* # allocated size */
  int *p ;      /* row pointers (size rows+1) OR row indices */
  int *i ;      /* col indices, size nzmax */
} csr_t ;
 
 
typedef struct { int a; int b; } int_pair;
 
 
#if defined(__cplusplus) && !defined (_MSC_VER)
extern "C" {
#endif
 
  size_t mzd_weight(const mzd_t *A);
 
  size_t mzd_weight_naive(const mzd_t *A);
 
  size_t mzd_weight_row(const mzd_t *A, rci_t i);
 
  static inline int nextelement(const word * const set1, const int m, const int pos){
    word setwd;
    int w;
    if (pos < 0){
      w = 0;
      setwd = set1[0];
    }
    else{
    w = SETWD(pos);
    if (w >= m) return -1;
    setwd = set1[w] & (m4ri_ffff<< SETBT(pos));
      }
 
    for (;;){
      if (setwd != 0)
    return TIMESWORDSIZE(w) + FIRSTBIT(setwd);
      if (++w == m)
    return -1;
      setwd = set1[w];
    }
  }
 
 
  rci_t mzd_gauss_naive(mzd_t *M, mzp_t *q, int full);
 
  int csr_max_row_wght(const csr_t * const p);
  
  csr_t * csr_transpose(csr_t *dst, const csr_t * const p);
 
  
  mzd_t *mzd_from_csr(mzd_t *dst, const csr_t *p);
 
  mzd_t *mzd_generator_from_csr(mzd_t *G, const csr_t * const H);
 
  mzd_t * csr_mzd_mul(mzd_t *C, const csr_t *S, const mzd_t *B, int clear);
 
  mzd_t * syndrome_vector(mzd_t *syndrome, mzd_t *row, csr_t *spaQ, int clear);
 
 
  size_t product_weight_csr_mzd(const csr_t *A, const mzd_t *B, int transpose);
 
  int rand_uniform(const int max);
 
  mzp_t * mzp_rand_len(mzp_t *q, rci_t length);
 
  static inline mzp_t * mzp_rand(mzp_t *q){
    if (q==NULL){
      printf("mzp_rand: permutation must be initialized!");
      exit(-1);
    }
    return mzp_rand_len(q,q->length);
  }
 
  void mzp_out(mzp_t const *p);
 
  mzp_t *perm_p(mzp_t *q, const mzp_t *p,rci_t start);
 
  mzp_t *perm_p_trans(mzp_t *q, const mzp_t *p,const rci_t start);
 
  csr_t *csr_free(csr_t *p);
 
  csr_t *csr_init(csr_t *mat, int rows, int cols, int nzmax);
 
  void csr_compress(csr_t *mat);
 
  csr_t * csr_from_pairs(csr_t *mat, const int nz, int_pair * const prs, const int nrows, const int ncols);
 
  void csr_out(const csr_t *mat);
 
  void csr_print(const csr_t * const smat, const char str[]);
 
  csr_t *csr_mm_read(char *fin, csr_t *mat, int transpose);
 
  csr_t *csr_apply_perm(csr_t *dst, const csr_t * const src, const mzp_t * const perm);
 
 
static inline void mzd_flip_bit(mzd_t * const M, rci_t const row, rci_t const col ) {
  word * const rawrow = mzd_row(M,row);
  __M4RI_FLIP_BIT(rawrow[col/m4ri_radix], col%m4ri_radix);
}
 
static inline int gauss_one(mzd_t *M, const int idx, const int begrow){
  rci_t startrow = begrow;
  rci_t pivots = 0;
  const rci_t i = idx;
  //  for (rci_t i = startcol; i < endcol ; ++i) {
  for(rci_t j = startrow ; j < M->nrows; ++j) {
    if (mzd_read_bit(M, j, i)) {
      mzd_row_swap(M, startrow, j);
      ++pivots;
      for(rci_t ii = 0 ;  ii < M->nrows; ++ii) {
        if (ii != startrow) {
          if (mzd_read_bit(M, ii, i)) {
            mzd_row_add_offset(M, ii, startrow,0);
          }
        }
      }
      startrow = startrow + 1;
      break;
    }
  }
  //  }
  return pivots; 
  // if one, need to update the current pivot list
}
 
static inline int sparse_syndrome_non_zero(const csr_t * const H, const int cnt, const int ee[]){
  for(int ir=0; ir < H->rows; ir++){
    int nz=0;
    for(int iL = H->p[ir], iE = 0; iL < H->p[ir+1]; iL++){
      int ic = H->i[iL];
      while((iE < cnt) && (ee[iE] < ic))
    iE++;
      if(iE >= cnt)
    break;
      if(ee[iE]==ic)
    nz ^= 1;      
    }
    if (nz)
      return 1;
  }
  return 0;
}
 
  int csr_csr_mul_non_zero(const csr_t * const A, const csr_t * const B);
  
  int syndrome_bit_count(const mzd_t * const row, const csr_t * const spaQ);
 
  int do_reduce(mzd_t *row, const mzd_t *matP0, const rci_t rankP0);
 
  void make_err(mzd_t *row, double p);
 
  int do_dist_clus(const csr_t * const P, const mzd_t * const G, int debug, int wmax, int start, const int rankG);
 
  csr_t * Lx_for_CSS_code(const csr_t * const Hx, const csr_t *const Hz);
  
#if defined(__cplusplus) && !defined (_MSC_VER)
}
#endif
  
#endif /* UTIL_M4RI_H */