# source:git/Singular/LIB/decodegb.lib@76afcd

spielwiese
Last change on this file since 76afcd was 76afcd, checked in by Stanislav Bulygin <bulygin@…>, 14 years ago
• Property mode set to 100644
File size: 52.1 KB
Line
1///////////////////////////////////////////////////////////////////////////////
2version="\$Id: decodegb.lib,v 1.10 2009-03-25 11:15:16 bulygin Exp \$";
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 original method of quadratic equations
12 is worked out here as well. We also (for comparison) enable to work with the
13 system of Fitzgerald-Lax. We provide some auxiliary functions for further
14 manipulations and decoding. For an overview of the methods mentioned above @ref{Decoding codes with GB}.
15 For the vanishing ideal computation the algorithm of Farr and Gao is
16 implemented.
17
18MAIN PROCEDURES:
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
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 int tim=rtimer;
726 matrix transf=inverse(transpose(h_full));
727
728 //------ expression matrix of check vectors w.r.t. the MDS basis -----------
729 tim=rtimer;
730 for (i=1; i<=red ; i++)
731 {
732  a[i]=transpose(submat(check,i..i,1..n));
733  a[i]=transf*a[i];
734 }
735
736 //----------- compute the structure constants ------------------------
737 tim=rtimer;
738 matrix te[n][1];
739 for (i=1; i<=n; i++)
740 {
741  for (j=1; j<=t+1; j++)
742  {
743   if ((j<i)&&(i<=t+1)) {c[i][j]=c[j][i];}
744   else
745   {
746    if (i+j<=n+1)
747    {
748     c[i][j]=te;
749     c[i][j][i+j-1,1]=1;
750    }
751    else
752    {
753     c[i][j]=star(h_full,i,j);
754     c[i][j]=transf*c[i][j];
755    }
756   }
757  }
758 }
759
760
761 tim=rtimer;
762 if (formal==0)
763 {
764  matrix s[red][1]=syndrome(check,y);
765  for (j=1; j<=red; j++)
766  {
767   sum1=0;
768   for (l=1; l<=n; l++)
769   {
770    sum1=sum1+a[j][l,1]*U(l);
771   }
772   result=result,sum1-s[j,1];
773  }
774 } else
775 {
776  for (j=1; j<=red; j++)
777  {
778   sum1=0;
779   for (l=1; l<=n; l++)
780   {
781    sum1=sum1+a[j][l,1]*U(l);
782   }
783   result=result,sum1-s(j);
784  }
785  for (j=1; j<=red; j++)
786  {
787     result=result,s(j)^(formal)-s(j);
788  }
789 }
790 if (fieldeq)
791 {
792  for (i=1; i<=n; i++)
793  {
794   result=result,U(i)^q-U(i);
795  }
796  for (j=1; j<=t; j++)
797  {
798     result=result,V(j)^q-V(j);
799  }
800 }
801
802 //----- form the quadratic equations according to the theory -----------
803 for (i=1; i<=n; i++)
804 {
805  sum1=0;
806  for (j=1; j<=t; j++)
807  {
808   sum2=0;
809   for (l=1; l<=n; l++)
810   {
811    sum2=sum2+c[i][j][l,1]*U(l);
812   }
813   sum1=sum1+sum2*V(j);
814  }
815  sum3=0;
816  for (l=1; l<=n; l++)
817  {
818   sum3=sum3+c[i][t+1][l,1]*U(l);
819  }
820  result=result,sum1-sum3;
821 }
822
823 result=simplify(result,2);
824
825 ideal qe=result;
826 export qe;
827 return(work);
828}
829example
830{
831     "EXAMPLE:";  echo = 2;
832     intvec v = option(get);
833
834     //correct 2 errors in [7,3] 8-ary code RS code
835     int t=2; int q=8; int n=7; int redun=4;
836     ring r=(q,a),x,dp;
837     matrix h_full=genMDSMat(n,a);
838     matrix h=submat(h_full,1..redun,1..n);
839     matrix g=dual_code(h);
840     matrix x[1][3]=0,0,1,0;
841     matrix y[1][7]=encode(x,g);
842
843     //disturb with 2 errors
844     matrix rec[1][7]=errorInsert(y,list(2,4),list(1,a));
845
846     //generate the system
847     def A=sysQE(h,rec,t);
848     setring A;
849     print(qe);
850
851     //let us decode
852     option(redSB);
853     ideal sys_qe=std(qe);
854     print(sys_qe);
855
856     option(set,v);
857}
858
859///////////////////////////////////////////////////////////////////////////////
860
861proc errorInsert(matrix y, list pos, list val)
862"USAGE:  errorInsert(y,pos,val); y is matrix, pos,val are list of int's
863@format
864        - y is a (code) word,
865        - pos = positions where errors occured,
866        - val = their corresponding values
867@end format
869EXAMPLE: example errorInsert; shows an example
870"
871{
872 matrix result[1][ncols(y)]=y;
873 if (size(pos)!=size(val)) {print("ERRORerror");}
874 for (int i=1; i<=size(pos); i++)
875 {
876  result[1,pos[i]]=y[1,pos[i]]+val[i];
877 }
878 return(result);
879}
880example
881{
882     "EXAMPLE:";  echo = 2;
883     //correct 2 errors in [7,3] 8-ary code RS code
884     int t=2; int q=8; int n=7; int redun=4;
885     ring r=(q,a),x,dp;
886     matrix h_full=genMDSMat(n,a);
887     matrix h=submat(h_full,1..redun,1..n);
888     matrix g=dual_code(h);
889     matrix x[1][3]=0,0,1,0;
890     matrix y[1][7]=encode(x,g);
891     print(y);
892
893     //disturb with 2 errors
894     matrix rec[1][7]=errorInsert(y,list(2,4),list(1,a));
895     print(rec);
896     print(rec-y);
897}
898
899///////////////////////////////////////////////////////////////////////////////
900
901proc errorRand(matrix y, int num, int e)
902"USAGE:    errorRand(y, num, e); y is matrix, num,e are int
903@format
904          - y is a (code) word,
905          - num is the number of errors,
906          - e is an extension degree (if one wants values to be from GF(p^e))
907@end format
909EXAMPLE:   example errorRand; shows an example
910"
911{
912 matrix result[1][ncols(y)]=y;
913 int i,j, flag, temp;
914 list pos, val;
915 matrix tempnum;
916
917 for (i=1; i<=num; i++)
918 {
919  while(1)
920  {
921   temp=random(1,ncols(y));
922   flag=1;
923   for (j=1; j<=size(pos); j++)
924   {
925    if (temp==pos[j]) {flag=0;}
926   }
927   if (flag) {pos[i]=temp;break;}
928  }
929 }
930
931 for (i=1; i<=num; i++)
932 {
933  flag=1;
934  while(flag)
935  {
936   tempnum=randomvector(1,e);
937   if (tempnum!=0) {flag=0;}
938  }
939  val[i]=tempnum;
940 }
941
942 for (i=1; i<=size(pos); i++)
943 {
944  result[1,pos[i]]=y[1,pos[i]]+val[i];
945 }
946 return(result);
947}
948example
949{
950  "EXAMPLE:";  echo = 2;
951     //correct 2 errors in [7,3] 8-ary code RS code
952     int t=2; int q=8; int n=7; int redun=4;
953     ring r=(q,a),x,dp;
954     matrix h_full=genMDSMat(n,a);
955     matrix h=submat(h_full,1..redun,1..n);
956     matrix g=dual_code(h);
957     matrix x[1][3]=0,0,1,0;
958     matrix y[1][7]=encode(x,g);
959
960     //disturb with 2 random errors
961     matrix rec[1][7]=errorRand(y,2,3);
962     print(rec);
963     print(rec-y);
964}
965
966///////////////////////////////////////////////////////////////////////////////
967
968proc randomCheck(int m, int n, int e)
969"USAGE:    randomCheck(m, n, e); m,n,e are int
970@format
971          - m x n are dimensions of the matrix,
972          - e is an extension degree (if one wants values to be from GF(p^e))
973@end format
974RETURN:    random check matrix
975EXAMPLE:   example randomCheck; shows an example
976"
977{
978 matrix result[m][n];
979 matrix rand[m][n-m];
980 int i,j;
981 matrix temp;
982 for (i=1; i<=m; i++)
983 {
984  temp=randomvector(n-m,e);
985  for (j=1; j<=n-m; j++)
986  {
987   rand[i,j]=temp[j,1];
988  }
989 }
990 result=concat(rand,unitmat(m));
991 return(result);
992}
993example
994{
995  "EXAMPLE:";  echo = 2;
996     int redun=5; int n=15;
997     ring r=2,x,dp;
998
999     //generate random check matrix for a [15,5] binary code
1000     matrix h=randomCheck(redun,n,1);
1001     print(h);
1002
1003     //corresponding generator matrix
1004     matrix g=dual_code(h);
1005     print(g);
1006}
1007
1008///////////////////////////////////////////////////////////////////////////////
1009
1010proc genMDSMat(int n, number a)
1011"USAGE:   genMDSMat(n, a); n is int, a is number
1012@format
1013        - n x n are dimensions of the MDS matrix,
1014        - a is a primitive element of the field.
1015@end format
1016NOTE:   An MDS matrix is constructed in the following way. We take 'a' to be a
1017        generator of the multiplicative group of the field. Then we construct
1018        the Vandermonde matrix with this 'a'.
1019ASSUME:   extension field should already be defined
1020RETURN:   a matrix with the MDS property. @xref{Decoding method based on quadratic equations}
1022EXAMPLE:  example genMDSMat; shows an example
1023"
1024{
1025 int i,j;
1026 matrix result[n][n];
1027 for (i=0; i<=n-1; i++)
1028 {
1029  for (j=0; j<=n-1; j++)
1030  {
1031   result[j+1,i+1]=(a^i)^j;
1032  }
1033 }
1034 return(result);
1035}
1036example
1037{
1038     "EXAMPLE:";  echo = 2;
1039     int q=16; int n=15;
1040     ring r=(q,a),x,dp;
1041
1042     //generate an MDS (Vandermonde) matrix
1043     matrix h_full=genMDSMat(n,a);
1044     print(h_full);
1045}
1046
1047///////////////////////////////////////////////////////////////////////////////
1048
1049
1050proc mindist (matrix check)
1051"USAGE:  mindist (check, q); check matrix, q int
1052@format
1053        - check is a check matrix,
1054        - q is the field size
1055@end format
1056RETURN:  minimum distance of the code together with the time needed for its
1057         calculation
1058EXAMPLE: example mindist; shows an example
1059"
1060{
1061 intvec vopt = option(get);
1062
1063 int n=ncols(check); int redun=nrows(check); int t=redun+1;
1064
1065 def br=basering;
1066 list rl=ringlist(br);
1067 int q=rl[1][1];
1068
1069 ring work=(q,a),(V(1..t),U(1..n)),dp;
1070 matrix check=imap(br,check);
1071
1072 ideal temp;
1073 int count=1;
1074 int flag=1;
1075 int flag2;
1076 int i, tim, timsolve;
1077 matrix z[1][n];
1078 option(redSB);
1079 def A=sysQE(check,z,count);
1080
1081 //proceed with solving the system w.r.t zero vector until some solutions
1082 //are found
1083 while (flag)
1084 {
1085    A=sysQE(check,z,count);
1086    setring A;
1087    ideal temp=qe;
1088    tim=rtimer;
1089    temp=std(temp);
1090    timsolve=timsolve+rtimer-tim;
1091    flag2=1;
1092    setring work;
1093    temp=imap(A,temp);
1094    for (i=1; i<=n; i++)
1095    {
1096      if
1097        (temp[i]!=U(n-i+1))
1098        {
1099          flag2=0;
1100        }
1101    }
1102    if (!flag2)
1103    {
1104      flag=0;
1105    }
1106    else
1107    {
1108      count++;
1109    }
1110 }
1111 list result=list(count,timsolve);
1112
1113 option(set,vopt);
1114 return(result);
1115}
1116example
1117{
1118     "EXAMPLE:";  echo = 2;
1119     //determine a minimum distance for a [7,3] binary code
1120     int q=8; int n=7; int redun=4; int t=redun+1;
1121     ring r=(q,a),x,dp;
1122
1123     //generate random check matrix
1124     matrix h=randomCheck(redun,n,1);
1125     print(h);
1126     list l=mindist(h);
1127     print(l[1]);
1128     //time for the comutation in secs
1129     print(l[2]);
1130}
1131
1132///////////////////////////////////////////////////////////////////////////////
1133
1134proc decode(matrix check, matrix rec)
1135"USAGE:    decode(check, rec, t); check, rec matrix, t int
1136@format
1137          - check is the check matrix of the code,
1138          - rec is a received word,
1139          - t is an upper bound for the number of errors one wants to correct
1140@end format
1141NOTE:     The method described in @ref{Decoding method based on quadratic equations}
1142          is used for decoding.
1143ASSUME:   Errors in rec should be correctable, otherwise the output is
1144          unpredictable
1145RETURN:   a codeword that is closest to rec
1146EXAMPLE:  example decode; shows an example
1147"
1148{
1149 intvec vopt = option(get);
1150
1151 def br=basering;
1152 int n=ncols(check);
1153
1154 int count=1;
1155 def A=sysQE(check,rec,count);
1156 while(1)
1157 {
1158  A=sysQE(check,rec,count);
1159  setring A;
1160  matrix h_full=genMDSMat(n,a);
1161  matrix rec=imap(br,rec);
1162  option(redSB);
1163  ideal qe_red=std(qe);
1164  if (qe_red[1]!=1)
1165  {
1166    break;
1167  }
1168  else
1169  {
1170    count++;
1171  }
1172  setring br;
1173 }
1174
1175 setring A;
1176
1177 //obtain a codeword
1178 //this works only if our code is indeed can correct these errors
1179 matrix syn[n][1];
1180 for (int i=1; i<=n; i++)
1181 {
1183 }
1184
1185 matrix real_syn=inverse(h_full)*syn;
1186 setring br;
1187 matrix real_syn=imap(A,real_syn);
1188
1189 option(set,vopt);
1190 return(rec-transpose(real_syn));
1191}
1192example
1193{
1194     "EXAMPLE:";  echo = 2;
1195     //correct 1 error in [15,7] binary code
1196     int t=1; int q=16; int n=15; int redun=10;
1197     ring r=(q,a),x,dp;
1198
1199     //generate random check matrix
1200     matrix h=randomCheck(redun,n,1);
1201     matrix g=dual_code(h);
1202     matrix x[1][n-redun]=0,0,1,0,1,0,1;
1203     matrix y[1][n]=encode(x,g);
1204     print(y);
1205
1206     // find out the minimum distance of the code
1207     list l=mindist(h);
1208
1209     //disturb with errors
1210     "Correct ",(l[1]-1) div 2," errors";
1211     matrix rec[1][n]=errorRand(y,(l[1]-1) div 2,1);
1212     print(rec);
1213
1214     //let us decode
1215     matrix dec_word=decode(h,rec);
1216     print(dec_word);
1217}
1218
1219///////////////////////////////////////////////////////////////////////////////
1220
1221
1222proc decodeRandom(int n, int redun, int ncodes, int ntrials, list #)
1223"USAGE:    decodeRandom(redun,q,ncodes,ntrials,[e]); all parameters int
1224@format
1225          - redun is a redundabcy of a (random) code,
1226          - q is the field size,
1227          - ncodes is the number of random codes to be processed,
1228          - ntrials is the number of received vectors per code to be corrected
1229          - If e is given it sets the correction capacity explicitly. It
1230          should be used in case one expects some lower bound,
1231          otherwise the procedure tries to compute the real minimum distance
1232          to find out the error-correction capacity
1233@end format
1234RETURN:    nothing;
1235EXAMPLE:   example decodeRandom; shows an example
1236"
1237{
1238 intvec vopt = option(get);
1239
1240 int i,j;
1241 matrix h;
1242 int dist, t, tim, tim2, tim3, timdist, timdec, timdist2, timdec2, timdec3;
1243 ideal sys;
1244 list tmp;
1245 int e;
1246 if (size(#)>0)
1247 {
1248  e=#[1];
1249 }
1250
1251 option(redSB);
1252 def br=basering;
1253 matrix h_full=genMDSMat(n,a);
1254 matrix z[1][ncols(h_full)];
1255
1256 //------------------ determine error-correction capacity -------------------
1257 for (i=1; i<=ncodes; i++)
1258 {
1259  setring br;
1260  h=randomCheck(redun,n,1);
1261  "check matrix:";
1262  print(h);
1263  if (e>0)
1264  {
1265     t=e;
1266  } else {
1267     tim=rtimer;
1268     tmp=mindist(h);
1269     timdist=timdist+rtimer-tim;
1270     timdist2=timdist2+tmp[2];
1271     dist=tmp[1];
1272     printf("d= %p",dist);
1273     t=(dist-1) div 2;
1274  }
1275  tim2=rtimer;
1276  tim3=rtimer;
1277
1278  //------------- generate the template system ----------------------
1279  def A=sysQE(h,z,t);
1280  setring A;
1281  matrix word,y,rec;
1282  ideal sys2,sys3;
1283  matrix h=imap(br,h);
1284  matrix g=dual_code(h);
1285  ideal sys=qe;
1286  print("The system is generated");
1287  timdec3=timdec3+rtimer-tim3;
1288
1289  //------ modify the template according to every received word --------------
1290  for (j=1; j<=ntrials; j++)
1291  {
1292   word=randomvector(n-redun,1);
1293   y=encode(transpose(word),g);
1294   rec=errorRand(y,t,1);
1296
1297   tim=rtimer;
1298   sys3=std(sys2);
1299   timdec=timdec+rtimer-tim;
1300  }
1301  timdec2=timdec2+rtimer-tim2;
1302  kill A;
1303  option(set,vopt);
1304 }
1305 printf("Time for mindist: %p", timdist);
1306 printf("Time for GB in mindist: %p", timdist);
1307 printf("Time for decoding: %p", timdec2);
1308 printf("Time for GB in decoding: %p", timdec);
1309 printf("Time for sysQE in decoding: %p", timdec3);
1310}
1311example
1312{
1313     "EXAMPLE:";  echo = 2;
1314     int q=32; int n=25; int redun=n-11; int t=redun+1;
1315     ring r=(q,a),x,dp;
1316
1317     // correct 2 errors in 5 random binary codes, 50 trials each
1318     decodeRandom(n,redun,5,50,2);
1319}
1320
1321///////////////////////////////////////////////////////////////////////////////
1322
1323
1324proc decodeCode(matrix check, int ntrials, list #)
1325"USAGE:     decodeCode(check, ntrials, [e]); check matrix, ntrials,e int
1326@format
1327           - check is a parity check matrix for the code,
1328           - ntrials is the number of received vectors per code to be
1329           corrected.
1330           - If e is given it sets the correction capacity explicitly. It
1331           should be used in case one expects some lower bound,
1332           otherwise the procedure tries to compute the real minimum distance
1333           to find out the error-correction capacity
1334@end format
1335RETURN:     nothing;
1336EXAMPLE:    example decodeCode; shows an example
1337"
1338{
1339 intvec vopt = option(get);
1340
1341 int n=ncols(check);
1342 int redun=nrows(check);
1343 int i,j;
1344 matrix h;
1345 int dist, t, tim, tim2, tim3, timdist, timdec, timdist2, timdec2, timdec3;
1346 ideal sys;
1347 list tmp;
1348 int e;
1349 if (size(#)>0)
1350 {
1351  e=#[1];
1352 }
1353
1354 option(redSB);
1355 def br=basering;
1356 matrix h_full=genMDSMat(n,a);
1357 matrix z[1][ncols(h_full)];
1358 setring br;
1359 h=check;
1360 "check matrix:";
1361 print(h);
1362
1363 //------------------ determine error-correction capacity -------------------
1364 if (e>0)
1365 {
1366    t=e;
1367 } else {
1368   tim=rtimer;
1369   tmp=mindist(h);
1370   timdist=timdist+rtimer-tim;
1371   timdist2=timdist2+tmp[2];
1372   dist=tmp[1];
1373   printf("d= %p",dist);
1374   t=(dist-1) div 2;
1375 }
1376 tim2=rtimer;
1377 tim3=rtimer;
1378
1379 //------------- generate the template system ----------------------
1380 def A=sysQE(h,z,t);
1381 setring A;
1382 matrix word,y,rec;
1383 ideal sys2,sys3;
1384 matrix h=imap(br,h);
1385 matrix g=dual_code(h);
1386 ideal sys=qe;
1387 print("The system is generated");
1388 timdec3=timdec3+rtimer-tim3;
1389
1390 //--- modify the template according to every received word ---------------
1391 for (j=1; j<=ntrials; j++)
1392 {
1393   word=randomvector(n-redun,1);
1394   y=encode(transpose(word),g);
1395   rec=errorRand(y,t,1);
1397
1398   tim=rtimer;
1399   sys3=std(sys2);
1400   timdec=timdec+rtimer-tim;
1401 }
1402 timdec2=timdec2+rtimer-tim2;
1403
1404 printf("Time for mindist: %p", timdist);
1405 printf("Time for GB in mindist: %p", timdist);
1406 printf("Time for decoding: %p", timdec2);
1407 printf("Time for GB in decoding: %p", timdec);
1408 printf("Time for sysQE in decoding: %p", timdec3);
1409
1410 option(set,vopt);
1411}
1412example
1413{
1414     "EXAMPLE:";  echo = 2;
1415     int q=32; int n=25; int redun=n-11; int t=redun+1;
1416     ring r=(q,a),x,dp;
1417     matrix check=randomCheck(redun,n,1);
1418
1419     // correct 2 errors in using the code above, 50 trials
1420     decodeCode(check,50,2);
1421}
1422
1423
1424///////////////////////////////////////////////////////////////////////////////
1425// adding syndrome values to the template system
1426static proc add_synd (matrix rec, matrix check, int redun, ideal sys)
1427{
1428     ideal result=sys;
1429     matrix s[redun][1]=syndrome(check,rec);
1430     for (int i=1; i<=redun; i++)
1431
1432     {
1433          result[i]=result[i]-s[i,1];
1434     }
1435     return(result);
1436}
1437
1438///////////////////////////////////////////////////////////////////////////////
1439// evaluate a polynomial at a given point
1440static proc ev (poly f, matrix p)
1441{
1442     if (ncols(p)>1) {ERROR("not a column vector");};
1443     int m=size(p);
1444     poly temp=f;
1445     for (int i=1; i<=m; i++)
1446     {
1447          temp=subst(temp,x(i),p[i,1]);
1448     }
1449     return(number(temp));
1450}
1451
1452///////////////////////////////////////////////////////////////////////////////
1453// return index of an element in the ideal where it does not vanish at the
1454//given point
1455static proc find_index (ideal G, matrix p)
1456{
1457     if (ncols(p)>1) {ERROR("not a column vector");};
1458     int i=1;
1459     int n=size(G);
1460     while(i<=n)
1461     {
1462          if (ev(G[i],p)!=0) {return(i);}
1463          i++;
1464     }
1465     return(-1);
1466}
1467
1468///////////////////////////////////////////////////////////////////////////////
1469// convert ideal to list
1470static proc ideal2list (ideal id)
1471{
1472     list l;
1473     for (int i=1; i<=size(id); i++)
1474     {
1475          l[i]=id[i];
1476     }
1477     return(l);
1478}
1479
1480///////////////////////////////////////////////////////////////////////////////
1481// convert list to ideal
1482static proc list2ideal (list l)
1483{
1484     ideal id;
1485     for (int i=1; i<=size(l); i++)
1486     {
1487          id[i]=l[i];
1488     }
1489     return(id);
1490}
1491
1492///////////////////////////////////////////////////////////////////////////////
1493// check whether given polynomial is divisible by some leading monomial of the
1494//ideal
1495static proc divisible (poly m, ideal G)
1496{
1497     for (int i=1; i<=size(G); i++)
1498     {
1500     }
1501     return(0);
1502}
1503
1504///////////////////////////////////////////////////////////////////////////////
1505
1506proc vanishId (list points)
1507"USAGE:  vanishId (points); point is a list of matrices
1508        'points' is a list of points for which the vanishing ideal is to be
1509        constructed
1510RETURN:  Vanishing ideal corresponding to the given set of points
1511EXAMPLE: example vanishId; shows an example
1512"
1513{
1514     int m=size(points[1]);
1515     int n=size(points);
1516
1517     ideal G=1;
1518     int i,k,j;
1519     list temp;
1520     poly h,cur;
1521
1522     //------------- proceed according to Farr-Gao algorithm ----------------
1523     for (k=1; k<=n; k++)
1524     {
1525          i=find_index(G,points[k]);
1526          cur=G[i];
1527          for(j=i+1; j<=size(G); j++)
1528          {
1529               G[j]=G[j]-ev(G[j],points[k])/ev(G[i],points[k])*G[i];
1530          }
1531          G=simplify(G,2);
1532          temp=ideal2list(G);
1533          temp=delete(temp,i);
1534          G=list2ideal(temp);
1535          for (j=1; j<=m; j++)
1536          {
1538               {
1539                    attrib(G,"isSB",1);
1540                    h=NF((x(j)-points[k][j,1])*cur,G);
1541                    temp=ideal2list(G);
1542                    temp=insert(temp,h);
1543                    G=list2ideal(temp);
1544                    G=sort(G)[1];
1545               }
1546          }
1547     }
1548     attrib(G,"isSB",1);
1549     return(G);
1550}
1551example
1552{
1553     "EXAMPLE:";  echo = 2;
1554      ring r=3,(x(1..3)),dp;
1555
1556     //generate all 3-vectors over GF(3)
1557     list points=pointsGen(3,1);
1558
1559     list points2=convPoints(points);
1560
1561     //grasps the first 11 points
1562     list p=graspList(points2,1,11);
1563     print(p);
1564
1565     //construct the vanishing ideal
1566     ideal id=vanishId(p);
1567     print(id);
1568}
1569
1570///////////////////////////////////////////////////////////////////////////////
1571// construct the list of all vectors of length m with elements in p^e, where p
1572//is characteristics
1573proc pointsGen (int m, int e)
1574{
1575     if (e>1)
1576     {
1577     list result;
1578     int count=1;
1579     int i,j;
1580     list l=ringlist(basering);
1581     int charac=l[1][1];
1582     number a=par(1);
1583     list tmp;
1584     for (i=1; i<=charac^(e*m); i++)
1585     {
1586          result[i]=tmp;
1587     }
1588     if (m==1)
1589     {
1590          result[count][m]=0;
1591          count++;
1592          for (j=1; j<=charac^(e)-1; j++)
1593          {
1594               result[count][m]=a^j;
1595               count++;
1596          }
1597          return(result);
1598     }
1599     list prev=pointsGen(m-1,e);
1600     for (i=1; i<=size(prev); i++)
1601     {
1602          result[count]=prev[i];
1603          result[count][m]=0;
1604          count++;
1605          for (j=1; j<=charac^(e)-1; j++)
1606          {
1607               result[count]=prev[i];
1608               result[count][m]=a^j;
1609               count++;
1610          }
1611     }
1612     return(result);
1613     }
1614
1615     if (e==1)
1616     {
1617     list result;
1618     int count=1;
1619     int i,j;
1620     list l=ringlist(basering);
1621     int charac=l[1][1];
1622     list tmp;
1623     for (i=1; i<=charac^m; i++)
1624     {
1625          result[i]=tmp;
1626     }
1627     if (m==1)
1628     {
1629          for (j=0; j<=charac-1; j++)
1630          {
1631               result[count][m]=number(j);
1632               count++;
1633          }
1634          return(result);
1635     }
1636     list prev=pointsGen(m-1,e);
1637     for (i=1; i<=size(prev); i++)
1638     {
1639          for (j=0; j<=charac-1; j++)
1640          {
1641               result[count]=prev[i];
1642               result[count][m]=number(j);
1643               count++;
1644          }
1645     }
1646     return(result);
1647     }
1648
1649}
1650
1651///////////////////////////////////////////////////////////////////////////////
1652// convert list to a column vector
1653static proc list2vec (list l)
1654{
1655     matrix m[size(l)][1];
1656     for (int i=1; i<=size(l); i++)
1657     {
1658          m[i,1]=l[i];
1659     }
1660     return(m);
1661}
1662
1663///////////////////////////////////////////////////////////////////////////////
1664// convert all the point in the list with list2vec
1665proc convPoints (list points)
1666{
1667     for (int i=1; i<=size(points); i++)
1668     {
1669          points[i]=list2vec(points[i]);
1670     }
1671     return(points);
1672}
1673
1674///////////////////////////////////////////////////////////////////////////////
1675// extracts elements from l in the range m..n
1676proc graspList (list l, int m, int n)
1677{
1678     list result;
1679     int count=1;
1680     for (int i=m; i<=n; i++)
1681     {
1682          result[count]=l[i];
1683          count++;
1684     }
1685     return(result);
1686}
1687
1688///////////////////////////////////////////////////////////////////////////////
1689// "characteristic" polynomial
1690static proc xi_gen (matrix p, int e, int s)
1691{
1692     poly prod=1;
1693     list rl=ringlist(basering);
1694     int charac=rl[1][1];
1695     int l;
1696     for (l=1; l<=s; l++)
1697     {
1698          prod=prod*(1-(x(l)-p[l,1])^(charac^e-1));
1699     }
1700     return(prod);
1701}
1702
1703///////////////////////////////////////////////////////////////////////////////
1704// generating polynomials in Fitzgerald-Lax construction
1705static proc gener_funcs (matrix check, list points, int e, ideal id, int s)
1706{
1707     int n=ncols(check);
1708     if (n!=size(points)) {ERROR("Incompatible sizes of check and points");}
1709     ideal xi;
1710     int i,j;
1711     for (i=1; i<=n; i++)
1712     {
1713          xi[i]=xi_gen(points[i],e,s);
1714     }
1715     ideal result;
1716     int m=nrows(check);
1717     poly sum;
1718     for (i=1; i<=m; i++)
1719     {
1720          sum=0;
1721          for (j=1; j<=n; j++)
1722          {
1723               sum=sum+check[i,j]*xi[j];
1724          }
1725          result[i]=NF(sum,id);
1726     }
1727     return(result);
1728}
1729
1730///////////////////////////////////////////////////////////////////////////////
1731
1732proc sysFL (matrix check, matrix y, int t, int e, int s)
1733"USAGE:    sysFL (check,y,t,e,s); check,y matrix, t,e,s int
1734@format
1735          - check is a parity check matrix of the code,
1736          - y is a received word,
1737          - t the number of errors to correct,
1738          - e is the extension degree,
1739          - s is the dimension of the point for the vanishing ideal
1740@end format
1741RETURN:  the system of Fitzgerald-Lax for the given decoding problem
1742THEORY:  Based on 'check' of the given linear code, the procedure constructs
1743         the corresponding ideal constructed with a generalization of
1744         Cooper's philosophy. For basics of the method @ref{Fitzgerald-Lax method}.
1746EXAMPLE:   example sysFL; shows an example
1747"
1748{
1749     list rl=ringlist(basering);
1750     int charac=rl[1][1];
1751     int n=ncols(check);
1752     int m=nrows(check);
1753     list points=pointsGen(s,e);
1754     list points2=convPoints(points);
1755     list p=graspList(points2,1,n);
1756     ideal id=vanishId(p,e);
1757     ideal funcs=gener_funcs(check,p,e,id,s);
1758
1759     ideal result;
1760     poly temp;
1761     int i,j,k;
1762
1763     //--------------- add vanishing realtions ---------------------
1764     for (i=1; i<=t; i++)
1765     {
1766          for (j=1; j<=size(id); j++)
1767          {
1768               temp=id[j];
1769               for (k=1; k<=s; k++)
1770               {
1771                    temp=subst(temp,x(k),x_var(i,k,s));
1772               }
1773               result=result,temp;
1774          }
1775     }
1776
1777     //--------------- add field equations --------------------
1778     for (i=1; i<=t; i++)
1779     {
1780          for (k=1; k<=s; k++)
1781          {
1782               result=result,x_var(i,k,s)^(charac^e)-x_var(i,k,s);
1783          }
1784     }
1785     for (i=1; i<=t; i++)
1786     {
1787          result=result,e(i)^(charac^e-1)-1;
1788     }
1789
1790     result=simplify(result,8);
1791
1792     //--------------- add check realtions --------------------
1793     poly sum;
1794     matrix syn[m][1]=syndrome(check,y);
1795     for (i=1; i<=size(funcs); i++)
1796     {
1797          sum=0;
1798          for (j=1; j<=t; j++)
1799          {
1800               temp=funcs[i];
1801               for (k=1; k<=s; k++)
1802               {
1803                    temp=subst(temp,x(k),x_var(j,k,s));
1804               }
1805               sum=sum+temp*e(j);
1806          }
1807          result=result,sum-syn[i,1];
1808     }
1809
1810     result=simplify(result,2);
1811
1812     points=points2;
1813     export points;
1814     return(result);
1815}
1816example
1817{
1818     "EXAMPLE:";  echo = 2;
1819     intvec vopt = option(get);
1820
1821     list l=FLpreprocess(3,1,11,2,"");
1822     def r=l[1];
1823     setring r;
1824     int s_work=l[2];
1825
1826     //the check matrix of [11,6,5] ternary code
1827     matrix h[5][11]=1,0,0,0,0,1,1,1,-1,-1,0,
1828          0,1,0,0,0,1,1,-1,1,0,-1,
1829          0,0,1,0,0,1,-1,1,0,1,-1,
1830          0,0,0,1,0,1,-1,0,1,-1,1,
1831          0,0,0,0,1,1,0,-1,-1,1,1;
1832     matrix g=dual_code(h);
1833     matrix x[1][6];
1834     matrix y[1][11]=encode(x,g);
1835     //disturb with 2 errors
1836     matrix rec[1][11]=errorInsert(y,list(2,4),list(1,-1));
1837
1838     //the Fitzgerald-Lax system
1839     ideal sys=sysFL(h,rec,2,1,s_work);
1840     print(sys);
1841     option(redSB);
1842     ideal red_sys=std(sys);
1843     red_sys;
1844     // read the solutions from this redGB
1845     // the points are (0,0,1) and (0,1,0) with error values 1 and -1 resp.
1846     // use list points to find error positions;
1847     points;
1848
1849     option(set,vopt);
1850}
1851
1852///////////////////////////////////////////////////////////////////////////////
1853// preprocessing steps for the Fitzgerald-Lax scheme
1854proc FLpreprocess (int p, int e, int n, int t, string minp)
1855{
1856     ring r1=p,x,dp;
1857     int s=1;
1858     while(p^(s*e)<n)
1859     {
1860          s++;
1861     }
1862     list var_ord;
1863     int i,j;
1864     int count=1;
1865     for (i=s; i>=1; i--)
1866     {
1867          var_ord[count]=string("x("+string(i)+")");
1868          count++;
1869     }
1870     for (i=t; i>=1; i--)
1871     {
1872          var_ord[count]=string("e("+string(i)+")");
1873          count++;
1874          for (j=s; j>=1; j--)
1875          {
1876               var_ord[count]=string("x1("+string(s*(i-1)+j)+")");
1877               count++;
1878          }
1879     }
1880
1881     list rl;
1882     list tmp;
1883
1884     if (e>1)
1885     {
1886          rl[1]=tmp;
1887          rl[1][1]=p;
1888          rl[1][2]=tmp;
1889          rl[1][2][1]=string("a");
1890          rl[1][3]=tmp;
1891          rl[1][3][1]=tmp;
1892          rl[1][3][1][1]=string("lp");
1893          rl[1][3][1][2]=1;
1894          rl[1][4]=ideal(0);
1895     } else {
1896          rl[1]=p;
1897     }
1898
1899     rl[2]=var_ord;
1900
1901     rl[3]=tmp;
1902     rl[3][1]=tmp;
1903     rl[3][1][1]=string("lp");
1904     intvec v=1;
1905     for (i=1; i<=size(var_ord)-1; i++)
1906     {
1907          v=v,1;
1908     }
1909     rl[3][1][2]=v;
1910     rl[3][2]=tmp;
1911     rl[3][2][1]=string("C");
1912     rl[3][2][2]=intvec(0);
1913
1914     rl[4]=ideal(0);
1915
1916     def r2=ring(rl);
1917     setring r2;
1918     list l=ringlist(r2);
1919     if (e>1)
1920     {
1921          execute(string("poly f="+minp));
1922          ideal id=f;
1923          l[1][4]=id;
1924     }
1925
1926     def r=ring(l);
1927     setring r;
1928
1929     return(list(r,s));
1930}
1931
1932///////////////////////////////////////////////////////////////////////////////
1933// imitating two indeces
1934static proc x_var (int i, int j, int s)
1935{
1936     return(x1(s*(i-1)+j));
1937}
1938
1939///////////////////////////////////////////////////////////////////////////////
1940// random vector of length n with entries from p^e, p the characteristic
1941static proc randomvector(int n, int e)
1942{
1943    int i;
1944    matrix result[n][1];
1945    for (i=1; i<=n; i++)
1946    {
1947        result[i,1]=asElement(random_prime_vector(e));
1948    }
1949    return(result);
1950}
1951
1952///////////////////////////////////////////////////////////////////////////////
1953// "convert" representation of an element from the field extension from vector
1954//to an elelemnt
1955static proc asElement(list l)
1956{
1957  number s;
1958  int i;
1959  number w=1;
1960  if (size(l)>1) {w=par(1);}
1961  for (i=0; i<=size(l)-1; i++)
1962  {
1963    s=s+w^i*l[i+1];
1964  }
1965  return(s);
1966}
1967
1968///////////////////////////////////////////////////////////////////////////////
1969// random vector of length n with entries from p, p the characteristic
1970static proc random_prime_vector (int n)
1971{
1972  list rl=ringlist(basering);
1973  int i, charac;
1974  for (i=2; i<=rl[1][1]; i++)
1975  {
1976    if (rl[1][1] mod i ==0)
1977    {
1978      break;
1979    }
1980  }
1981  charac=i;
1982
1983  list l;
1984
1985  for (i=1; i<=n; i++)
1986  {
1987    l=l+list(random(0,charac-1));
1988  }
1989  return(l);
1990}
1991
1992///////////////////////////////////////////////////////////////////////////////
1993
1994proc decodeRandomFL(int n, int redun, int p, int e, int t, int ncodes, int ntrials, string minpol)
1995"USAGE:    decodeRandomFL(redun,p,e,n,t,ncodes,ntrials,minpol);
1996@format
1997          - n is length of codes generated,
1998          - redun = redundancy of codes generated,
1999          - p is characteristics,
2000          - e is the extension degree,
2001          - t is the number of errors to correct,
2002          - ncodes is the number of random codes to be processed,
2003          - ntrials is the number of received vectors per code to be corrected,
2004          - minpol: due to some pecularities of SINGULAR one needs to provide
2005          minimal polynomial for the extension explicitly
2006@end format
2007RETURN:    nothing
2008EXAMPLE:   example decodeRandomFL; shows an example
2009"
2010{
2011 intvec vopt = option(get);
2012
2013 list l=FLpreprocess(p,e,n,t,minpol);
2014
2015 def r=l[1];
2016 int s_work=l[2];
2017 export(s_work);
2018 setring r;
2019
2020 int i,j;
2021 matrix h, g, word, y, rec;
2022 int dist, tim, tim2, tim3, timdist, timdec, timdist2, timdec2, timdec3;
2023 ideal sys, sys2, sys3;
2024 list tmp;
2025
2026 option(redSB);
2027 matrix z[1][n];
2028
2029 for (i=1; i<=ncodes; i++)
2030 {
2031     h=randomCheck(redun,n,e);
2032     g=dual_code(h);
2033     tim2=rtimer;
2034     tim3=rtimer;
2035
2036     //---------------- generate the template system -----------------------
2037     sys=sysFL(h,z,t,e,s_work);
2038     timdec3=timdec3+rtimer-tim3;
2039
2040     //------ modifying the template according to the received word ---------
2041     for (j=1; j<=ntrials; j++)
2042     {
2043          word=randomvector(n-redun,1);
2044          y=encode(transpose(word),g);
2045          rec=errorRand(y,t,e);
2047          tim=rtimer;
2048          sys3=std(sys2);
2049          timdec=timdec+rtimer-tim;
2050     }
2051     timdec2=timdec2+rtimer-tim2;
2052 }
2053
2054 printf("Time for decoding: %p", timdec2);
2055 printf("Time for GB in decoding: %p", timdec);
2056 printf("Time for generating Fitzgerald-Lax system during decoding: %p", timdec3);
2057
2058 option(set,vopt);
2059}
2060example
2061{
2062     "EXAMPLE:";  echo = 2;
2063
2064     // correcting one error for one random binary code of length 25,
2065     //redundancy 14; 300 words are processed
2066     decodeRandomFL(25,14,2,1,1,1,300,"");
2067}
2068
2069///////////////////////////////////////////////////////////////////////////////
2070// add syndrome values to the template system in FL
2071static proc LF_add_synd (matrix rec, matrix check, ideal sys)
2072{
2073     int redun=nrows(check);
2074     ideal result=sys;
2075     matrix s[redun][1]=syndrome(check,rec);
2076     for (int i=size(sys)-redun+1; i<=size(sys); i++)
2077     {
2078          result[i]=result[i]-s[i-size(sys)+redun,1];
2079     }
2080     return(result);
2081}
2082
2083
2084/*
2085//////////////     SOME RELATIVELY EASY EXAMPLES    //////////////
2086///////////////////  THAT RUN AROUND ONE MINUTE  ////////////////
2087
2088"EXAMPLE:";  echo = 2;
2089int q=128; int n=120; int redun=n-30;
2090ring r=(q,a),x,dp;
2091decodeRandom(n,redun,1,1,6);
2092
2093int q=128; int n=120; int redun=n-20;
2094ring r=(q,a),x,dp;
2095decodeRandom(n,redun,1,1,9);
2096
2097int q=128; int n=120; int redun=n-10;
2098ring r=(q,a),x,dp;
2099decodeRandom(n,redun,1,1,19);
2100
2101int q=256; int n=150; int redun=n-10;
2102ring r=(q,a),x,dp;
2103decodeRandom(n,redun,1,1,22);
2104
2105//////////////     SOME HARD EXAMPLES    //////////////////////
2106//////      THAT MAYBE WILL BE DOABLE LATER     ///////////////
2107
21081.) These random instances are not doable in <=1000 sec.
2109
2110"EXAMPLE:";  echo = 2;
2111int q=128; int n=120; int redun=n-40;
2112ring r=(q,a),x,dp;
2113decodeRandom(n,redun,1,1,6);
2114
2115redun=n-30;
2116decodeRandom(n,redun,1,1,8);
2117
2118redun=n-20;
2119decodeRandom(n,redun,1,1,12);
2120
2121redun=n-10;
2122decodeRandom(n,redun,1,1,24);
2123
2124int q=256; int n=150; int redun=n-10;
2125ring r=(q,a),x,dp;
2126decodeRandom(n,redun,1,1,26);
2127
2128
21292.) Generic decoding is hard!
2130
2131int q=32; int n=31; int redun=n-16; int t=3;
2132ring r=(q,a),(V(1..n),U(n..1),s(redun..1)),(dp(n),lp(n),dp(redun));
2133matrix check[redun][n]= 1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,
21340,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,
21350,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
21360,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,
21370,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,
21380,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
21390,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,
21400,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,
21410,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,
21420,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,
21430,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,
21441,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,
21450,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,
21461,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,
21471,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,
21480,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,
21490,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
21500,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,
21510,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,
21520,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,
21530,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,
21540,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,
21550,1,0,1,0,0,1,0,0,1;
2156matrix rec[1][n];
2157
2158def A=sysQE(check,rec,t,1,2);
2159setring A;
2160print(qe);
2161ideal red_qe=stdfglm(qe);
2162
2163*/
Note: See TracBrowser for help on using the repository browser.