source: git/Singular/LIB/decodegb.lib @ 6e118cf

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

OVERVIEW:
 In this library we generate several systems used for decoding cyclic codes and finding 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. 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) (*).
 For the vanishing ideal computation the algorithm of Farr and Gao is implemented.

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

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

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

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

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


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

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

///////////////////////////////////////////////////////////////////////////////
// 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, int#)
"USAGE:    sysBin (v, Q, n, #);  v a number if errors, Q is a generating set of the code, n the length, # 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
RETURN:    keeps the ring with the resulting system, which ideal is called 'bin'
EXAMPLE:   example sysBin; shows an example
"
{
 //ring r=2,(sigma(1..v),S(1..n)),(lp(v),dp(n));
 ring r=2,(S(1..n),sigma(1..v)),lp;
 list cyclot;
 ideal result;
 int i,j,k,s;
 list comb;
 poly sum_, mon;
 int count1, count2, upper, coef_, flag, gener;
 list Q_update;
 if (#==1)
 {
  for (i=1; i<=n; i++)
  {
   cyclot[i]=0;
  }
  for (i=1; i<=size(Q); i++)
  {
   flag=1;
   gener=Q[i];
   while(flag)
   {
    cyclot[gener]=1;
    gener=2*gener mod n;
    if (gener == Q[i]) {flag=0;}
   }
  }
  for (i=1; i<=n; i++)
  {
   if ((cyclot[i] == 1)&&(i mod 2 == 1)) {Q_update[size(Q_update)+1]=i;}
  }
 }
 else
 {
  Q_update=Q;
 }
 
 //-------------- 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, int fieldeq, int formal)
"USAGE:   sysQE(check, y, t, fieldeq, formal); check is the check matrix of the code   
          y is a received word, t the number of errors to be corrected, 
          if fieldeq=1, then field equations are added; if formal=0, 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
RETURN:   the ring to work with together with the resulting system
EXAMPLE:  example sysQE; shows an example
"
{
 def br=basering; 
 list rl=ringlist(br);
 
 int red=nrows(check);
 int n=ncols(check);
 int q=rl[1][1];
 
 if (formal==0)
 { 
  ring work=(q,a),(V(1..t),U(1..n)),dp;
 } else
 {
  ring work=(q,a),(V(1..t),U(1..n),s(1..red)),(dp(t),lp(n),dp(red));
 }
 
 matrix check=imap(br,check);
 matrix y=imap(br,y);
 
 matrix h_full=genMDSMat(n,a);
 matrix h=submat(h_full,1..red,1..n);
 if (nrows(y)!=1) {print("ERROR1Pell");}
 if (ncols(y)!=n) {print("ERROR2Pell");}
  
 ideal result;
 
 list c;
 list a;
 list tmp,tmp2;
 int i,j,l,k;
 number sum,prod,sig;
        poly sum1,sum2,sum3;
 for (i=1; i<=n; i++)
 {
  c[i]=tmp;
 }
 
 int tim=rtimer;
 matrix transf=inverse(transpose(h_full));
  
 //----------- 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);
 //exportto(Top,h_full);   
} example 
{
     "EXAMPLE:";  echo = 2;
     //correct 2 errors in [7,3] 8-ary code RS code
     int t=2; int q=8; int n=7; int redun=4;
     ring r=(q,a),x,dp;
     matrix h_full=genMDSMat(n,a);
     matrix h=submat(h_full,1..redun,1..n);
     matrix g=dual_code(h);
     matrix x[1][3]=0,0,1,0;
     matrix y[1][7]=encode(x,g);
     //disturb with 2 errors
     matrix rec[1][7]=errorInsert(y,list(2,4),list(1,a));
     //generate the system
     def A=sysQE(h,rec,t,0,0);
     setring A;
     print(qe);
     //let us decode
     option(redSB);
     ideal sys_qe=std(qe);
     print(sys_qe);     
}

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

proc errorInsert(matrix y, list pos, list val)
"USAGE:  errorInsert(y, pos, val); y is a (code) word, pos = positions where errors occured, val = their corresponding values
RETURN:  corresponding received word
EXAMPLE: example error; shows an example
"
{
 matrix result[1][ncols(y)]=y;
 if (size(pos)!=size(val)) {print("ERRORerror");}
 for (int i=1; i<=size(pos); i++)
 {
  result[1,pos[i]]=y[1,pos[i]]+val[i];
 }
 return(result);
} example 
{
     "EXAMPLE:";  echo = 2;
     //correct 2 errors in [7,3] 8-ary code RS code
     int t=2; int q=8; int n=7; int redun=4;
     ring r=(q,a),x,dp;
     matrix h_full=genMDSMat(n,a);
     matrix h=submat(h_full,1..redun,1..n);
     matrix g=dual_code(h);
     matrix x[1][3]=0,0,1,0;
     matrix y[1][7]=encode(x,g);
     print(y);
     
     //disturb with 2 errors
     matrix rec[1][7]=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 a (code) word, num is the number of errors, e is an extension degree (if one wants values to
           be from GF(p^e))
RETURN:    corresponding received word
EXAMPLE:   example errorRand; shows an example
"
{
 matrix result[1][ncols(y)]=y;
 int i,j, flag, temp;
 list pos, val;
 matrix tempnum;

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

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

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

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

proc genMDSMat(int n, number a)
"USAGE:   genMDSMat(n, a); n x n are dimensions of the matrix, a is a primitive element of the field.          
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
EXAMPLE:  example genMDSMat; shows an example  
"
{
 int i,j;
 matrix result[n][n];
 for (i=0; i<=n-1; i++)
 {
  for (j=0; j<=n-1; j++)
  {
   result[j+1,i+1]=(a^i)^j;
  } 
 }
 return(result);
} example 
{
     "EXAMPLE:";  echo = 2;
     int q=16; int n=15;
     ring r=(q,a),x,dp;
     
     //generate an MDS (Vandermonde) matrix
     matrix h_full=genMDSMat(n,a);
     print(h_full);     
}

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


proc mindist (matrix check)
"USAGE:  mindist (check, q); check is a check matrix, q = field size
RETURN:  minimum distance of the code together with the time needed for its calculation
EXAMPLE: example mindist; shows an example 
"
{
 int n=ncols(check); int redun=nrows(check); int t=redun+1;
 
 def br=basering; 
 list rl=ringlist(br); 
 int q=rl[1][1];
 
 ring work=(q,a),(V(1..t),U(1..n)),dp;  
 matrix check=imap(br,check); 
 
 ideal temp;
 int count=1;
 int flag=1;
 int flag2;
 int i, tim, timsolve;
 matrix z[1][n]; 
 option(redSB);
 def A=sysQE(check,z,count,0,0);

 //----------- proceed with solving the system w.r.t zero vector until some solutions are found -------------------- 
 while (flag)
 {
    A=sysQE(check,z,count,0,0);
    setring A;
    ideal temp=qe;
    tim=rtimer;
    temp=std(temp);
    timsolve=timsolve+rtimer-tim;  
    flag2=1;
    setring work;
    temp=imap(A,temp);    
    for (i=1; i<=n; i++)
    {
      if 
        (temp[i]!=U(n-i+1)) 
        {
          flag2=0;
        }
    }
    if (!flag2) 
    {
      flag=0;
    }
    else 
    {
      count++;
    }
 }
 list result=list(count,timsolve); 
 return(result); 
} example 
{
     "EXAMPLE:";  echo = 2;
     //determine a minimum distance for a [7,3] binary code 
     int q=8; int n=7; int redun=4; int t=redun+1; 
     ring r=(q,a),x,dp;     
               
     //generate random check matrix
     matrix h=randomCheck(redun,n,1);
     print(h);
     list l=mindist(h);     
     print(l[1]);  
     //time for the comutation in secs   
     print(l[2]);     
}

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

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

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


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

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

 //------------------ determine error-correction capacity -------------------
 for (i=1; i<=ncodes; i++)
 {
  setring br;
  h=randomCheck(redun,n,1);
  print(h);
  if (#>0)
  {
     t=#;
  } else {
     tim=rtimer;
     tmp=mindist(h);
     timdist=timdist+rtimer-tim;
     timdist2=timdist2+tmp[2];
     dist=tmp[1];
     printf("d= %p",dist);
     t=(dist-1) div 2; 
  }     
  tim2=rtimer;  
  tim3=rtimer;

  //------------- generate the template system ----------------------
  def A=sysQE(h,z,t,0,0);
  setring A;
  matrix word,y,rec;
  ideal sys2,sys3;
  matrix h=imap(br,h);
  matrix g=dual_code(h);  
  ideal sys=qe;
  print("The system is generated");
  timdec3=timdec3+rtimer-tim3;  

  //------------- 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;
 } 
 printf("Time for mindist: %p", timdist);
 printf("Time for GB in mindist: %p", timdist);
 printf("Time for decoding: %p", timdec2);
 printf("Time for GB in decoding: %p", timdec);
 printf("Time for sysQE in decoding: %p", timdec3); 
} example 
{
     "EXAMPLE:";  echo = 2;
     int q=32; int n=25; int redun=n-11; int t=redun+1;
     ring r=(q,a),x,dp;
     
     // correct 2 errors in 5 random binary codes, 50 trials each
     solveForRandom(n,redun,5,50,2); 
}

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


proc solveForCode(matrix check, int ntrials, int #)
"USAGE:     solveForCode(check, ntrials); 
            check is a check matrix for the code, ntrials = number of received vectors per code to be corrected.
            If # is given it sets the correction 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
RETURN:     nothing;  
EXAMPLE:    example solveForCode; shows an example 
"
{
 int n=ncols(check);
 int redun=nrows(check);
 int i,j;
 matrix h;
 int dist, t, tim, tim2, tim3, timdist, timdec, timdist2, timdec2, timdec3;
 ideal sys;
 list tmp;

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

 //------------------ determine error-correction capacity -------------------
 if (#>0)
 {
    t=#;
 } else {
   tim=rtimer;
   tmp=mindist(h);
   timdist=timdist+rtimer-tim;
   timdist2=timdist2+tmp[2];
   dist=tmp[1];
   printf("d= %p",dist);
   t=(dist-1) div 2; 
 }     
 tim2=rtimer;  
 tim3=rtimer;

 //------------- generate the template system ----------------------
 def A=sysQE(h,z,t,0,0);
 setring A;
 matrix word,y,rec;
 ideal sys2,sys3;
 matrix h=imap(br,h);
 matrix g=dual_code(h);  
 ideal sys=qe;
 print("The system is generated");
 timdec3=timdec3+rtimer-tim3; 

 //------------- 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); 
} example 
{
     "EXAMPLE:";  echo = 2;
     int q=32; int n=25; int redun=n-11; int t=redun+1;
     ring r=(q,a),x,dp;
     matrix check=randomCheck(redun,n,1);
     
     // correct 2 errors in using the code above, 50 trials 
     solveForCode(check,50,2); 
}


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

////////////////////////////////////////////////////////////////////////////////////
// checl 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,e); 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 is a check matrix of the code, y is a received word, t the number of errors to correct,
           e is the extension degree, s is the dimension of the point for the vanishing ideal
RETURN:    the system of Fitzgerald-Lax for the given decoding problem
EXAMPLE:   example sysFL; shows an example
"
{
     list rl=ringlist(basering);        
     int charac=rl[1][1];     
     int n=ncols(check);
     int m=nrows(check);          
     list points=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;
     
     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;     
}

///////////////////////////////////////////////////////////////////////////////
// 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("dp("+string((t-1)*(s+1)+s)+"),lp("+string(s+1)+")");
          rl[1][3][1][1]=string("lp");
          rl[1][3][1][2]=1;
          rl[1][4]=ideal(0);
     } else {
          rl[1]=p;
     }
     
     rl[2]=var_ord;
     
     rl[3]=tmp;
     rl[3][1]=tmp;     
     //rl[3][1][1]=string("dp("+string((t-1)*(s+1)+s)+"),lp("+string(s+1)+")");
     rl[3][1][1]=string("lp");
     intvec v=1;
     for (i=1; i<=size(var_ord)-1; i++)
     {
          v=v,1;
     }
     rl[3][1][2]=v; 
     rl[3][2]=tmp;
     rl[3][2][1]=string("C");    
     rl[3][2][2]=intvec(0);
     
     rl[4]=ideal(0);     
     
     def r2=ring(rl);
     setring r2;
     list l=ringlist(r2);     
     if (e>1)
     {
          execute(string("poly f="+minp));     
          ideal id=f;     
          l[1][4]=id;
     }
     
     def r=ring(l);
     setring r;      
     
     return(list(r,s));
}

///////////////////////////////////////////////////////////////////////////////
// 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 #)
{     
    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, int #)
{
  if (#==1)
  {
     list rl=ringlist(basering);     
     int charac=rl[1][1];
  } else {
     int charac=2;
  }
  list l;
  int i;
  for (i=1; i<=n; i++)
  {
    l=l+list(random(0,charac-1));
  }
  return(l);
}

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

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

 option(redSB); 
 matrix z[1][n]; 
 
 for (i=1; i<=ncodes; i++)
 {
     h=randomCheck(redun,n,e,1);            
     g=dual_code(h);
     tim2=rtimer;       
     tim3=rtimer;

     //------------ 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); 
} example
{
     "EXAMPLE:";  echo = 2;
     
     // correcting one error for one random binary code of length 25, redundancy 14; 300 words are processed
     solveForRandomFL(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;
solveForRandom(n,redun,1,1,6);

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

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

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

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

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

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

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


2.) Generic decoding is hard!

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

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

*/