Index: ngular/LIB/decode.lib
===================================================================
--- Singular/LIB/decode.lib	(revision 0c97d628cec893e481da6229e10b8066c8a40642)
+++ 	(revision )
@@ -1,1868 +1,0 @@
-///////////////////////////////////////////////////////////////////////////////
-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);
-}
