Changeset fba86f8 in git for Singular/LIB/decodegb.lib


Ignore:
Timestamp:
Oct 8, 2008, 6:52:24 PM (16 years ago)
Author:
Stanislav Bulygin <bulygin@…>
Branches:
(u'spielwiese', '17f1d200f27c5bd38f5dfc6e8a0879242279d1d8')
Children:
7d56875cf7720719f5b6acd193768a4d5b57848f
Parents:
02efaee061610b6730796d2e862e2ab63663c50a
Message:
indentation, names, formatting,#,option


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

Legend:

Unmodified
Added
Removed
  • Singular/LIB/decodegb.lib

    r02efae rfba86f8  
    11///////////////////////////////////////////////////////////////////////////////
    2 version="$Id: decodegb.lib,v 1.5 2008-10-07 14:14:56 bulygin Exp $";
     2version="$Id: decodegb.lib,v 1.6 2008-10-08 16:52:24 bulygin Exp $";
    33category="Coding theory";
    44info="
    5 LIBRARY:   decodedistGB.lib     Generating and solving systems of polynomial equations for decoding and finding the minimum distance of linear codes
    6 AUTHORS:   Stanislav Bulygin,   bulygin@mathematik.uni-kl.de                     
     5LIBRARY: decodegb.lib         Decoding and min distance of linear codes with GB
     6AUTHOR:  Stanislav Bulygin,   bulygin@mathematik.uni-kl.de                     
    77
    88OVERVIEW:
    9  In this library we generate several systems used for decoding cyclic codes and finding their minimum distance.
    10  Namely, we work with the Cooper's philosophy and generalized Newton identities.
    11  The original method of quadratic equations is worked out here as well.
    12  We also (for comparison) enable to work with the system of Fitzgerald-Lax.
    13  We provide some auxiliary functions for further manipulations and decoding.
    14  For an overview of the methods mentioned above, cf. Stanislav Bulygin, Ruud Pellikaan: 'Decoding and finding the minimum distance with
    15  Groebner bases: history and new insights', in 'Selected Topics in Information and Coding Theory', World Scientific (2008) (to appear) (*).
    16  For the vanishing ideal computation the algorithm of Farr and Gao is implemented.
     9 In this library we generate several systems used for decoding cyclic codes and
     10 finding their minimum distance. Namely, we work with the Cooper's philosophy
     11 and generalized Newton identities. The original method of quadratic equations
     12 is worked out here as well. We also (for comparison) enable to work with the
     13 system of Fitzgerald-Lax. We provide some auxiliary functions for further
     14 manipulations and decoding. For an overview of the methods mentioned above,
     15 cf.
     16 @format
     17  Stanislav Bulygin, Ruud Pellikaan: 'Decoding and finding the minimum
     18  distance with Groebner bases: history and new insights', in 'Selected Topics
     19  in Information and Coding Theory', World Scientific (2008) (to appear) (*).
     20 @end format
     21 For the vanishing ideal computation the algorithm of Farr and Gao is
     22 implemented.
    1723
    1824MAIN PROCEDURES:
    19  sysCRHT(n,defset,e,q,m,#); generates the CRHT-ideal that follows Cooper's philosophy
    20  sysCRHTMindistBinary(n,defset,w); generates the CRHT-ideal to find the minimum distance in the binary case
    21  sysNewton(n,defset,t,q,m,#); generates the ideal with the Generalized Newton identities
    22  sysBin(v,Q,n,#);  generates Bin system as in the work of Augot et.al, cf. [*] for the reference
    23  encode(x,g);  encodes given message with the given generator matrix
    24  syndrome(h, c);  computes a syndrome w.r.t. the given check matrix
    25  sysQE(check,y,t,fieldeq,formal); generates the system of quadratic equations as in [*]
    26  errorInsert(y,pos,val);  inserts errors in a word
    27  errorRand(y,num,e);  inserts random errors in a word
    28  randomCheck(m,n,e);  generates a random check matrix
    29  genMDSMat(n,a);  generates an MDS (actually an RS) matrix
    30  mindist(check);   computes the minimum distance of a code
    31  decode(rec);  decoding of a word using the system of quadratic equations
    32  solveForRandom(redun,ncodes,ntrials,#); a procedure for manipulation with random codes
    33  solveForCode(check,ntrials,#); a procedure for manipulation with the given codes
    34  vanishId(points); computes the vanishing ideal for the given set of points
    35  sysFL(check,y,t,e,s); generates the Fitzgerald-Lax system
    36  solveForRandomFL(n,redun,p,e,t,ncodes,ntrials,minpol); a procedure for manipulation with random codes via Fitzgerald-Lax
     25 sysCRHT(..);        generates the CRHT-ideal as in Cooper's philosophy
     26 sysCRHTMindist(..); CRHT-ideal to find the minimum distance in the binary case
     27 sysNewton(..);      generates the ideal with the generalized Newton identities
     28 sysBin(..);         generates Bin system using Waring function
     29 encode(x,g);        encodes given message x with the given generator matrix g
     30 syndrome(h,c);      computes a syndrome w.r.t. the given check matrix
     31 sysQE(..);          generates the system of quadratic equations for decoding
     32 errorInsert(..);    inserts errors in a word
     33 errorRand(y,num,e); inserts random errors in a word
     34 randomCheck(m,n,e); generates a random check matrix
     35 genMDSMat(n,a);     generates an MDS (actually an RS) matrix
     36 mindist(check);     computes the minimum distance of a code
     37 decode(rec);        decoding of a word using the system of quadratic equations
     38 decodeRandom(..); a procedure for manipulation with random codes
     39 decodeCode(..);   a procedure for manipulation with the given code
     40 vanishId(points);   computes the vanishing ideal for the given set of points
     41 sysFL(..);          generates the Fitzgerald-Lax system
     42 decodeRandomFL(..); manipulation with random codes via Fitzgerald-Lax
    3743 
    3844
     
    7682
    7783///////////////////////////////////////////////////////////////////////////////
    78 // the poolynomial for Sala's restrictions
     84// the polynomial for Sala's restrictions
    7985static proc p_poly(int n, int a, int b)
    8086{
     
    8995///////////////////////////////////////////////////////////////////////////////
    9096
    91 proc sysCRHT (int n, list defset, int e, int q, int m, int #)
    92 "USAGE:   sysCRHT(n,defset,e,q,m,#);  n length of the cyclic code, defset is a list representing the defining set,
    93           e the error-correcting capacity, m degree extension of the splitting field, if #>0 additional equations
    94           representing the fact that every two error positions are either different or at least one of them is zero
    95 RETURN:   the ring to work with the CRHT-ideal (with Sala's additions), the ideal itself is exported with the name 'crht'
    96 EXAMPLE:  example sysCRHT; shows an example
     97proc sysCRHT (int n, list defset, int e, int q, int m, list #)
     98"USAGE:   sysCRHT(n,defset,e,q,m,[k]); n,e,q,m,k int, defset list of int's
     99@format
     100         - n length of the cyclic code,
     101         - defset is a list representing the defining set,
     102         - e the error-correcting capacity,
     103         - q field size
     104         - m degree extension of the splitting field,
     105         - if k>0 additional equations representing the fact that every two
     106         error positions are either different or at least one of them is zero
     107@end format
     108RETURN: 
     109@format
     110        the ring to work with the CRHT-ideal (with Sala's additions),
     111        containig an ideal with name 'crht'
     112@end format
     113EXAMPLE: example sysCRHT; shows an example
    97114"
    98115{
     
    102119  int i,j;
    103120  poly sum;
     121  int k;
     122  if ( size(#) > 0)
     123  {
     124    k = #[1];
     125  }
    104126 
    105   //------------ add check equations -----------------
     127  //---------------------- add check equations --------------------------
    106128  for (i=1; i<=r; i++)
    107129  {
     
    114136  }
    115137 
    116   //------------ field equations on syndromes -----------
     138  //--------------------- field equations on syndromes ------------------
    117139  for (i=1; i<=r; i++)
    118140  {
     
    120142  }
    121143 
    122   //------------ restrictions on error-locations: n-th roots of unity --------------
     144  //------ restrictions on error-locations: n-th roots of unity ----------
    123145  for (i=1; i<=e; i++)
    124146  {
     
    131153  } 
    132154 
    133   //------------ add Sala's additional conditions if necessary -----------------
    134   if (#)
     155  //--------- add Sala's additional conditions if necessary --------------
     156  if ( k > 0 )
     157 
    135158  {
    136159    for (i=1; i<=e; i++)
     
    144167  export crht; 
    145168  return(@crht);   
    146 } example
    147 {
    148   "EXAMPLE:"; echo=2;
     169}
     170example
     171{ "EXAMPLE:"; echo=2;
    149172  // binary cyclic [15,7,5] code with defining set (1,3)
    150  
    151   list defset=1,3; // defining set
    152  
    153   int n=15; // length
    154   int e=2; // error-correcting capacity
    155   int q=2; // basefield size
    156   int m=4; // degree extension of the splitting field
    157   int sala=1; // indicator to add additional equations
     173  intvec v = option(get);
     174
     175  list defset=1,3;           // defining set 
     176  int n=15;                  // length
     177  int e=2;                   // error-correcting capacity
     178  int q=2;                   // basefield size
     179  int m=4;                   // degree extension of the splitting field
     180  int sala=1;                // indicator to add additional equations
    158181   
    159182  def A=sysCRHT(n,defset,e,q,m);
    160183  setring A;
    161   A; // shows the ring we are working in
    162   print(crht); // the CRHT-ideal
     184  A;                         // shows the ring we are working in
     185  print(crht);               // the CRHT-ideal
    163186  option(redSB);
    164   ideal red_crht=std(crht);
    165   // reduced Groebner basis
     187  ideal red_crht=std(crht);  // reduced Groebner basis
    166188  print(red_crht);
    167189 
     
    169191  A=sysCRHT(n,defset,e,q,m,sala);
    170192  setring A; 
    171   print(crht); // the CRHT-ideal with additional equations from Sala
     193  print(crht);                // CRHT-ideal with additional equations from Sala
    172194  option(redSB);
    173   ideal red_crht=std(crht);
    174   // reduced Groebner basis
    175   print(red_crht);
    176   // general error-locator polynomial for this code
    177   red_crht[5];
    178 }
    179 
    180 ///////////////////////////////////////////////////////////////////////////////
    181 
    182 
    183 proc sysCRHTMindistBinary (int n, list defset, int w)
    184 "USAGE:  sysCRHTMindistBinary(n,defset,w);  n length of the cyclic code, defset is a list representing the defining set,
    185          w is a candidate for the minimum distance
    186 RETURN:  the ring to work with the Sala's ideal for the minimum distance, the ideal itself is exported with the name 'crht_md'
    187 EXAMPLE: example sysCRHTMindistBinary; shows an example
     195  ideal red_crht=std(crht);   // reduced Groebner basis
     196  print(red_crht);           
     197  red_crht[5];                // general error-locator polynomial for this code
     198  option(set,v);
     199}
     200
     201///////////////////////////////////////////////////////////////////////////////
     202
     203
     204proc sysCRHTMindist (int n, list defset, int w)
     205"USAGE:  sysCRHTMindist(n,defset,w);  n,w are int, defset is list of int's 
     206@format
     207        - n length of the cyclic code,
     208        - defset is a list representing the defining set,
     209        - w is a candidate for the minimum distance
     210@end format
     211RETURN:
     212@format
     213         the ring to work with the Sala's ideal for the minimum distance
     214         containing the ideal with name 'crht_md'
     215@end format
     216EXAMPLE: example sysCRHTMindist; shows an example
    188217"
    189218{
     
    223252  export crht_md; 
    224253  return(@crht_md);   
    225 } example
     254}
     255example
    226256{
    227257  "EXAMPLE:"; echo=2;
     258  intvec v = option(get);
    228259  // binary cyclic [15,7,5] code with defining set (1,3)
    229260 
    230   list defset=1,3; // defining set
     261  list defset=1,3;             // defining set 
     262  int n=15;                    // length
     263  int d=5;                     // candidate for the minimum distance 
     264     
     265  def A=sysCRHTMindist(n,defset,d);
     266  setring A;
     267  A;                           // shows the ring we are working in
     268  print(crht_md);              // the Sala's ideal for mindist
     269  option(redSB);
     270  ideal red_crht_md=std(crht_md);     
     271  print(red_crht_md);          // reduced Groebner basis
    231272 
    232   int n=15; // length
    233   int d=5; // candidate for the minimum distance 
    234      
    235   def A=sysCRHTMindistBinary(n,defset,d);
    236   setring A;
    237   A; // shows the ring we are working in
    238   print(crht_md); // the Sala's ideal for mindist
    239   option(redSB);
    240   ideal red_crht_md=std(crht_md);
    241   // reduced Groebner basis
    242   print(red_crht_md); 
     273  option(set,v);
    243274}
    244275
     
    253284///////////////////////////////////////////////////////////////////////////////
    254285
    255 proc sysNewton (int n, list defset, int t, int q, int m, int #)
    256 "USAGE:   sysNewton (n, defset, t, q, m, #);  n is length, defset is the defining set,
    257           t is the number of errors, q is basefield size, m is degree extension of the splitting field
    258           if triangular>0 it indicates that Newton identities in triangular form should be constructed
    259 RETURN:   the ring to work with the generalized Newton identities (in triangular form if applicable),
    260           the ideal itself is exported with the name 'newton'
     286proc sysNewton (int n, list defset, int t, int q, int m, list #)
     287"USAGE:   sysNewton (n,defset,t,q,m,[tr]); n,t,q,m,tr int, defset list int's
     288@format
     289         - n is length,
     290         - defset is the defining set,
     291         - t is the number of errors,
     292         - q is basefield size,
     293         - m is degree extension of the splitting field,
     294         - if tr>0 it indicates that Newton identities in triangular
     295           form should be constructed
     296@end format
     297RETURN:   
     298@format
     299        the ring to work with the generalized Newton identities (in triangular
     300        form if applicable) containing the ideal with name 'newton'
     301@end format
    261302EXAMPLE:  example sysNewton; shows an example
    262303"
     
    265306 int i,j;
    266307 int flag;
     308 int tr;
     309 
     310 if (size(#)>0)
     311 {
     312  tr=#[1];
     313 }
     314 
    267315 for (i=n; i>=1; i--)
    268316 {
     
    295343 
    296344 //------------ generate generalized Newton identities ----------
    297  if (#)
     345 if (tr)
    298346 {
    299347  for (i=1; i<=t; i++)
     
    342390 export newton;
    343391 return(@newton);
    344 } example
     392}
     393example
    345394{
    346395     "EXAMPLE:";  echo = 2;
    347      // Newton identities for a binary 3-error-correcting cyclic code of length 31 with defining set (1,5,7)
    348      
    349      int n=31; // length
     396     // Newton identities for a binary 3-error-correcting cyclic code of
     397     //length 31 with defining set (1,5,7)
     398     
     399     int n=31;          // length
    350400     list defset=1,5,7; //defining set
    351      int t=3; // number of errors
    352      int q=2; // basefield size
    353      int m=5; // degree extension of the splitting field
    354      int triangular=1; // indicator of triangular form of Newton identities       
     401     int t=3;           // number of errors
     402     int q=2;           // basefield size
     403     int m=5;           // degree extension of the splitting field
     404     int tr=1;          // indicator of triangular form of Newton identities
     405     
    355406     
    356407     def A=sysNewton(n,defset,t,q,m);
    357408     setring A;
    358      A; // shows the ring we are working in
    359      print(newton); // generalized Newton identities
     409     A;                 // shows the ring we are working in
     410     print(newton);     // generalized Newton identities
    360411     
    361412     //===============================
    362      A=sysNewton(n,defset,t,q,m,triangular);
     413     A=sysNewton(n,defset,t,q,m,tr);
    363414     setring A;     
    364      print(newton); // generalized Newton identities in triangular form
    365 }
    366 
    367 ///////////////////////////////////////////////////////////////////////////////
    368 // forms a list of special combinations needed for computation of Waring's function
     415     print(newton);     // generalized Newton identities in triangular form
     416}
     417
     418///////////////////////////////////////////////////////////////////////////////
     419// forms a list of special combinations needed for computation of Waring's
     420//function
    369421static proc combinations_sum (int m, int n)
    370422{
     
    417469
    418470
    419 proc sysBin (int v, list Q, int n, int#)
    420 "USAGE:    sysBin (v, Q, n, #);  v a number if errors, Q is a generating set of the code, n the length, # is an additional parameter: if 
    421            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
    422 RETURN:    keeps the ring with the resulting system, which ideal is called 'bin'
     471proc sysBin (int v, list Q, int n, list #)
     472"USAGE:    sysBin (v,Q,n,[odd]);  v,n,odd are int, Q is list of int's
     473@format
     474          - v a number if errors,
     475          - Q is a generating set of the code,
     476          - n the length,
     477          - odd is an additional parameter: if 
     478           set to 1, then the generating set is enlarged by odd elements,
     479           which are 2^(some power)*(some elment in the gen.set) mod n
     480@end format
     481THEORY:   
     482@format
     483          The system using Waring's function due to Augot et.al.
     484          is built as in [*].
     485@end format
     486RETURN:    the ring with the resulting system called 'bin'
    423487EXAMPLE:   example sysBin; shows an example
    424488"
    425489{
     490 int odd;
     491 if (size(#)>0)
     492 {
     493  odd=#[1];
     494 }
     495
    426496 //ring r=2,(sigma(1..v),S(1..n)),(lp(v),dp(n));
    427497 ring r=2,(S(1..n),sigma(1..v)),lp;
     
    433503 int count1, count2, upper, coef_, flag, gener;
    434504 list Q_update;
    435  if (#==1)
     505 if (odd==1)
    436506 {
    437507  for (i=1; i<=n; i++)
     
    460530 }
    461531 
    462  //-------------- form polynomials for the Bin system via Waring's function -------------- 
     532 //---- form polynomials for the Bin system via Waring's function --------- 
    463533 for (i=1; i<=size(Q_update); i++)
    464534 {
     
    583653///////////////////////////////////////////////////////////////////////////////
    584654
    585 proc sysQE(matrix check, matrix y, int t, int fieldeq, int formal)
    586 "USAGE:   sysQE(check, y, t, fieldeq, formal); check is the check matrix of the code   
    587           y is a received word, t the number of errors to be corrected,
    588           if fieldeq=1, then field equations are added; if formal=0, field equations on (known) syndrome variables
    589           are not added, in order to add them (note that the exponent should be as a number of elements in the INITIAL alphabet) one
     655proc sysQE(matrix check, matrix y, int t, list #)
     656"USAGE:   sysQE(check,y,t,[fieldeq,formal]);check,y matrix;t,fieldeq,formal int
     657@format
     658        - check is the check matrix of the code   
     659        - y is a received word,
     660        - t the number of errors to be corrected,
     661        - if fieldeq=1, then field equations are added,
     662        - if formal=0, field equations on (known) syndrome variables
     663          are not added, in order to add them (note that the exponent should
     664          be as a number of elements in the INITIAL alphabet) one
    590665          needs to set formal>0 for the exponent
    591 RETURN:   the ring to work with together with the resulting system
     666@end format
     667THEORY:   The system of quadratic equations is built as in [*].
     668RETURN:   the ring to work with together with the resulting system called 'qe'
    592669EXAMPLE:  example sysQE; shows an example
    593670"
    594671{
     672 int fieldeq;
     673 int formal;
     674 if (size(#)>0)
     675 {
     676  fieldeq=#[1];
     677 }
     678 if (size(#)>1)
     679 {
     680  formal=#[2];
     681 }
     682 
    595683 def br=basering;
    596684 list rl=ringlist(br);
     
    632720 matrix transf=inverse(transpose(h_full));
    633721 
    634  //----------- expression matrix of check vectors w.r.t. the MDS basis --------------
     722 //------ expression matrix of check vectors w.r.t. the MDS basis -----------
    635723 tim=rtimer;
    636724 for (i=1; i<=red ; i++)
     
    706794 }
    707795 
    708  //--------- form the quadratic equations according to the theory ---------------
     796 //----- form the quadratic equations according to the theory -----------
    709797 for (i=1; i<=n; i++)
    710798 {
     
    731819 ideal qe=result;
    732820 export qe;
    733  return(work);
    734  //exportto(Top,h_full);   
     821 return(work);
    735822} example
    736823{
    737824     "EXAMPLE:";  echo = 2;
     825     intvec v = option(get);
     826     
    738827     //correct 2 errors in [7,3] 8-ary code RS code
    739828     int t=2; int q=8; int n=7; int redun=4;
     
    744833     matrix x[1][3]=0,0,1,0;
    745834     matrix y[1][7]=encode(x,g);
     835     
    746836     //disturb with 2 errors
    747      matrix rec[1][7]=errorInsert(y,list(2,4),list(1,a));
     837     matrix rec[1][7]=errorInsert(y,list(2,4),list(1,a));
     838     
    748839     //generate the system
    749      def A=sysQE(h,rec,t,0,0);
     840     def A=sysQE(h,rec,t);
    750841     setring A;
    751842     print(qe);
     843     
    752844     //let us decode
    753845     option(redSB);
    754846     ideal sys_qe=std(qe);
    755847     print(sys_qe);     
     848     
     849     option(set,v);
    756850}
    757851
     
    759853
    760854proc errorInsert(matrix y, list pos, list val)
    761 "USAGE:  errorInsert(y, pos, val); y is a (code) word, pos = positions where errors occured, val = their corresponding values
     855"USAGE:  errorInsert(y,pos,val); y is matrix, pos,val list of int's
     856@format
     857        - y is a (code) word,
     858        - pos = positions where errors occured,
     859        - val = their corresponding values
     860@end format
    762861RETURN:  corresponding received word
    763 EXAMPLE: example error; shows an example
     862EXAMPLE: example errorInsert; shows an example
    764863"
    765864{
     
    793892
    794893proc errorRand(matrix y, int num, int e)
    795 "USAGE:    errorRand(y, num, e); y is a (code) word, num is the number of errors, e is an extension degree (if one wants values to
    796            be from GF(p^e))
     894"USAGE:    errorRand(y, num, e); y matrix, num,e int
     895@format
     896          - y is a (code) word,
     897          - num is the number of errors,
     898          - e is an extension degree (if one wants values to be from GF(p^e))
     899@end format
    797900RETURN:    corresponding received word
    798901EXAMPLE:   example errorRand; shows an example
     
    854957///////////////////////////////////////////////////////////////////////////////
    855958
    856 proc randomCheck(int m, int n, int e, int #)
    857 "USAGE:    randomCheck(m, n, e); m x n are dimensions of the matrix, e is an extension degree (if one wants values to
    858            be from GF(p^e))
     959proc randomCheck(int m, int n, int e)
     960"USAGE:    randomCheck(m, n, e); m,n,e are int
     961@format
     962          - m x n are dimensions of the matrix,
     963          - e is an extension degree (if one wants values to be from GF(p^e))
     964@end format
    859965RETURN:    random check matrix
    860966EXAMPLE:   example randomCheck; shows an example
     
    867973 for (i=1; i<=m; i++)
    868974 {
    869   temp=randomvector(n-m,e,#); 
     975  temp=randomvector(n-m,e);
    870976  for (j=1; j<=n-m; j++)
    871977  {
     
    893999
    8941000proc genMDSMat(int n, number a)
    895 "USAGE:   genMDSMat(n, a); n x n are dimensions of the matrix, a is a primitive element of the field.         
    896 NOTE:     An MDS matrix is constructed in the following way. We take a to be a generator of the multiplicative group of the field.
    897           Then we construct the Vandermonde matrix with this a.
     1001"USAGE:   genMDSMat(n, a); n is int, a is number
     1002@format
     1003        - n x n are dimensions of the MDS matrix,
     1004        - a is a primitive element of the field.
     1005@end format   
     1006NOTE:     
     1007@format
     1008        An MDS matrix is constructed in the following way. We take a to be a
     1009        generator of the multiplicative group of the field. Then we construct
     1010        the Vandermonde matrix with this a.
     1011@end format
    8981012ASSUME:   extension field should already be defined
    8991013RETURN:   a matrix with the MDS property
     
    9261040
    9271041proc mindist (matrix check)
    928 "USAGE:  mindist (check, q); check is a check matrix, q = field size
    929 RETURN:  minimum distance of the code together with the time needed for its calculation
     1042"USAGE:  mindist (check, q); check matrix, q int
     1043@format
     1044        - check is a check matrix,
     1045        - q is the field size
     1046@end format
     1047RETURN:
     1048@format
     1049        minimum distance of the code together with the time needed for its
     1050        calculation
     1051@end format
    9301052EXAMPLE: example mindist; shows an example
    9311053"
    9321054{
     1055 intvec vopt = option(get);
     1056
    9331057 int n=ncols(check); int redun=nrows(check); int t=redun+1;
    9341058 
     
    9471071 matrix z[1][n];
    9481072 option(redSB);
    949  def A=sysQE(check,z,count,0,0);
    950 
    951  //----------- proceed with solving the system w.r.t zero vector until some solutions are found --------------------
     1073 def A=sysQE(check,z,count);
     1074
     1075 //proceed with solving the system w.r.t zero vector until some solutions
     1076 //are found
    9521077 while (flag)
    9531078 {
    954     A=sysQE(check,z,count,0,0);
     1079    A=sysQE(check,z,count);
    9551080    setring A;
    9561081    ideal temp=qe;
     
    9791104 }
    9801105 list result=list(count,timsolve);
     1106 
     1107 option(set,vopt);
    9811108 return(result);
    982 } example
     1109}
     1110example
    9831111{
    9841112     "EXAMPLE:";  echo = 2;
     
    9991127
    10001128proc decode(matrix check, matrix rec)
    1001 "USAGE:    decode(check, rec, t); check is the check matrix of the code
    1002            rec is a received word, t is an upper bound for the number of errors one wants to correct
    1003 ASSUME:    Errors in rec should be correctable, otherwise the output is unpredictable
     1129"USAGE:    decode(check, rec, t); check, rec matrix, t int
     1130@format
     1131          - check is the check matrix of the code,
     1132          - rec is a received word,
     1133          - t is an upper bound for the number of errors one wants to correct
     1134@end format
     1135ASSUME:
     1136@format
     1137          Errors in rec should be correctable, otherwise the output is
     1138          unpredictable
     1139@end format
    10041140RETURN:    a codeword that is closest to rec
    10051141EXAMPLE:   example decode; shows an example
    10061142"
    1007 {
     1143{
     1144 intvec vopt = option(get);
     1145 
    10081146 def br=basering;
    10091147 int n=ncols(check);
    10101148 
    10111149 int count=1;
    1012  def A=sysQE(check,rec,count,0,0);
     1150 def A=sysQE(check,rec,count);
    10131151 while(1)
    10141152 {
    1015   A=sysQE(check,rec,count,0,0);
     1153  A=sysQE(check,rec,count);
    10161154  setring A;     
    10171155  matrix h_full=genMDSMat(n,a);
     
    10431181 setring br;
    10441182 matrix real_syn=imap(A,real_syn);
     1183 
     1184 option(set,vopt);
    10451185 return(rec-transpose(real_syn));     
    1046 } example
     1186}
     1187example
    10471188{
    10481189     "EXAMPLE:";  echo = 2;     
     
    10741215
    10751216
    1076 proc solveForRandom(int n, int redun, int ncodes, int ntrials, int #)
    1077 "USAGE:    solveForRandom(redun, q, ncodes, ntrials); redun is a redundabcy of a (random) code,
    1078            q = field size, ncodes = number of random codes to be processed,
    1079            ntrials = number of received vectors per code to be corrected.
    1080            If # is given it sets the correction capacity explicitly. It should be used in case one expexts some lower bound,
    1081            otherwise the procedure tries to compute the real minimum distance to find out the error-correction capacity
     1217proc decodeRandom(int n, int redun, int ncodes, int ntrials, list #)
     1218"USAGE:    decodeRandom(redun,q,ncodes,ntrials,[e]); all parameters int
     1219@format
     1220          - redun is a redundabcy of a (random) code,
     1221          - q is the field size,
     1222          - ncodes is the number of random codes to be processed,
     1223          - ntrials is the number of received vectors per code to be corrected
     1224          - If e is given it sets the correction capacity explicitly. It
     1225          should be used in case one expects some lower bound,
     1226          otherwise the procedure tries to compute the real minimum distance
     1227          to find out the error-correction capacity
     1228@end format
    10821229RETURN:    nothing;
    1083 EXAMPLE:   example solveForRandom; shows an example
     1230EXAMPLE:   example decodeRandom; shows an example
    10841231"
    10851232{
     1233 intvec vopt = option(get);
     1234
    10861235 int i,j;
    10871236 matrix h;
     
    10891238 ideal sys;
    10901239 list tmp;
     1240 int e;
     1241 if (size(#)>0)
     1242 {
     1243  e=#[1];
     1244 }
    10911245
    10921246 option(redSB);
     
    11021256  "check matrix:";
    11031257  print(h);
    1104   if (#>0)
    1105   {
    1106      t=#;
     1258  if (e>0)
     1259  {
     1260     t=e;
    11071261  } else {
    11081262     tim=rtimer;
     
    11181272
    11191273  //------------- generate the template system ----------------------
    1120   def A=sysQE(h,z,t,0,0);
     1274  def A=sysQE(h,z,t);
    11211275  setring A;
    11221276  matrix word,y,rec;
     
    11281282  timdec3=timdec3+rtimer-tim3; 
    11291283
    1130   //------------- modify the template according to every received word -------------------
     1284  //------ modify the template according to every received word --------------
    11311285  for (j=1; j<=ntrials; j++)
    11321286  {
     
    11421296  timdec2=timdec2+rtimer-tim2;
    11431297  kill A;
     1298  option(set,vopt);
    11441299 }
    11451300 printf("Time for mindist: %p", timdist);
     
    11481303 printf("Time for GB in decoding: %p", timdec);
    11491304 printf("Time for sysQE in decoding: %p", timdec3);
    1150 } example
     1305}
     1306example
    11511307{
    11521308     "EXAMPLE:";  echo = 2;
     
    11551311     
    11561312     // correct 2 errors in 5 random binary codes, 50 trials each
    1157      solveForRandom(n,redun,5,50,2);
    1158 }
    1159 
    1160 ///////////////////////////////////////////////////////////////////////////////
    1161 
    1162 
    1163 proc solveForCode(matrix check, int ntrials, int #)
    1164 "USAGE:     solveForCode(check, ntrials);
    1165             check is a check matrix for the code, ntrials = number of received vectors per code to be corrected.
    1166             If # is given it sets the correction capacity explicitly. It should be used in case one expects some lower bound,
    1167             otherwise the procedure tries to compute the real minimum distance to find out the error-correction capacity
     1313     decodeRandom(n,redun,5,50,2);
     1314}
     1315
     1316///////////////////////////////////////////////////////////////////////////////
     1317
     1318
     1319proc decodeCode(matrix check, int ntrials, list #)
     1320"USAGE:     decodeCode(check, ntrials, [e]); check matrix, ntrials,e int
     1321@format           
     1322           - check is a check matrix for the code,
     1323           - ntrials is the number of received vectors per code to be
     1324           corrected.
     1325           - If e is given it sets the correction capacity explicitly. It
     1326           should be used in case one expects some lower bound,
     1327           otherwise the procedure tries to compute the real minimum distance
     1328           to find out the error-correction capacity
     1329@end format
    11681330RETURN:     nothing; 
    1169 EXAMPLE:    example solveForCode; shows an example
     1331EXAMPLE:    example decodeCode; shows an example
    11701332"
    11711333{
     1334 intvec vopt = option(get);
     1335
    11721336 int n=ncols(check);
    11731337 int redun=nrows(check);
     
    11771341 ideal sys;
    11781342 list tmp;
     1343 int e;
     1344 if (size(#)>0)
     1345 {
     1346  e=#[1];
     1347 }
    11791348
    11801349 option(redSB);
     
    11881357
    11891358 //------------------ determine error-correction capacity -------------------
    1190  if (#>0)
    1191  {
    1192     t=#;
     1359 if (e>0)
     1360 {
     1361    t=e;
    11931362 } else {
    11941363   tim=rtimer;
     
    12041373
    12051374 //------------- generate the template system ----------------------
    1206  def A=sysQE(h,z,t,0,0);
     1375 def A=sysQE(h,z,t);
    12071376 setring A;
    12081377 matrix word,y,rec;
     
    12141383 timdec3=timdec3+rtimer-tim3;
    12151384
    1216  //------------- modify the template according to every received word -------------------
     1385 //--- modify the template according to every received word ---------------
    12171386 for (j=1; j<=ntrials; j++)
    12181387 {
     
    12331402 printf("Time for GB in decoding: %p", timdec);
    12341403 printf("Time for sysQE in decoding: %p", timdec3);
    1235 } example
     1404 
     1405 option(set,vopt);
     1406}
     1407example
    12361408{
    12371409     "EXAMPLE:";  echo = 2;
     
    12411413     
    12421414     // correct 2 errors in using the code above, 50 trials
    1243      solveForCode(check,50,2);
     1415     decodeCode(check,50,2);
    12441416}
    12451417
     
    12731445}
    12741446
    1275 //////////////////////////////////////////////////////////////////////////////////////
    1276 // return index of an element in the ideal where it does not vanish at the given point
     1447///////////////////////////////////////////////////////////////////////////////
     1448// return index of an element in the ideal where it does not vanish at the
     1449//given point
    12771450static proc find_index (ideal G, matrix p)
    12781451{
     
    13121485}
    13131486
    1314 ////////////////////////////////////////////////////////////////////////////////////
    1315 // checl whether given polynomial is divisible by some leading monomial of the ideal
     1487///////////////////////////////////////////////////////////////////////////////
     1488// check whether given polynomial is divisible by some leading monomial of the
     1489//ideal
    13161490static proc divisible (poly m, ideal G)
    13171491{
     
    13261500
    13271501proc vanishId (list points)
    1328 "USAGE:  vanishId (points,e); points is a list of points, for which the vanishing ideal is to be constructed.
     1502"USAGE:  vanishId (points); point is a list of matrices
     1503@format
     1504        - points is a list of points for which the vanishing ideal is to be
     1505        constructed       
     1506@end format
    13291507RETURN:  Vanishing ideal corresponding to the given set of points
    13301508EXAMPLE: example vanishId; shows an example
     
    13861564}
    13871565
    1388 //////////////////////////////////////////////////////////////////////////////////////////////////
    1389 // construct the list of all vectors of length m with elements in p^e, where p is characteristics
     1566///////////////////////////////////////////////////////////////////////////////
     1567// construct the list of all vectors of length m with elements in p^e, where p
     1568//is characteristics
    13901569proc pointsGen (int m, int e)
    13911570{
     
    15481727
    15491728proc sysFL (matrix check, matrix y, int t, int e, int s)
    1550 "USAGE:    sysFL (check,y,t,e,s); check is a check matrix of the code, y is a received word, t the number of errors to correct,
    1551            e is the extension degree, s is the dimension of the point for the vanishing ideal
     1729"USAGE:    sysFL (check,y,t,e,s); check,y matrix, t,e,s int
     1730@format
     1731          - check is a check matrix of the code,
     1732          - y is a received word,
     1733          - t the number of errors to correct,
     1734          - e is the extension degree,
     1735          - s is the dimension of the point for the vanishing ideal
     1736@end format
    15521737RETURN:    the system of Fitzgerald-Lax for the given decoding problem
    15531738EXAMPLE:   example sysFL; shows an example
     
    16201805     export points;     
    16211806     return(result);
    1622 } example
     1807}
     1808example
    16231809{
    16241810     "EXAMPLE:";  echo = 2;
     1811     intvec vopt = option(get);
    16251812     
    16261813     list l=FLpreprocess(3,1,11,2,"");
     
    16461833     option(redSB);
    16471834     ideal red_sys=std(sys);
    1648      red_sys; // read the solutions from this redGB
     1835     red_sys;
     1836     // read the solutions from this redGB
    16491837     // the points are (0,0,1) and (0,1,0) with error values 1 and -1 resp.
    16501838     // use list points to find error positions;
    1651      points;     
     1839     points;
     1840     
     1841     option(set,vopt);
    16521842}
    16531843
     
    16911881          rl[1][2][1]=string("a");
    16921882          rl[1][3]=tmp;
    1693           rl[1][3][1]=tmp;
    1694           //rl[1][3][1][1]=string("dp("+string((t-1)*(s+1)+s)+"),lp("+string(s+1)+")");
     1883          rl[1][3][1]=tmp;         
    16951884          rl[1][3][1][1]=string("lp");
    16961885          rl[1][3][1][2]=1;
     
    17041893     rl[3]=tmp;
    17051894     rl[3][1]=tmp;     
    1706      //rl[3][1][1]=string("dp("+string((t-1)*(s+1)+s)+"),lp("+string(s+1)+")");
    17071895     rl[3][1][1]=string("lp");
    17081896     intvec v=1;
     
    17431931///////////////////////////////////////////////////////////////////////////////
    17441932// random vector of length n with entries from p^e, p the characteristic
    1745 static proc randomvector(int n, int e, int #)
     1933static proc randomvector(int n, int e)
    17461934{     
    17471935    int i;
     
    17491937    for (i=1; i<=n; i++)
    17501938    {
    1751         result[i,1]=asElement(random_prime_vector(e,#));
     1939        result[i,1]=asElement(random_prime_vector(e));
    17521940    }
    17531941    return(result);       
    17541942}
    17551943
    1756 ////////////////////////////////////////////////////////////////////////////////////////////
    1757 // "convert" representation of an element from the field extension from vector to an elelemnt
     1944///////////////////////////////////////////////////////////////////////////////
     1945// "convert" representation of an element from the field extension from vector
     1946//to an elelemnt
    17581947static proc asElement(list l)
    17591948{
     
    17711960///////////////////////////////////////////////////////////////////////////////
    17721961// random vector of length n with entries from p, p the characteristic
    1773 static proc random_prime_vector (int n, int #)
    1774 {
    1775   if (#==1)
    1776   {
    1777      list rl=ringlist(basering);     
    1778      int charac=rl[1][1];
    1779   } else {
    1780      int charac=2;
    1781   }
     1962static proc random_prime_vector (int n)
     1963
     1964  list rl=ringlist(basering);
     1965  int i, charac;
     1966  for (i=2; i<=rl[1][1]; i++)
     1967  {
     1968    if (rl[1][1] mod i ==0)
     1969    {     
     1970      break;
     1971    }
     1972  }
     1973  charac=i; 
     1974 
    17821975  list l;
    1783   int i;
     1976 
    17841977  for (i=1; i<=n; i++)
    17851978  {
     
    17911984///////////////////////////////////////////////////////////////////////////////
    17921985
    1793 proc solveForRandomFL(int n, int redun, int p, int e, int t, int ncodes, int ntrials, string minpol)
    1794 "USAGE:    solveForRandomFL(redun,p,e,n,t,ncodes,ntrials,minpol); n = length of codes generated, redun = redundancy of codes generated,
    1795            p = characteristics, e is the extension degree,
    1796            t = number of errors to correct, ncodes = number of random codes to be processed
    1797            ntrials = number of received vectors per code to be corrected
    1798            due to some pecularities of SINGULAR one needs to provide minimal polynomial for the extension explicitly
     1986proc decodeRandomFL(int n, int redun, int p, int e, int t, int ncodes, int ntrials, string minpol)
     1987"USAGE:    decodeRandomFL(redun,p,e,n,t,ncodes,ntrials,minpol);
     1988@format
     1989          - n is length of codes generated, redun = redundancy of codes
     1990          generated,
     1991          - p is characteristics,
     1992          - e is the extension degree,
     1993          - t is the number of errors to correct,
     1994          - ncodes is the number of random codes to be processed,
     1995          - ntrials is the number of received vectors per code to be corrected,
     1996          - minpol: due to some pecularities of SINGULAR one needs to provide
     1997          minimal polynomial for the extension explicitly
     1998@end format
    17991999RETURN:    nothing
    1800 EXAMPLE:   example solveForRandomFL; shows an example
    1801 {
     2000EXAMPLE:   example decodeRandomFL; shows an example
     2001"
     2002{
     2003 intvec vopt = option(get);
     2004
    18022005 list l=FLpreprocess(p,e,n,t,minpol);
    18032006 
     
    18182021 for (i=1; i<=ncodes; i++)
    18192022 {
    1820      h=randomCheck(redun,n,e,1);           
     2023     h=randomCheck(redun,n,e);           
    18212024     g=dual_code(h);
    18222025     tim2=rtimer;       
    18232026     tim3=rtimer;
    18242027
    1825      //------------ generate the template system ------------------              
     2028     //---------------- generate the template system ----------------------- 
    18262029     sys=sysFL(h,z,t,e,s_work);     
    18272030     timdec3=timdec3+rtimer-tim3;
    18282031   
    1829      //------------- modifying the template according to the received word ---------------
     2032     //------ modifying the template according to the received word ---------
    18302033     for (j=1; j<=ntrials; j++)
    18312034     {
     
    18442047 printf("Time for GB in decoding: %p", timdec);
    18452048 printf("Time for generating Fitzgerald-Lax system during decoding: %p", timdec3);
    1846 } example
     2049 
     2050 option(set,vopt);
     2051}
     2052example
    18472053{
    18482054     "EXAMPLE:";  echo = 2;
    18492055     
    1850      // correcting one error for one random binary code of length 25, redundancy 14; 300 words are processed
    1851      solveForRandomFL(25,14,2,1,1,1,300,"");
     2056     // correcting one error for one random binary code of length 25,
     2057     //redundancy 14; 300 words are processed
     2058     decodeRandomFL(25,14,2,1,1,1,300,"");
    18522059}
    18532060
     
    18742081int q=128; int n=120; int redun=n-30;
    18752082ring r=(q,a),x,dp;
    1876 solveForRandom(n,redun,1,1,6);
     2083decodeRandom(n,redun,1,1,6);
    18772084
    18782085int q=128; int n=120; int redun=n-20;
    18792086ring r=(q,a),x,dp;
    1880 solveForRandom(n,redun,1,1,9);
     2087decodeRandom(n,redun,1,1,9);
    18812088
    18822089int q=128; int n=120; int redun=n-10;
    18832090ring r=(q,a),x,dp;
    1884 solveForRandom(n,redun,1,1,19);
     2091decodeRandom(n,redun,1,1,19);
    18852092
    18862093int q=256; int n=150; int redun=n-10;
    18872094ring r=(q,a),x,dp;
    1888 solveForRandom(n,redun,1,1,22);
     2095decodeRandom(n,redun,1,1,22);
    18892096
    18902097//////////////     SOME HARD EXAMPLES    //////////////////////
     
    18962103int q=128; int n=120; int redun=n-40;
    18972104ring r=(q,a),x,dp;
    1898 solveForRandom(n,redun,1,1,6);
     2105decodeRandom(n,redun,1,1,6);
    18992106
    19002107redun=n-30;
    1901 solveForRandom(n,redun,1,1,8);
     2108decodeRandom(n,redun,1,1,8);
    19022109
    19032110redun=n-20;
    1904 solveForRandom(n,redun,1,1,12);
     2111decodeRandom(n,redun,1,1,12);
    19052112
    19062113redun=n-10;
    1907 solveForRandom(n,redun,1,1,24);
     2114decodeRandom(n,redun,1,1,24);
    19082115
    19092116int q=256; int n=150; int redun=n-10;
    19102117ring r=(q,a),x,dp;
    1911 solveForRandom(n,redun,1,1,26);
     2118decodeRandom(n,redun,1,1,26);
    19122119
    19132120
Note: See TracChangeset for help on using the changeset viewer.