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

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