source: git/Singular/LIB/decodegb.lib @ 3e0148

///////////////////////////////////////////////////////////////////////////////
version="$Id: decodegb.lib,v 1.11 2009-04-06 13:30:27 Singular 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 @ref{Decoding codes with GB}.
 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 are 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: the ring to work with the CRHT-ideal (with Sala's additions),
        containig an ideal with name 'crht'
THEORY:  Based on 'defset' of the given cyclic code, the procedure constructs
         the corresponding Cooper-Reed-Heleseth-Truong ideal 'crht'. With its
         help one can solve the decoding problem. For basics of the method @ref{Cooper philosophy}.
SEE ALSO: sysNewton, sysBin
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:  the ring to work with the Sala's ideal for the minimum distance
         containing the ideal with name 'crht_md'
THEORY:  Based on 'defset' of the given cyclic code, the procedure constructs
         the corresponding Cooper-Reed-Heleseth-Truong ideal 'crht_md'. With
         its help one can find minimum distance of the code in the binary
         case. For basics of the method @ref{Cooper philosophy}.
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 is 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:  the ring to work with the generalized Newton identities (in
         triangular form if applicable) containing the ideal with name 'newton'
THEORY:  Based on 'defset' of the given cyclic code, the procedure constructs
         the corresponding ideal 'newton' with the generalized Newton
         identities. With its help one can solve the decoding problem. For
         basics of the method @ref{Generalized Newton identities}.
SEE ALSO: sysCRHT, sysBin
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 defining set of the code,
          - n the length,
          - odd is an additional parameter: if
           set to 1, then the defining set is enlarged by odd elements,
           which are 2^(some power)*(some elment in the def.set) mod n
@end format
RETURN:    the ring with the resulting system called 'bin'
THEORY:  Based on Q of the given cyclic code, the procedure constructs
         the corresponding ideal 'bin' with the use of the Waring function.
         With its help one can solve the decoding problem.
         For basics of the method @ref{Generalized Newton identities}.
SEE ALSO: sysNewton, sysCRHT
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 a parity 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 equal to the number of elements in the INITIAL alphabet) one
          needs to set formal>0 for the exponent
@end format
RETURN:   the ring to work with together with the resulting system called 'qe'
THEORY:  Based on 'check' of the given linear code, the procedure constructs
         the corresponding ideal that gives an opportunity to compute
         unknown syndrome of the received word y. After computing the unknown
         syndromes one is able to solve the decoding problem.
         For basics of the method @ref{Decoding method based on quadratic equations}.
SEE ALSO: sysFL
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 are 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 is matrix, num,e are 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:   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.
SEE ALSO: Decoding method based on quadratic equations
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:  minimum distance of the code together with the time needed for its
         calculation
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
NOTE:     The method described in @ref{Decoding method based on quadratic equations}
          is used for decoding.
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
"
{
 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 parity 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
        'points' is a list of points for which the vanishing ideal is to be
        constructed
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 parity 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
THEORY:  Based on 'check' of the given linear code, the procedure constructs
         the corresponding ideal constructed with a generalization of
         Cooper's philosophy. For basics of the method @ref{Fitzgerald-Lax method}.
SEE ALSO: sysQE
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);

*/