diff --git a/src/dist_cc.c b/src/dist_cc.c index 2fe9864..c42c8af 100644 --- a/src/dist_cc.c +++ b/src/dist_cc.c @@ -168,10 +168,10 @@ static inline void one_ordered_pos_del(one_vec_t * const err, _maybe_unused cons #ifndef NDEBUG if ((pos<0) || (pos >= err->wei) || (err->wei == 0) || (err->vec[pos] != val)) ERROR("this should not happen!"); -#endif +#endif + err->wei --; for(int i=pos; i < err->wei; i++) err->vec[i] = err->vec[i+1]; - err->wei --; } /** @brief recursively construct codewords diff --git a/src/dist_m4ri.c b/src/dist_m4ri.c index 2e41e47..798c0c1 100644 --- a/src/dist_m4ri.c +++ b/src/dist_m4ri.c @@ -27,7 +27,7 @@ #include "dist_m4ri.h" /** @brief prepare an ordered pivot-skip list of length `n-rank` */ -mzp_t * do_skip_pivs(const int rank, const mzp_t * const pivs){ +mzp_t * do_skip_pivs(mzp_t * ans, const int rank, const mzp_t * const pivs){ const rci_t n=pivs->length; int *arr=calloc(rank,sizeof(int)); if(!arr) @@ -35,7 +35,13 @@ mzp_t * do_skip_pivs(const int rank, const mzp_t * const pivs){ for(int i=0; i< rank; i++) arr[i] = pivs->values[i]; qsort(arr, rank, sizeof(arr[0]), cmp_rci_t); - mzp_t *ans=mzp_init(n-rank); + if((ans) && (ans->length != n-rank)){ + mzp_free(ans); + ans=NULL; + } + if (!ans) + ans=mzp_init(n-rank); + rci_t j=0, j1=0; /** position to insert the next result */ for(rci_t i1=0; i1 < rank; i1++, j++){ const rci_t piv=arr[i1]; @@ -87,6 +93,8 @@ int do_RW_dist(const csr_t * const spaH0, const csr_t * const spaL0, steps, wmin, classical, nvar); mzd_t * mH = mzd_from_csr(NULL, spaH0); + mzd_t *mHT = NULL; + mzp_t * skip_pivs = NULL; /** actual `vector` in sparse form */ rci_t *ee = malloc(nvar*sizeof(rci_t)); @@ -94,7 +102,6 @@ int do_RW_dist(const csr_t * const spaH0, const csr_t * const spaL0, ERROR("memory allocation failed!\n"); /** 1. Construct random column permutation P */ - mzp_t * perm=mzp_init(nvar); /** identity column permutation */ mzp_t * pivs=mzp_init(nvar); /** list of pivot columns */ if((!pivs) || (!perm)) @@ -115,10 +122,14 @@ int do_RW_dist(const csr_t * const spaH0, const csr_t * const spaL0, } /** construct skip-pivot permutation */ - mzp_t * skip_pivs = do_skip_pivs(rank, pivs); - - /** TODO: would it be faster to transpose `mH` first? */ - + skip_pivs = do_skip_pivs(skip_pivs, rank, pivs); +#ifndef NEW +# define NEW 1 +#endif +#if (NEW==1) + /** it is a bit faster to transpose `mH` first. */ + mHT = mzd_transpose(mHT,mH); +#endif /** calculate sparse version of each vector (list of positions) * `p` `p``p` # pivot columns marked with `p` * [1 a1 b1 ] -> [a1 1 a2 a3 0 ] @@ -129,13 +140,43 @@ int do_RW_dist(const csr_t * const spaH0, const csr_t * const spaL0, for (int ir=0; ir< k; ir++){ /** each row in the dual matrix */ int cnt=0; /** how many non-zero elements */ const int col = ee[cnt++] = skip_pivs->values[ir]; +#if (NEW==0) /** older version going over columns of `H` */ for(int ix=0; ixvalues[ix]; if (cnt >= minW) /** `cw` of no interest */ break; } - if(cnt < minW){ +#elif (NEW==2) /** + function `mzd_find_pivot()` walks over columns + one-by-one to find a non-zero bit. Returns `1` if + a non-zero bit was found. + WARNING: this is the slowest option!!! + */ + 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 >= minW) /** `cw` of no interest */ + break; + } + else + break; + } +#else /** NEW, use transposed `H` -- the `fastest` version of the code*/ + word * rawrow = mHT->rows[col]; + rci_t j=-1; + while(cnt < minW){/** `cw` of no interest */ + j=nextelement(rawrow,mHT->width,j); + if(j==-1) // empty line after simplification + break; + // printf("cnt=%d j=%d\n",cnt,j); + ee[cnt++] = pivs->values[j++]; + } +#endif /* NEW */ + if (cnt < minW){ /** sort the column indices */ qsort(ee, cnt, sizeof(rci_t), cmp_rci_t); #ifndef NDEBUG @@ -158,33 +199,35 @@ int do_RW_dist(const csr_t * const spaH0, const csr_t * const spaL0, /** TODO: try local search to `lerr` (if 2 or larger) */ /** at this point we have `cnt` codeword indices in `ee` */ if(debug&16){ - printf("# minW=%d found cw of weight %d: [",minW,cnt); - int max = ((cnt<50) || (debug&2048)) ? cnt : 50 ; - for(int i=0; i< max;i++) - printf("%d%s",ee[i], i+1!=max?" ": (cnt==max ? "]\n" : "...]\n")); + printf("# step=%d row=%d minW=%d found cw of W=%d: [",ii,ir,minW,cnt); + 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")); } minW=cnt; if (minW <= wmin){ /** early termination condition */ minW = - minW; /** this distance value is of little interest; */ } } - } + } } /** end of the dual matrix rows loop */ if(debug&8){ if(ii%1000==999) printf("# round=%d of %d minW=%d\n", ii+1, steps, minW); } - mzp_free(skip_pivs); }/** end of `steps` random window */ //alldone: /** early termination label */ /** clean up */ + if(skip_pivs) + mzp_free(skip_pivs); mzp_free(perm); mzp_free(pivs); free(ee); - + if(mHT) + mzd_free(mHT); mzd_free(mH); return minW; diff --git a/src/util_io.c b/src/util_io.c index 13c3733..7023891 100644 --- a/src/util_io.c +++ b/src/util_io.c @@ -1,5 +1,5 @@ +#include #include "util_io.h" -// #include "util.h" params_t prm={ .debug=3, @@ -137,11 +137,10 @@ void var_init(int argc, char **argv, params_t * const p){ if (p->debug&4) printf("# read %s, seed=%d\n",argv[i],p->seed); if (p->seed<=0){ - p->seed=time(NULL)+1000*p->seed; + p->seed=time(NULL)+1000*p->seed+10*getpid(); if(p->debug&4) printf("# initializing rng from time(NULL), seed=%d\n",p->seed); } - srand(p->seed); } else{ /* unrecognized option */ printf("# unrecognized parameter \"%s\" at position %d\n",argv[i],i); @@ -230,6 +229,8 @@ void var_init(int argc, char **argv, params_t * const p){ p->spaG=NULL; } + srand(p->seed); + rci_t n = (p->spaH)-> cols; if ((!p->classical) && ((p->spaG) && (n != (p->spaG) -> cols))) ERROR("Column count mismatch in H and G matrices: %d != %d", diff --git a/src/util_m4ri.h b/src/util_m4ri.h index 99bfd53..5ef8c23 100644 --- a/src/util_m4ri.h +++ b/src/util_m4ri.h @@ -145,21 +145,20 @@ extern "C" { static inline int nextelement(const word * const set1, const int m, const int pos){ word setwd; int w; -#if 1 if (pos < 0){ w = 0; setwd = set1[0]; } - else -#endif - // { - w = SETWD(pos); - setwd = set1[w] & (m4ri_ffff<< SETBT(pos)); - // } + else{ + w = SETWD(pos); + setwd = set1[w] & (m4ri_ffff<< SETBT(pos)); + } for (;;){ - if (setwd != 0) return TIMESWORDSIZE(w) + FIRSTBIT(setwd); - if (++w == m) return -1; + if (setwd != 0) + return TIMESWORDSIZE(w) + FIRSTBIT(setwd); + if (++w == m) + return -1; setwd = set1[w]; } }