# source:git/Singular/LIB/decodegb.lib@bf42f0

jengelh-datetimespielwiese
Last change on this file since bf42f0 was bf42f0, checked in by Hans Schoenemann <hannes@…>, 12 years ago
• Property mode set to `100644`
File size: 50.9 KB
Line
1///////////////////////////////////////////////////////////////////////////////
2version="\$Id\$";
3category="Coding theory";
4info="
5LIBRARY: decodegb.lib         Decoding and min distance of linear codes with GB
6AUTHOR:  Stanislav Bulygin,   bulygin@mathematik.uni-kl.de
7
8OVERVIEW:
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 origindeal 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 @ref{Decoding codes with Groebner bases}.
15 For the vanishing ideal computation the algorithm of Farr and Gao is
16 implemented.
17
18PROCEDURES:
19 sysCRHT(..);        generates the CRHT-ideal as in Cooper's philosophy
20 sysCRHTMindist(..); CRHT-ideal to find the minimum distance in the binary case
21 sysNewton(..);      generates the ideal with the generalized Newton identities
22 sysBin(..);         generates Bin system using Waring function
23 encode(x,g);        encodes given message x with the given generator matrix g
24 syndrome(h,c);      computes a syndrome w.r.t. the given check matrix
25 sysQE(..);          generates the system of quadratic equations for decoding
26 errorInsert(..);    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 decodeRandom(..); a procedure for manipulation with random codes
33 decodeCode(..);   a procedure for manipulation with the given code
34 vanishId(points);   computes the vanishing ideal for the given set of points
35 sysFL(..);          generates the Fitzgerald-Lax system
36 decodeRandomFL(..); manipulation with random codes via Fitzgerald-Lax
37
38
39KEYWORDS:  Cyclic code; Linear code; Decoding;
40           Minimum distance; Groebner bases, decodeGB
41";
42
43LIB "linalg.lib";
44LIB "brnoeth.lib";
45
46///////////////////////////////////////////////////////////////////////////////
47// creates a list result, where result[i]=i, 1<=i<=n
48static proc lis (int n)
49{
50 list result;
51 if (n<=0) {print("ERRORlis");}
52 for (int i=1; i<=n; i++)
53 {
54  result=result+list(i);
55 }
56 return(result);
57}
58
59///////////////////////////////////////////////////////////////////////////////
60// creates a list of all combinations without repititions of m objects out of n
61static proc combinations (int m, int n)
62{
63 list result;
64 if (m>n) {print("ERRORcombinations");}
65 if (m==n) {result[size(result)+1]=lis(m);return(result);}
66 if (m==0) {result[size(result)+1]=list();return(result);}
67 list temp=combinations(m-1,n-1);
68 for (int i=1; i<=size(temp); i++)
69 {
70  temp[i]=temp[i]+list(n);
71 }
72 result=combinations(m,n-1)+temp;
73 return(result);
74}
75
76
77///////////////////////////////////////////////////////////////////////////////
78// the polynomial for Sala's restrictions
79static proc p_poly(int n, int a, int b)
80{
81  poly f;
82  for (int i=0; i<=n-1; i++)
83  {
84    f=f+Z(a)^i*Z(b)^(n-1-i);
85  }
86  return(f);
87}
88
89///////////////////////////////////////////////////////////////////////////////
90
91proc sysCRHT (int n, list defset, int e, int q, int m, list #)
92"USAGE:   sysCRHT(n,defset,e,q,m,[k]); n,e,q,m,k are int, defset list of int's
93@format
94         - n length of the cyclic code,
95         - defset is a list representing the defining set,
96         - e the error-correcting capacity,
97         - q field size
98         - m degree extension of the splitting field,
99         - if k>0 additional equations representing the fact that every two
100         error positions are either different or at least one of them is zero
101@end format
102RETURN: the ring to work with the CRHT-ideal (with Sala's additions),
103        containig an ideal with name 'crht'
104THEORY:  Based on 'defset' of the given cyclic code, the procedure constructs
105         the corresponding Cooper-Reed-Heleseth-Truong ideal 'crht'. With its
106         help one can solve the decoding problem. For basics of the method @ref{Cooper philosophy}.
108EXAMPLE: example sysCRHT; shows an example
109"
110{
111  int r=size(defset);
112  ring @crht=(q,a),(Y(e..1),Z(1..e),X(r..1)),lp;
113  ideal crht;
114  int i,j;
115  poly sum;
116  int k;
117  if ( size(#) > 0)
118  {
119    k = #[1];
120  }
121
122  //---------------------- add check equations --------------------------
123  for (i=1; i<=r; i++)
124  {
125    sum=0;
126    for (j=1; j<=e; j++)
127    {
128      sum=sum+Y(j)*Z(j)^defset[i];
129    }
130    crht[i]=sum-X(i);
131  }
132
133  //--------------------- field equations on syndromes ------------------
134  for (i=1; i<=r; i++)
135  {
136    crht=crht,X(i)^(q^m)-X(i);
137  }
138
139  //------ restrictions on error-locations: n-th roots of unity ----------
140  for (i=1; i<=e; i++)
141  {
142    crht=crht,Z(i)^(n+1)-Z(i);
143  }
144
145  for (i=1; i<=e; i++)
146  {
147    crht=crht,Y(i)^(q-1)-1;
148  }
149
151  if ( k > 0 )
152
153  {
154    for (i=1; i<=e; i++)
155    {
156      for (j=i+1; j<=e; j++)
157      {
158        crht=crht,Z(i)*Z(j)*p_poly(n,i,j);
159      }
160    }
161  }
162  export crht;
163  return(@crht);
164}
165example
166{ "EXAMPLE:"; echo=2;
167  // binary cyclic [15,7,5] code with defining set (1,3)
168  intvec v = option(get);
169
170  list defset=1,3;           // defining set
171  int n=15;                  // length
172  int e=2;                   // error-correcting capacity
173  int q=2;                   // basefield size
174  int m=4;                   // degree extension of the splitting field
176
177  def A=sysCRHT(n,defset,e,q,m);
178  setring A;
179  A;                         // shows the ring we are working in
180  print(crht);               // the CRHT-ideal
181  option(redSB);
182  ideal red_crht=std(crht);  // reduced Groebner basis
183  print(red_crht);
184
185  //============================
186  A=sysCRHT(n,defset,e,q,m,sala);
187  setring A;
188  print(crht);                // CRHT-ideal with additional equations from Sala
189  option(redSB);
190  ideal red_crht=std(crht);   // reduced Groebner basis
191  print(red_crht);
192  red_crht[5];                // general error-locator polynomial for this code
193  option(set,v);
194}
195
196///////////////////////////////////////////////////////////////////////////////
197
198
199proc sysCRHTMindist (int n, list defset, int w)
200"USAGE:  sysCRHTMindist(n,defset,w);  n,w are int, defset is list of int's
201@format
202        - n length of the cyclic code,
203        - defset is a list representing the defining set,
204        - w is a candidate for the minimum distance
205@end format
206RETURN:  the ring to work with the Sala's ideal for the minimum distance
207         containing the ideal with name 'crht_md'
208THEORY:  Based on 'defset' of the given cyclic code, the procedure constructs
209         the corresponding Cooper-Reed-Heleseth-Truong ideal 'crht_md'. With
210         its help one can find minimum distance of the code in the binary
211         case. For basics of the method @ref{Cooper philosophy}.
212EXAMPLE: example sysCRHTMindist; shows an example
213"
214{
215  int r=size(defset);
216  ring @crht_md=2,Z(1..w),lp;
217  ideal crht_md;
218  int i,j;
219  poly sum;
220
221  //------------ add check equations --------------
222  for (i=1; i<=r; i++)
223  {
224    sum=0;
225    for (j=1; j<=w; j++)
226    {
227      sum=sum+Z(j)^defset[i];
228    }
229    crht_md[i]=sum;
230  }
231
232
233  //----------- locations are n-th roots of unity ------------
234  for (i=1; i<=w; i++)
235  {
236    crht_md=crht_md,Z(i)^n-1;
237  }
238
239  //------------ adding conditions on locations being different ------------
240  for (i=1; i<=w; i++)
241  {
242    for (j=i+1; j<=w; j++)
243    {
244      crht_md=crht_md,Z(i)*Z(j)*p_poly(n,i,j);
245    }
246  }
247
248  export crht_md;
249  return(@crht_md);
250}
251example
252{
253  "EXAMPLE:"; echo=2;
254  intvec v = option(get);
255  // binary cyclic [15,7,5] code with defining set (1,3)
256
257  list defset=1,3;             // defining set
258  int n=15;                    // length
259  int d=5;                     // candidate for the minimum distance
260
261  def A=sysCRHTMindist(n,defset,d);
262  setring A;
263  A;                           // shows the ring we are working in
264  print(crht_md);              // the Sala's ideal for mindist
265  option(redSB);
266  ideal red_crht_md=std(crht_md);
267  print(red_crht_md);          // reduced Groebner basis
268
269  option(set,v);
270}
271
272///////////////////////////////////////////////////////////////////////////////
273// slightly modified mod function
274static proc mod_ (int n, int m)
275{
276  if (n mod m==0) {return(m);}
277  if (n mod m!=0) {return(n mod m);}
278}
279
280///////////////////////////////////////////////////////////////////////////////
281
282proc sysNewton (int n, list defset, int t, int q, int m, list #)
283"USAGE:   sysNewton (n,defset,t,q,m,[tr]); n,t,q,m,tr int, defset is list int's
284@format
285         - n is length,
286         - defset is the defining set,
287         - t is the number of errors,
288         - q is basefield size,
289         - m is degree extension of the splitting field,
290         - if tr>0 it indicates that Newton identities in triangular
291           form should be constructed
292@end format
293RETURN:  the ring to work with the generalized Newton identities (in
294         triangular form if applicable) containing the ideal with name 'newton'
295THEORY:  Based on 'defset' of the given cyclic code, the procedure constructs
296         the corresponding ideal 'newton' with the generalized Newton
297         identities. With its help one can solve the decoding problem. For
298         basics of the method @ref{Generalized Newton identities}.
300EXAMPLE:  example sysNewton; shows an example
301"
302{
303 string s="ring @newton=("+string(q)+",a),(";
304 int i,j;
305 int flag;
306 int tr;
307
308 if (size(#)>0)
309 {
310  tr=#[1];
311 }
312
313 for (i=n; i>=1; i--)
314 {
315  for (j=1; j<=size(defset); j++)
316  {
317    flag=1;
318    if (i==defset[j])
319    {
320      flag=0;
321      break;
322    }
323  }
324  if (flag)
325  {
326    s=s+"S("+string(i)+"),";
327  }
328 }
329 s=s+"sigma(1.."+string(t)+"),";
330 for (i=size(defset); i>=2; i--)
331 {
332  s=s+"S("+string(defset[i])+"),";
333 }
334 s=s+"S("+string(defset[1])+")),lp;";
335
336 execute(s);
337
338 ideal newton;
339 poly sum;
340
341
342 //------------ generate generalized Newton identities ----------
343 if (tr)
344 {
345  for (i=1; i<=t; i++)
346  {
347    sum=0;
348    for (j=1; j<=i-1; j++)
349    {
350      sum=sum+sigma(j)*S(i-j);
351    }
352    newton=newton,S(i)+sum+number(i)*sigma(i);
353  }
354 } else
355 {
356  for (i=1; i<=t; i++)
357  {
358    sum=0;
359    for (j=1; j<=t; j++)
360    {
361      sum=sum+sigma(j)*S(mod_(i-j,n));
362    }
363    newton=newton,S(i)+sum;
364  }
365 }
366 for (i=1; i<=n-t; i++)
367 {
368  sum=0;
369  for (j=1; j<=t; j++)
370  {
371    sum=sum+sigma(j)*S(t+i-j);
372  }
373  newton=newton,S(t+i)+sum;
374 }
375
376 //----------- add field equations on sigma's --------------
377 for (i=1; i<=t; i++)
378 {
379  newton=newton,sigma(i)^(q^m)-sigma(i);
380 }
381
382 //----------- add conjugacy relations ------------------
383 for (i=1; i<=n; i++)
384 {
385  newton=newton,S(i)^q-S(mod_(q*i,n));
386 }
387 newton=simplify(newton,2);
388 export newton;
389 return(@newton);
390}
391example
392{
393     "EXAMPLE:";  echo = 2;
394     // Newton identities for a binary 3-error-correcting cyclic code of
395     //length 31 with defining set (1,5,7)
396
397     int n=31;          // length
398     list defset=1,5,7; //defining set
399     int t=3;           // number of errors
400     int q=2;           // basefield size
401     int m=5;           // degree extension of the splitting field
402     int tr=1;          // indicator of triangular form of Newton identities
403
404
405     def A=sysNewton(n,defset,t,q,m);
406     setring A;
407     A;                 // shows the ring we are working in
408     print(newton);     // generalized Newton identities
409
410     //===============================
411     A=sysNewton(n,defset,t,q,m,tr);
412     setring A;
413     print(newton);     // generalized Newton identities in triangular form
414}
415
416///////////////////////////////////////////////////////////////////////////////
417// forms a list of special combinations needed for computation of Waring's
418//function
419static proc combinations_sum (int m, int n)
420{
421 list result;
422 list comb=combinations(m-1,n+m-1);
423 int i,j,flag,count;
424 list interm=comb;
425 for (i=1; i<=size(comb); i++)
426 {
427  interm[i][1]=comb[i][1]-1;
428  for (j=2; j<=m-1; j++)
429  {
430   interm[i][j]=comb[i][j]-comb[i][j-1]-1;
431  }
432  interm[i][m]=n+m-comb[i][m-1]-1;
433  flag=1;
434  count=2;
435  while ((flag)&&(count<=m))
436  {
437   if (interm[i][count] mod count != 0) {flag=0;}
438   count++;
439  }
440  if (flag)
441  {
442   for (j=2; j<=m; j++)
443   {
444    interm[i][j]=interm[i][j] div j;
445   }
446   result[size(result)+1]=interm[i];
447  }
448 }
449 return(result);
450}
451
452///////////////////////////////////////////////////////////////////////////////
453//if n=q^e*m, m and n are coprime, then return e
454static proc exp_count (int n, int q)
455{
456 int flag=1;
457 int result=0;
458 while(flag)
459 {
460  if (n mod q != 0) {flag=0;}
461   else {n=n div q; result++;}
462 }
463 return(result);
464}
465
466///////////////////////////////////////////////////////////////////////////////
467
468
469proc sysBin (int v, list Q, int n, list #)
470"USAGE:    sysBin (v,Q,n,[odd]);  v,n,odd are int, Q is list of int's
471@format
472          - v a number if errors,
473          - Q is a defining set of the code,
474          - n the length,
475          - odd is an additional parameter: if
476           set to 1, then the defining set is enlarged by odd elements,
477           which are 2^(some power)*(some elment in the def.set) mod n
478@end format
479RETURN:    the ring with the resulting system called 'bin'
480THEORY:  Based on Q of the given cyclic code, the procedure constructs
481         the corresponding ideal 'bin' with the use of the Waring function.
482         With its help one can solve the decoding problem.
483         For basics of the method @ref{Generalized Newton identities}.
485EXAMPLE:   example sysBin; shows an example
486"
487{
488 int odd;
489 if (size(#)>0)
490 {
491  odd=#[1];
492 }
493
494 //ring r=2,(sigma(1..v),S(1..n)),(lp(v),dp(n));
495 ring r=2,(S(1..n),sigma(1..v)),lp;
496 list cyclot;
497 ideal result;
498 int i,j,k,s;
499 list comb;
500 poly sum_, mon;
501 int count1, count2, upper, coef_, flag, gener;
502 list Q_update;
503 if (odd==1)
504 {
505  for (i=1; i<=n; i++)
506  {
507   cyclot[i]=0;
508  }
509  for (i=1; i<=size(Q); i++)
510  {
511   flag=1;
512   gener=Q[i];
513   while(flag)
514   {
515    cyclot[gener]=1;
516    gener=2*gener mod n;
517    if (gener == Q[i]) {flag=0;}
518   }
519  }
520  for (i=1; i<=n; i++)
521  {
522   if ((cyclot[i] == 1)&&(i mod 2 == 1)) {Q_update[size(Q_update)+1]=i;}
523  }
524 }
525 else
526 {
527  Q_update=Q;
528 }
529
530 //---- form polynomials for the Bin system via Waring's function ---------
531 for (i=1; i<=size(Q_update); i++)
532 {
533  comb=combinations_sum(v,Q_update[i]);
534  sum_=0;
535  for (k=1; k<=size(comb); k++)
536  {
537   upper=0;
538   for (j=1; j<=v; j++)
539   {
540    upper=upper+comb[k][j];
541   }
542   count1=0;
543   for (j=2; j<=upper-1; j++)
544   {
545    count1=count1+exp_count(j,2);
546   }
547   count1=count1+exp_count(Q_update[i],2);
548   count2=0;
549   for (j=1; j<=v; j++)
550   {
551    for (s=2; s<=comb[k][j]; s++)
552    {
553     count2=count2+exp_count(s,2);
554    }
555   }
556   if (count1<count2) {print("ERRORsysBin");}
557   if (count1>count2) {coef_=0;}
558   if (count1 == count2) {coef_=1;}
559   mon=1;
560   for (j=1; j<=v; j++)
561   {
562    mon=mon*sigma(j)^(comb[k][j]);
563   }
564   sum_=sum_+coef_*mon;
565  }
566  result=result,S(Q_update[i])-sum_;
567 }
568 ideal bin=simplify(result,2);
569 export bin;
570 return(r);
571}
572example
573{
574     "EXAMPLE:";  echo = 2;
575     // [31,16,7] quadratic residue code
576     list l=1,5,7,9,19,25;
577     // we do not need even synromes here
578     def A=sysBin(3,l,31);
579     setring A;
580     print(bin);
581}
582
583///////////////////////////////////////////////////////////////////////////////
584
585proc encode (matrix x, matrix g)
586"USAGE:  encode (x, g);  x a row vector (message), and g a generator matrix
587RETURN:  corresponding codeword
588EXAMPLE: example encode; shows an example
589"
590{
591 if (nrows(x)>1) {print("ERRORencode1!");}
592 if (ncols(x)!=nrows(g)) {print("ERRORencode2!");}
593 return(x*g);
594}
595example
596{
597     "EXAMPLE:";  echo = 2;
598     ring r=2,x,dp;
599     matrix x[1][4]=1,0,1,0;
600     matrix g[4][7]=1,0,0,0,0,1,1,
601                    0,1,0,0,1,0,1,
602                    0,0,1,0,1,1,1,
603                    0,0,0,1,1,1,0;
604     //encode x with the generator matrix g
605     print(encode(x,g));
606}
607
608///////////////////////////////////////////////////////////////////////////////
609
610proc syndrome (matrix h, matrix c)
611"USAGE:  syndrome (h, c);  h a check matrix, c a row vector (codeword)
612RETURN:  corresponding syndrome
613EXAMPLE: example syndrome; shows an example
614"
615{
616 if (nrows(c)>1) {print("ERRORsyndrome1!");}
617 if (ncols(c)!=ncols(h)) {print("ERRORsyndrome2!");}
618 return(h*transpose(c));
619}
620example
621{
622     "EXAMPLE:";  echo = 2;
623     ring r=2,x,dp;
624     matrix x[1][4]=1,0,1,0;
625     matrix g[4][7]=1,0,0,0,0,1,1,
626                    0,1,0,0,1,0,1,
627                    0,0,1,0,1,1,1,
628                    0,0,0,1,1,1,0;
629     //encode x with the generator matrix g
630     matrix c=encode(x,g);
631     // disturb
632     c[1,3]=0;
633     //compute syndrome
634     //corresponding check matrix
635     matrix check[3][7]=1,0,0,1,1,0,1,0,1,0,1,0,1,1,0,0,1,0,1,1,1;
636     print(syndrome(check,c));
637     c[1,3]=1;
638     //now c is a codeword
639     print(syndrome(check,c));
640}
641
642///////////////////////////////////////////////////////////////////////////////
643// (coordinatewise) star product of two vectors
644static proc star(matrix m, int i, int j)
645{
646 matrix result[ncols(m)][1];
647 for (int k=1; k<=ncols(m); k++)
648 {
649  result[k,1]=m[i,k]*m[j,k];
650 }
651 return(result);
652}
653
654///////////////////////////////////////////////////////////////////////////////
655
656proc sysQE(matrix check, matrix y, int t, list #)
657"USAGE:   sysQE(check,y,t,[fieldeq,formal]);check,y matrix;t,fieldeq,formal int
658@format
659        - check is a parity check matrix of the code
660        - y is a received word,
661        - t the number of errors to be corrected,
662        - if fieldeq=1, then field equations are added,
663        - if formal=0, field equations on (known) syndrome variables
664          are not added, in order to add them (note that the exponent should
665          be equal to the number of elements in the INITIAL alphabet) one
666          needs to set formal>0 for the exponent
667@end format
668RETURN:   the ring to work with together with the resulting system called 'qe'
669THEORY:  Based on 'check' of the given linear code, the procedure constructs
670         the corresponding ideal that gives an opportunity to compute
671         unknown syndrome of the received word y. After computing the unknown
672         syndromes one is able to solve the decoding problem.
673         For basics of the method @ref{Decoding method based on quadratic equations}.
675EXAMPLE:  example sysQE; shows an example
676"
677{
678 int fieldeq;
679 int formal;
680 if (size(#)>0)
681 {
682  fieldeq=#[1];
683 }
684 if (size(#)>1)
685 {
686  formal=#[2];
687 }
688
689 def br=basering;
690 list rl=ringlist(br);
691
692 int red=nrows(check);
693 int n=ncols(check);
694 int q=rl[1][1];
695
696 if (formal==0)
697 {
698  ring work=(q,a),(V(1..t),U(1..n)),dp;
699 } else
700 {
701  ring work=(q,a),(V(1..t),U(1..n),s(1..red)),(dp(t),lp(n),dp(red));
702 }
703
704 matrix check=imap(br,check);
705 matrix y=imap(br,y);
706
707 matrix h_full=genMDSMat(n,a);
708 matrix h=submat(h_full,1..red,1..n);
709 if (nrows(y)!=1) {print("ERROR1Pell");}
710 if (ncols(y)!=n) {print("ERROR2Pell");}
711
712 ideal result;
713
714 list c;
715 list a;
716 list tmp,tmp2;
717 int i,j,l,k;
718 number sum,prod,sig;
719        poly sum1,sum2,sum3;
720 for (i=1; i<=n; i++)
721 {
722  c[i]=tmp;
723 }
724
725 matrix transf=inverse(transpose(h_full));
726
727 //------ expression matrix of check vectors w.r.t. the MDS basis -----------
728 for (i=1; i<=red ; i++)
729 {
730  a[i]=transpose(submat(check,i..i,1..n));
731  a[i]=transf*a[i];
732 }
733
734 //----------- compute the structure constants ------------------------
735 matrix te[n][1];
736 for (i=1; i<=n; i++)
737 {
738  for (j=1; j<=t+1; j++)
739  {
740   if ((j<i)&&(i<=t+1)) {c[i][j]=c[j][i];}
741   else
742   {
743    if (i+j<=n+1)
744    {
745     c[i][j]=te;
746     c[i][j][i+j-1,1]=1;
747    }
748    else
749    {
750     c[i][j]=star(h_full,i,j);
751     c[i][j]=transf*c[i][j];
752    }
753   }
754  }
755 }
756
757
758 if (formal==0)
759 {
760  matrix s[red][1]=syndrome(check,y);
761  for (j=1; j<=red; j++)
762  {
763   sum1=0;
764   for (l=1; l<=n; l++)
765   {
766    sum1=sum1+a[j][l,1]*U(l);
767   }
768   result=result,sum1-s[j,1];
769  }
770 } else
771 {
772  for (j=1; j<=red; j++)
773  {
774   sum1=0;
775   for (l=1; l<=n; l++)
776   {
777    sum1=sum1+a[j][l,1]*U(l);
778   }
779   result=result,sum1-s(j);
780  }
781  for (j=1; j<=red; j++)
782  {
783     result=result,s(j)^(formal)-s(j);
784  }
785 }
786 if (fieldeq)
787 {
788  for (i=1; i<=n; i++)
789  {
790   result=result,U(i)^q-U(i);
791  }
792  for (j=1; j<=t; j++)
793  {
794     result=result,V(j)^q-V(j);
795  }
796 }
797
798 //----- form the quadratic equations according to the theory -----------
799 for (i=1; i<=n; i++)
800 {
801  sum1=0;
802  for (j=1; j<=t; j++)
803  {
804   sum2=0;
805   for (l=1; l<=n; l++)
806   {
807    sum2=sum2+c[i][j][l,1]*U(l);
808   }
809   sum1=sum1+sum2*V(j);
810  }
811  sum3=0;
812  for (l=1; l<=n; l++)
813  {
814   sum3=sum3+c[i][t+1][l,1]*U(l);
815  }
816  result=result,sum1-sum3;
817 }
818
819 result=simplify(result,2);
820
821 ideal qe=result;
822 export qe;
823 return(work);
824}
825example
826{
827     "EXAMPLE:";  echo = 2;
828     intvec v = option(get);
829
830     //correct 2 errors in [7,3] 8-ary code RS code
831     int t=2; int q=8; int n=7; int redun=4;
832     ring r=(q,a),x,dp;
833     matrix h_full=genMDSMat(n,a);
834     matrix h=submat(h_full,1..redun,1..n);
835     matrix g=dual_code(h);
836     matrix x[1][3]=0,0,1,0;
837     matrix y[1][7]=encode(x,g);
838
839     //disturb with 2 errors
840     matrix rec[1][7]=errorInsert(y,list(2,4),list(1,a));
841
842     //generate the system
843     def A=sysQE(h,rec,t);
844     setring A;
845     print(qe);
846
847     //let us decode
848     option(redSB);
849     ideal sys_qe=std(qe);
850     print(sys_qe);
851
852     option(set,v);
853}
854
855///////////////////////////////////////////////////////////////////////////////
856
857proc errorInsert(matrix y, list pos, list val)
858"USAGE:  errorInsert(y,pos,val); y is matrix, pos,val are list of int's
859@format
860        - y is a (code) word,
861        - pos = positions where errors occured,
862        - val = their corresponding values
863@end format
865EXAMPLE: example errorInsert; shows an example
866"
867{
868 matrix result[1][ncols(y)]=y;
869 if (size(pos)!=size(val)) {print("ERRORerror");}
870 for (int i=1; i<=size(pos); i++)
871 {
872  result[1,pos[i]]=y[1,pos[i]]+val[i];
873 }
874 return(result);
875}
876example
877{
878     "EXAMPLE:";  echo = 2;
879     //correct 2 errors in [7,3] 8-ary code RS code
880     int t=2; int q=8; int n=7; int redun=4;
881     ring r=(q,a),x,dp;
882     matrix h_full=genMDSMat(n,a);
883     matrix h=submat(h_full,1..redun,1..n);
884     matrix g=dual_code(h);
885     matrix x[1][3]=0,0,1,0;
886     matrix y[1][7]=encode(x,g);
887     print(y);
888
889     //disturb with 2 errors
890     matrix rec[1][7]=errorInsert(y,list(2,4),list(1,a));
891     print(rec);
892     print(rec-y);
893}
894
895///////////////////////////////////////////////////////////////////////////////
896
897proc errorRand(matrix y, int num, int e)
898"USAGE:    errorRand(y, num, e); y is matrix, num,e are int
899@format
900          - y is a (code) word,
901          - num is the number of errors,
902          - e is an extension degree (if one wants values to be from GF(p^e))
903@end format
905EXAMPLE:   example errorRand; shows an example
906"
907{
908 matrix result[1][ncols(y)]=y;
909 int i,j, flag, temp;
910 list pos, val;
911 matrix tempnum;
912
913 for (i=1; i<=num; i++)
914 {
915  while(1)
916  {
917   temp=random(1,ncols(y));
918   flag=1;
919   for (j=1; j<=size(pos); j++)
920   {
921    if (temp==pos[j]) {flag=0;}
922   }
923   if (flag) {pos[i]=temp;break;}
924  }
925 }
926
927 for (i=1; i<=num; i++)
928 {
929  flag=1;
930  while(flag)
931  {
932   tempnum=randomvector(1,e);
933   if (tempnum!=0) {flag=0;}
934  }
935  val[i]=tempnum;
936 }
937
938 for (i=1; i<=size(pos); i++)
939 {
940  result[1,pos[i]]=y[1,pos[i]]+val[i];
941 }
942 return(result);
943}
944example
945{
946  "EXAMPLE:";  echo = 2;
947     //correct 2 errors in [7,3] 8-ary code RS code
948     int t=2; int q=8; int n=7; int redun=4;
949     ring r=(q,a),x,dp;
950     matrix h_full=genMDSMat(n,a);
951     matrix h=submat(h_full,1..redun,1..n);
952     matrix g=dual_code(h);
953     matrix x[1][3]=0,0,1,0;
954     matrix y[1][7]=encode(x,g);
955
956     //disturb with 2 random errors
957     matrix rec[1][7]=errorRand(y,2,3);
958     print(rec);
959     print(rec-y);
960}
961
962///////////////////////////////////////////////////////////////////////////////
963
964proc randomCheck(int m, int n, int e)
965"USAGE:    randomCheck(m, n, e); m,n,e are int
966@format
967          - m x n are dimensions of the matrix,
968          - e is an extension degree (if one wants values to be from GF(p^e))
969@end format
970RETURN:    random check matrix
971EXAMPLE:   example randomCheck; shows an example
972"
973{
974 matrix result[m][n];
975 matrix rand[m][n-m];
976 int i,j;
977 matrix temp;
978 for (i=1; i<=m; i++)
979 {
980  temp=randomvector(n-m,e);
981  for (j=1; j<=n-m; j++)
982  {
983   rand[i,j]=temp[j,1];
984  }
985 }
986 result=concat(rand,unitmat(m));
987 return(result);
988}
989example
990{
991  "EXAMPLE:";  echo = 2;
992     int redun=5; int n=15;
993     ring r=2,x,dp;
994
995     //generate random check matrix for a [15,5] binary code
996     matrix h=randomCheck(redun,n,1);
997     print(h);
998
999     //corresponding generator matrix
1000     matrix g=dual_code(h);
1001     print(g);
1002}
1003
1004///////////////////////////////////////////////////////////////////////////////
1005
1006proc genMDSMat(int n, number a)
1007"USAGE:   genMDSMat(n, a); n is int, a is number
1008@format
1009        - n x n are dimensions of the MDS matrix,
1010        - a is a primitive element of the field.
1011@end format
1012NOTE:   An MDS matrix is constructed in the following way. We take 'a' to be a
1013        generator of the multiplicative group of the field. Then we construct
1014        the Vandermonde matrix with this 'a'.
1015ASSUME:   extension field should already be defined
1016RETURN:   a matrix with the MDS property.
1018EXAMPLE:  example genMDSMat; shows an example
1019"
1020{
1021 int i,j;
1022 matrix result[n][n];
1023 for (i=0; i<=n-1; i++)
1024 {
1025  for (j=0; j<=n-1; j++)
1026  {
1027   result[j+1,i+1]=(a^i)^j;
1028  }
1029 }
1030 return(result);
1031}
1032example
1033{
1034     "EXAMPLE:";  echo = 2;
1035     int q=16; int n=15;
1036     ring r=(q,a),x,dp;
1037
1038     //generate an MDS (Vandermonde) matrix
1039     matrix h_full=genMDSMat(n,a);
1040     print(h_full);
1041}
1042
1043///////////////////////////////////////////////////////////////////////////////
1044
1045
1046proc mindist (matrix check)
1047"USAGE:  mindist (check, q); check matrix, q int
1048@format
1049        - check is a check matrix,
1050        - q is the field size
1051@end format
1052RETURN:  minimum distance of the code
1053EXAMPLE: example mindist; shows an example
1054"
1055{
1056 intvec vopt = option(get);
1057
1058 int n=ncols(check); int redun=nrows(check); int t=redun+1;
1059
1060 def br=basering;
1061 list rl=ringlist(br);
1062 int q=rl[1][1];
1063
1064 ring work=(q,a),(V(1..t),U(1..n)),dp;
1065 matrix check=imap(br,check);
1066
1067 ideal temp;
1068 int count=1;
1069 int flag=1;
1070 int flag2;
1071 int i;
1072 matrix z[1][n];
1073 option(redSB);
1074 def A=sysQE(check,z,count);
1075
1076 //proceed with solving the system w.r.t zero vector until some solutions
1077 //are found
1078 while (flag)
1079 {
1080    A=sysQE(check,z,count);
1081    setring A;
1082    ideal temp=qe;
1083    temp=std(temp);
1084    flag2=1;
1085    setring work;
1086    temp=imap(A,temp);
1087    for (i=1; i<=n; i++)
1088    {
1089      if
1090        (temp[i]!=U(n-i+1))
1091        {
1092          flag2=0;
1093        }
1094    }
1095    if (!flag2)
1096    {
1097      flag=0;
1098    }
1099    else
1100    {
1101      count++;
1102    }
1103 }
1104 int result=count;
1105
1106 option(set,vopt);
1107 return(result);
1108}
1109example
1110{
1111     "EXAMPLE:";  echo = 2;
1112     //determine a minimum distance for a [7,3] binary code
1113     int q=8; int n=7; int redun=4; int t=redun+1;
1114     ring r=(q,a),x,dp;
1115
1116     //generate random check matrix
1117     matrix h=randomCheck(redun,n,1);
1118     print(h);
1119     int l=mindist(h);
1120     l;
1121}
1122
1123///////////////////////////////////////////////////////////////////////////////
1124
1125proc decode(matrix check, matrix rec)
1126"USAGE:    decode(check, rec, t); check, rec matrix, t int
1127@format
1128          - check is the check matrix of the code,
1129          - rec is a received word,
1130          - t is an upper bound for the number of errors one wants to correct
1131@end format
1132NOTE:     The method described in @ref{Decoding method based on quadratic equations}
1133          is used for decoding.
1134ASSUME:   Errors in rec should be correctable, otherwise the output is
1135          unpredictable
1136RETURN:   a codeword that is closest to rec
1137EXAMPLE:  example decode; shows an example
1138"
1139{
1140 intvec vopt = option(get);
1141
1142 def br=basering;
1143 int n=ncols(check);
1144
1145 int count=1;
1146 def A=sysQE(check,rec,count);
1147 while(1)
1148 {
1149  A=sysQE(check,rec,count);
1150  setring A;
1151  matrix h_full=genMDSMat(n,a);
1152  matrix rec=imap(br,rec);
1153  option(redSB);
1154  ideal qe_red=std(qe);
1155  if (qe_red[1]!=1)
1156  {
1157    break;
1158  }
1159  else
1160  {
1161    count++;
1162  }
1163  setring br;
1164 }
1165
1166 setring A;
1167
1168 //obtain a codeword
1169 //this works only if our code is indeed can correct these errors
1170 matrix syn[n][1];
1171 for (int i=1; i<=n; i++)
1172 {
1174 }
1175
1176 matrix real_syn=inverse(h_full)*syn;
1177 setring br;
1178 matrix real_syn=imap(A,real_syn);
1179
1180 option(set,vopt);
1181 return(rec-transpose(real_syn));
1182}
1183example
1184{
1185     "EXAMPLE:";  echo = 2;
1186     //correct 1 error in [15,7] binary code
1187     int t=1; int q=16; int n=15; int redun=10;
1188     ring r=(q,a),x,dp;
1189
1190     //generate random check matrix
1191     matrix h=randomCheck(redun,n,1);
1192     matrix g=dual_code(h);
1193     matrix x[1][n-redun]=0,0,1,0,1,0,1;
1194     matrix y[1][n]=encode(x,g);
1195     print(y);
1196
1197     // find out the minimum distance of the code
1198     list l=mindist(h);
1199
1200     //disturb with errors
1201     "Correct ",(l[1]-1) div 2," errors";
1202     matrix rec[1][n]=errorRand(y,(l[1]-1) div 2,1);
1203     print(rec);
1204
1205     //let us decode
1206     matrix dec_word=decode(h,rec);
1207     print(dec_word);
1208}
1209
1210///////////////////////////////////////////////////////////////////////////////
1211
1212
1213proc decodeRandom(int n, int redun, int ncodes, int ntrials, list #)
1214"USAGE:    decodeRandom(redun,q,ncodes,ntrials,[e]); all parameters int
1215@format
1216          - redun is a redundabcy of a (random) code,
1217          - q is the field size,
1218          - ncodes is the number of random codes to be processed,
1219          - ntrials is the number of received vectors per code to be corrected
1220          - If e is given it sets the correction capacity explicitly. It
1221          should be used in case one expects some lower bound,
1222          otherwise the procedure tries to compute the real minimum distance
1223          to find out the error-correction capacity
1224@end format
1225RETURN:    nothing;
1226EXAMPLE:   example decodeRandom; shows an example
1227"
1228{
1229 intvec vopt = option(get);
1230
1231 int i,j;
1232 matrix h;
1233 int dist, t;
1234 ideal sys;
1235 int tmp;
1236 int e;
1237 if (size(#)>0)
1238 {
1239  e=#[1];
1240 }
1241
1242 option(redSB);
1243 def br=basering;
1244 matrix h_full=genMDSMat(n,a);
1245 matrix z[1][ncols(h_full)];
1246
1247 //------------------ determine error-correction capacity -------------------
1248 for (i=1; i<=ncodes; i++)
1249 {
1250  setring br;
1251  h=randomCheck(redun,n,1);
1252  "check matrix:";
1253  print(h);
1254  if (e>0)
1255  {
1256     t=e;
1257  } else {
1258     tmp=mindist(h);
1259     dist=tmp;
1260     printf("d= %p",dist);
1261     t=(dist-1) div 2;
1262  }
1263
1264  //------------- generate the template system ----------------------
1265  def A=sysQE(h,z,t);
1266  setring A;
1267  matrix word,y,rec;
1268  ideal sys2,sys3;
1269  matrix h=imap(br,h);
1270  matrix g=dual_code(h);
1271  ideal sys=qe;
1272  print("The system is generated");
1273
1274  //------ modify the template according to every received word --------------
1275  for (j=1; j<=ntrials; j++)
1276  {
1277   word=randomvector(n-redun,1);
1278   y=encode(transpose(word),g);
1279   print("Codeword:");
1280   print(y);
1281   rec=errorRand(y,t,1);
1283   print(rec);
1285   option(redSB);
1286   sys3=std(sys2);
1287   print("The Groebenr basis of the QE system:");
1288   print(sys3);
1289  }
1290  kill A;
1291  option(set,vopt);
1292 }
1293}
1294example
1295{
1296     "EXAMPLE:";  echo = 2;
1297     int q=32; int n=25; int redun=n-11; int t=redun+1;
1298     ring r=(q,a),x,dp;
1299
1300     // correct 2 errors in 2 random binary codes, 3 trials each
1301     decodeRandom(n,redun,2,3,2);
1302}
1303
1304///////////////////////////////////////////////////////////////////////////////
1305
1306
1307proc decodeCode(matrix check, int ntrials, list #)
1308"USAGE:     decodeCode(check, ntrials, [e]); check matrix, ntrials,e int
1309@format
1310           - check is a parity check matrix for the code,
1311           - ntrials is the number of received vectors per code to be
1312           corrected.
1313           - If e is given it sets the correction capacity explicitly. It
1314           should be used in case one expects some lower bound,
1315           otherwise the procedure tries to compute the real minimum distance
1316           to find out the error-correction capacity
1317@end format
1318RETURN:     nothing;
1319EXAMPLE:    example decodeCode; shows an example
1320"
1321{
1322 intvec vopt = option(get);
1323
1324 int n=ncols(check);
1325 int redun=nrows(check);
1326 int i,j;
1327 matrix h;
1328 int dist, t;
1329 ideal sys;
1330 int tmp;
1331 int e;
1332 if (size(#)>0)
1333 {
1334  e=#[1];
1335 }
1336
1337 option(redSB);
1338 def br=basering;
1339 matrix h_full=genMDSMat(n,a);
1340 matrix z[1][ncols(h_full)];
1341 setring br;
1342 h=check;
1343 "check matrix:";
1344 print(h);
1345
1346 //------------------ determine error-correction capacity -------------------
1347 if (e>0)
1348 {
1349    t=e;
1350 } else {
1351   tmp=mindist(h);
1352   dist=tmp;
1353   printf("d= %p",dist);
1354   t=(dist-1) div 2;
1355 }
1356
1357 //------------- generate the template system ----------------------
1358 def A=sysQE(h,z,t);
1359 setring A;
1360 matrix word,y,rec;
1361 ideal sys2,sys3;
1362 matrix h=imap(br,h);
1363 matrix g=dual_code(h);
1364 ideal sys=qe;
1365 print("The system is generated");
1366
1367 //--- modify the template according to every received word ---------------
1368 for (j=1; j<=ntrials; j++)
1369 {
1370   word=randomvector(n-redun,1);
1371   y=encode(transpose(word),g);
1372   print("Codeword:");
1373   print(y);
1374   rec=errorRand(y,t,1);
1376   print(rec);
1378   option(redSB);
1379   sys3=std(sys2);
1380   print("Groebner basis of the QE system:");
1381   print(sys3);
1382 }
1383
1384 option(set,vopt);
1385}
1386example
1387{
1388     "EXAMPLE:";  echo = 2;
1389     int q=32; int n=25; int redun=n-11; int t=redun+1;
1390     ring r=(q,a),x,dp;
1391     matrix check=randomCheck(redun,n,1);
1392
1393     // correct 2 errors in using the code above, 3 trials
1394     decodeCode(check,3,2);
1395}
1396
1397
1398///////////////////////////////////////////////////////////////////////////////
1399// adding syndrome values to the template system
1400static proc add_synd (matrix rec, matrix check, int redun, ideal sys)
1401{
1402     ideal result=sys;
1403     matrix s[redun][1]=syndrome(check,rec);
1404     for (int i=1; i<=redun; i++)
1405
1406     {
1407          result[i]=result[i]-s[i,1];
1408     }
1409     return(result);
1410}
1411
1412///////////////////////////////////////////////////////////////////////////////
1413// evaluate a polynomial at a given point
1414static proc ev (poly f, matrix p)
1415{
1416     if (ncols(p)>1) {ERROR("not a column vector");};
1417     int m=size(p);
1418     poly temp=f;
1419     for (int i=1; i<=m; i++)
1420     {
1421          temp=subst(temp,var(i),p[i,1]);
1422     }
1423     return(number(temp));
1424}
1425
1426///////////////////////////////////////////////////////////////////////////////
1427// return index of an element in the ideal where it does not vanish at the
1428//given point
1429static proc find_index (ideal G, matrix p)
1430{
1431     if (ncols(p)>1) {ERROR("not a column vector");};
1432     int i=1;
1433     int n=size(G);
1434     while(i<=n)
1435     {
1436          if (ev(G[i],p)!=0) {return(i);}
1437          i++;
1438     }
1439     return(-1);
1440}
1441
1442///////////////////////////////////////////////////////////////////////////////
1443// convert ideal to list
1444static proc ideal2list (ideal id)
1445{
1446     list l;
1447     for (int i=1; i<=size(id); i++)
1448     {
1449          l[i]=id[i];
1450     }
1451     return(l);
1452}
1453
1454///////////////////////////////////////////////////////////////////////////////
1455// convert list to ideal
1456static proc list2ideal (list l)
1457{
1458     ideal id;
1459     for (int i=1; i<=size(l); i++)
1460     {
1461          id[i]=l[i];
1462     }
1463     return(id);
1464}
1465
1466///////////////////////////////////////////////////////////////////////////////
1467// check whether given polynomial is divisible by some leading monomial of the
1468//ideal
1469static proc divisible (poly m, ideal G)
1470{
1471     for (int i=1; i<=size(G); i++)
1472     {
1474     }
1475     return(0);
1476}
1477
1478///////////////////////////////////////////////////////////////////////////////
1479
1480proc vanishId (list points)
1481"USAGE:  vanishId (points); point is a list of matrices
1482        'points' is a list of points for which the vanishing ideal is to be
1483        constructed
1484RETURN:  Vanishing ideal corresponding to the given set of points
1485EXAMPLE: example vanishId; shows an example
1486"
1487{
1488     int m=size(points[1]);
1489     int n=size(points);
1490
1491     ideal G=1;
1492     int i,k,j;
1493     list temp;
1494     poly h,cur;
1495
1496     //------------- proceed according to Farr-Gao algorithm ----------------
1497     for (k=1; k<=n; k++)
1498     {
1499          i=find_index(G,points[k]);
1500          cur=G[i];
1501          for(j=i+1; j<=size(G); j++)
1502          {
1503               G[j]=G[j]-ev(G[j],points[k])/ev(G[i],points[k])*G[i];
1504          }
1505          G=simplify(G,2);
1506          temp=ideal2list(G);
1507          temp=delete(temp,i);
1508          G=list2ideal(temp);
1509          for (j=1; j<=m; j++)
1510          {
1512               {
1513                    attrib(G,"isSB",1);
1514                    h=NF((var(j)-points[k][j,1])*cur,G);
1515                    temp=ideal2list(G);
1516                    temp=insert(temp,h);
1517                    G=list2ideal(temp);
1518                    G=sort(G)[1];
1519               }
1520          }
1521     }
1522     attrib(G,"isSB",1);
1523     return(G);
1524}
1525example
1526{
1527     "EXAMPLE:";  echo = 2;
1528      ring r=3,(x(1..3)),dp;
1529
1530     //generate all 3-vectors over GF(3)
1531     list points=pointsGen(3,1);
1532
1533     list points2=convPoints(points);
1534
1535     //grasps the first 11 points
1536     list p=graspList(points2,1,11);
1537     print(p);
1538
1539     //construct the vanishing ideal
1540     ideal id=vanishId(p);
1541     print(id);
1542}
1543
1544///////////////////////////////////////////////////////////////////////////////
1545// construct the list of all vectors of length m with elements in p^e, where p
1546// is theharacteristic
1547proc pointsGen (int m, int e)
1548{
1549     if (e>1)
1550     {
1551     list result;
1552     int count=1;
1553     int i,j;
1554     list l=ringlist(basering);
1555     int charac=l[1][1];
1556     number a=par(1);
1557     list tmp;
1558     for (i=1; i<=charac^(e*m); i++)
1559     {
1560          result[i]=tmp;
1561     }
1562     if (m==1)
1563     {
1564          result[count][m]=0;
1565          count++;
1566          for (j=1; j<=charac^(e)-1; j++)
1567          {
1568               result[count][m]=a^j;
1569               count++;
1570          }
1571          return(result);
1572     }
1573     list prev=pointsGen(m-1,e);
1574     for (i=1; i<=size(prev); i++)
1575     {
1576          result[count]=prev[i];
1577          result[count][m]=0;
1578          count++;
1579          for (j=1; j<=charac^(e)-1; j++)
1580          {
1581               result[count]=prev[i];
1582               result[count][m]=a^j;
1583               count++;
1584          }
1585     }
1586     return(result);
1587     }
1588
1589     if (e==1)
1590     {
1591     list result;
1592     int count=1;
1593     int i,j;
1594     list l=ringlist(basering);
1595     int charac=l[1][1];
1596     list tmp;
1597     for (i=1; i<=charac^m; i++)
1598     {
1599          result[i]=tmp;
1600     }
1601     if (m==1)
1602     {
1603          for (j=0; j<=charac-1; j++)
1604          {
1605               result[count][m]=number(j);
1606               count++;
1607          }
1608          return(result);
1609     }
1610     list prev=pointsGen(m-1,e);
1611     for (i=1; i<=size(prev); i++)
1612     {
1613          for (j=0; j<=charac-1; j++)
1614          {
1615               result[count]=prev[i];
1616               result[count][m]=number(j);
1617               count++;
1618          }
1619     }
1620     return(result);
1621     }
1622
1623}
1624
1625///////////////////////////////////////////////////////////////////////////////
1626// convert list to a column vector
1627static proc list2vec (list l)
1628{
1629     matrix m[size(l)][1];
1630     for (int i=1; i<=size(l); i++)
1631     {
1632          m[i,1]=l[i];
1633     }
1634     return(m);
1635}
1636
1637///////////////////////////////////////////////////////////////////////////////
1638// convert all the point in the list with list2vec
1639proc convPoints (list points)
1640{
1641     for (int i=1; i<=size(points); i++)
1642     {
1643          points[i]=list2vec(points[i]);
1644     }
1645     return(points);
1646}
1647
1648///////////////////////////////////////////////////////////////////////////////
1649// extracts elements from l in the range m..n
1650proc graspList (list l, int m, int n)
1651{
1652     list result;
1653     int count=1;
1654     for (int i=m; i<=n; i++)
1655     {
1656          result[count]=l[i];
1657          count++;
1658     }
1659     return(result);
1660}
1661
1662///////////////////////////////////////////////////////////////////////////////
1663// "characteristic" polynomial
1664static proc xi_gen (matrix p, int e, int s)
1665{
1666     poly prod=1;
1667     list rl=ringlist(basering);
1668     int charac=rl[1][1];
1669     int l;
1670     for (l=1; l<=s; l++)
1671     {
1672          prod=prod*(1-(var(l)-p[l,1])^(charac^e-1));
1673     }
1674     return(prod);
1675}
1676
1677///////////////////////////////////////////////////////////////////////////////
1678// generating polynomials in Fitzgerald-Lax construction
1679static proc gener_funcs (matrix check, list points, int e, ideal id, int s)
1680{
1681     int n=ncols(check);
1682     if (n!=size(points)) {ERROR("Incompatible sizes of check and points");}
1683     ideal xi;
1684     int i,j;
1685     for (i=1; i<=n; i++)
1686     {
1687          xi[i]=xi_gen(points[i],e,s);
1688     }
1689     ideal result;
1690     int m=nrows(check);
1691     poly sum;
1692     for (i=1; i<=m; i++)
1693     {
1694          sum=0;
1695          for (j=1; j<=n; j++)
1696          {
1697               sum=sum+check[i,j]*xi[j];
1698          }
1699          result[i]=NF(sum,id);
1700     }
1701     return(result);
1702}
1703
1704///////////////////////////////////////////////////////////////////////////////
1705
1706proc sysFL (matrix check, matrix y, int t, int e, int s)
1707"USAGE:    sysFL (check,y,t,e,s); check,y matrix, t,e,s int
1708@format
1709          - check is a parity check matrix of the code,
1710          - y is a received word,
1711          - t the number of errors to correct,
1712          - e is the extension degree,
1713          - s is the dimension of the point for the vanishing ideal
1714@end format
1715RETURN:  the system of Fitzgerald-Lax for the given decoding problem
1716THEORY:  Based on 'check' of the given linear code, the procedure constructs
1717         the corresponding ideal constructed with a generalization of
1718         Cooper's philosophy. For basics of the method @ref{Fitzgerald-Lax method}.
1720EXAMPLE:   example sysFL; shows an example
1721"
1722{
1723     list rl=ringlist(basering);
1724     int charac=rl[1][1];
1725     int n=ncols(check);
1726     int m=nrows(check);
1727     list points=pointsGen(s,e);
1728     list points2=convPoints(points);
1729     list p=graspList(points2,1,n);
1730     ideal id=vanishId(p,e);
1731     ideal funcs=gener_funcs(check,p,e,id,s);
1732
1733     ideal result;
1734     poly temp;
1735     int i,j,k;
1736
1737     //--------------- add vanishing realtions ---------------------
1738     for (i=1; i<=t; i++)
1739     {
1740          for (j=1; j<=size(id); j++)
1741          {
1742               temp=id[j];
1743               for (k=1; k<=s; k++)
1744               {
1745                    temp=subst(temp,var(k),x_var(i,k,s));
1746               }
1747               result=result,temp;
1748          }
1749     }
1750
1751     //--------------- add field equations --------------------
1752     for (i=1; i<=t; i++)
1753     {
1754          for (k=1; k<=s; k++)
1755          {
1756               result=result,x_var(i,k,s)^(charac^e)-x_var(i,k,s);
1757          }
1758     }
1759     for (i=1; i<=t; i++)
1760     {
1761          result=result,e(i)^(charac^e-1)-1;
1762     }
1763
1764     result=simplify(result,8);
1765
1766     //--------------- add check realtions --------------------
1767     poly sum;
1768     matrix syn[m][1]=syndrome(check,y);
1769     for (i=1; i<=size(funcs); i++)
1770     {
1771          sum=0;
1772          for (j=1; j<=t; j++)
1773          {
1774               temp=funcs[i];
1775               for (k=1; k<=s; k++)
1776               {
1777                    temp=subst(temp,var(k),x_var(j,k,s));
1778               }
1779               sum=sum+temp*e(j);
1780          }
1781          result=result,sum-syn[i,1];
1782     }
1783
1784     result=simplify(result,2);
1785
1786     points=points2;
1787     export points;
1788     return(result);
1789}
1790example
1791{
1792     "EXAMPLE:";  echo = 2;
1793     intvec vopt = option(get);
1794
1795     list l=FLpreprocess(3,1,11,2,"");
1796     def r=l[1];
1797     setring r;
1798     int s_work=l[2];
1799
1800     //the check matrix of [11,6,5] ternary code
1801     matrix h[5][11]=1,0,0,0,0,1,1,1,-1,-1,0,
1802          0,1,0,0,0,1,1,-1,1,0,-1,
1803          0,0,1,0,0,1,-1,1,0,1,-1,
1804          0,0,0,1,0,1,-1,0,1,-1,1,
1805          0,0,0,0,1,1,0,-1,-1,1,1;
1806     matrix g=dual_code(h);
1807     matrix x[1][6];
1808     matrix y[1][11]=encode(x,g);
1809     //disturb with 2 errors
1810     matrix rec[1][11]=errorInsert(y,list(2,4),list(1,-1));
1811
1812     //the Fitzgerald-Lax system
1813     ideal sys=sysFL(h,rec,2,1,s_work);
1814     print(sys);
1815     option(redSB);
1816     ideal red_sys=std(sys);
1817     red_sys;
1818     // read the solutions from this redGB
1819     // the points are (0,0,1) and (0,1,0) with error values 1 and -1 resp.
1820     // use list points to find error positions;
1821     points;
1822
1823     option(set,vopt);
1824}
1825
1826///////////////////////////////////////////////////////////////////////////////
1827// preprocessing steps for the Fitzgerald-Lax scheme
1828proc FLpreprocess (int p, int e, int n, int t, string minp)
1829{
1830     ring r1=p,x,dp;
1831     int s=1;
1832     while(p^(s*e)<n)
1833     {
1834          s++;
1835     }
1836     list var_ord;
1837     int i,j;
1838     int count=1;
1839     for (i=s; i>=1; i--)
1840     {
1841          var_ord[count]=string("x("+string(i)+")");
1842          count++;
1843     }
1844     for (i=t; i>=1; i--)
1845     {
1846          var_ord[count]=string("e("+string(i)+")");
1847          count++;
1848          for (j=s; j>=1; j--)
1849          {
1850               var_ord[count]=string("x1("+string(s*(i-1)+j)+")");
1851               count++;
1852          }
1853     }
1854
1855     list rl;
1856     list tmp;
1857
1858     if (e>1)
1859     {
1860          rl[1]=tmp;
1861          rl[1][1]=p;
1862          rl[1][2]=tmp;
1863          rl[1][2][1]=string("a");
1864          rl[1][3]=tmp;
1865          rl[1][3][1]=tmp;
1866          rl[1][3][1][1]=string("lp");
1867          rl[1][3][1][2]=1;
1868          rl[1][4]=ideal(0);
1869     } else {
1870          rl[1]=p;
1871     }
1872
1873     rl[2]=var_ord;
1874
1875     rl[3]=tmp;
1876     rl[3][1]=tmp;
1877     rl[3][1][1]=string("lp");
1878     intvec v=1;
1879     for (i=1; i<=size(var_ord)-1; i++)
1880     {
1881          v=v,1;
1882     }
1883     rl[3][1][2]=v;
1884     rl[3][2]=tmp;
1885     rl[3][2][1]=string("C");
1886     rl[3][2][2]=intvec(0);
1887
1888     rl[4]=ideal(0);
1889
1890     def r2=ring(rl);
1891     setring r2;
1892     list l=ringlist(r2);
1893     if (e>1)
1894     {
1895          execute(string("poly f="+minp));
1896          ideal id=f;
1897          l[1][4]=id;
1898     }
1899
1900     def r=ring(l);
1901     setring r;
1902
1903     return(list(r,s));
1904}
1905
1906///////////////////////////////////////////////////////////////////////////////
1907// imitating two indeces
1908static proc x_var (int i, int j, int s)
1909{
1910     return(x1(s*(i-1)+j));
1911}
1912
1913///////////////////////////////////////////////////////////////////////////////
1914// random vector of length n with entries from p^e, p the characteristic
1915static proc randomvector(int n, int e)
1916{
1917    int i;
1918    matrix result[n][1];
1919    for (i=1; i<=n; i++)
1920    {
1921        result[i,1]=asElement(random_prime_vector(e));
1922    }
1923    return(result);
1924}
1925
1926///////////////////////////////////////////////////////////////////////////////
1927// "convert" representation of an element from the field extension from vector
1928//to an elelemnt
1929static proc asElement(list l)
1930{
1931  number s;
1932  int i;
1933  number w=1;
1934  if (size(l)>1) {w=par(1);}
1935  for (i=0; i<=size(l)-1; i++)
1936  {
1937    s=s+w^i*l[i+1];
1938  }
1939  return(s);
1940}
1941
1942///////////////////////////////////////////////////////////////////////////////
1943// random vector of length n with entries from p, p the characteristic
1944static proc random_prime_vector (int n)
1945{
1946  list rl=ringlist(basering);
1947  int i, charac;
1948  for (i=2; i<=rl[1][1]; i++)
1949  {
1950    if (rl[1][1] mod i ==0)
1951    {
1952      break;
1953    }
1954  }
1955  charac=i;
1956
1957  list l;
1958
1959  for (i=1; i<=n; i++)
1960  {
1961    l=l+list(random(0,charac-1));
1962  }
1963  return(l);
1964}
1965
1966///////////////////////////////////////////////////////////////////////////////
1967
1968proc decodeRandomFL(int n, int redun, int p, int e, int t, int ncodes, int ntrials, string minpol)
1969"USAGE:    decodeRandomFL(redun,p,e,n,t,ncodes,ntrials,minpol);
1970@format
1971          - n is length of codes generated,
1972          - redun = redundancy of codes generated,
1973          - p is the characteristic,
1974          - e is the extension degree,
1975          - t is the number of errors to correct,
1976          - ncodes is the number of random codes to be processed,
1977          - ntrials is the number of received vectors per code to be corrected,
1978          - minpol: due to some pecularities of SINGULAR one needs to provide
1979          minimal polynomial for the extension explicitly
1980@end format
1981RETURN:    nothing
1982EXAMPLE:   example decodeRandomFL; shows an example
1983"
1984{
1985 intvec vopt = option(get);
1986
1987 list l=FLpreprocess(p,e,n,t,minpol);
1988
1989 def r=l[1];
1990 int s_work=l[2];
1991 export(s_work);
1992 setring r;
1993
1994 int i,j;
1995 matrix h, g, word, y, rec;
1996 ideal sys, sys2, sys3;
1997
1998
1999 option(redSB);
2000 matrix z[1][n];
2001
2002 for (i=1; i<=ncodes; i++)
2003 {
2004     h=randomCheck(redun,n,e);
2005     g=dual_code(h);
2006
2007     //---------------- generate the template system -----------------------
2008     sys=sysFL(h,z,t,e,s_work);
2009
2010     //------ modifying the template according to the received word ---------
2011     for (j=1; j<=ntrials; j++)
2012     {
2013          word=randomvector(n-redun,1);
2014          y=encode(transpose(word),g);
2015          print("Codeword:");
2016          print(y);
2017          rec=errorRand(y,t,e);
2019          print(rec);
2021          sys3=std(sys2);
2022          print("Groebner basis of the FL system:");
2023          print(sys3);
2024     }
2025 }
2026
2027 option(set,vopt);
2028}
2029example
2030{
2031     "EXAMPLE:";  echo = 2;
2032
2033     // correcting one error for one random binary code of length 25,
2034     // redundancy 14; 10 words are processed
2035     decodeRandomFL(25,14,2,1,1,1,10,"");
2036}
2037
2038///////////////////////////////////////////////////////////////////////////////
2039// add syndrome values to the template system in FL
2040static proc LF_add_synd (matrix rec, matrix check, ideal sys)
2041{
2042     int redun=nrows(check);
2043     ideal result=sys;
2044     matrix s[redun][1]=syndrome(check,rec);
2045     for (int i=size(sys)-redun+1; i<=size(sys); i++)
2046     {
2047          result[i]=result[i]-s[i-size(sys)+redun,1];
2048     }
2049     return(result);
2050}
2051
2052
2053/*
2054//////////////     SOME RELATIVELY EASY EXAMPLES    //////////////
2055///////////////////  THAT RUN AROUND ONE MINUTE  ////////////////
2056
2057"EXAMPLE:";  echo = 2;
2058int q=128; int n=120; int redun=n-30;
2059ring r=(q,a),x,dp;
2060decodeRandom(n,redun,1,1,6);
2061
2062int q=128; int n=120; int redun=n-20;
2063ring r=(q,a),x,dp;
2064decodeRandom(n,redun,1,1,9);
2065
2066int q=128; int n=120; int redun=n-10;
2067ring r=(q,a),x,dp;
2068decodeRandom(n,redun,1,1,19);
2069
2070int q=256; int n=150; int redun=n-10;
2071ring r=(q,a),x,dp;
2072decodeRandom(n,redun,1,1,22);
2073
2074//////////////     SOME HARD EXAMPLES    //////////////////////
2075//////      THAT MAYBE WILL BE DOABLE LATER     ///////////////
2076
20771.) These random instances are not doable in <=1000 sec.
2078
2079"EXAMPLE:";  echo = 2;
2080int q=128; int n=120; int redun=n-40;
2081ring r=(q,a),x,dp;
2082decodeRandom(n,redun,1,1,6);
2083
2084redun=n-30;
2085decodeRandom(n,redun,1,1,8);
2086
2087redun=n-20;
2088decodeRandom(n,redun,1,1,12);
2089
2090redun=n-10;
2091decodeRandom(n,redun,1,1,24);
2092
2093int q=256; int n=150; int redun=n-10;
2094ring r=(q,a),x,dp;
2095decodeRandom(n,redun,1,1,26);
2096
2097
20982.) Generic decoding is hard!
2099
2100int q=32; int n=31; int redun=n-16; int t=3;
2101ring r=(q,a),(V(1..n),U(n..1),s(redun..1)),(dp(n),lp(n),dp(redun));
2102matrix check[redun][n]= 1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,
21030,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,
21040,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
21050,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,
21060,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,
21070,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
21080,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,
21090,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,
21100,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,
21110,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,
21120,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,
21131,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,
21140,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,
21151,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,
21161,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,
21170,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,
21180,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
21190,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,
21200,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,
21210,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,
21220,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,
21230,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,
21240,1,0,1,0,0,1,0,0,1;
2125matrix rec[1][n];
2126
2127def A=sysQE(check,rec,t,1,2);
2128setring A;
2129print(qe);
2130ideal red_qe=stdfglm(qe);
2131
2132*/
Note: See TracBrowser for help on using the repository browser.