dist_rw.c

Go to the documentation of this file.

 
// #include <m4ri/config.h>
#include <inttypes.h>
#include <strings.h>
#include <stdlib.h>
#include <time.h>
#include <m4ri/m4ri.h>
 
#include "mmio.h"
#include "util_m4ri.h"
#include "util_io.h"
#include "dist_m4ri.h"
 
 
int do_RW_dist(params_t * const p){
  const csr_t * const spaH0 = p->spaH;
  const csr_t * const spaL0 = p->spaL;
  const int steps = p->steps;
  const int wmin = p->wmin;
  const int wmax = p->wmax;
  const int classical = p->classical;
  const int debug = p->debug;
  const int nvar = spaH0->cols;
  if(((!classical)&&(spaL0==NULL)) ||
     ((classical)&&(spaL0!=NULL))){
    printf("L0 %s NULL classical=%d\n",spaL0==NULL ? "=" : "!=", classical);       
    ERROR("L0 should be non-NULL only for classical code!\n");
  }
 
  int minW = nvar + 1;
  if (p->dmax > 0 && p->dmax < minW) {
    minW = p->dmax;
  }
  if (p->min_w != INT_MAX && p->min_w < minW) {
    minW = p->min_w;
  }
 
  if(debug&2)
    printf("# running do_RW_dist() with steps=%d wmin=%d wmax=%d classical=%d nvar=%d\n",
       steps, wmin, wmax, classical, nvar);
  
  mzd_t * mH = mzd_from_csr(NULL, spaH0);
  mzd_t *mHT = NULL;
  rci_t *ee = malloc(nvar*sizeof(rci_t)); 
  
  if((!mH) || (!ee))
    ERROR("memory allocation failed!\n");
 
  mzp_t * perm=mzp_init(nvar); 
  mzp_t * pivs=mzp_init(nvar); 
  mzp_t * pivs_srtd=mzp_init(nvar); 
  mzp_t * skip_pivs=mzp_init(nvar); 
  if((!pivs) || (!perm))
    ERROR("memory allocation failed!\n");
 
  for (int ii=0; ii< steps; ii++){
    pivs=mzp_rand(pivs); 
    mzp_set_ui(perm,1);
    perm=perm_p_trans(perm,pivs,0); 
    int rank=0;
    for(int i=0; i< nvar; i++){
      int col=perm->values[i];
      int ret=gauss_one(mH, col, rank);
      if(ret)
        pivs->values[rank++]=col;
    }
 
    pivs_srtd = mzp_copy(pivs_srtd,pivs);
    qsort(pivs_srtd->values, rank, sizeof(pivs->values[0]), cmp_rci_t);
    int end=-1, num=0;
    for(int i=0; i < rank; i++){
      int beg = end + 1;
      end = pivs_srtd->values[i];
      for(int j = beg; j < end; j++)
    skip_pivs->values[num++] = j;
    }
    for(int j = end + 1 ; j < nvar; j++)
      skip_pivs->values[num++] = j;
    
#ifndef NDEBUG
    if (num + rank != nvar)
      ERROR("mismatch: rank=%d and num=%d do not add to nvar=%d\n",rank,num,nvar);
#endif
    skip_pivs->length = num;
 
#ifndef NEW
# define NEW 1
#endif 
#if (NEW==1)
    mHT = mzd_transpose(mHT,mH);
#endif     
    int k = nvar - rank;
    for (int ir=0; ir< k; ir++){ 
      int cnt=0; 
      const int col = ee[cnt++] = skip_pivs->values[ir];
      int limit = nvar + 1;
      int cur_d = (minW <= nvar) ? minW : 0;
      if (cur_d > 0) {
        if ((p->outC || p->maxC || p->dW > 0) && p->dW >= 0) {
          limit = minint(limit, cur_d + p->dW + 1);
        } else {
          limit = minint(limit, cur_d);
        }
      }
#if (NEW==0) 
      for(int ix=0; ix<rank; ix++){
        if(mzd_read_bit(mH,ix,col))
          ee[cnt++] = pivs->values[ix];
    if (cnt >= limit) 
      break;
      }
#elif (NEW==2) 
      rci_t ic=col, ix=0;
      while (ix < rank){
    int res = mzd_find_pivot(mH, ix, col, &ix, &ic);
    if((res)&&(ic==col)){
      ee[cnt++] = pivs->values[ix++];
      //      printf("cnt=%d j=%d\n",cnt,ix); 
      if (cnt >= limit) 
        break;
    }
    else
      break;
      }
#else 
      word * rawrow = mzd_row(mHT,col);  
      rci_t j=-1;
      while(cnt < limit){
    j=nextelement(rawrow,mHT->width,j);
    if(j==-1) // empty line after simplification
      break; 
    ee[cnt++] = pivs->values[j++];
      }
#endif /* NEW */              
      if (cnt < limit){
    qsort(ee, cnt, sizeof(rci_t), cmp_rci_t);
#ifndef NDEBUG
    if(sparse_syndrome_non_zero(spaH0, cnt, ee)){
      printf("# cw of weight %d: [",cnt);
      for(int i=0; i<cnt;i++)
        printf("%d%s",ee[i],i+1==cnt?" ":"]\n");
      ERROR("this should not happen: cw not orthogonal to H");
    }
#endif /* NDEBUG */
      
    int nz;
    if (classical)
      nz=1; 
    else
      nz = sparse_syndrome_non_zero(spaL0, cnt, ee);    
    if(nz){ 
          p->codewords = codeword_add_maybe(p, ee, cnt);
      if(debug&16){
        printf("# step=%d row=%d minW=%d found cw of W=%d: [",ii,ir,minW,cnt);
        const int max = ((cnt<25) || (debug&2048)) ?  cnt : 25 ;
        for(int i=0; i< max; i++)
          printf("%d%s", ee[i], i+1!=max?" ": (cnt==max ? "]\n" : "...]\n"));
      }
          if (cnt < minW) {
            minW = cnt;
          }
          if (p->maxC && p->num_cws >= p->maxC) {
            goto alldone;
          }
      if (minW <= wmin){ 
        minW = - minW;   
        goto alldone; 
      }
    }
      }      
    } 
    if(debug&8){
      if(ii%1000==999)
    printf("# round=%d of %d minW=%d\n", ii+1, steps, minW);
    }
    
  }
 alldone: 
  if(skip_pivs)
    mzp_free(skip_pivs);
  if(pivs_srtd)
    mzp_free(pivs_srtd);
  mzp_free(perm);
  mzp_free(pivs);
  free(ee);
  if(mHT)
    mzd_free(mHT);
  mzd_free(mH);
  
  if (minW < 0) {
    return minW;
  }
  if (p->min_w == INT_MAX) {
    return 0;
  }
  if (wmax > 0 && p->min_w > wmax) {
    return 0;
  }
  return p->min_w;
}
 
 
#ifdef STANDALONE
 
int do_CC_dist(params_t * const p);
 
 
int main(int argc, char **argv){
  params_t * const p = &prm;
 
  var_init(argc,argv,p);
 
  if (p->finC) {
    nzlist_read(p->finC, p);
  }
 
  //  const int n=p->nvar;
 
  if (prm.method & 1){ /* RW method */
    
    prm.dist_max=do_RW_dist(p);
 
    if (prm.debug&1){
      if (prm.dist_max != 0)
        printf("### RW upper bound on the distance: %d\n",prm.dist_max);
      else
        printf("### RW: no upper bound found up to wmax = %d\n", prm.wmax);
      if(prm.dist_max <0)
        printf("### negative distance due to wmin=%d set (early termination)\n",prm.wmin);
      else if (prm.dist_max ==0)
        printf("### no codewords of weight <= %d found\n",prm.wmax);
    }
    if (prm.dist_max < 0) {
      if(prm.debug) {
        if (prm.method == 1) 
          printf("RW algorithm upper bound for the distance d=%d\n", prm.dist_max);
        else
          printf("RW algorithm upper bound for the distance d=%d (early termination due to wmin=%d)\n", prm.dist_max, prm.wmin);
      }
      printf("%d\n",prm.dist_max);
      prm.dist_min = prm.dist_max;
      goto end_all;
    }
    if (prm.method==1){ 
      if(prm.debug) {
        if (prm.dist_max > 0)
          printf("RW algorithm upper bound for the distance d=%d\n", prm.dist_max);
        else
          printf("RW algorithm: no upper bound found up to wmax = %d\n", prm.wmax);
      }
      printf("%d\n",prm.dist_max);
    }
    else{
      if (prm.wmax==0) {
        if (prm.dist_max == 0) {
          ERROR("RW found no codewords, cannot run CC with wmax=0. Please specify wmax > 0.\n");
        }
        prm.wmax=abs(prm.dist_max)-1;
      }
      else if (prm.dist_max != 0) {
        prm.wmax=minint(prm.wmax, abs(prm.dist_max)-1);   
      }
      if (prm.wmax == 0) {
        prm.dist_min = 1;
        prm.dist_max = 1;
        if (prm.debug & 1) {
          printf("success (distance is 1) d=1\n");
        }
        printf("1\n");
        goto end_all;
      }
    }
  }
  
  if (prm.method & 2){ /* cluster method */
    int dmin=do_CC_dist(p);
 
    if (dmin>0){ 
      if (prm.debug&1)
    printf("### Cluster (actual min-weight codeword found): d=%d\n",dmin);
      printf("%d\n",dmin);
      prm.dist_min = dmin; /* actual distance found */
      prm.dist_max = dmin;
      goto end_all;
    }
    else if (dmin<0){
      if (prm.debug&1)
    printf("### Cluster dmin=%d  (no codewords of weight up to %d)\n",dmin,-dmin);
      if (-dmin==abs(prm.dist_max)-1){
    prm.dist_min=abs(prm.dist_max); /* OK */
        if (prm.debug&1)
          printf("success (two distance bounds coincide) d=%d\n",prm.dist_min);
        printf("%d\n",prm.dist_min);       
        goto end_all;
      }
      else{
    prm.dist_min=-dmin;
        if(prm.debug){
          if (prm.dist_max>prm.dist_min)
            printf("# distance in the interval (inclusive) %d to %d\n", prm.dist_min,prm.dist_max);
          else
            printf("# cluster algorithm failed to find a codeword up to wmax=%d\n",-dmin);
        }
        printf("%d\n",dmin);
      }
    }
    else
      ERROR("unexpected dmin=0\n");
  }
 end_all:
    if (p->outC) {
      char comment[256];
      sprintf(comment, "generated by dist_m4ri");
      nzlist_write(p->outC, comment, p);
    }
    if (p->debug & 32) {
      cw_vec_t *cw;
      for(cw = p->codewords; cw != NULL; cw = (cw_vec_t *)(cw->hh.next)){
        printf("# cw: [ ");
        for(int i=0; i<cw->weight; i++) printf("%d ", 1 + cw->arr[i]);
        printf("] cnt=%d\n", cw->cnt);
      }
    }
    var_kill(p);
    
    return 0;
  }
 
#endif /* STANDALONE */