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
}

#ifdef __M4RI_MAX_MZD_BLOCKSIZE
#define MZD_FIRST_ROW(M) mzd_row((M), 0)
#else
#define MZD_FIRST_ROW(M) ((M)->data)
#endif
 

#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;

typedef struct {
  int rows;   
  int cols;   
  int kw;     
  word *mask; 
} colmask_t;

typedef struct {
  int m;       
  int *rmap;   
  int *rlist;  
  word *bits;  
  size_t cap;  
} cols_ws_t;
 
#define SORT_NAME rci
#define SORT_TYPE rci_t
#define SORT_CMP(x, y) ((x) - (y))
#define SORT_DEF static inline
#include "sort.h"
 
 
#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);
}
 
#ifdef __M4RI_MAX_MZD_BLOCKSIZE
int gauss_one_rows_blocks(mzd_t *M, const int idx, const int begrow, const rci_t nrows_eff);
#endif

static inline int gauss_one_rows(mzd_t *M, const int idx, const int begrow,
                                 const rci_t nrows_eff) {
  if (begrow >= nrows_eff) return 0;
#ifdef __M4RI_MAX_MZD_BLOCKSIZE
  /* The rows i are at mzd_row(M, 0) + i * rowstride, see MZD_FIRST_ROW(); with M4RI versions before 2024, unless the
   * matrix spans several blocks */
  if (__builtin_expect((uintptr_t)mzd_row(M, nrows_eff - 1) - (uintptr_t)mzd_row(M, 0) !=
                       (uintptr_t)(nrows_eff - 1) * (uintptr_t)M->rowstride * sizeof(word), 0)) {
    return gauss_one_rows_blocks(M, idx, begrow, nrows_eff);
  }
#endif
  const wi_t word_idx = idx >> 6;
  const word bit_mask = (word)1 << (idx & 63);
  const wi_t width = M->width;
  const wi_t rowstride = M->rowstride;
  word *row_j = mzd_row(M, begrow);
  for (rci_t j = begrow; j < nrows_eff; ++j, row_j += rowstride) {
    if (row_j[word_idx] & bit_mask) {
      if (j != begrow) {
        mzd_row_swap(M, begrow, j);
      }
      const word * const __restrict piv_row = mzd_row(M, begrow);
      word * __restrict dst_row = MZD_FIRST_ROW(M);
      for (rci_t ii = 0; ii < begrow; ++ii, dst_row += rowstride) {
        if (dst_row[word_idx] & bit_mask) {
          for (wi_t w = 0; w < width; ++w) {
            dst_row[w] ^= piv_row[w];
          }
        }
      }
      dst_row += rowstride;
      for (rci_t ii = begrow + 1; ii < nrows_eff; ++ii, dst_row += rowstride) {
        if (dst_row[word_idx] & bit_mask) {
          for (wi_t w = 0; w < width; ++w) {
            dst_row[w] ^= piv_row[w];
          }
        }
      }
      return 1;
    }
  }
  return 0;
}

static inline int gauss_one(mzd_t *M, const int idx, const int begrow) {
  return gauss_one_rows(M, idx, begrow, M->nrows);
}

 
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;
}

static inline int colmask_syndrome_non_zero(const colmask_t * const cm, const int cnt, const int ee[]) {
  const int kw = cm->kw;
  const word * const mask = cm->mask;
  if (kw == 1) {
    word s = 0;
    for (int i = 0; i < cnt; i++)
      s ^= mask[ee[i]];
    return s != 0;
  }
  for (int w = 0; w < kw; w++) {
    word s = 0;
    for (int i = 0; i < cnt; i++)
      s ^= mask[(size_t)ee[i] * kw + w];
    if (s)
      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);

  mzd_t * mzd_nullspace(const csr_t * const H);

  colmask_t * colmask_from_csr(const csr_t * const A);

  colmask_t * colmask_free(colmask_t *cm);

  cols_ws_t * cols_ws_init(const int m);

  cols_ws_t * cols_ws_free(cols_ws_t *ws);

  int csr_cols_rank(const csr_t * const HT, const int cnt, const int cols[], cols_ws_t * const ws);

  int csr_cols_dependency(const csr_t * const HT, const int cnt, const int cols[], cols_ws_t * const ws,
                          int sub[]);

  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);

  void refresh_nullspace_basis(mzd_t *N, const int *cw_arr, int cw_wt, uint64_t *rng_state);

  int csr_all_ones_in_rowspan(const csr_t * const H);
  
#if defined(__cplusplus) && !defined (_MSC_VER)
}
#endif
  
#endif /* UTIL_M4RI_H */