source: git/Singular/LIB/decodegb.lib @ 5a9fa5

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