///////////////////////////////////////////////////////////////////////////////
version="$Id: decode.lib,v 1.1 2008-08-08 11:58:19 Singular Exp $";
info="
LIBRARY:   decode.lib  generating and solving systems of polynomial equations for different problems in coding theory

AUTHORS:   Stanislav Bulygin,   bulygin@mathematik.uni-kl.de,
@* Ruud Pellikaan,  g.r.pellikaan@tue.nl           

STATUS: Complete

OVERVIEW:
@* In this library we generate several systems used for decoding cyclic codes. Namely
@* SynSym and SynSymPos which represent Cooper's philosophy, Newton based on Newton identities, and Bin
@* which uses Waring's function, cf. Augot et. al. INRIA Report no. 4652, 2002. 
@* System by Pellikaan-Bulygin is generated. Also (for comparison) the system by Fitzgerald-Lax.
@* We provide also some auxiliary functions for further manipulations and decoding

MAIN PROCEDURES:
 
 SynSymPos(v,Q,n); generates SynSymPos system
 NewtonPos(v,n);  generates NewtonPos system
 Bin(v,Q,n,#);  generates Bin system
 enc(x,g);  encodes given message with a given generator matrix
 syndrome(h, c);  computes a syndroem w.r.t. a given check matrix
 PellSys(check,y,t,q,fieldeq,formal);  generates a Pellikaan system
 error(y,pos,val);  inserts errors in a word
 error_rand(y,num,e);  inserts random errors in a word
 random_check(m,n,e);  generates a random check matrix
 MDSmat(n,a);  generates an MDS (actually an RS) matrix
 mindist(check,q);   computes the minimum distance of the code
 decode(syn_full);  decoding of a word using calculated unknown syndromes 
 gauss_preprocess(work,linsyn);  substitution procedure for Gaussian elimination
 decode_preprocess(sys);  extracts values of variables after GB coputation
 solve_for_random(redun,q,ncodes,ntrials,#);  a procedure for manipulation with random codes
 solve_for_code(check,q,ntrials,#); a procedure for manipulation with the given codes
 vanish_id(points); computes the vanishing ideal for the given set of points
 FLsystem(check,y,t,e,s); generates the Fitzgerald-Lax system
 asElement(l); returns an element s of GF(2^n), represented by a vector (l[1],...,l[n]) over GF(2), n=size(l)
 FL_solve_for_random(n,redun,p,e,t,ncodes,ntrials,minpol); a procedure for manipulation with random codes with Fitzgerald-Lax
 
AUXILIARY PROCEDURES:
 determ_id(m);
 ev(f,p);
 find_index(G,p);
 ideal2list(id);
 list2ideal(l);
 divisible(m,G);
 points_gen(m,e);
 list2vec(l);
 conv_points(points);
 grasp_list(l,m,n);
 xi_gen(p,e,s); 
 gener_funcs(check,points,e,id,s); 
 FLpreprocess(p,e,n,t,minp);
 x_var(i,j);
 randomvector(n,e,#);
 random_prime_vector(n,#);
 LF_add_synd(rec,check,sys);
 lis(n);  
 combinations(m,n);
 combinsert(temp,i);
 combinations_sum(m,n);
 exp_count(n,q);
 star(m,i,j);
 add_synd(rec,check,redun,sys);
 list2intvec(l);
KEYWORDS:  cyclic code, linear code, decoding, Groebner bases
";

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

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

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

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

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 CRHT (int n, list defset, int e, int q, int m, int #)
"USAGE:  CRHT(n,defset,e,q,m,#);  n length of the cyclic code, defset is a list representing the defining set,
@ e the error-correcting capacity, m degree extension of the splitting field, if #>0 additional equations
@ representing the fact that every two error positions are either different or at least one of them is zero
RETURN:  a ring to work with the CRHT-ideal (with Sala's additions), the ideal itself is exported with the name "crht"
EXAMPLE: example CRHT; shows an example
{
  int r=size(defset);  
  ring @crht=(q,a),(Y(e..1),Z(1..e),X(r..1)),lp;
  ideal crht;
  int i,j;
  poly sum;
  
  // check equations
  for (i=1; i<=r; i++)
  {
    sum=0;
    for (j=1; j<=e; j++)
    {
      sum=sum+Y(j)*Z(j)^defset[i];
    }
    crht[i]=sum-X(i);
  }
  
  // restrictions on syndromes
  for (i=1; i<=r; i++)
  {
    crht=crht,X(i)^(q^m)-X(i);
  }
  
  // n-th roots of unity
  for (i=1; i<=e; i++)
  {
    crht=crht,Z(i)^(n+1)-Z(i);
  }
  
  for (i=1; i<=e; i++)
  {
    crht=crht,Y(i)^(q-1)-1;
  }  
  
  if (#)
  {
    for (i=1; i<=e; i++)
    {
      for (j=i+1; j<=e; j++)
      {
        crht=crht,Z(i)*Z(j)*p_poly(n,i,j);
      }
    }
  }  
  export crht;  
  return(@crht);    
} example
{
  "EXAMPLE:"; echo=2;
  // binary cyclic [15,7,5] code with defining set (1,3)
  
  list defset=1,3; // defining set
  
  int n=15; // length
  int e=2; // error-correcting capacity
  int q=2; // basefield size
  int m=4; // degree extension of the splitting field
  int sala=1; // indicator to add additional equations
    
  def A=CRHT(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=CRHT(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 CRHT_mindist_binary (int n, list defset, int w)
"USAGE:  CRHT(n,defset,w);  n length of the cyclic code, defset is a list representing the defining set,
@ w is a candidate for the minimum distance
RETURN:  a ring to work with the Sala's ideal for mindist, the ideal itself is exported with the name "crht_md"
EXAMPLE: example CRHT_mindist_binary; shows an example
{
  int r=size(defset);  
  ring @crht_md=2,Z(1..w),lp;
  ideal crht_md;
  int i,j;
  poly sum;
  
  // check equations
  for (i=1; i<=r; i++)
  {
    sum=0;
    for (j=1; j<=w; j++)
    {
      sum=sum+Z(j)^defset[i];
    }
    crht_md[i]=sum;
  }  
  
  
  // n-th roots of unity
  for (i=1; i<=w; i++)
  {
    crht_md=crht_md,Z(i)^n-1;
  }  
  
  
  for (i=1; i<=w; i++)
  {
    for (j=i+1; j<=w; j++)
    {
      crht_md=crht_md,Z(i)*Z(j)*p_poly(n,i,j);
    }
  }
    
  export crht_md;  
  return(@crht_md);    
} example
{
  "EXAMPLE:"; echo=2;
  // binary cyclic [15,7,5] code with defining set (1,3)
  
  list defset=1,3; // defining set
  
  int n=15; // length
  int d=5; // candidate for the minimum distance  
      
  def A=CRHT_mindist_binary(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);  
}

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

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

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

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 Bin (int v, list Q, int n, int#)
"USAGE:  Bin (v, Q, n, #);  v a number if errors, Q is a generating set of the code, n the length, is additional parameter is 
@* set to 1, then the generating set is enlarged by odd elements, which are 2^(some power)*(some elment in the gen.set) mod n
RETURN:  keeps a ring with the resulting system
EXAMPLE: example Bin; shows an example
{
 //ring r=2,(sigma(1..v),S(1..n)),(lp(v),dp(n));
 ring r=2,(sigma(1..v),S(1..n)),dp;
 list cyclot;
 ideal result;
 int i,j,k,s;
 list comb;
 poly sum_, mon;
 int count1, count2, upper, coef_, flag, gener;
 list Q_update;
 if (#==1)
 {
  for (i=1; i<=n; i++)
  {
   cyclot[i]=0;
  }
  for (i=1; i<=size(Q); i++)
  {
   flag=1;
   gener=Q[i];
   while(flag)
   {
    cyclot[gener]=1;
    gener=2*gener mod n;
    if (gener == Q[i]) {flag=0;}
   }
  }
  for (i=1; i<=n; i++)
  {
   if ((cyclot[i] == 1)&&(i mod 2 == 1)) {Q_update[size(Q_update)+1]=i;}
  }
 }
 else
 {
  Q_update=Q;
 }
   
 for (i=1; i<=size(Q_update); i++)
 {
  comb=combinations_sum(v,Q_update[i]);
  sum_=0;
  for (k=1; k<=size(comb); k++)
  {
   upper=0;
   for (j=1; j<=v; j++)
   {
    upper=upper+comb[k][j];
   }
   count1=0;
   for (j=2; j<=upper-1; j++) 
   {
    count1=count1+exp_count(j,2);
   }
   count1=count1+exp_count(Q_update[i],2);
   count2=0;
   for (j=1; j<=v; j++)
   {
    for (s=2; s<=comb[k][j]; s++)
    {
     count2=count2+exp_count(s,2);
    }
   }
   if (count1<count2) {print("ERRORBin");}
   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_;  
 }
 export(result);
 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
     Bin(3,l,31);
     print(result);
}

proc enc (matrix x, matrix g)
"USAGE:  enc (x, g);  x a row vector (message), and g a generator matrix
RETURN:  corresponding codeword
EXAMPLE: example enc; shows an example
{
 if (nrows(x)>1) {print("ERRORenc1!");}
 if (ncols(x)!=nrows(g)) {print("ERRORenc2!");}
 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(enc(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=enc(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));
}

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 PellSys(matrix check, matrix y, int t, int q_orig, int fieldeq, int formal)
"USAGE:  PellSys(check, y, t, q, fieldeq, formal); check is the check matrix of the code   
@* y is a received word, t the number of errors to be corrected, q the size of an (extended) alphabet,
@* if fieldeq=1, then field equations are added; if formal=0, fields equations on (known) syndrome variables (they have to be defined previously)
@* 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
ASSUME:  working ring should already be defined
RETURN:  the resulting system
EXAMPLE: example PellSys; shows an example
{
 int red=nrows(check);
 int n=ncols(check);
 
 matrix h_full=MDSmat(n,a);
 matrix h=submat(h_full,1..red,1..n);
 if (nrows(y)!=1) {print("ERROR1Pell");}
 if (ncols(y)!=n) {print("ERROR2Pell");}
  
 ideal result;
 
 list c;
 list a;
 list tmp,tmp2;
 int i,j,l,k;
 number sum,prod,sig;
        poly sum1,sum2,sum3;
 for (i=1; i<=n; i++)
 {
  c[i]=tmp;
 }
 
 int tim=rtimer;
 matrix transf=inverse(transpose(h_full));
  
 tim=rtimer;
 for (i=1; i<=red ; i++)
 {
  a[i]=transpose(submat(check,i..i,1..n));
  a[i]=transf*a[i];
 }
 
 tim=rtimer;
 matrix te[n][1]; 
 for (i=1; i<=n; i++)
 {
  for (j=1; j<=t+1; j++)
  {
   if ((j<i)&&(i<=t+1)) {c[i][j]=c[j][i];}
   else 
   {
    if (i+j<=n+1) 
    {
     c[i][j]=te;
     c[i][j][i+j-1,1]=1; 
    }
    else
    {
     c[i][j]=star(h_full,i,j);
     c[i][j]=transf*c[i][j];
    }   
   }
  }
 }
 
 
 tim=rtimer; 
 if (formal==0)
 {
  matrix s[red][1]=syndrome(check,y);
  for (j=1; j<=red; j++)
  {
   sum1=0;
   for (l=1; l<=n; l++)
   {
    sum1=sum1+a[j][l,1]*U(l);
   }
   result=result,sum1-s[j,1];
  }
 } else
 {
  for (j=1; j<=red; j++)
  {
   sum1=0;
   for (l=1; l<=n; l++)
   {
    sum1=sum1+a[j][l,1]*U(l);
   }
   result=result,sum1-s(j);
  }
  for (j=1; j<=red; j++)
  {
     result=result,s(j)^(formal)-s(j);
  }
 }
 if (fieldeq)
 {
  for (i=1; i<=n; i++)
  {
   result=result,U(i)^q-U(i);
  }
  for (j=1; j<=t; j++)
  {
     result=result,V(j)^q-V(j);
  }
 } 
 for (i=1; i<=n; i++)
 {
  sum1=0;
  for (j=1; j<=t; j++)
  {
   sum2=0;
   for (l=1; l<=n; l++)
   {
    sum2=sum2+c[i][j][l,1]*U(l);
   }
   sum1=sum1+sum2*V(j);
  }
  sum3=0;
  for (l=1; l<=n; l++)
  {
   sum3=sum3+c[i][t+1][l,1]*U(l);
  }
  result=result,sum1-sum3;
 } 
 
 result=simplify(result,2);
 exportto(Top,h_full);
 return(result);   
} example 
{
     "EXAMPLE:";  echo = 2;
     LIB "brnoeth.lib";
     //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),(V(1..t),U(n..1)),dp;
     matrix h_full=MDSmat(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]=enc(x,g);
     //disturb with 2 errors
     matrix rec[1][7]=error(y,list(2,4),list(1,a));
     //generate the system
     ideal sys=PellSys(h,rec,t,q,0,0);
     print(sys);
     //let us decode
     option(redSB);
     ideal sys_red=std(sys);
     print(sys_red);     
}

proc error(matrix y, list pos, list val)
"USAGE:  error(y, pos, val); y is a (code) word, pos = positions where errors occured, val = their corresponding values
RETURN:  corresponding received word
EXAMPLE: example error; shows an example
{
 matrix result[1][ncols(y)]=y;
 if (size(pos)!=size(val)) {print("ERRORerror");}
 for (int i=1; i<=size(pos); i++)
 {
  result[1,pos[i]]=y[1,pos[i]]+val[i];
 }
 return(result);
} example 
{
     "EXAMPLE:";  echo = 2;
     LIB "brnoeth.lib";
     //correct 2 errors in [7,3] 8-ary RS code
     int t=2; int q=8; int n=7; int redun=4;
     ring r=(q,a),(V(1..t),U(n..1)),dp;
     matrix h_full=MDSmat(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]=enc(x,g);
     
     //disturb with 2 errors
     matrix rec[1][7]=error(y,list(2,4),list(1,a));
     
     //generate the system
     ideal sys=PellSys(h,rec,t,q,0,0);
     print(sys);
     //let us decode
     option(redSB);
     ideal sys_red=std(sys);
     print(sys_red);     
}

proc error_rand(matrix y, int num, int e)
"USAGE:  error_rand(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 error_rand; 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;
     LIB "brnoeth.lib";
     //correct 2 errors in [7,3] 8-ary RS code
     int t=2; int q=8; int n=7; int redun=4;
     ring r=(q,a),(V(1..t),U(n..1)),dp;
     matrix h_full=MDSmat(n,a);
     matrix h=submat(h_full,1..redun,1..n);
     print(h);
     matrix g=dual_code(h);
     matrix x[1][n-redun]=0,0,1,0;
     matrix y[1][n]=enc(x,g);
     print(y);
     
     //disturb with 2 random errors
     matrix rec[1][7]=error_rand(y,2,3);
     print(rec);     
     
     //generate the system
     ideal sys=PellSys(h,rec,t,q,0,0);
     print(sys);
     //let us decode
     option(redSB);
     ideal sys_red=std(sys);
     print(sys_red);     
}

proc random_check(int m, int n, int e, int #)
"USAGE:  random_check(m, n, e); mxn 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 random_check; 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;
     LIB "brnoeth.lib";
     //correct 1 error in [7,3] binary code 
     int t=1; int q=8; int n=7; int redun=4;
     ring r=(q,a),(V(1..t),U(n..1)),dp;
          
     //generate random check matrix
     matrix h=random_check(redun,n,1);
     print(h);
     
     matrix g=dual_code(h);
     matrix x[1][n-redun]=0,0,1,0;
     matrix y[1][n]=enc(x,g);
     //disturb with 1 random error
     matrix rec[1][7]=error_rand(y,1,1);     
     //generate the system
     ideal sys=PellSys(h,rec,t,q,0,0);
     print(sys);
     //let us decode
     option(redSB);
     ideal sys_red=std(sys);
     print(sys_red);     
}

proc MDSmat(int n, number a)
"USAGE:  MDSmat(n, a); nxn are dimensions of the matrix, a is a primitive element of the field
RETURN:  a matrix with the MDS property
EXAMPLE: example random_check; 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;
     LIB "brnoeth.lib";
     //correct 1 error in [7,3] binary code 
     int t=1; int q=8; int n=7; int redun=4;
     ring r=(q,a),(V(1..t),U(n..1)),dp;
     
     //generate an MDS (RS) matrix
     matrix h_full=MDSmat(n,a);
     print(h_full);
          
     //generate random check matrix
     matrix h=random_check(redun,n,1);     
     matrix g=dual_code(h);
     matrix x[1][n-redun]=0,0,1,0;
     matrix y[1][n]=enc(x,g);
     //disturb with 1 random error
     matrix rec[1][n]=error_rand(y,1,1);     
     //generate the system
     ideal sys=PellSys(h,rec,t,q,0,0);
     print(sys);
     //let us decode
     option(redSB);
     ideal sys_red=std(sys);
     print(sys_red);     
}


proc mindist (matrix check, int q)
"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;
 
 ideal temp;
 int count=1;
 int flag=1;
 int flag2;
 int i, tim, timsolve;
 matrix z[1][n]; 
 option(redSB);
 while (flag)
 {
  temp=PellSys(check,z,count,q,0,0);  
  tim=rtimer;
  temp=std(temp);
  timsolve=timsolve+rtimer-tim;
  
  flag2=1;
  for (i=1; i<=ncols(check); i++)
  {
   if (temp[i]!=U(i)) {flag2=0;}
  }
  if (!flag2) {flag=0;}
  else {count++;}
 }
 return(list(count,timsolve)); 
} 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),(V(1..t),U(n..1)),dp;     
               
     //generate random check matrix
     matrix h=random_check(redun,n,1);
     print(h);
     list l=mindist(h,q);
     print(l[1]);  
     //time for the comutation in secs   
     print(l[2]);     
}

proc decode(matrix syn_full)
"USAGE:  decode(h_full, syn_full); syn_full is a full syndrome
RETURN:  corresponding codeword
EXAMPLE: example decode; shows an example 
{
 int n=nrows(syn_full);
 matrix h_full=MDSmat(n,a);
 return(inverse(h_full)*syn_full); 
} example 
{
     "EXAMPLE:";  echo = 2;
     LIB "brnoeth.lib";
     //correct 1 error in [7,3] binary code 
     int t=1; int q=8; int n=7; int redun=4;
     ring r=(q,a),(V(1..t),U(n..1)),dp;     
             
     //generate random check matrix
     matrix h=random_check(redun,n,1);     
     matrix g=dual_code(h);
     matrix x[1][n-redun]=0,0,1,0;
     matrix y[1][n]=enc(x,g);
     print(y);
     //disturb with 1 random error
     matrix rec[1][n]=error_rand(y,1,1);
     print(rec);     
     //generate the system
     ideal sys=PellSys(h,rec,t,q,0,0);
     print(sys);
     //let us decode
     option(redSB);
     ideal sys_red=std(sys);
     print(sys_red);
     
     //obtain a codeword
     //this works only if our code is indeed can correct one error
     matrix syn=decode_preprocess(ideal(sys_red[1..7]));
     print(syn);
     print(rec-transpose(decode(syn)));     
}

proc gauss_preprocess(ideal work, matrix linsyn)
"USAGE:  gauss_preprocess(work, linsyn); work is an ideal, linsyn is a matrix in row-reduced form
@* We then express variables from linsyn and substitute them in work
RETURN:  an ideal with substituted variables
NOTE:  becuase of use of subst, it might take more time, than one could expect
EXAMPLE: example gauss_preprocess; shows an example 
{
 int m=size(work);
 int i,j,k;
 poly temp;
 list substit;
 list substvar;
 for (i=1; i<=nrows(linsyn); i++)
 { 
  temp=linsyn[i,1];
  j=ncols(linsyn);
  while(linsyn[i,j]==0)
  {
   j--;
  }
  for (k=2; k<=j-1; k++)
  {
   temp=temp+linsyn[i,k]*U(k-1);
  }
  substit[i]=-temp;
  substvar[i]=U(j-1);
 }
 for (i=1; i<=m; i++)
 {
  for (j=1; j<=size(substit); j++)
  {
   work[i]=subst(work[i],substvar[j],substit[j]);
  }
 }
 export(substvar);
 return(work);
}

proc decode_preprocess(ideal sys)
"USAGE:  decode_preprocess(sys); sys is a reduced Groebner basis of PellSys system
ASSUME: it is of the form: (var - val.of.var.);
RETURN:  correpsponding values
EXAMPLE: example decode_preprocess; shows an example 
{
 matrix result[size(sys)][1];
 for (int i=1; i<=size(sys); i++)
 {
  result[i,1]=-sys[i]+lead(sys[i]);
 }
 return(result);
} example 
{
     "EXAMPLE:";  echo = 2;
     LIB "brnoeth.lib";
     //correct 1 error in [7,3] binary code 
     int t=1; int q=8; int n=7; int redun=4;
     ring r=(q,a),(V(1..t),U(n..1)),dp;     
             
     //generate random check matrix
     matrix h=random_check(redun,n,1);     
     matrix g=dual_code(h);
     matrix x[1][n-redun]=0,0,1,0;
     matrix y[1][n]=enc(x,g);
     print(y);
     //disturb with 1 random error
     matrix rec[1][n]=error_rand(y,1,1);     
     //generate the system
     ideal sys=PellSys(h,rec,t,q,0,0);
     print(sys);
     //let us decode
     option(redSB);
     ideal sys_red=std(sys);
     print(sys_red);
     
     //obtain a codeword
     //this works only if our code is indeed can correct one error
     matrix syn=decode_preprocess(ideal(sys_red[1..7]));
     print(syn);
     print(rec-transpose(decode(syn)));     
}

proc solve_for_random(int n, int redun, int q, int ncodes, int ntrials, int #)
"USAGE:  solve_for_random(redun, q, ncodes, ntrials); redun is a redundabcy of the (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 min.dist explicitly. It should be used in case one expexts some lower bound
RETURN:  
EXAMPLE: example solve_for_random; shows an example 
{
 int i,j;
 matrix h, g, word, y, rec;
 int dist, t, tim, tim2, tim3, timdist, timdec, timdist2, timdec2, timdec3;
 ideal sys, sys2, sys3;
 list tmp;

 option(redSB);
 matrix h_full=MDSmat(n,a); 
 matrix z[1][ncols(h_full)];
 int n=ncols(h_full);
 for (i=1; i<=ncodes; i++)
 {
  h=random_check(redun,n,1);
  if (#>0)
  {
     dist=#;
  } else {
     tim=rtimer;
     tmp=mindist(h,q);
     timdist=timdist+rtimer-tim;
     timdist2=timdist2+tmp[2];
     dist=tmp[1];
     printf("d= %p",dist);
  }  
  t=(dist-1) div 2;  
  tim2=rtimer;
  g=dual_code(h);  
  tim3=rtimer;
  sys=PellSys(h,z,t,q,0,0);
  print("system is generated");
  timdec3=timdec3+rtimer-tim3;  
  for (j=1; j<=ntrials; j++)
  {
   word=randomvector(n-redun,1);
   y=enc(transpose(word),g);
   rec=error_rand(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 PellSys in decoding: %p", timdec3);
 export(sys2);
 export(sys3);
 return("Done"); 
} example 
{
     "EXAMPLE:";  echo = 2;
     int q=32; int n=25; int redun=n-11; int t=redun+1;
     ring r=(q,a),(V(1..t),U(n..1)),dp;
     
     solve_for_random(n,redun,q,5,50); 
}

proc solve_for_code(matrix check, int q, int ntrials, int #)
"USAGE:  solve_for_code(check, q, ntrials); 
@* check is a check matrix for the code, q = field size, 
@* ntrials = number of received vectors per code to be corrected
@* if # is given it sets min.dist explicitly. It should be used in case one expexts some lower bound
RETURN:  
EXAMPLE: example solve_for_code; shows an example 
{
 int redun=nrows(check);
 int i,j;
 matrix h, g, word, y, rec;
 int dist, t, tim, tim2, tim3, timdist, timdec, timdist2, timdec2, timdec3;
 ideal sys, sys2, sys3;
 list tmp;

 option(redSB); 
 int n=ncols(check);
 matrix h_full=MDSmat(n,a);
 matrix z[1][ncols(h_full)];
 
  h=check;
  if (#>0)
  {
     dist=#;
  } else {
     tim=rtimer;
     tmp=mindist(h,q);
     timdist=timdist+rtimer-tim;
     timdist2=timdist2+tmp[2];
     dist=tmp[1];
     printf("d= %p",dist);
  }  
  t=(dist-1) div 2;  
  tim2=rtimer;
  g=dual_code(h);  
  tim3=rtimer;
  sys=PellSys(h,z,t,q,0,0);
  print("system is generated");
  timdec3=timdec3+rtimer-tim3;  
  for (j=1; j<=ntrials; j++)
  {
   word=randomvector(n-redun,1);
   y=enc(transpose(word),g);
   rec=error_rand(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 PellSys in decoding: %p", timdec3);
 export(sys2);
 export(sys3);
 return("Done"); 
} example 
{
     "EXAMPLE:";  echo = 2;
     int q=32; int n=25; int redun=n-11; int t=redun+1;
     ring r=(q,a),(V(1..t),U(n..1)),dp;
     
     //generate random check matrix
     matrix check=random_check(redun,n,1);
     print(check);
     
     solve_for_code(check,q,50); 
}

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

proc determ_id (matrix m)
/*"USAGE:  determ_id (m); m a matrix with variable entries
RETURN:  an ideal I(t,U,V) with t=ncols(m)-1*/
{
 int t=ncols(m)-1;
 int n=nrows(m);
 list comb1, comb2;
 comb1=combinations(t,n);
 comb2=combinations(t,t+1);
 int i,j;
 ideal result; 
 for (i=1; i<=size(comb1); i++)
 {
      for (j=1; j<=size(comb2); j++)
      {
          result=result,det(submat(m,list2intvec(comb1[i]),list2intvec(comb2[j])));
      }
 }
 return(result);    
} 

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

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

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

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

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

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

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

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

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

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

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

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 FLsystem (matrix check, matrix y, int t, int e, int s)
"USAGE:  FLsystem (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 FLsystem; shows an example
{
     list rl=ringlist(basering);        
     int charac=rl[1][1];     
     int n=ncols(check);
     int m=nrows(check);          
     list points=points_gen(s,e);
     list points2=conv_points(points);
     list p=grasp_list(points2,1,n);     
     ideal id=vanish_id(p,e);      
     ideal funcs=gener_funcs(check,p,e,id,s);          
     
     ideal result;
     poly temp;
     int i,j,k;
     
     //vanishing realtions     
     for (i=1; i<=t; i++)
     {          
          for (j=1; j<=size(id); j++)
          {
               temp=id[j];               
               for (k=1; k<=s; k++)
               {
                    temp=subst(temp,x(k),x_var(i,k,s));
               }               
               result=result,temp;
          }
     }          
     
     //field equations
     for (i=1; i<=t; i++)
     {
          for (k=1; k<=s; k++)
          {
               result=result,x_var(i,k,s)^(charac^e)-x_var(i,k,s);
          }
     }
     for (i=1; i<=t; i++)
     {
          result=result,e(i)^(charac^e-1)-1;
     }
     
     result=simplify(result,8);
     
     //check realtions
     poly sum;
     matrix syn[m][1]=syndrome(check,y);          
     for (i=1; i<=size(funcs); i++)
     {          
          sum=0;
          for (j=1; j<=t; j++)
          {
               temp=funcs[i];
               for (k=1; k<=s; k++)
               {
                    temp=subst(temp,x(k),x_var(j,k,s));
               }
               sum=sum+temp*e(j);         
          }
          result=result,sum-syn[i,1];
     }
     
     result=simplify(result,2);
     
     points=points2;
     export points;     
     return(result);
} example
{
     "EXAMPLE:";  echo = 2;
     
     list l=FLpreprocess(3,1,11,2,"");
     def r=l[1];
     setring r;
     int s_work=l[2];     
     
     //the check matrix of [11,6,5] ternary code
     matrix h[5][11]=1,0,0,0,0,1,1,1,-1,-1,0,
          0,1,0,0,0,1,1,-1,1,0,-1,
          0,0,1,0,0,1,-1,1,0,1,-1,
          0,0,0,1,0,1,-1,0,1,-1,1,
          0,0,0,0,1,1,0,-1,-1,1,1;
     matrix g=dual_code(h);
     matrix x[1][6];
     matrix y[1][11]=enc(x,g);
     //disturb with 2 errors
     matrix rec[1][11]=error(y,list(2,4),list(1,-1));

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

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

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

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

proc asElement(list l)
"USAGE:  asElement(l);  l a list
ASSUME:  the extension ring is defined with the explicit minpoly
RETURN:  number
NOTE:    returns an element s of GF(2^n), represented by a vector (l[1],...,l[n]) over GF(2)
EXAMPLE: example asElement; shows an example
"{
  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);
}
example
{
  "EXAMPLE:";  echo = 2;
  ring r  = (2,a),x(1..5),dp;
  minpoly = a5+a4+a3+a+1;
  asElement(list(0,1,1,0,0));
}

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 FL_solve_for_random(int n, int redun, int p, int e, int t, int ncodes, int ntrials, string minpol)
"USAGE:  FL_solve_for_random(redun,p,e,n,t,ncodes,ntrials,minpol); n = length of codes generated, redun = redundancy of codes generated, 
@* p = characteristics, e is the extension degree, 
@* q = number of errors to correct, ncodes = number of random codes to be processed
@* ntrials = number of received vectors per code to be corrected
@* due to some pecularities of SINGULAR one needs to provide minimal polynomial for the extension explicitly
RETURN:  the system of Fitzgerald-Lax for the given decoding problem
EXAMPLE: example FLsystem; 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=random_check(redun,n,e,1);            
     g=dual_code(h);
     tim2=rtimer;       
     tim3=rtimer;               
     sys=FLsystem(h,z,t,e,s_work);     
     timdec3=timdec3+rtimer-tim3;
   
     for (j=1; j<=ntrials; j++)
     {
          word=randomvector(n-redun,1);
          y=enc(transpose(word),g);
          rec=error_rand(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 LF in decoding: %p", timdec3);
 
 return("Done"); 
} example
{
     "EXAMPLE:";  echo = 2;
     
     // decoding for one random binary code of length 25, redundancy 14; 300 words are processed
     FL_solve_for_random(25,14,2,1,1,1,300,"");
}

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