Context Navigation

← Previous Change
Next Change →

decodegb.lib

Timestamp:

Jan 14, 2009, 5:07:05 PM (15 years ago)

Author:

Hans Schönemann <hannes@…>

Branches:

(u'spielwiese', 'd0474371d8c5d8068ab70bfb42719c97936b18a6')

Children:

0721816437af5ddabc83aa203a12d9b58b42a33c

Parents:

95edd5641377e851d4a1d4e986853687991d0895

Message:

*hannes: format


git-svn-id: file:///usr/local/Singular/svn/trunk@11306 2c84dea3-7e68-4137-9b89-c4e89433aadc

File:

: 1 edited

Singular/LIB/decodegb.lib (modified) (113 diffs)

Legend:

: Unmodified
: Added
: Removed

Singular/LIB/decodegb.lib

-                      r95edd5
+                      r7f3ad4
 ///////////////////////////////////////////////////////////////////////////////
 version="$Id: decodegb.lib,v 1.8 2008-10-22 13:31:48 bulygin Exp $";
+version="$Id: decodegb.lib,v 1.9 2009-01-14 16:07:03 Singular Exp $";
 category="Coding theory";
 info="
 LIBRARY: decodegb.lib         Decoding and min distance of linear codes with GB
 AUTHOR:  Stanislav Bulygin,   bulygin@mathematik.uni-kl.de
+AUTHOR:  Stanislav Bulygin,   bulygin@mathematik.uni-kl.de
 OVERVIEW:
  In this library we generate several systems used for decoding cyclic codes and
  finding their minimum distance. Namely, we work with the Cooper's philosophy
  and generalized Newton identities. The original method of quadratic equations
  is worked out here as well. We also (for comparison) enable to work with the
  system of Fitzgerald-Lax. We provide some auxiliary functions for further
  manipulations and decoding. For an overview of the methods mentioned above @ref{Decoding codes with GB}
  For the vanishing ideal computation the algorithm of Farr and Gao is
+ finding their minimum distance. Namely, we work with the Cooper's philosophy
+ and generalized Newton identities. The original method of quadratic equations
+ is worked out here as well. We also (for comparison) enable to work with the
+ system of Fitzgerald-Lax. We provide some auxiliary functions for further
+ manipulations and decoding. For an overview of the methods mentioned above @ref{Decoding codes with GB}
+ For the vanishing ideal computation the algorithm of Farr and Gao is
  implemented.
 MAIN PROCEDURES:
+MAIN PROCEDURES:
  sysCRHT(..);        generates the CRHT-ideal as in Cooper's philosophy
  sysCRHTMindist(..); CRHT-ideal to find the minimum distance in the binary case
 …
  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
+ decode(rec);        decoding of a word using the system of quadratic equations
  decodeRandom(..); a procedure for manipulation with random codes
  decodeCode(..);   a procedure for manipulation with the given code
 …
  sysFL(..);          generates the Fitzgerald-Lax system
  decodeRandomFL(..); manipulation with random codes via Fitzgerald-Lax
 KEYWORDS:  Cyclic code; Linear code; Decoding;
+KEYWORDS:  Cyclic code; Linear code; Decoding;
            Minimum distance; Groebner bases
 ";
 …
 ///////////////////////////////////////////////////////////////////////////////
 // the polynomial for Sala's restrictions
+// the polynomial for Sala's restrictions
 static proc p_poly(int n, int a, int b)
+{
 …
 "USAGE:   sysCRHT(n,defset,e,q,m,[k]); n,e,q,m,k int, defset list of int's
 @format
          - n length of the cyclic code,
+         - n length of the cyclic code,
          - defset is a list representing the defining set,
          - e the error-correcting capacity,
+         - e the error-correcting capacity,
          - q field size
          - m degree extension of the splitting field,
          - if k>0 additional equations representing the fact that every two
+         - m degree extension of the splitting field,
+         - if k>0 additional equations representing the fact that every two
          error positions are either different or at least one of them is zero
 @end format
 RETURN: the ring to work with the CRHT-ideal (with Sala's additions),
+RETURN: the ring to work with the CRHT-ideal (with Sala's additions),
         containig an ideal with name 'crht'
 THEORY:  Based on 'defset' of the given cyclic code, the procedure constructs
          the corresponding Cooper-Reed-Heleseth-Truong ideal 'crht'. With its
+THEORY:  Based on 'defset' of the given cyclic code, the procedure constructs
+         the corresponding Cooper-Reed-Heleseth-Truong ideal 'crht'. With its
          help one can solve the decoding problem. For basics of the method @ref{Cooper philosophy}.
 SEE ALSO: sysNewton, sysBin
 …
+"
+{
   int r=size(defset);
+  int r=size(defset);
   ring @crht=(q,a),(Y(e..1),Z(1..e),X(r..1)),lp;
   ideal crht;
 …
   poly sum;
   int k;
   if ( size(#) > 0)
+  {
     k = #[1];
+  }
+  if ( size(#) > 0)
+  {
+    k = #[1];
+  }
   //---------------------- add check equations --------------------------
   for (i=1; i<=r; i++)
 …
     crht[i]=sum-X(i);
+  }
   //--------------------- field equations on syndromes ------------------
   for (i=1; i<=r; i++)
 …
     crht=crht,X(i)^(q^m)-X(i);
+  }
   //------ restrictions on error-locations: n-th roots of unity ----------
   for (i=1; i<=e; i++)
 …
     crht=crht,Z(i)^(n+1)-Z(i);
+  }
   for (i=1; i<=e; i++)
+  {
     crht=crht,Y(i)^(q-1)-1;
+  }
+  }
   //--------- add Sala's additional conditions if necessary --------------
   if ( k > 0 )
+  {
     for (i=1; i<=e; i++)
 …
+      }
+    }
+  }
   export crht;
   return(@crht);
+}
+  }
+  export crht;
+  return(@crht);
+}
 example
 { "EXAMPLE:"; echo=2;
 …
   intvec v = option(get);
   list defset=1,3;           // defining set
+  list defset=1,3;           // defining set
   int n=15;                  // length
   int e=2;                   // error-correcting capacity
 …
   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);
+  def A=sysCRHT(n,defset,e,q,m);
   setring A;
   A;                         // shows the ring we are working in
   print(crht);               // the CRHT-ideal
+  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;
+  A=sysCRHT(n,defset,e,q,m,sala);
+  setring A;
   print(crht);                // CRHT-ideal with additional equations from Sala
   option(redSB);
   ideal red_crht=std(crht);   // reduced Groebner basis
   print(red_crht);
+  print(red_crht);
   red_crht[5];                // general error-locator polynomial for this code
   option(set,v);
 …
 proc sysCRHTMindist (int n, list defset, int w)
 "USAGE:  sysCRHTMindist(n,defset,w);  n,w are int, defset is list of int's
+"USAGE:  sysCRHTMindist(n,defset,w);  n,w are int, defset is list of int's
 @format
         - n length of the cyclic code,
+        - n length of the cyclic code,
         - defset is a list representing the defining set,
         - w is a candidate for the minimum distance
 @end format
 RETURN:  the ring to work with the Sala's ideal for the minimum distance
+RETURN:  the ring to work with the Sala's ideal for the minimum distance
          containing the ideal with name 'crht_md'
 THEORY:  Based on 'defset' of the given cyclic code, the procedure constructs
          the corresponding Cooper-Reed-Heleseth-Truong ideal 'crht_md'. With
          its help one can find minimum distance of the code in the binary
+THEORY:  Based on 'defset' of the given cyclic code, the procedure constructs
+         the corresponding Cooper-Reed-Heleseth-Truong ideal 'crht_md'. With
+         its help one can find minimum distance of the code in the binary
          case. For basics of the method @ref{Cooper philosophy}.
 EXAMPLE: example sysCRHTMindist; shows an example
+"
+{
   int r=size(defset);
+  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++)
 …
+    }
     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++)
 …
+    }
+  }
   export crht_md;
   return(@crht_md);
+}
+  export crht_md;
+  return(@crht_md);
+}
 example
+{
 …
   intvec v = option(get);
   // binary cyclic [15,7,5] code with defining set (1,3)
   list defset=1,3;             // defining set
   int n=15;                    // length
   int d=5;                     // candidate for the minimum distance
   def A=sysCRHTMindist(n,defset,d);
+  list defset=1,3;             // defining set
+  int n=15;                    // length
+  int d=5;                     // candidate for the minimum distance
+  def A=sysCRHTMindist(n,defset,d);
   setring A;
   A;                           // shows the ring we are working in
   print(crht_md);              // the Sala's ideal for mindist
   option(redSB);
   ideal red_crht_md=std(crht_md);
+  ideal red_crht_md=std(crht_md);
   print(red_crht_md);          // reduced Groebner basis
   option(set,v);
+}
 …
 proc sysNewton (int n, list defset, int t, int q, int m, list #)
 "USAGE:   sysNewton (n,defset,t,q,m,[tr]); n,t,q,m,tr int, defset list int's
+"USAGE:   sysNewton (n,defset,t,q,m,[tr]); n,t,q,m,tr int, defset list int's
 @format
          - n is length,
+         - n is length,
          - defset is the defining set,
          - t is the number of errors,
          - q is basefield size,
+         - t is the number of errors,
+         - q is basefield size,
          - m is degree extension of the splitting field,
          - if tr>0 it indicates that Newton identities in triangular
+         - if tr>0 it indicates that Newton identities in triangular
            form should be constructed
 @end format
 RETURN:  the ring to work with the generalized Newton identities (in
+RETURN:  the ring to work with the generalized Newton identities (in
          triangular form if applicable) containing the ideal with name 'newton'
 THEORY:  Based on 'defset' of the given cyclic code, the procedure constructs
          the corresponding ideal 'newton' with the generalized Newton
          identities. With its help one can solve the decoding problem. For
+THEORY:  Based on 'defset' of the given cyclic code, the procedure constructs
+         the corresponding ideal 'newton' with the generalized Newton
+         identities. With its help one can solve the decoding problem. For
          basics of the method @ref{Cooper philosophy}.
 SEE ALSO: sysCRHT, sysBin
 EXAMPLE:  example sysNewton; shows an example
+"
+{
+{
  string s="ring @newton=("+string(q)+",a),(";
  int i,j;
  int flag;
  int tr;
  if (size(#)>0)
+ {
   tr=#[1];
+ }
  for (i=n; i>=1; i--)
+ {
 …
+    {
       flag=0;
       break;
+      break;
+    }
+  }
 …
   s=s+"S("+string(defset[i])+"),";
+ }
  s=s+"S("+string(defset[1])+")),lp;";
+ s=s+"S("+string(defset[1])+")),lp;";
  execute(s);
  ideal newton;
+ ideal newton;
  poly sum;
  //------------ generate generalized Newton identities ----------
  if (tr)
 …
+  }
  } else
+ {
+ {
   for (i=1; i<=t; i++)
+  {
 …
+  }
   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=simplify(newton,2);
  export newton;
  return(@newton);
+}
 example
+ return(@newton);
+}
+example
+{
      "EXAMPLE:";  echo = 2;
      // Newton identities for a binary 3-error-correcting cyclic code of
+     // 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 m=5;           // degree extension of the splitting field
      int tr=1;          // indicator of triangular form of Newton identities
      def A=sysNewton(n,defset,t,q,m);
      setring A;
      A;                 // shows the ring we are working in
      print(newton);     // generalized Newton identities
      //===============================
      A=sysNewton(n,defset,t,q,m,tr);
      setring A;
+     setring A;
      print(newton);     // generalized Newton identities in triangular form
+}
 ///////////////////////////////////////////////////////////////////////////////
 // forms a list of special combinations needed for computation of Waring's
+// forms a list of special combinations needed for computation of Waring's
 //function
 static proc combinations_sum (int m, int n)
 …
    count++;
+  }
   if (flag)
+  if (flag)
+  {
    for (j=2; j<=m; j++)
 …
     interm[i][j]=interm[i][j] div j;
+   }
    result[size(result)+1]=interm[i];
+   result[size(result)+1]=interm[i];
+  }
+ }
 …
 "USAGE:    sysBin (v,Q,n,[odd]);  v,n,odd are int, Q is list of int's
 @format
           - v a number if errors,
           - Q is a generating set of the code,
           - n the length,
           - odd is an additional parameter: if
            set to 1, then the generating set is enlarged by odd elements,
+          - v a number if errors,
+          - Q is a generating set of the code,
+          - n the length,
+          - odd is an additional parameter: if
+           set to 1, then the generating set is enlarged by odd elements,
            which are 2^(some power)*(some elment in the gen.set) mod n
 @end format
 RETURN:    the ring with the resulting system called 'bin'
 THEORY:  Based on Q of the given cyclic code, the procedure constructs
          the corresponding ideal 'bin' with the use of Waring function.
          With its help one can solve the decoding problem.
+THEORY:  Based on Q of the given cyclic code, the procedure constructs
+         the corresponding ideal 'bin' with the use of Waring function.
+         With its help one can solve the decoding problem.
          For basics of the method @ref{Generalized Newton identities}.
 SEE ALSO: sysNewton, sysCRHT
 …
   Q_update=Q;
+ }
  //---- form polynomials for the Bin system via Waring's function ---------
+ //---- form polynomials for the Bin system via Waring's function ---------
  for (i=1; i<=size(Q_update); i++)
+ {
 …
+   }
    count1=0;
    for (j=2; j<=upper-1; j++)
+   for (j=2; j<=upper-1; j++)
+   {
     count1=count1+exp_count(j,2);
 …
    sum_=sum_+coef_*mon;
+  }
   result=result,S(Q_update[i])-sum_;
+  result=result,S(Q_update[i])-sum_;
+ }
  ideal bin=simplify(result,2);
  export bin;
  return(r);
+}
 example
+}
+example
+{
      "EXAMPLE:";  echo = 2;
 …
  if (ncols(x)!=nrows(g)) {print("ERRORencode2!");}
  return(x*g);
+}
 example
+}
+example
+{
      "EXAMPLE:";  echo = 2;
 …
 ,0,0,1,1,1,0;
      //encode x with the generator matrix g
      print(encode(x,g));
+     print(encode(x,g));
+}
 …
  if (nrows(c)>1) {print("ERRORsyndrome1!");}
  if (ncols(c)!=ncols(h)) {print("ERRORsyndrome2!");}
  return(h*transpose(c));
+}
 example
+ return(h*transpose(c));
+}
+example
+{
      "EXAMPLE:";  echo = 2;
 …
      print(syndrome(check,c));
      c[1,3]=1;
      //now c is a codeword
+     //now c is a codeword
      print(syndrome(check,c));
+}
 …
 "USAGE:   sysQE(check,y,t,[fieldeq,formal]);check,y matrix;t,fieldeq,formal int
 @format
         - check is the check matrix of the code
         - y is a received word,
         - t the number of errors to be corrected,
         - if fieldeq=1, then field equations are added,
         - if formal=0, field equations on (known) syndrome variables
           are not added, in order to add them (note that the exponent should
+        - check is the check matrix of the code
+        - y is a received word,
+        - t the number of errors to be corrected,
+        - if fieldeq=1, then field equations are added,
+        - if formal=0, field equations on (known) syndrome variables
+          are not added, in order to add them (note that the exponent should
           be as a number of elements in the INITIAL alphabet) one
           needs to set formal>0 for the exponent
 @end format
 RETURN:   the ring to work with together with the resulting system called 'qe'
 THEORY:  Based on 'check' of the given linear code, the procedure constructs
+THEORY:  Based on 'check' of the given linear code, the procedure constructs
          the corresponding ideal that gives an opportunity to compute
          unknown syndrome of the received word y. Further
          one is able to solve the decoding problem.
+         one is able to solve the decoding problem.
          For basics of the method @ref{Decoding method based on quadratic equations}.
 SEE ALSO: sysFL
 …
   formal=#[2];
+ }
  def br=basering;
+ 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;
 …
   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;
 …
   a[i]=transf*a[i];
+ }
  //----------- compute the structure constants ------------------------
  tim=rtimer;
  matrix te[n][1];
+ matrix te[n][1];
  for (i=1; i<=n; i++)
+ {
 …
+  {
    if ((j<i)&&(i<=t+1)) {c[i][j]=c[j][i];}
    else
+   else
+   {
     if (i+j<=n+1)
+    if (i+j<=n+1)
+    {
      c[i][j]=te;
      c[i][j][i+j-1,1]=1;
+     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;
+ tim=rtimer;
  if (formal==0)
+ {
 …
      result=result,V(j)^q-V(j);
+  }
+ }
+ }
  //----- form the quadratic equations according to the theory -----------
  for (i=1; i<=n; i++)
 …
+  }
   result=result,sum1-sum3;
+ }
+ }
  result=simplify(result,2);
  ideal qe=result;
  export qe;
  return(work);
+}
 example
+ return(work);
+}
+example
+{
      "EXAMPLE:";  echo = 2;
      intvec v = option(get);
      //correct 2 errors in [7,3] 8-ary code RS code
      int t=2; int q=8; int n=7; int redun=4;
 …
      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));
+     matrix rec[1][7]=errorInsert(y,list(2,4),list(1,a));
      //generate the system
      def A=sysQE(h,rec,t);
      setring A;
      print(qe);
      //let us decode
      option(redSB);
      ideal sys_qe=std(qe);
      print(sys_qe);
+     print(sys_qe);
      option(set,v);
+}
 …
 "USAGE:  errorInsert(y,pos,val); y is matrix, pos,val list of int's
 @format
         - y is a (code) word,
         - pos = positions where errors occured,
+        - y is a (code) word,
+        - pos = positions where errors occured,
         - val = their corresponding values
 @end format
 …
+ }
  return(result);
+}
 example
+}
+example
+{
      "EXAMPLE:";  echo = 2;
 …
      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));
 …
 "USAGE:    errorRand(y, num, e); y matrix, num,e int
 @format
           - y is a (code) word,
           - num is the number of errors,
+          - y is a (code) word,
+          - num is the number of errors,
           - e is an extension degree (if one wants values to be from GF(p^e))
 @end format
 …
+   }
    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;
+ }
+  {
+   tempnum=randomvector(1,e);
+   if (tempnum!=0) {flag=0;}
+  }
+  val[i]=tempnum;
+ }
  for (i=1; i<=size(pos); i++)
+ {
 …
+ }
  return(result);
+}
 example
+}
+example
+{
   "EXAMPLE:";  echo = 2;
 …
      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);
+     print(rec-y);
+}
 …
 "USAGE:    randomCheck(m, n, e); m,n,e are int
 @format
           - m x n are dimensions of the matrix,
+          - m x n are dimensions of the matrix,
           - e is an extension degree (if one wants values to be from GF(p^e))
 @end format
 …
+  {
    rand[i,j]=temp[j,1];
+  }
+  }
+ }
  result=concat(rand,unitmat(m));
  return(result);
+}
 example
+{
   "EXAMPLE:";  echo = 2;
+}
+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);
+     print(g);
+}
 …
 "USAGE:   genMDSMat(n, a); n is int, a is number
 @format
         - n x n are dimensions of the MDS matrix,
+        - n x n are dimensions of the MDS matrix,
         - a is a primitive element of the field.
 @end format
 NOTE:   An MDS matrix is constructed in the following way. We take a to be a
         generator of the multiplicative group of the field. Then we construct
+@end format
+NOTE:   An MDS matrix is constructed in the following way. We take a to be a
+        generator of the multiplicative group of the field. Then we construct
         the Vandermonde matrix with this a.
 ASSUME:   extension field should already be defined
+ASSUME:   extension field should already be defined
 RETURN:   a matrix with the MDS property
 EXAMPLE:  example genMDSMat; shows an example
+EXAMPLE:  example genMDSMat; shows an example
+"
+{
 …
+  {
    result[j+1,i+1]=(a^i)^j;
+  }
+  }
+ }
  return(result);
+}
 example
+}
+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);
+     print(h_full);
+}
 …
 "USAGE:  mindist (check, q); check matrix, q int
 @format
         - check is a check matrix,
+        - check is a check matrix,
         - q is the field size
 @end format
 RETURN:  minimum distance of the code together with the time needed for its
+RETURN:  minimum distance of the code together with the time needed for its
          calculation
 EXAMPLE: example mindist; shows an example
+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);
+ 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);
+ ring work=(q,a),(V(1..t),U(1..n)),dp;
+ matrix check=imap(br,check);
  ideal temp;
  int count=1;
 …
  int flag2;
  int i, tim, timsolve;
  matrix z[1][n];
+ matrix z[1][n];
  option(redSB);
  def A=sysQE(check,z,count);
  //proceed with solving the system w.r.t zero vector until some solutions
+ //proceed with solving the system w.r.t zero vector until some solutions
  //are found
  while (flag)
 …
     tim=rtimer;
     temp=std(temp);
     timsolve=timsolve+rtimer-tim;
+    timsolve=timsolve+rtimer-tim;
     flag2=1;
     setring work;
     temp=imap(A,temp);
+    temp=imap(A,temp);
     for (i=1; i<=n; i++)
+    {
       if
         (temp[i]!=U(n-i+1))
+      if
+        (temp[i]!=U(n-i+1))
+        {
           flag2=0;
+        }
+    }
     if (!flag2)
+    if (!flag2)
+    {
       flag=0;
+    }
     else
+    else
+    {
       count++;
+    }
+ }
  list result=list(count,timsolve);
+ list result=list(count,timsolve);
  option(set,vopt);
  return(result);
+}
 example
+ 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;
+     //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]);
+     list l=mindist(h);
+     print(l[1]);
+     //time for the comutation in secs
+     print(l[2]);
+}
 …
 @format
           - check is the check matrix of the code,
           - rec is a received word,
+          - rec is a received word,
           - t is an upper bound for the number of errors one wants to correct
 @end format
 ASSUME:   Errors in rec should be correctable, otherwise the output is
+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
+EXAMPLE:  example decode; shows an example
+"
+{
  intvec vopt = option(get);
+ intvec vopt = option(get);
  def br=basering;
  int n=ncols(check);
  int count=1;
+ int n=ncols(check);
+ int count=1;
  def A=sysQE(check,rec,count);
  while(1)
+ {
   A=sysQE(check,rec,count);
   setring A;
+  setring A;
   matrix h_full=genMDSMat(n,a);
   matrix rec=imap(br,rec);
   option(redSB);
   ideal qe_red=std(qe);
+  ideal qe_red=std(qe);
   if (qe_red[1]!=1)
+  {
     break;
+  }
   else
+  else
+  {
     count++;
 …
   setring br;
+ }
  setring A;
+ setring A;
  //obtain a codeword
  //this works only if our code is indeed can correct these errors
 …
+ {
   syn[i,1]=-qe_red[n-i+1]+lead(qe_red[n-i+1]);
+ }
  matrix real_syn=inverse(h_full)*syn;
+ }
+ matrix real_syn=inverse(h_full)*syn;
  setring br;
  matrix real_syn=imap(A,real_syn);
  option(set,vopt);
  return(rec-transpose(real_syn));
+}
 example
+{
      "EXAMPLE:";  echo = 2;
      //correct 1 error in [15,7] binary code
+ 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;
+     ring r=(q,a),x,dp;
      //generate random check matrix
      matrix h=randomCheck(redun,n,1);
+     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);
+     print(rec);
      //let us decode
      matrix dec_word=decode(h,rec);
      print(dec_word);
+     print(dec_word);
+}
 …
 @format
           - redun is a redundabcy of a (random) code,
           - q is the field size,
+          - q is the field size,
           - ncodes is the number of random codes to be processed,
           - ntrials is the number of received vectors per code to be corrected
           - If e is given it sets the correction capacity explicitly. It
+          - If e is given it sets the correction capacity explicitly. It
           should be used in case one expects some lower bound,
           otherwise the procedure tries to compute the real minimum distance
+          otherwise the procedure tries to compute the real minimum distance
           to find out the error-correction capacity
 @end format
 RETURN:    nothing;
 EXAMPLE:   example decodeRandom; shows an example
+RETURN:    nothing;
+EXAMPLE:   example decodeRandom; shows an example
+"
+{
 …
  option(redSB);
  def br=basering;
  matrix h_full=genMDSMat(n,a);
+ matrix h_full=genMDSMat(n,a);
  matrix z[1][ncols(h_full)];
 …
      dist=tmp[1];
      printf("d= %p",dist);
      t=(dist-1) div 2;
+  }
   tim2=rtimer;
+     t=(dist-1) div 2;
+  }
+  tim2=rtimer;
   tim3=rtimer;
 …
   ideal sys2,sys3;
   matrix h=imap(br,h);
   matrix g=dual_code(h);
+  matrix g=dual_code(h);
   ideal sys=qe;
   print("The system is generated");
   timdec3=timdec3+rtimer-tim3;
+  timdec3=timdec3+rtimer-tim3;
   //------ modify the template according to every received word --------------
   for (j=1; j<=ntrials; j++)
+  {
    word=randomvector(n-redun,1);
+   word=randomvector(n-redun,1);
    y=encode(transpose(word),g);
    rec=errorRand(y,t,1);
    sys2=add_synd(rec,h,redun,sys);
+   rec=errorRand(y,t,1);
+   sys2=add_synd(rec,h,redun,sys);
    tim=rtimer;
    sys3=std(sys2);
    timdec=timdec+rtimer-tim;
+  }
+   sys3=std(sys2);
+   timdec=timdec+rtimer-tim;
+  }
   timdec2=timdec2+rtimer-tim2;
   kill A;
   option(set,vopt);
+ }
+ }
  printf("Time for mindist: %p", timdist);
  printf("Time for GB in mindist: %p", timdist);
  printf("Time for decoding: %p", timdec2);
  printf("Time for GB in decoding: %p", timdec);
  printf("Time for sysQE in decoding: %p", timdec3);
+}
 example
+ printf("Time for sysQE in decoding: %p", timdec3);
+}
+example
+{
      "EXAMPLE:";  echo = 2;
      int q=32; int n=25; int redun=n-11; int t=redun+1;
      ring r=(q,a),x,dp;
      // correct 2 errors in 5 random binary codes, 50 trials each
      decodeRandom(n,redun,5,50,2);
+     decodeRandom(n,redun,5,50,2);
+}
 …
 proc decodeCode(matrix check, int ntrials, list #)
 "USAGE:     decodeCode(check, ntrials, [e]); check matrix, ntrials,e int
 @format
            - check is a check matrix for the code,
            - ntrials is the number of received vectors per code to be
+"USAGE:     decodeCode(check, ntrials, [e]); check matrix, ntrials,e int
+@format
+           - check is a check matrix for the code,
+           - ntrials is the number of received vectors per code to be
            corrected.
            - If e is given it sets the correction capacity explicitly. It
+           - If e is given it sets the correction capacity explicitly. It
            should be used in case one expects some lower bound,
            otherwise the procedure tries to compute the real minimum distance
+           otherwise the procedure tries to compute the real minimum distance
            to find out the error-correction capacity
 @end format
 RETURN:     nothing;
 EXAMPLE:    example decodeCode; shows an example
+RETURN:     nothing;
+EXAMPLE:    example decodeCode; shows an example
+"
+{
 …
  option(redSB);
  def br=basering;
  matrix h_full=genMDSMat(n,a);
+ matrix h_full=genMDSMat(n,a);
  matrix z[1][ncols(h_full)];
  setring br;
 …
    dist=tmp[1];
    printf("d= %p",dist);
    t=(dist-1) div 2;
+ }
  tim2=rtimer;
+   t=(dist-1) div 2;
+ }
+ tim2=rtimer;
  tim3=rtimer;
 …
  ideal sys2,sys3;
  matrix h=imap(br,h);
  matrix g=dual_code(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 ---------------
+ timdec3=timdec3+rtimer-tim3;
+ //--- modify the template according to every received word ---------------
  for (j=1; j<=ntrials; j++)
+ {
    word=randomvector(n-redun,1);
+   word=randomvector(n-redun,1);
    y=encode(transpose(word),g);
    rec=errorRand(y,t,1);
    sys2=add_synd(rec,h,redun,sys);
+   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;
+   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);
+ printf("Time for sysQE in decoding: %p", timdec3);
  option(set,vopt);
+}
 example
+}
+example
+{
      "EXAMPLE:";  echo = 2;
 …
      ring r=(q,a),x,dp;
      matrix check=randomCheck(redun,n,1);
      // correct 2 errors in using the code above, 50 trials
      decodeCode(check,50,2);
+     // correct 2 errors in using the code above, 50 trials
+     decodeCode(check,50,2);
+}
 …
      matrix s[redun][1]=syndrome(check,rec);
      for (int i=1; i<=redun; i++)
+     {
           result[i]=result[i]-s[i,1];
+     {
+          result[i]=result[i]-s[i,1];
+     }
      return(result);
 …
 ///////////////////////////////////////////////////////////////////////////////
 // return index of an element in the ideal where it does not vanish at the
+// 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)
 …
 ///////////////////////////////////////////////////////////////////////////////
 // check whether given polynomial is divisible by some leading monomial of the
+// check whether given polynomial is divisible by some leading monomial of the
 //ideal
 static proc divisible (poly m, ideal G)
 …
      for (int i=1; i<=size(G); i++)
+     {
           if (m/leadmonom(G[i])!=0) {return(1);}
+          if (m/leadmonom(G[i])!=0) {return(1);}
+     }
      return(0);
 …
 proc vanishId (list points)
 "USAGE:  vanishId (points); point is a list of matrices
         'points' is a list of points for which the vanishing ideal is to be
         constructed
+        '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
+EXAMPLE: example vanishId; shows an example
+"
+{
      int m=size(points[1]);
      int n=size(points);
      ideal G=1;
      int i,k,j;
 …
      //------------- proceed according to Farr-Gao algorithm ----------------
      for (k=1; k<=n; k++)
+     {
+     {
           i=find_index(G,points[k]);
           cur=G[i];
+          cur=G[i];
           for(j=i+1; j<=size(G); j++)
+          {
 …
           temp=ideal2list(G);
           temp=delete(temp,i);
           G=list2ideal(temp);
+          G=list2ideal(temp);
           for (j=1; j<=m; j++)
+          {
 …
+               {
                     attrib(G,"isSB",1);
                     h=NF((x(j)-points[k][j,1])*cur,G);
+                    h=NF((x(j)-points[k][j,1])*cur,G);
                     temp=ideal2list(G);
                     temp=insert(temp,h);
                     G=list2ideal(temp);
                     G=sort(G)[1];
+                    G=sort(G)[1];
+               }
+          }
 …
      attrib(G,"isSB",1);
      return(G);
+}
 example
+}
+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);
 …
 ///////////////////////////////////////////////////////////////////////////////
 // construct the list of all vectors of length m with elements in p^e, where p
+// construct the list of all vectors of length m with elements in p^e, where p
 //is characteristics
 proc pointsGen (int m, int e)
 …
           count++;
           for (j=1; j<=charac^(e)-1; j++)
+          {
+          {
                result[count][m]=a^j;
                count++;
 …
           return(result);
+     }
      list prev=pointsGen(m-1,e);
+     list prev=pointsGen(m-1,e);
      for (i=1; i<=size(prev); i++)
+     {
 …
      return(result);
+     }
      if (e==1)
+     {
 …
      int i,j;
      list l=ringlist(basering);
      int charac=l[1][1];
+     int charac=l[1][1];
      list tmp;
      for (i=1; i<=charac^m; i++)
 …
+     {
           for (j=0; j<=charac-1; j++)
+          {
+          {
                result[count][m]=number(j);
                count++;
 …
           return(result);
+     }
      list prev=pointsGen(m-1,e);
+     list prev=pointsGen(m-1,e);
      for (i=1; i<=size(prev); i++)
+     {
 …
      return(result);
+     }
+}
 …
+{
      poly prod=1;
      list rl=ringlist(basering);
+     list rl=ringlist(basering);
      int charac=rl[1][1];
      int l;
 …
 "USAGE:    sysFL (check,y,t,e,s); check,y matrix, t,e,s int
 @format
           - check is a check matrix of the code,
           - y is a received word,
+          - 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,
+          - e is the extension degree,
           - s is the dimension of the point for the vanishing ideal
 @end format
 RETURN:  the system of Fitzgerald-Lax for the given decoding problem
 THEORY:  Based on 'check' of the given linear code, the procedure constructs
+THEORY:  Based on 'check' of the given linear code, the procedure constructs
          the corresponding ideal constructed with a generalization of
          Cooper's philosophy. For basics of the method @ref{Fitzgerald-Lax method}.
 …
+"
+{
      list rl=ringlist(basering);
      int charac=rl[1][1];
+     list rl=ringlist(basering);
+     int charac=rl[1][1];
      int n=ncols(check);
      int m=nrows(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);
+     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];
+               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++)
 …
           result=result,e(i)^(charac^e-1)-1;
+     }
      result=simplify(result,8);
      //--------------- add check realtions --------------------
      poly sum;
      matrix syn[m][1]=syndrome(check,y);
+     matrix syn[m][1]=syndrome(check,y);
      for (i=1; i<=size(funcs); i++)
+     {
+     {
           sum=0;
           for (j=1; j<=t; j++)
 …
                     temp=subst(temp,x(k),x_var(j,k,s));
+               }
                sum=sum+temp*e(j);
+               sum=sum+temp*e(j);
+          }
           result=result,sum-syn[i,1];
+     }
      result=simplify(result,2);
      points=points2;
      export points;
+     export points;
      return(result);
+}
+}
 example
+{
      "EXAMPLE:";  echo = 2;
      intvec vopt = option(get);
      list l=FLpreprocess(3,1,11,2,"");
      def r=l[1];
      setring r;
      int s_work=l[2];
+     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,
 …
      matrix rec[1][11]=errorInsert(y,list(2,4),list(1,-1));
      //the Fitzgerald-Lax system
+     //the Fitzgerald-Lax system
      ideal sys=sysFL(h,rec,2,1,s_work);
      print(sys);
      option(redSB);
      ideal red_sys=std(sys);
      red_sys;
+     red_sys;
      // read the solutions from this redGB
      // the points are (0,0,1) and (0,1,0) with error values 1 and -1 resp.
      // use list points to find error positions;
      points;
      option(set,vopt);
+}
 …
+     {
           s++;
+     }
+     }
      list var_ord;
      int i,j;
 …
                count++;
+          }
+     }
+     }
      list rl;
      list tmp;
      if (e>1)
+     {
 …
           rl[1][2][1]=string("a");
           rl[1][3]=tmp;
           rl[1][3][1]=tmp;
+          rl[1][3][1]=tmp;
           rl[1][3][1][1]=string("lp");
           rl[1][3][1][2]=1;
 …
           rl[1]=p;
+     }
      rl[2]=var_ord;
      rl[3]=tmp;
      rl[3][1]=tmp;
+     rl[3][1]=tmp;
      rl[3][1][1]=string("lp");
      intvec v=1;
 …
           v=v,1;
+     }
      rl[3][1][2]=v;
+     rl[3][1][2]=v;
      rl[3][2]=tmp;
      rl[3][2][1]=string("C");
+     rl[3][2][1]=string("C");
      rl[3][2][2]=intvec(0);
      rl[4]=ideal(0);
+     rl[4]=ideal(0);
      def r2=ring(rl);
      setring r2;
      list l=ringlist(r2);
+     list l=ringlist(r2);
      if (e>1)
+     {
           execute(string("poly f="+minp));
           ideal id=f;
+          execute(string("poly f="+minp));
+          ideal id=f;
           l[1][4]=id;
+     }
      def r=ring(l);
      setring r;
+     setring r;
      return(list(r,s));
+}
 …
 // imitating two indeces
 static proc x_var (int i, int j, int s)
+{
+{
      return(x1(s*(i-1)+j));
+}
 …
 // random vector of length n with entries from p^e, p the characteristic
 static proc randomvector(int n, int e)
+{
+{
     int i;
     matrix result[n][1];
 …
         result[i,1]=asElement(random_prime_vector(e));
+    }
     return(result);
+}
 ///////////////////////////////////////////////////////////////////////////////
 // "convert" representation of an element from the field extension from vector
+    return(result);
+}
+///////////////////////////////////////////////////////////////////////////////
+// "convert" representation of an element from the field extension from vector
 //to an elelemnt
 static proc asElement(list l)
 …
   int i;
   number w=1;
   if (size(l)>1) {w=par(1);}
+  if (size(l)>1) {w=par(1);}
   for (i=0; i<=size(l)-1; i++)
+  {
 …
 // random vector of length n with entries from p, p the characteristic
 static proc random_prime_vector (int n)
+{
+{
   list rl=ringlist(basering);
   int i, charac;
 …
+  {
     if (rl[1][1] mod i ==0)
+    {
+    {
       break;
+    }
+  }
   charac=i;
+  charac=i;
   list l;
   for (i=1; i<=n; i++)
+  {
 …
 proc decodeRandomFL(int n, int redun, int p, int e, int t, int ncodes, int ntrials, string minpol)
 "USAGE:    decodeRandomFL(redun,p,e,n,t,ncodes,ntrials,minpol);
+"USAGE:    decodeRandomFL(redun,p,e,n,t,ncodes,ntrials,minpol);
 @format
           - n is length of codes generated, redun = redundancy of codes
           generated,
           - p is characteristics,
           - e is the extension degree,
           - t is the number of errors to correct,
+          - n is length of codes generated, redun = redundancy of codes
+          generated,
+          - p is characteristics,
+          - e is the extension degree,
+          - t is the number of errors to correct,
           - ncodes is the number of random codes to be processed,
           - ntrials is the number of received vectors per code to be corrected,
           - minpol: due to some pecularities of SINGULAR one needs to provide
+          - minpol: due to some pecularities of SINGULAR one needs to provide
           minimal polynomial for the extension explicitly
 @end format
 …
  intvec vopt = option(get);
  list l=FLpreprocess(p,e,n,t,minpol);
+ list l=FLpreprocess(p,e,n,t,minpol);
  def r=l[1];
  int s_work=l[2];
+ int s_work=l[2];
  export(s_work);
  setring r;
  int i,j;
  matrix h, g, word, y, rec;
  int dist, tim, tim2, tim3, timdist, timdec, timdist2, timdec2, timdec3;
  ideal sys, sys2, sys3;
  list tmp;
  option(redSB);
  matrix z[1][n];
+ list tmp;
+ option(redSB);
+ matrix z[1][n];
  for (i=1; i<=ncodes; i++)
+ {
      h=randomCheck(redun,n,e);
+     h=randomCheck(redun,n,e);
      g=dual_code(h);
      tim2=rtimer;
+     tim2=rtimer;
      tim3=rtimer;
      //---------------- generate the template system -----------------------
      sys=sysFL(h,z,t,e,s_work);
+     //---------------- 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++)
 …
           y=encode(transpose(word),g);
           rec=errorRand(y,t,e);
           sys2=LF_add_synd(rec,h,sys);
+          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);
+ printf("Time for generating Fitzgerald-Lax system during decoding: %p", timdec3);
  option(set,vopt);
+}
+}
 example
+{
      "EXAMPLE:";  echo = 2;
      // correcting one error for one random binary code of length 25,
+     // correcting one error for one random binary code of length 25,
      //redundancy 14; 300 words are processed
      decodeRandomFL(25,14,2,1,1,1,300,"");
 …
 "EXAMPLE:";  echo = 2;
 int q=128; int n=120; int redun=n-30;
+int q=128; int n=120; int redun=n-30;
 ring r=(q,a),x,dp;
 decodeRandom(n,redun,1,1,6);
 int q=128; int n=120; int redun=n-20;
+int q=128; int n=120; int redun=n-20;
 ring r=(q,a),x,dp;
 decodeRandom(n,redun,1,1,9);
 …
 "EXAMPLE:";  echo = 2;
 int q=128; int n=120; int redun=n-40;
+int q=128; int n=120; int redun=n-40;
 ring r=(q,a),x,dp;
 decodeRandom(n,redun,1,1,6);
 …
 decodeRandom(n,redun,1,1,24);
 int q=256; int n=150; int redun=n-10;
+int q=256; int n=150; int redun=n-10;
 ring r=(q,a),x,dp;
 decodeRandom(n,redun,1,1,26);
 …
 setring A;
 print(qe);
 ideal red_qe=stdfglm(qe);
+ideal red_qe=stdfglm(qe);
 */

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset 7f3ad4 in git for Singular/LIB/decodegb.lib

Legend:

Singular/LIB/decodegb.lib

Download in other formats: