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

///////////////////////////////////////////////////////////////////////////////
version="$Id: decodegb.lib,v 1.6 2008-10-08 16:52:24 bulygin Exp $";
category="Coding theory";
info="
LIBRARY: decodegb.lib         Decoding and min distance of linear codes with GB
AUTHOR:  Stanislav Bulygin,   bulygin@mathematik.uni-kl.de                     

OVERVIEW:
 In this library we generate several systems used for decoding cyclic codes and
 finding their 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 some auxiliary functions for further 
 manipulations and decoding. For an overview of the methods mentioned above, 
 cf. 
 @format
  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) (*).
 @end format
 For the vanishing ideal computation the algorithm of Farr and Gao is 
 implemented.

MAIN PROCEDURES: 
 sysCRHT(..);        generates the CRHT-ideal as in Cooper's philosophy
 sysCRHTMindist(..); CRHT-ideal to find the minimum distance in the binary case
 sysNewton(..);      generates the ideal with the generalized Newton identities
 sysBin(..);         generates Bin system using Waring function
 encode(x,g);        encodes given message x with the given generator matrix g
 syndrome(h,c);      computes a syndrome w.r.t. the given check matrix
 sysQE(..);          generates the system of quadratic equations for decoding
 errorInsert(..);    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 a code
 decode(rec);        decoding of a word using the system of quadratic equations 
 decodeRandom(..); a procedure for manipulation with random codes
 decodeCode(..);   a procedure for manipulation with the given code
 vanishId(points);   computes the vanishing ideal for the given set of points
 sysFL(..);          generates the Fitzgerald-Lax system
 decodeRandomFL(..); manipulation with random codes via Fitzgerald-Lax
 

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

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

///////////////////////////////////////////////////////////////////////////////
// creates a list result, where result[i]=i, 1<=i<=n
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);
}

///////////////////////////////////////////////////////////////////////////////
// creates a list of all combinations without repititions of m objects out of n
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);
}


///////////////////////////////////////////////////////////////////////////////
// the polynomial for Sala's restrictions 
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, list #)
"USAGE:   sysCRHT(n,defset,e,q,m,[k]); n,e,q,m,k int, defset list of int's
@format
         - n length of the cyclic code, 
         - defset is a list representing the defining set,
         - e the error-correcting capacity, 
         - q field size
         - m degree extension of the splitting field, 
         - if k>0 additional equations representing the fact that every two 
         error positions are either different or at least one of them is zero
@end format
RETURN:  
@format
        the ring to work with the CRHT-ideal (with Sala's additions), 
        containig an ideal with name 'crht'
@end format
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;
  int k;
  if ( size(#) > 0) 
  { 
    k = #[1]; 
  }
  
  //---------------------- add 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);
  }
  
  //--------------------- field equations on syndromes ------------------
  for (i=1; i<=r; i++)
  {
    crht=crht,X(i)^(q^m)-X(i);
  }
  
  //------ restrictions on error-locations: 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;
  }  
  
  //--------- add Sala's additional conditions if necessary --------------
  if ( k > 0 )
  
  {
    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)
  intvec v = option(get);

  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);                // CRHT-ideal with additional equations from Sala
  option(redSB);
  ideal red_crht=std(crht);   // reduced Groebner basis
  print(red_crht);            
  red_crht[5];                // general error-locator polynomial for this code
  option(set,v);
}

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


proc sysCRHTMindist (int n, list defset, int w)
"USAGE:  sysCRHTMindist(n,defset,w);  n,w are int, defset is list of int's  
@format
        - n length of the cyclic code, 
        - defset is a list representing the defining set,
        - w is a candidate for the minimum distance
@end format
RETURN:
@format
         the ring to work with the Sala's ideal for the minimum distance 
         containing the ideal with name 'crht_md'
@end format
EXAMPLE: example sysCRHTMindist; shows an example
"
{
  int r=size(defset);  
  ring @crht_md=2,Z(1..w),lp;
  ideal crht_md;
  int i,j;
  poly sum;
  
  //------------ add 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;
  }  
  
  
  //----------- locations are n-th roots of unity ------------
  for (i=1; i<=w; i++)
  {
    crht_md=crht_md,Z(i)^n-1;
  }  
  
  //------------ adding conditions on locations being different ------------
  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;
  intvec v = option(get);
  // 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=sysCRHTMindist(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);     
  print(red_crht_md);          // reduced Groebner basis
  
  option(set,v);
}

///////////////////////////////////////////////////////////////////////////////
// slightly modified mod function
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, list #)
"USAGE:   sysNewton (n,defset,t,q,m,[tr]); n,t,q,m,tr int, defset list int's 
@format
         - 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 tr>0 it indicates that Newton identities in triangular 
           form should be constructed
@end format
RETURN:   
@format
        the ring to work with the generalized Newton identities (in triangular 
        form if applicable) containing the ideal with name 'newton'
@end format
EXAMPLE:  example sysNewton; shows an example
"
{  
 string s="ring @newton=("+string(q)+",a),(";
 int i,j;
 int flag;
 int tr;
 
 if (size(#)>0)
 {
  tr=#[1];
 }
 
 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 (tr)
 {
  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;
 } 
 
 //----------- add field equations on sigma's --------------
 for (i=1; i<=t; i++)
 {
  newton=newton,sigma(i)^(q^m)-sigma(i);
 }
 
 //----------- add 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 tr=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,tr);
     setring A;     
     print(newton);     // generalized Newton identities in triangular form
}

///////////////////////////////////////////////////////////////////////////////
// forms a list of special combinations needed for computation of Waring's 
//function
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);
}

///////////////////////////////////////////////////////////////////////////////
//if n=q^e*m, m and n are coprime, then return e
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, list #)
"USAGE:    sysBin (v,Q,n,[odd]);  v,n,odd are int, Q is list of int's
@format
          - v a number if errors, 
          - Q is a generating set of the code, 
          - n the length, 
          - odd is an additional parameter: if  
           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
@end format
THEORY:    
@format
          The system using Waring's function due to Augot et.al. 
          is built as in [*].
@end format
RETURN:    the ring with the resulting system called 'bin'
EXAMPLE:   example sysBin; shows an example
"
{
 int odd;
 if (size(#)>0)
 {
  odd=#[1];
 }

 //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 (odd==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;
 }
 
 //---- form polynomials for the Bin system via Waring's function ---------  
 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));
}

///////////////////////////////////////////////////////////////////////////////
// (coordinatewise) star product of two vectors
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, list #)
"USAGE:   sysQE(check,y,t,[fieldeq,formal]);check,y matrix;t,fieldeq,formal int
@format
        - 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, field 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
@end format
THEORY:   The system of quadratic equations is built as in [*].
RETURN:   the ring to work with together with the resulting system called 'qe'
EXAMPLE:  example sysQE; shows an example
"
{
 int fieldeq;
 int formal;
 if (size(#)>0)
 {
  fieldeq=#[1];
 }
 if (size(#)>1)
 {
  formal=#[2];
 }
 
 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));
  
 //------ expression matrix of check vectors w.r.t. the MDS basis -----------
 tim=rtimer;
 for (i=1; i<=red ; i++)
 {
  a[i]=transpose(submat(check,i..i,1..n));
  a[i]=transf*a[i];
 }
 
 //----------- compute the structure constants ------------------------
 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);
  }
 } 
 
 //----- form the quadratic equations according to the theory -----------
 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); 
} example 
{
     "EXAMPLE:";  echo = 2;
     intvec v = option(get);
     
     //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]=errorInsert(y,list(2,4),list(1,a)); 
     
     //generate the system
     def A=sysQE(h,rec,t);
     setring A;
     print(qe);
     
     //let us decode
     option(redSB);
     ideal sys_qe=std(qe);
     print(sys_qe);     
     
     option(set,v);
}

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

proc errorInsert(matrix y, list pos, list val)
"USAGE:  errorInsert(y,pos,val); y is matrix, pos,val list of int's
@format
        - y is a (code) word, 
        - pos = positions where errors occured, 
        - val = their corresponding values
@end format
RETURN:  corresponding received word
EXAMPLE: example errorInsert; 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]=errorInsert(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 matrix, num,e int
@format
          - 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))
@end format
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)
"USAGE:    randomCheck(m, n, e); m,n,e are int
@format
          - m x n are dimensions of the matrix, 
          - e is an extension degree (if one wants values to be from GF(p^e))
@end format
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 is int, a is number
@format
        - n x n are dimensions of the MDS matrix, 
        - a is a primitive element of the field.
@end format    
NOTE:     
@format
        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.
@end format
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 matrix, q int
@format
        - check is a check matrix, 
        - q is the field size
@end format
RETURN:
@format
        minimum distance of the code together with the time needed for its 
        calculation
@end format
EXAMPLE: example mindist; shows an example 
"
{
 intvec vopt = option(get);

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

 //proceed with solving the system w.r.t zero vector until some solutions 
 //are found
 while (flag)
 {
    A=sysQE(check,z,count);
    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); 
 
 option(set,vopt);
 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, rec matrix, t int
@format
          - 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
@end format
ASSUME:
@format
          Errors in rec should be correctable, otherwise the output is 
          unpredictable
@end format
RETURN:    a codeword that is closest to rec
EXAMPLE:   example decode; shows an example 
"
{
 intvec vopt = option(get); 
 
 def br=basering;
 int n=ncols(check); 
 
 int count=1; 
 def A=sysQE(check,rec,count);
 while(1)
 {
  A=sysQE(check,rec,count);
  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);
 
 option(set,vopt);
 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 decodeRandom(int n, int redun, int ncodes, int ntrials, list #)
"USAGE:    decodeRandom(redun,q,ncodes,ntrials,[e]); all parameters int
@format
          - redun is a redundabcy of a (random) code,
          - q is the field size, 
          - ncodes is the number of random codes to be processed,
          - ntrials is the number of received vectors per code to be corrected
          - If e is given it sets the correction capacity explicitly. It 
          should be used in case one expects some lower bound,
          otherwise the procedure tries to compute the real minimum distance 
          to find out the error-correction capacity
@end format
RETURN:    nothing; 
EXAMPLE:   example decodeRandom; shows an example 
"
{
 intvec vopt = option(get);

 int i,j;
 matrix h;
 int dist, t, tim, tim2, tim3, timdist, timdec, timdist2, timdec2, timdec3;
 ideal sys;
 list tmp;
 int e;
 if (size(#)>0)
 {
  e=#[1];
 }

 option(redSB);
 def br=basering;
 matrix h_full=genMDSMat(n,a); 
 matrix z[1][ncols(h_full)];

 //------------------ determine error-correction capacity -------------------
 for (i=1; i<=ncodes; i++)
 {
  setring br;
  h=randomCheck(redun,n,1);
  "check matrix:";
  print(h);
  if (e>0)
  {
     t=e;
  } 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;

  //------------- generate the template system ----------------------
  def A=sysQE(h,z,t);
  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;  

  //------ modify the template according to every received word --------------
  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);   
   timdec=timdec+rtimer-tim;   
  } 
  timdec2=timdec2+rtimer-tim2;
  kill A;
  option(set,vopt);
 } 
 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
     decodeRandom(n,redun,5,50,2); 
}

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


proc decodeCode(matrix check, int ntrials, list #)
"USAGE:     decodeCode(check, ntrials, [e]); check matrix, ntrials,e int 
@format            
           - check is a check matrix for the code, 
           - ntrials is the number of received vectors per code to be 
           corrected.
           - If e is given it sets the correction capacity explicitly. It 
           should be used in case one expects some lower bound,
           otherwise the procedure tries to compute the real minimum distance 
           to find out the error-correction capacity
@end format
RETURN:     nothing;  
EXAMPLE:    example decodeCode; shows an example 
"
{
 intvec vopt = option(get);

 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;
 int e;
 if (size(#)>0)
 {
  e=#[1];
 }

 option(redSB);
 def br=basering;
 matrix h_full=genMDSMat(n,a); 
 matrix z[1][ncols(h_full)];
 setring br;
 h=check;
 "check matrix:";
 print(h);

 //------------------ determine error-correction capacity -------------------
 if (e>0)
 {
    t=e;
 } 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;

 //------------- generate the template system ----------------------
 def A=sysQE(h,z,t);
 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; 

 //--- modify the template according to every received word --------------- 
 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);   
   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); 
 
 option(set,vopt);
} 
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 
     decodeCode(check,50,2); 
}


///////////////////////////////////////////////////////////////////////////////
// adding syndrome values to the template system
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);
}

///////////////////////////////////////////////////////////////////////////////
// evaluate a polynomial at a given point
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));
}

///////////////////////////////////////////////////////////////////////////////
// return index of an element in the ideal where it does not vanish at the 
//given point
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);
}

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

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

///////////////////////////////////////////////////////////////////////////////
// check whether given polynomial is divisible by some leading monomial of the 
//ideal
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); point is a list of matrices
@format
        - points is a list of points for which the vanishing ideal is to be 
        constructed        
@end format
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;

     //------------- proceed according to Farr-Gao algorithm ----------------
     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=pointsGen(3,1);
     
     list points2=convPoints(points);
     
     //grasps the first 11 points
     list p=graspList(points2,1,11);
     print(p);
     
     //construct the vanishing ideal
     ideal id=vanishId(p);
     print(id);
}

///////////////////////////////////////////////////////////////////////////////
// construct the list of all vectors of length m with elements in p^e, where p 
//is characteristics
proc pointsGen (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=pointsGen(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=pointsGen(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);
     }
     
}

///////////////////////////////////////////////////////////////////////////////
// convert list to a column vector
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);
}

///////////////////////////////////////////////////////////////////////////////
// convert all the point in the list with list2vec
proc convPoints (list points)
{
     for (int i=1; i<=size(points); i++)
     {
          points[i]=list2vec(points[i]);
     }
     return(points);
}

///////////////////////////////////////////////////////////////////////////////
// extracts elements from l in the range m..n
proc graspList (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);
}

///////////////////////////////////////////////////////////////////////////////
// "characteristic" polynomial
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);
}

///////////////////////////////////////////////////////////////////////////////
// generating polynomials in Fitzgerald-Lax construction
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,y matrix, t,e,s int
@format
          - 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
@end format
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=pointsGen(s,e);
     list points2=convPoints(points);
     list p=graspList(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;
     
     //--------------- add 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;
          }
     }          
     
     //--------------- add 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);
     
     //--------------- add 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;
     intvec vopt = option(get);
     
     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]=errorInsert(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;
     
     option(set,vopt);
}

///////////////////////////////////////////////////////////////////////////////
// preprocessing steps for the Fitzgerald-Lax scheme
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("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("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));
}

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

///////////////////////////////////////////////////////////////////////////////
// random vector of length n with entries from p^e, p the characteristic
static proc randomvector(int n, int e)
{     
    int i;
    matrix result[n][1];
    for (i=1; i<=n; i++)
    {
        result[i,1]=asElement(random_prime_vector(e));
    }
    return(result);        
}

///////////////////////////////////////////////////////////////////////////////
// "convert" representation of an element from the field extension from vector 
//to an elelemnt
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);
}

///////////////////////////////////////////////////////////////////////////////
// random vector of length n with entries from p, p the characteristic
static proc random_prime_vector (int n)
{  
  list rl=ringlist(basering);
  int i, charac;
  for (i=2; i<=rl[1][1]; i++)
  {
    if (rl[1][1] mod i ==0)
    {      
      break;
    }
  }
  charac=i;  
  
  list l;
  
  for (i=1; i<=n; i++)
  {
    l=l+list(random(0,charac-1));
  }
  return(l);
}

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

proc decodeRandomFL(int n, int redun, int p, int e, int t, int ncodes, int ntrials, string minpol)
"USAGE:    decodeRandomFL(redun,p,e,n,t,ncodes,ntrials,minpol); 
@format
          - n is length of codes generated, redun = redundancy of codes 
          generated, 
          - p is characteristics, 
          - e is the extension degree, 
          - t is the number of errors to correct, 
          - ncodes is the number of random codes to be processed,
          - ntrials is the number of received vectors per code to be corrected,
          - minpol: due to some pecularities of SINGULAR one needs to provide 
          minimal polynomial for the extension explicitly
@end format
RETURN:    nothing
EXAMPLE:   example decodeRandomFL; shows an example
"
{
 intvec vopt = option(get);

 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;
 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);            
     g=dual_code(h);
     tim2=rtimer;       
     tim3=rtimer;

     //---------------- generate the template system -----------------------  
     sys=sysFL(h,z,t,e,s_work);     
     timdec3=timdec3+rtimer-tim3;
   
     //------ modifying the template according to the received word ---------
     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); 
 
 option(set,vopt);
} 
example
{
     "EXAMPLE:";  echo = 2;
     
     // correcting one error for one random binary code of length 25, 
     //redundancy 14; 300 words are processed
     decodeRandomFL(25,14,2,1,1,1,300,"");
}

///////////////////////////////////////////////////////////////////////////////
// add syndrome values to the template system in FL
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 RELATIVELY EASY EXAMPLES    //////////////
///////////////////  THAT RUN AROUND ONE MINUTE  ////////////////

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

int q=128; int n=120; int redun=n-20; 
ring r=(q,a),x,dp;
decodeRandom(n,redun,1,1,9);

int q=128; int n=120; int redun=n-10;
ring r=(q,a),x,dp;
decodeRandom(n,redun,1,1,19);

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

//////////////     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;
decodeRandom(n,redun,1,1,6);

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

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

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

int q=256; int n=150; int redun=n-10; 
ring r=(q,a),x,dp;
decodeRandom(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); 

*/