source: git/Singular/LIB/decodegb.lib @ 288382

///////////////////////////////////////////////////////////////////////////////
version="$Id: decodegb.lib,v 1.1 2008-08-13 15:44:20 bulygin Exp $";
category="Coding theory";
info="
LIBRARY:   decodedistGB.lib     Generating and solving systems of polynomial equations for decoding and finding the minimum distance of linear codes
AUTHORS:   Stanislav Bulygin,   bulygin@mathematik.uni-kl.de                     

OVERVIEW:
 In this library we generate several systems used for decoding cyclic codes. And finding the minimum distance
 Namely, we work with the Cooper's philosophy and generalized Newton identities.
 The original method of quadratic equations is worked out here as well. 
 We also (for comparison) enable to work with the system of Fitzgerald-Lax.
 We provide also some auxiliary functions for further manipulations and decoding.
 For an overview of the methods mentioned above, cf. Stanislav Bulygin, Ruud Pellikaan: 'Decoding and finding the minimum distance with 
 Groebner bases: history and new insights', in 'Selected Topics in Information and Coding Theory', World Scientific (2008) (to appear) (*).

MAIN PROCEDURES: 
 sysCRHT(n,defset,e,q,m,#); generates the CRHT-ideal that follows Cooper's philosophy, Sala's extentions are available
 sysCRHTMindistBinary(n,defset,w); generates the ideal from Cooper's philosophy to find the minimum distance in the binary case
 sysNewton(n,defset,t,q,m,#); generates the ideal with the Generalized Newton identities
 sysBin(v,Q,n,#);  generates Bin system as in the work of Augot et.al, cf. [*] for a reference
 encode(x,g);  encodes given message with a given generator matrix
 syndrome(h, c);  computes a syndroem w.r.t. a given check matrix
 sysQE(check,y,t,fieldeq,formal); generates the system of quadratic equations as in the method of Pellikaan and Bulygin
 error(y,pos,val);  inserts errors in a word
 errorRand(y,num,e);  inserts random errors in a word
 randomCheck(m,n,e);  generates a random check matrix
 genMDSMat(n,a);  generates an MDS (actually an RS) matrix
 mindist(check);   computes the minimum distance of the code via solving systems of quadratic equations
 decode(rec);  decoding of a word using the systems of quadratic equations 
 solveForRandom(redun,ncodes,ntrials,#);  a procedure for manipulation with random codes
 solveForCode(check,ntrials,#); a procedure for manipulation with a given codes
 vanishId(points); computes the vanishing ideal for the given set of points. The algorithm of Farr and Gao is implemented
 sysFL(check,y,t,e,s); generates the Fitzgerald-Lax system
 FLSolveForRandom(n,redun,p,e,t,ncodes,ntrials,minpol); a procedure for manipulation with random codes via Fitzgerald-Lax
 

KEYWORDS:  Cyclic code; Linear code; Decoding; 
           Minimum distance; Groebner bases
";

LIB "linalg.lib";
LIB "brnoeth.lib";

///////////////////////////////////////////////////////////////////////////////

static proc lis (int n)
{
 list result;
 if (n<=0) {print("ERRORlis");}
 for (int i=1; i<=n; i++)
 {
  result=result+list(i);
 }
 return(result);
}

static proc combinations (int m, int n)
{
 list result;
 if (m>n) {print("ERRORcombinations");}
 if (m==n) {result[size(result)+1]=lis(m);return(result);}
 if (m==0) {result[size(result)+1]=list();return(result);}
 list temp=combinations(m-1,n-1);
 for (int i=1; i<=size(temp); i++)
 {
  temp[i]=temp[i]+list(n);
 }
 result=combinations(m,n-1)+temp;
 return(result);
}

static proc combinsert (list temp, int i)
{
 list result;
 list tmp;
 int j,k;
 for (j=1; j<=size(temp); j++)
 {
  result[j]=tmp;
 }
 for (j=1; j<=size(temp); j++)
 {
  for (k=1; k<=size(temp[j]); k++)
  {
   if (temp[j][k]<i) {result[j][k]=temp[j][k];}
   else {result[j][k]=temp[j][k]+1;}
  }
 }
 return(result);   
}

static proc p_poly(int n, int a, int b)
{
  poly f;
  for (int i=0; i<=n-1; i++)
  {
    f=f+Z(a)^i*Z(b)^(n-1-i);
  }
  return(f);
}

proc sysCRHT (int n, list defset, int e, int q, int m, int #)
"USAGE:   sysCRHT(n,defset,e,q,m,#);  n length of the cyclic code, defset is a list representing the defining set,
          e the error-correcting capacity, m degree extension of the splitting field, if #>0 additional equations
          representing the fact that every two error positions are either different or at least one of them is zero
RETURN:   a ring to work with the CRHT-ideal (with Sala's additions), the ideal itself is exported with the name 'crht'
EXAMPLE:  example sysCRHT; shows an example
"
{
  int r=size(defset);  
  ring @crht=(q,a),(Y(e..1),Z(1..e),X(r..1)),lp;
  ideal crht;
  int i,j;
  poly sum;
  
  // check equations
  for (i=1; i<=r; i++)
  {
    sum=0;
    for (j=1; j<=e; j++)
    {
      sum=sum+Y(j)*Z(j)^defset[i];
    }
    crht[i]=sum-X(i);
  }
  
  // restrictions on syndromes
  for (i=1; i<=r; i++)
  {
    crht=crht,X(i)^(q^m)-X(i);
  }
  
  // n-th roots of unity
  for (i=1; i<=e; i++)
  {
    crht=crht,Z(i)^(n+1)-Z(i);
  }
  
  for (i=1; i<=e; i++)
  {
    crht=crht,Y(i)^(q-1)-1;
  }  
  
  if (#)
  {
    for (i=1; i<=e; i++)
    {
      for (j=i+1; j<=e; j++)
      {
        crht=crht,Z(i)*Z(j)*p_poly(n,i,j);
      }
    }
  }  
  export crht;  
  return(@crht);    
} example
{
  "EXAMPLE:"; echo=2;
  // binary cyclic [15,7,5] code with defining set (1,3)
  
  list defset=1,3; // defining set
  
  int n=15; // length
  int e=2; // error-correcting capacity
  int q=2; // basefield size
  int m=4; // degree extension of the splitting field
  int sala=1; // indicator to add additional equations
    
  def A=sysCRHT(n,defset,e,q,m); 
  setring A;
  A; // shows the ring we are working in
  print(crht); // the CRHT-ideal 
  option(redSB);
  ideal red_crht=std(crht);
  // reduced Groebner basis
  print(red_crht);
  
  //============================
  A=sysCRHT(n,defset,e,q,m,sala); 
  setring A;  
  print(crht); // the CRHT-ideal with additional equations from Sala
  option(redSB);
  ideal red_crht=std(crht);
  // reduced Groebner basis
  print(red_crht);
  // general error-locator polynomial for this code
  red_crht[5];
}


proc sysCRHTMindistBinary (int n, list defset, int w)
"USAGE:  sysCRHTMindistBinary(n,defset,w);  n length of the cyclic code, defset is a list representing the defining set,
         w is a candidate for the minimum distance
RETURN:  a ring to work with the Sala's ideal for mindist, the ideal itself is exported with the name 'crht_md'
EXAMPLE: example sysCRHTMindistBinary; shows an example
"
{
  int r=size(defset);  
  ring @crht_md=2,Z(1..w),lp;
  ideal crht_md;
  int i,j;
  poly sum;
  
  // check equations
  for (i=1; i<=r; i++)
  {
    sum=0;
    for (j=1; j<=w; j++)
    {
      sum=sum+Z(j)^defset[i];
    }
    crht_md[i]=sum;
  }  
  
  
  // n-th roots of unity
  for (i=1; i<=w; i++)
  {
    crht_md=crht_md,Z(i)^n-1;
  }  
  
  
  for (i=1; i<=w; i++)
  {
    for (j=i+1; j<=w; j++)
    {
      crht_md=crht_md,Z(i)*Z(j)*p_poly(n,i,j);
    }
  }
    
  export crht_md;  
  return(@crht_md);    
} example
{
  "EXAMPLE:"; echo=2;
  // binary cyclic [15,7,5] code with defining set (1,3)
  
  list defset=1,3; // defining set
  
  int n=15; // length
  int d=5; // candidate for the minimum distance  
      
  def A=sysCRHTMindistBinary(n,defset,d); 
  setring A;
  A; // shows the ring we are working in
  print(crht_md); // the Sala's ideal for mindist
  option(redSB);
  ideal red_crht_md=std(crht_md);
  // reduced Groebner basis
  print(red_crht_md);  
}

static proc mod_ (int n, int m)
{
  if (n mod m==0) {return(m);}
  if (n mod m!=0) {return(n mod m);}
}

proc sysNewton (int n, list defset, int t, int q, int m, int #)
"USAGE:   sysNewton (n, defset, t, q, m, #);  n is length, defset is the defining set,
          t is the number of errors, q is basefield size, m is degree extension of the splitting field
          if triangular>0 it indicates that Newton identities in triangular form should be constructed
RETURN:   a ring to work with the generalized Newton identities (in triangular form if applicable), 
          the ideal itself is exported with the name 'newton'
EXAMPLE:  example sysNewton; shows an example
"
{  
 string s="ring @newton=("+string(q)+",a),(";
 int i,j;
 int flag;
 for (i=n; i>=1; i--)
 {
  for (j=1; j<=size(defset); j++)
  {
    flag=1;
    if (i==defset[j])
    {
      flag=0;
      break;      
    }
  }
  if (flag)
  {
    s=s+"S("+string(i)+"),";
  }
 }
 s=s+"sigma(1.."+string(t)+"),";
 for (i=size(defset); i>=2; i--)
 {
  s=s+"S("+string(defset[i])+"),";
 }
 s=s+"S("+string(defset[1])+")),lp;"; 
  
 execute(s);
 
 ideal newton; 
 poly sum;
 
 
 // generate generalized Newton identities
 if (#)
 {
  for (i=1; i<=t; i++)
  {
    sum=0;
    for (j=1; j<=i-1; j++)
    {
      sum=sum+sigma(j)*S(i-j);
    }
    newton=newton,S(i)+sum+number(i)*sigma(i);
  }
 } else
 { 
  for (i=1; i<=t; i++)
  {
    sum=0;
    for (j=1; j<=t; j++)
    {
      sum=sum+sigma(j)*S(mod_(i-j,n));
    }
    newton=newton,S(i)+sum;
  }
 }
 for (i=1; i<=n-t; i++)
 {
  sum=0;
  for (j=1; j<=t; j++)
  {
    sum=sum+sigma(j)*S(t+i-j);
  }
  newton=newton,S(t+i)+sum;
 } 
 
 // field equations on sigma's
 for (i=1; i<=t; i++)
 {
  newton=newton,sigma(i)^(q^m)-sigma(i);
 }
 
 // conjugacy relations
 for (i=1; i<=n; i++)
 {
  newton=newton,S(i)^q-S(mod_(q*i,n));
 }
 newton=simplify(newton,2);
 export newton;
 return(@newton); 
} example 
{
     "EXAMPLE:";  echo = 2;
     // Newton identities for a binary 3-error-correcting cyclic code of length 31 with defining set (1,5,7)
     
     int n=31; // length
     list defset=1,5,7; //defining set
     int t=3; // number of errors
     int q=2; // basefield size
     int m=5; // degree extension of the splitting field
     int triangular=1; // indicator of triangular form of Newton identities        
     
     def A=sysNewton(n,defset,t,q,m);
     setring A;
     A; // shows the ring we are working in
     print(newton); // generalized Newton identities
     
     //===============================
     A=sysNewton(n,defset,t,q,m,triangular);
     setring A;     
     print(newton); // generalized Newton identities in triangular form
}

static proc combinations_sum (int m, int n)
{
 list result;
 list comb=combinations(m-1,n+m-1);
 int i,j,flag,count;
 list interm=comb;
 for (i=1; i<=size(comb); i++)
 {
  interm[i][1]=comb[i][1]-1;
  for (j=2; j<=m-1; j++)
  {
   interm[i][j]=comb[i][j]-comb[i][j-1]-1;
  }
  interm[i][m]=n+m-comb[i][m-1]-1;
  flag=1;
  count=2;
  while ((flag)&&(count<=m))
  {
   if (interm[i][count] mod count != 0) {flag=0;}
   count++;
  }
  if (flag) 
  {
   for (j=2; j<=m; j++)
   {
    interm[i][j]=interm[i][j] div j;
   }
   result[size(result)+1]=interm[i];  
  }
 }
 return(result);
}

static proc exp_count (int n, int q)
{
 int flag=1;
 int result=0;
 while(flag)
 {
  if (n mod q != 0) {flag=0;}
   else {n=n div q; result++;}
 }
 return(result);
}


proc sysBin (int v, list Q, int n, int#)
"USAGE:    sysBin (v, Q, n, #);  v a number if errors, Q is a generating set of the code, n the length, # is additional parameter is 
           set to 1, then the generating set is enlarged by odd elements, which are 2^(some power)*(some elment in the gen.set) mod n
RETURN:    keeps a ring with the resulting system, which ideal is called 'bin'
EXAMPLE:   example sysBin; shows an example
"
{
 //ring r=2,(sigma(1..v),S(1..n)),(lp(v),dp(n));
 ring r=2,(S(1..n),sigma(1..v)),lp;
 list cyclot;
 ideal result;
 int i,j,k,s;
 list comb;
 poly sum_, mon;
 int count1, count2, upper, coef_, flag, gener;
 list Q_update;
 if (#==1)
 {
  for (i=1; i<=n; i++)
  {
   cyclot[i]=0;
  }
  for (i=1; i<=size(Q); i++)
  {
   flag=1;
   gener=Q[i];
   while(flag)
   {
    cyclot[gener]=1;
    gener=2*gener mod n;
    if (gener == Q[i]) {flag=0;}
   }
  }
  for (i=1; i<=n; i++)
  {
   if ((cyclot[i] == 1)&&(i mod 2 == 1)) {Q_update[size(Q_update)+1]=i;}
  }
 }
 else
 {
  Q_update=Q;
 }
   
 for (i=1; i<=size(Q_update); i++)
 {
  comb=combinations_sum(v,Q_update[i]);
  sum_=0;
  for (k=1; k<=size(comb); k++)
  {
   upper=0;
   for (j=1; j<=v; j++)
   {
    upper=upper+comb[k][j];
   }
   count1=0;
   for (j=2; j<=upper-1; j++) 
   {
    count1=count1+exp_count(j,2);
   }
   count1=count1+exp_count(Q_update[i],2);
   count2=0;
   for (j=1; j<=v; j++)
   {
    for (s=2; s<=comb[k][j]; s++)
    {
     count2=count2+exp_count(s,2);
    }
   }
   if (count1<count2) {print("ERRORsysBin");}
   if (count1>count2) {coef_=0;}
   if (count1 == count2) {coef_=1;}
   mon=1;
   for (j=1; j<=v; j++)
   {
    mon=mon*sigma(j)^(comb[k][j]);
   }
   sum_=sum_+coef_*mon;
  }
  result=result,S(Q_update[i])-sum_;  
 }
 ideal bin=simplify(result,2);
 export bin;
 return(r);
} example 
{
     "EXAMPLE:";  echo = 2;
     // [31,16,7] quadratic residue code
     list l=1,5,7,9,19,25;
     // we do not need even synromes here
     def A=sysBin(3,l,31);
     setring A;
     print(bin);
}

proc encode (matrix x, matrix g)
"USAGE:  encode (x, g);  x a row vector (message), and g a generator matrix
RETURN:  corresponding codeword
EXAMPLE: example encode; shows an example
"
{
 if (nrows(x)>1) {print("ERRORencode1!");}
 if (ncols(x)!=nrows(g)) {print("ERRORencode2!");}
 return(x*g);
} example 
{
     "EXAMPLE:";  echo = 2;
     ring r=2,x,dp;
     matrix x[1][4]=1,0,1,0;
     matrix g[4][7]=1,0,0,0,0,1,1,
                    0,1,0,0,1,0,1,
                    0,0,1,0,1,1,1,
                    0,0,0,1,1,1,0;
     //encode x with the generator matrix g
     print(encode(x,g)); 
}

proc syndrome (matrix h, matrix c)
"USAGE:  syndrome (h, c);  h a check matrix, c a row vector (codeword)
RETURN:  corresponding syndrome
EXAMPLE: example syndrome; shows an example
"
{
 if (nrows(c)>1) {print("ERRORsyndrome1!");}
 if (ncols(c)!=ncols(h)) {print("ERRORsyndrome2!");}
 return(h*transpose(c));  
} example 
{
     "EXAMPLE:";  echo = 2;
     ring r=2,x,dp;
     matrix x[1][4]=1,0,1,0;
     matrix g[4][7]=1,0,0,0,0,1,1,
                    0,1,0,0,1,0,1,
                    0,0,1,0,1,1,1,
                    0,0,0,1,1,1,0;
     //encode x with the generator matrix g
     matrix c=encode(x,g);
     // disturb
     c[1,3]=0;
     //compute syndrome
     //corresponding check matrix
     matrix check[3][7]=1,0,0,1,1,0,1,0,1,0,1,0,1,1,0,0,1,0,1,1,1;
     print(syndrome(check,c));
     c[1,3]=1;
     //now c is a codeword 
     print(syndrome(check,c));
}

static proc star(matrix m, int i, int j)
{
 matrix result[ncols(m)][1];
 for (int k=1; k<=ncols(m); k++)
 {
  result[k,1]=m[i,k]*m[j,k];
 }
 return(result);
}

proc sysQE(matrix check, matrix y, int t, int fieldeq, int formal)
"USAGE:   sysQE(check, y, t, fieldeq, formal); check is the check matrix of the code   
          y is a received word, t the number of errors to be corrected, 
          if fieldeq=1, then field equations are added; if formal=0, fields equations on (known) syndrome variables 
          are not added, in order to add them (note that the exponent should be as a number of elements in the INITIAL alphabet) one
          needs to set formal>0 for the exponent
RETURN:   the ring to work with together with the resulting system
EXAMPLE:  example sysQE; shows an example
"
{
 def br=basering; 
 list rl=ringlist(br);
 
 int red=nrows(check);
 int n=ncols(check);
 int q=rl[1][1];
 
 if (formal==0)
 { 
  ring work=(q,a),(V(1..t),U(1..n)),dp;
 } else
 {
  ring work=(q,a),(V(1..t),U(1..n),s(1..red)),(dp(t),lp(n),dp(red));
 }
 
 matrix check=imap(br,check);
 matrix y=imap(br,y);
 
 matrix h_full=genMDSMat(n,a);
 matrix h=submat(h_full,1..red,1..n);
 if (nrows(y)!=1) {print("ERROR1Pell");}
 if (ncols(y)!=n) {print("ERROR2Pell");}
  
 ideal result;
 
 list c;
 list a;
 list tmp,tmp2;
 int i,j,l,k;
 number sum,prod,sig;
        poly sum1,sum2,sum3;
 for (i=1; i<=n; i++)
 {
  c[i]=tmp;
 }
 
 int tim=rtimer;
 matrix transf=inverse(transpose(h_full));
  
 tim=rtimer;
 for (i=1; i<=red ; i++)
 {
  a[i]=transpose(submat(check,i..i,1..n));
  a[i]=transf*a[i];
 }
 
 tim=rtimer;
 matrix te[n][1]; 
 for (i=1; i<=n; i++)
 {
  for (j=1; j<=t+1; j++)
  {
   if ((j<i)&&(i<=t+1)) {c[i][j]=c[j][i];}
   else 
   {
    if (i+j<=n+1) 
    {
     c[i][j]=te;
     c[i][j][i+j-1,1]=1; 
    }
    else
    {
     c[i][j]=star(h_full,i,j);
     c[i][j]=transf*c[i][j];
    }   
   }
  }
 }
 
 
 tim=rtimer; 
 if (formal==0)
 {
  matrix s[red][1]=syndrome(check,y);
  for (j=1; j<=red; j++)
  {
   sum1=0;
   for (l=1; l<=n; l++)
   {
    sum1=sum1+a[j][l,1]*U(l);
   }
   result=result,sum1-s[j,1];
  }
 } else
 {
  for (j=1; j<=red; j++)
  {
   sum1=0;
   for (l=1; l<=n; l++)
   {
    sum1=sum1+a[j][l,1]*U(l);
   }
   result=result,sum1-s(j);
  }
  for (j=1; j<=red; j++)
  {
     result=result,s(j)^(formal)-s(j);
  }
 }
 if (fieldeq)
 {
  for (i=1; i<=n; i++)
  {
   result=result,U(i)^q-U(i);
  }
  for (j=1; j<=t; j++)
  {
     result=result,V(j)^q-V(j);
  }
 } 
 for (i=1; i<=n; i++)
 {
  sum1=0;
  for (j=1; j<=t; j++)
  {
   sum2=0;
   for (l=1; l<=n; l++)
   {
    sum2=sum2+c[i][j][l,1]*U(l);
   }
   sum1=sum1+sum2*V(j);
  }
  sum3=0;
  for (l=1; l<=n; l++)
  {
   sum3=sum3+c[i][t+1][l,1]*U(l);
  }
  result=result,sum1-sum3;
 } 
 
 result=simplify(result,2);
 
 ideal qe=result;
 export qe;
 return(work);
 //exportto(Top,h_full);   
} example 
{
     "EXAMPLE:";  echo = 2;
     //correct 2 errors in [7,3] 8-ary code RS code
     int t=2; int q=8; int n=7; int redun=4;
     ring r=(q,a),x,dp;
     matrix h_full=genMDSMat(n,a);
     matrix h=submat(h_full,1..redun,1..n);
     matrix g=dual_code(h);
     matrix x[1][3]=0,0,1,0;
     matrix y[1][7]=encode(x,g);
     //disturb with 2 errors
     matrix rec[1][7]=error(y,list(2,4),list(1,a));
     //generate the system
     def A=sysQE(h,rec,t,0,0);
     setring A;
     print(qe);
     //let us decode
     option(redSB);
     ideal sys_qe=std(qe);
     print(sys_qe);     
}

proc error(matrix y, list pos, list val)
"USAGE:  error(y, pos, val); y is a (code) word, pos = positions where errors occured, val = their corresponding values
RETURN:  corresponding received word
EXAMPLE: example error; shows an example
"
{
 matrix result[1][ncols(y)]=y;
 if (size(pos)!=size(val)) {print("ERRORerror");}
 for (int i=1; i<=size(pos); i++)
 {
  result[1,pos[i]]=y[1,pos[i]]+val[i];
 }
 return(result);
} example 
{
     "EXAMPLE:";  echo = 2;
     //correct 2 errors in [7,3] 8-ary code RS code
     int t=2; int q=8; int n=7; int redun=4;
     ring r=(q,a),x,dp;
     matrix h_full=genMDSMat(n,a);
     matrix h=submat(h_full,1..redun,1..n);
     matrix g=dual_code(h);
     matrix x[1][3]=0,0,1,0;
     matrix y[1][7]=encode(x,g);
     print(y);
     
     //disturb with 2 errors
     matrix rec[1][7]=error(y,list(2,4),list(1,a));
     print(rec);
     print(rec-y);
}

proc errorRand(matrix y, int num, int e)
"USAGE:    errorRand(y, num, e); y is a (code) word, num is the number of errors, e is an extension degree (if one wants values to
           be from GF(p^e)
RETURN:    corresponding received word
EXAMPLE:   example errorRand; shows an example
"
{
 matrix result[1][ncols(y)]=y;
 int i,j, flag, temp;
 list pos, val;
 matrix tempnum;

 for (i=1; i<=num; i++)
 {
  while(1)
  {
   temp=random(1,ncols(y));
   flag=1;
   for (j=1; j<=size(pos); j++)
   {
    if (temp==pos[j]) {flag=0;}
   }
   if (flag) {pos[i]=temp;break;}
  } 
 }
 
 for (i=1; i<=num; i++)
 {
  flag=1;
  while(flag)
  {   
   tempnum=randomvector(1,e);   
   if (tempnum!=0) {flag=0;}   
  }
  val[i]=tempnum; 
 }
  
 for (i=1; i<=size(pos); i++)
 {
  result[1,pos[i]]=y[1,pos[i]]+val[i];
 }
 return(result);
} example 
{
  "EXAMPLE:";  echo = 2;
     //correct 2 errors in [7,3] 8-ary code RS code
     int t=2; int q=8; int n=7; int redun=4;
     ring r=(q,a),x,dp;
     matrix h_full=genMDSMat(n,a);
     matrix h=submat(h_full,1..redun,1..n);
     matrix g=dual_code(h);
     matrix x[1][3]=0,0,1,0;
     matrix y[1][7]=encode(x,g);
     
     //disturb with 2 random errors
     matrix rec[1][7]=errorRand(y,2,3);
     print(rec);
     print(rec-y);     
}

proc randomCheck(int m, int n, int e, int #)
"USAGE:    randomCheck(m, n, e); m x n are dimensions of the matrix, e is an extension degree (if one wants values to
           be from GF(p^e)
RETURN:    random check matrix
EXAMPLE:   example randomCheck; shows an example
"
{
 matrix result[m][n];
 matrix rand[m][n-m];
 int i,j;
 matrix temp;
 for (i=1; i<=m; i++)
 {
  temp=randomvector(n-m,e,#);  
  for (j=1; j<=n-m; j++)
  {
   rand[i,j]=temp[j,1];
  } 
 }
 result=concat(rand,unitmat(m));
 return(result);
} example 
{
  "EXAMPLE:";  echo = 2;     
     int redun=5; int n=15;
     ring r=2,x,dp;
     
     //generate random check matrix for a [15,5] binary code
     matrix h=randomCheck(redun,n,1);
     print(h);
     
     //corresponding generator matrix
     matrix g=dual_code(h);
     print(g);     
}

proc genMDSMat(int n, number a)
"USAGE:   genMDSMat(n, a); n x n are dimensions of the matrix, a is a primitive element of the field
          An MDS matrix is constructed in the following way. We take a to be a generator of the multiplicative group of the field.
          Then we construct the Vandermonde matrix with this a.
ASSUME:   extension field should already be defined 
RETURN:   a matrix with the MDS property
EXAMPLE:  example genMDSMat; shows an example  
"
{
 int i,j;
 matrix result[n][n];
 for (i=0; i<=n-1; i++)
 {
  for (j=0; j<=n-1; j++)
  {
   result[j+1,i+1]=(a^i)^j;
  } 
 }
 return(result);
} example 
{
     "EXAMPLE:";  echo = 2;
     int q=16; int n=15;
     ring r=(q,a),x,dp;
     
     //generate an MDS (Vandermonde) matrix
     matrix h_full=genMDSMat(n,a);
     print(h_full);     
}


proc mindist (matrix check)
"USAGE:  mindist (check, q); check is a check matrix, q = field size
RETURN:  minimum distance of the code together with the time needed for its calculation
EXAMPLE: example mindist; shows an example 
"
{
 int n=ncols(check); int redun=nrows(check); int t=redun+1;
 
 def br=basering; 
 list rl=ringlist(br); 
 int q=rl[1][1];
 
 ring work=(q,a),(V(1..t),U(1..n)),dp;  
 matrix check=imap(br,check); 
 
 ideal temp;
 int count=1;
 int flag=1;
 int flag2;
 int i, tim, timsolve;
 matrix z[1][n]; 
 option(redSB);
 def A=sysQE(check,z,count,0,0);
 while (flag)
 {
    A=sysQE(check,z,count,0,0);
    setring A;
    ideal temp=qe;
    tim=rtimer;
    temp=std(temp);
    timsolve=timsolve+rtimer-tim;  
    flag2=1;
    setring work;
    temp=imap(A,temp);    
    for (i=1; i<=n; i++)
    {
      if 
        (temp[i]!=U(n-i+1)) 
        {
          flag2=0;
        }
    }
    if (!flag2) 
    {
      flag=0;
    }
    else 
    {
      count++;
    }
 }
 list result=list(count,timsolve); 
 return(result); 
} example 
{
     "EXAMPLE:";  echo = 2;
     //determine a minimum distance for a [7,3] binary code 
     int q=8; int n=7; int redun=4; int t=redun+1; 
     ring r=(q,a),x,dp;     
               
     //generate random check matrix
     matrix h=randomCheck(redun,n,1);
     print(h);
     list l=mindist(h);     
     print(l[1]);  
     //time for the comutation in secs   
     print(l[2]);     
}

proc decode(matrix check, matrix rec)
"USAGE:    decode(check, rec, t); check is the check matrix of the code
           rec is a received word, t is an upper bound for the number of errors one wants to correct
ASSUME:    Errors in rec should be correctable, otherwise the output is unpredictable
RETURN:    a codeword that is closest to rec
EXAMPLE:   example decode; shows an example 
"
{ 
 def br=basering;
 int n=ncols(check); 
 
 int count=1; 
 def A=sysQE(check,rec,count,0,0);
 while(1)
 {
  A=sysQE(check,rec,count,0,0);
  setring A;     
  matrix h_full=genMDSMat(n,a);
  matrix rec=imap(br,rec);
  option(redSB);
  ideal qe_red=std(qe);   
  if (qe_red[1]!=1)
  {
    break;
  }
  else 
  {
    count++;
  }
  setring br;
 }
 
 setring A; 
     
 //obtain a codeword
 //this works only if our code is indeed can correct these errors
 matrix syn[n][1];
 for (int i=1; i<=n; i++)
 {
  syn[i,1]=-qe_red[n-i+1]+lead(qe_red[n-i+1]);
 }   
 
 matrix real_syn=inverse(h_full)*syn;  
 setring br;
 matrix real_syn=imap(A,real_syn);
 return(rec-transpose(real_syn));     
} example 
{
     "EXAMPLE:";  echo = 2;     
     //correct 1 error in [15,7] binary code 
     int t=1; int q=16; int n=15; int redun=10;
     ring r=(q,a),x,dp;     
             
     //generate random check matrix
     matrix h=randomCheck(redun,n,1);     
     matrix g=dual_code(h);
     matrix x[1][n-redun]=0,0,1,0,1,0,1;
     matrix y[1][n]=encode(x,g);
     print(y);
     
     // find out the minimum distance of the code
     list l=mindist(h);
     
     //disturb with errors
     "Correct ",(l[1]-1) div 2," errors";
     matrix rec[1][n]=errorRand(y,(l[1]-1) div 2,1);
     print(rec);     
     
     //let us decode
     matrix dec_word=decode(h,rec);
     print(dec_word);     
}


proc solveForRandom(int n, int redun, int ncodes, int ntrials, int #)
"USAGE:    solveForRandom(redun, q, ncodes, ntrials); redun is a redundabcy of a (random) code,
           q = field size, ncodes = number of random codes to be processed
           ntrials = number of received vectors per code to be corrected
           if # is given it sets the correction capacity explicitly. It should be used in case one expexts some lower bound,
           otherwise the procedure tries to compute the real minimum distance to find out the error-correction capacity
RETURN:    nothing; 
EXAMPLE:   example solveForRandom; shows an example 
"
{
 int i,j;
 matrix h;
 int dist, t, tim, tim2, tim3, timdist, timdec, timdist2, timdec2, timdec3;
 ideal sys;
 list tmp;

 option(redSB);
 def br=basering;
 matrix h_full=genMDSMat(n,a); 
 matrix z[1][ncols(h_full)];
 int n=ncols(h_full);
 for (i=1; i<=ncodes; i++)
 {
  setring br;
  h=randomCheck(redun,n,1);
  print(h);
  if (#>0)
  {
     t=#;
  } else {
     tim=rtimer;
     tmp=mindist(h);
     timdist=timdist+rtimer-tim;
     timdist2=timdist2+tmp[2];
     dist=tmp[1];
     printf("d= %p",dist);
     t=(dist-1) div 2; 
  }     
  tim2=rtimer;  
  tim3=rtimer;
  def A=sysQE(h,z,t,0,0);
  setring A;
  matrix word,y,rec;
  ideal sys2,sys3;
  matrix h=imap(br,h);
  matrix g=dual_code(h);  
  ideal sys=qe;
  print("The system is generated");
  timdec3=timdec3+rtimer-tim3;  
  for (j=1; j<=ntrials; j++)
  {
   word=randomvector(n-redun,1);   
   y=encode(transpose(word),g);
   rec=errorRand(y,t,1);   
   sys2=add_synd(rec,h,redun,sys); 
  
   tim=rtimer;
   sys3=std(sys2);          
   sys3;
   timdec=timdec+rtimer-tim;   
  } 
  timdec2=timdec2+rtimer-tim2;  
 } 
 printf("Time for mindist: %p", timdist);
 printf("Time for GB in mindist: %p", timdist);
 printf("Time for decoding: %p", timdec2);
 printf("Time for GB in decoding: %p", timdec);
 printf("Time for sysQE in decoding: %p", timdec3); 
} example 
{
     "EXAMPLE:";  echo = 2;
     int q=32; int n=25; int redun=n-11; int t=redun+1;
     ring r=(q,a),x,dp;
     
     // correct 2 errors in 5 random binary codes, 50 trials each
     solveForRandom(n,redun,5,50,2); 
}


proc solveForCode(matrix check, int ntrials, int #)
"USAGE:     solveForCode(check, ntrials); 
            check is a check matrix for the code, ntrials = number of received vectors per code to be corrected
            if # is given it sets the correction cpacity explicitly. It should be used in case one expexts some lower bound
            otherwise the procedure tries to compute the real minimum distance to find out the error-correction capacity
RETURN:     nothing;  
EXAMPLE:    example solveForCode; shows an example 
"
{
 int n=ncols(check);
 int redun=nrows(check);
 int i,j;
 matrix h;
 int dist, t, tim, tim2, tim3, timdist, timdec, timdist2, timdec2, timdec3;
 ideal sys;
 list tmp;

 option(redSB);
 def br=basering;
 matrix h_full=genMDSMat(n,a); 
 matrix z[1][ncols(h_full)];
 int n=ncols(h_full);
 setring br;
 h=check;
 print(h);
 if (#>0)
 {
    t=#;
 } else {
   tim=rtimer;
   tmp=mindist(h);
   timdist=timdist+rtimer-tim;
   timdist2=timdist2+tmp[2];
   dist=tmp[1];
   printf("d= %p",dist);
   t=(dist-1) div 2; 
 }     
 tim2=rtimer;  
 tim3=rtimer;
 def A=sysQE(h,z,t,0,0);
 setring A;
 matrix word,y,rec;
 ideal sys2,sys3;
 matrix h=imap(br,h);
 matrix g=dual_code(h);  
 ideal sys=qe;
 print("The system is generated");
 timdec3=timdec3+rtimer-tim3;  
 for (j=1; j<=ntrials; j++)
 {
   word=randomvector(n-redun,1);   
   y=encode(transpose(word),g);
   rec=errorRand(y,t,1);   
   sys2=add_synd(rec,h,redun,sys); 
   
   tim=rtimer;
   sys3=std(sys2);          
   sys3;
   timdec=timdec+rtimer-tim;   
 } 
 timdec2=timdec2+rtimer-tim2;  
  
 printf("Time for mindist: %p", timdist);
 printf("Time for GB in mindist: %p", timdist);
 printf("Time for decoding: %p", timdec2);
 printf("Time for GB in decoding: %p", timdec);
 printf("Time for sysQE in decoding: %p", timdec3); 
} example 
{
     "EXAMPLE:";  echo = 2;
     int q=32; int n=25; int redun=n-11; int t=redun+1;
     ring r=(q,a),x,dp;
     matrix check=randomCheck(redun,n,1);
     
     // correct 2 errors in using the code above, 50 trials 
     solveForCode(check,50,2); 
}


static proc list2intvec (list l)
{
 intvec result;
 for (int i=1; i<=size(l); i++)
 {
  result[i]=l[i];
 }
 return(result);
}
   

static proc add_synd (matrix rec, matrix check, int redun, ideal sys)
{
     ideal result=sys;
     matrix s[redun][1]=syndrome(check,rec);
     for (int i=1; i<=redun; i++)
     
     {
          result[i]=result[i]-s[i,1];          
     }
     return(result);
}

static proc ev (poly f, matrix p)
{
     if (ncols(p)>1) {ERROR("not a column vector");};
     int m=size(p);
     poly temp=f;
     for (int i=1; i<=m; i++)
     {
          temp=subst(temp,x(i),p[i,1]);
     }
     return(number(temp));
}

static proc find_index (ideal G, matrix p)
{
     if (ncols(p)>1) {ERROR("not a column vector");};
     int i=1;
     int n=size(G);
     while(i<=n)
     {
          if (ev(G[i],p)!=0) {return(i);}
          i++;
     }
     return(-1);
}

static proc ideal2list (ideal id)
{
     list l;
     for (int i=1; i<=size(id); i++)
     {
          l[i]=id[i];
     }
     return(l);
}

static proc list2ideal (list l)
{
     ideal id;
     for (int i=1; i<=size(l); i++)
     {
          id[i]=l[i];
     }
     return(id);
}

static proc divisible (poly m, ideal G)
{
     for (int i=1; i<=size(G); i++)
     {
          if (m/leadmonom(G[i])!=0) {return(1);}          
     }
     return(0);
}

proc vanishId (list points)
"USAGE:  vanishId (points,e); points is a list of point that define, where polynomials from the vanishing ideal will vanish,
RETURN:  Vanishing ideal corresponding to the given set of points
EXAMPLE: example vanishId; shows an example 
"
{
     int m=size(points[1]);
     int n=size(points);
     
     ideal G=1;
     int i,k,j;
     list temp;
     poly h,cur;
     for (k=1; k<=n; k++)
     {          
          i=find_index(G,points[k]);
          cur=G[i];          
          for(j=i+1; j<=size(G); j++)
          {
               G[j]=G[j]-ev(G[j],points[k])/ev(G[i],points[k])*G[i];
          }
          G=simplify(G,2);
          temp=ideal2list(G);
          temp=delete(temp,i);
          G=list2ideal(temp);          
          for (j=1; j<=m; j++)
          {
               if (!divisible(x(j)*leadmonom(cur),G))
               {
                    attrib(G,"isSB",1);
                    h=NF((x(j)-points[k][j,1])*cur,G);                    
                    temp=ideal2list(G);
                    temp=insert(temp,h);
                    G=list2ideal(temp);
                    G=sort(G)[1];     
               }
          }
     }
     attrib(G,"isSB",1);
     return(G);
} example 
{
     "EXAMPLE:";  echo = 2;
      ring r=3,(x(1..3)),dp;
     
     //generate all 3-vectors over GF(3)
     list points=points_gen(3,1);
     
     list points2=conv_points(points);
     
     //grasps the first 11 points
     list p=grasp_list(points2,1,11);
     print(p);
     
     //construct the vanishing ideal
     ideal id=vanishId(p);
     print(id);
}

proc points_gen (int m, int e)
{
     if (e>1)
     {
     list result;
     int count=1;
     int i,j;
     list l=ringlist(basering);
     int charac=l[1][1];
     number a=par(1);
     list tmp;
     for (i=1; i<=charac^(e*m); i++)
     {
          result[i]=tmp;
     }
     if (m==1)
     {
          result[count][m]=0;
          count++;
          for (j=1; j<=charac^(e)-1; j++)
          {               
               result[count][m]=a^j;
               count++;
          }
          return(result);
     }
     list prev=points_gen(m-1,e);     
     for (i=1; i<=size(prev); i++)
     {
          result[count]=prev[i];
          result[count][m]=0;
          count++;
          for (j=1; j<=charac^(e)-1; j++)
          {
               result[count]=prev[i];
               result[count][m]=a^j;
               count++;
          }
     }
     return(result);
     }
     
     if (e==1)
     {
     list result;
     int count=1;
     int i,j;
     list l=ringlist(basering);
     int charac=l[1][1];     
     list tmp;
     for (i=1; i<=charac^m; i++)
     {
          result[i]=tmp;
     }
     if (m==1)
     {
          for (j=0; j<=charac-1; j++)
          {               
               result[count][m]=number(j);
               count++;
          }
          return(result);
     }
     list prev=points_gen(m-1,e);     
     for (i=1; i<=size(prev); i++)
     {
          for (j=0; j<=charac-1; j++)
          {
               result[count]=prev[i];
               result[count][m]=number(j);
               count++;
          }
     }
     return(result);
     }
     
}

static proc list2vec (list l)
{
     matrix m[size(l)][1];
     for (int i=1; i<=size(l); i++)
     {
          m[i,1]=l[i];
     }
     return(m);
}

proc conv_points (list points)
{
     for (int i=1; i<=size(points); i++)
     {
          points[i]=list2vec(points[i]);
     }
     return(points);
}

proc grasp_list (list l, int m, int n)
{
     list result;
     int count=1;
     for (int i=m; i<=n; i++)
     {
          result[count]=l[i];
          count++;
     }
     return(result);
}

static proc xi_gen (matrix p, int e, int s)
{
     poly prod=1;
     list rl=ringlist(basering);     
     int charac=rl[1][1];
     int l;
     for (l=1; l<=s; l++)
     {
          prod=prod*(1-(x(l)-p[l,1])^(charac^e-1));
     }
     return(prod);
}

static proc gener_funcs (matrix check, list points, int e, ideal id, int s)
{
     int n=ncols(check);
     if (n!=size(points)) {ERROR("Incompatible sizes of check and points");}
     ideal xi;
     int i,j;
     for (i=1; i<=n; i++)
     {
          xi[i]=xi_gen(points[i],e,s);
     }
     ideal result;
     int m=nrows(check);
     poly sum;
     for (i=1; i<=m; i++)
     {
          sum=0;
          for (j=1; j<=n; j++)
          {
               sum=sum+check[i,j]*xi[j];
          }
          result[i]=NF(sum,id);
     }
     return(result);
}

proc sysFL (matrix check, matrix y, int t, int e, int s)
"USAGE:    sysFL (check,y,t,e,s); check is a check matrix of the code, y is a received word, t the number of errors to correct,
           e is the extension degree, s is the dimension of the point for the vanishing ideal
RETURN:    the system of Fitzgerald-Lax for the given decoding problem
EXAMPLE:   example sysFL; shows an example
"
{
     list rl=ringlist(basering);        
     int charac=rl[1][1];     
     int n=ncols(check);
     int m=nrows(check);          
     list points=points_gen(s,e);
     list points2=conv_points(points);
     list p=grasp_list(points2,1,n);     
     ideal id=vanishId(p,e);      
     ideal funcs=gener_funcs(check,p,e,id,s);          
     
     ideal result;
     poly temp;
     int i,j,k;
     
     //vanishing realtions     
     for (i=1; i<=t; i++)
     {          
          for (j=1; j<=size(id); j++)
          {
               temp=id[j];               
               for (k=1; k<=s; k++)
               {
                    temp=subst(temp,x(k),x_var(i,k,s));
               }               
               result=result,temp;
          }
     }          
     
     //field equations
     for (i=1; i<=t; i++)
     {
          for (k=1; k<=s; k++)
          {
               result=result,x_var(i,k,s)^(charac^e)-x_var(i,k,s);
          }
     }
     for (i=1; i<=t; i++)
     {
          result=result,e(i)^(charac^e-1)-1;
     }
     
     result=simplify(result,8);
     
     //check realtions
     poly sum;
     matrix syn[m][1]=syndrome(check,y);          
     for (i=1; i<=size(funcs); i++)
     {          
          sum=0;
          for (j=1; j<=t; j++)
          {
               temp=funcs[i];
               for (k=1; k<=s; k++)
               {
                    temp=subst(temp,x(k),x_var(j,k,s));
               }
               sum=sum+temp*e(j);         
          }
          result=result,sum-syn[i,1];
     }
     
     result=simplify(result,2);
     
     points=points2;
     export points;     
     return(result);
} example
{
     "EXAMPLE:";  echo = 2;
     
     list l=FLpreprocess(3,1,11,2,"");
     def r=l[1];
     setring r;
     int s_work=l[2];     
     
     //the check matrix of [11,6,5] ternary code
     matrix h[5][11]=1,0,0,0,0,1,1,1,-1,-1,0,
          0,1,0,0,0,1,1,-1,1,0,-1,
          0,0,1,0,0,1,-1,1,0,1,-1,
          0,0,0,1,0,1,-1,0,1,-1,1,
          0,0,0,0,1,1,0,-1,-1,1,1;
     matrix g=dual_code(h);
     matrix x[1][6];
     matrix y[1][11]=encode(x,g);
     //disturb with 2 errors
     matrix rec[1][11]=error(y,list(2,4),list(1,-1));

     //the Fitzgerald-Lax system     
     ideal sys=sysFL(h,rec,2,1,s_work);
     print(sys);
     option(redSB);
     ideal red_sys=std(sys);
     red_sys; // read the solutions from this redGB
     // the points are (0,0,1) and (0,1,0) with error values 1 and -1 resp.
     // use list points to find error positions;
     points;     
}

proc FLpreprocess (int p, int e, int n, int t, string minp)
{
     ring r1=p,x,dp;
     int s=1;
     while(p^(s*e)<n)
     {
          s++;
     }     
     list var_ord;
     int i,j;
     int count=1;
     for (i=s; i>=1; i--)
     {
          var_ord[count]=string("x("+string(i)+")");
          count++;
     }
     for (i=t; i>=1; i--)
     {
          var_ord[count]=string("e("+string(i)+")");
          count++;
          for (j=s; j>=1; j--)
          {
               var_ord[count]=string("x1("+string(s*(i-1)+j)+")");
               count++;
          }
     }     
     
     list rl;
     list tmp;
     
     if (e>1)
     {
          rl[1]=tmp;
          rl[1][1]=p;
          rl[1][2]=tmp;
          rl[1][2][1]=string("a");
          rl[1][3]=tmp;
          rl[1][3][1]=tmp;
          //rl[1][3][1][1]=string("dp("+string((t-1)*(s+1)+s)+"),lp("+string(s+1)+")");
          rl[1][3][1][1]=string("lp");
          rl[1][3][1][2]=1;
          rl[1][4]=ideal(0);
     } else {
          rl[1]=p;
     }
     
     rl[2]=var_ord;
     
     rl[3]=tmp;
     rl[3][1]=tmp;     
     //rl[3][1][1]=string("dp("+string((t-1)*(s+1)+s)+"),lp("+string(s+1)+")");
     rl[3][1][1]=string("lp");
     intvec v=1;
     for (i=1; i<=size(var_ord)-1; i++)
     {
          v=v,1;
     }
     rl[3][1][2]=v; 
     rl[3][2]=tmp;
     rl[3][2][1]=string("C");    
     rl[3][2][2]=intvec(0);
     
     rl[4]=ideal(0);     
     
     def r2=ring(rl);
     setring r2;
     list l=ringlist(r2);     
     if (e>1)
     {
          execute(string("poly f="+minp));     
          ideal id=f;     
          l[1][4]=id;
     }
     
     def r=ring(l);
     setring r;      
     
     return(list(r,s));
}

static proc x_var (int i, int j, int s)
{     
     return(x1(s*(i-1)+j));
}

static proc randomvector(int n, int e, int #)
{     
    int i;
    matrix result[n][1];
    for (i=1; i<=n; i++)
    {
        result[i,1]=asElement(random_prime_vector(e,#));
    }
    return(result);        
}

static proc asElement(list l)
{
  number s;
  int i;
  number w=1;
  if (size(l)>1) {w=par(1);}  
  for (i=0; i<=size(l)-1; i++)
  {
    s=s+w^i*l[i+1];
  }
  return(s);
}

proc random_prime_vector (int n, int #)
{
  if (#==1)
  {
     list rl=ringlist(basering);     
     int charac=rl[1][1];
  } else {
     int charac=2;
  }
  list l;
  int i;
  for (i=1; i<=n; i++)
  {
    l=l+list(random(0,charac-1));
  }
  return(l);
}

proc FLSolveForRandom(int n, int redun, int p, int e, int t, int ncodes, int ntrials, string minpol)
"USAGE:    FLSolveForRandom(redun,p,e,n,t,ncodes,ntrials,minpol); n = length of codes generated, redun = redundancy of codes generated, 
           p = characteristics, e is the extension degree, 
           q = number of errors to correct, ncodes = number of random codes to be processed
           ntrials = number of received vectors per code to be corrected
           due to some pecularities of SINGULAR one needs to provide minimal polynomial for the extension explicitly
RETURN:    nothing
EXAMPLE:   example FLSolveForRandom; shows an example
{
 list l=FLpreprocess(p,e,n,t,minpol); 
 
 def r=l[1];
 int s_work=l[2]; 
 export(s_work);
 setring r;
 
 int i,j;
 matrix h, g, word, y, rec;
 list l;
 int dist, tim, tim2, tim3, timdist, timdec, timdist2, timdec2, timdec3;
 ideal sys, sys2, sys3;
 list tmp; 

 option(redSB); 
 matrix z[1][n]; 
 
 for (i=1; i<=ncodes; i++)
 {
     h=randomCheck(redun,n,e,1);            
     g=dual_code(h);
     tim2=rtimer;       
     tim3=rtimer;               
     sys=sysFL(h,z,t,e,s_work);     
     timdec3=timdec3+rtimer-tim3;
   
     for (j=1; j<=ntrials; j++)
     {
          word=randomvector(n-redun,1);
          y=encode(transpose(word),g);
          rec=errorRand(y,t,e);
          sys2=LF_add_synd(rec,h,sys);   
          tim=rtimer;
          sys3=std(sys2);
          timdec=timdec+rtimer-tim;
     }            
     timdec2=timdec2+rtimer-tim2;
 } 
 
 printf("Time for decoding: %p", timdec2);
 printf("Time for GB in decoding: %p", timdec);
 printf("Time for generating Fitzgerald-Lax system during decoding: %p", timdec3); 
} example
{
     "EXAMPLE:";  echo = 2;
     
     // decoding for one random binary code of length 25, redundancy 14; 300 words are processed
     FLSolveForRandom(25,14,2,1,1,1,300,"");
}

static proc LF_add_synd (matrix rec, matrix check, ideal sys)
{
     int redun=nrows(check);
     ideal result=sys;
     matrix s[redun][1]=syndrome(check,rec);
     for (int i=size(sys)-redun+1; i<=size(sys); i++)
     {
          result[i]=result[i]-s[i-size(sys)+redun,1];
     }
     return(result);
}


/*
//////////////     SOME HARD EXAMPLES    //////////////////////
//////      THAT MAYBE WILL BE DOABLE LATER     ///////////////

1.) These random instances are not doable in <=1000 sec.

"EXAMPLE:";  echo = 2;
int q=128; int n=120; int redun=n-40; 
ring r=(q,a),x,dp;
solveForRandom(n,redun,1,1,6);

redun=n-30;
solveForRandom(n,redun,1,1,8);

redun=n-20;
solveForRandom(n,redun,1,1,12);

redun=n-10;
solveForRandom(n,redun,1,1,24);

int q=256; int n=150; int redun=n-10; 
ring r=(q,a),x,dp;
solveForRandom(n,redun,1,1,26);


2.) Generic decoding is hard!

int q=32; int n=31; int redun=n-16; int t=3;
ring r=(q,a),(V(1..n),U(n..1),s(redun..1)),(dp(n),lp(n),dp(redun));
matrix check[redun][n]= 1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,
0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,
0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,
0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,
0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,
1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,
1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,
1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,
0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,
0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,
0,1,0,1,0,0,1,0,0,1;
matrix rec[1][n];

def A=sysQE(check,rec,t,1,2);
setring A;
print(qe);
ideal red_qe=stdfglm(qe); 

*/