/////////////////////////////////////////////////////////////////////////////// version="$Id$"; category="Coding theory"; info=" LIBRARY: decodegb.lib Decoding and min distance of linear codes with GB AUTHOR: Stanislav Bulygin, bulygin@mathematik.uni-kl.de OVERVIEW: In this library we generate several systems used for decoding cyclic codes and finding their minimum distance. Namely, we work with the Cooper's philosophy and generalized Newton identities. The origindeal method of quadratic equations is worked out here as well. We also (for comparison) enable to work with the system of Fitzgerald-Lax. We provide some auxiliary functions for further manipulations and decoding. For an overview of the methods mentioned above @ref{Decoding codes with Groebner bases}. For the vanishing ideal computation the algorithm of Farr and Gao is implemented. PROCEDURES: sysCRHT(..); generates the CRHT-ideal as in Cooper's philosophy sysCRHTMindist(..); CRHT-ideal to find the minimum distance in the binary case sysNewton(..); generates the ideal with the generalized Newton identities sysBin(..); generates Bin system using Waring function encode(x,g); encodes given message x with the given generator matrix g syndrome(h,c); computes a syndrome w.r.t. the given check matrix sysQE(..); generates the system of quadratic equations for decoding errorInsert(..); inserts errors in a word errorRand(y,num,e); inserts random errors in a word randomCheck(m,n,e); generates a random check matrix genMDSMat(n,a); generates an MDS (actually an RS) matrix mindist(check); computes the minimum distance of a code decode(rec); decoding of a word using the system of quadratic equations decodeRandom(..); a procedure for manipulation with random codes decodeCode(..); a procedure for manipulation with the given code vanishId(points); computes the vanishing ideal for the given set of points sysFL(..); generates the Fitzgerald-Lax system decodeRandomFL(..); manipulation with random codes via Fitzgerald-Lax KEYWORDS: Cyclic code; Linear code; Decoding; Minimum distance; Groebner bases, decodeGB "; LIB "linalg.lib"; LIB "brnoeth.lib"; /////////////////////////////////////////////////////////////////////////////// // creates a list result, where result[i]=i, 1<=i<=n static proc lis (int n) { list result; if (n<=0) {print("ERRORlis");} for (int i=1; i<=n; i++) { result=result+list(i); } return(result); } /////////////////////////////////////////////////////////////////////////////// // creates a list of all combinations without repititions of m objects out of n static proc combinations (int m, int n) { list result; if (m>n) {print("ERRORcombinations");} if (m==n) {result[size(result)+1]=lis(m);return(result);} if (m==0) {result[size(result)+1]=list();return(result);} list temp=combinations(m-1,n-1); for (int i=1; i<=size(temp); i++) { temp[i]=temp[i]+list(n); } result=combinations(m,n-1)+temp; return(result); } /////////////////////////////////////////////////////////////////////////////// // the polynomial for Sala's restrictions static proc p_poly(int n, int a, int b) { poly f; for (int i=0; i<=n-1; i++) { f=f+Z(a)^i*Z(b)^(n-1-i); } return(f); } /////////////////////////////////////////////////////////////////////////////// proc sysCRHT (int n, list defset, int e, int q, int m, list #) "USAGE: sysCRHT(n,defset,e,q,m,[k]); n,e,q,m,k are int, defset list of int's @format - n length of the cyclic code, - defset is a list representing the defining set, - e the error-correcting capacity, - q field size - m degree extension of the splitting field, - if k>0 additional equations representing the fact that every two error positions are either different or at least one of them is zero @end format RETURN: the ring to work with the CRHT-ideal (with Sala's additions), containig an ideal with name 'crht' THEORY: Based on 'defset' of the given cyclic code, the procedure constructs the corresponding Cooper-Reed-Heleseth-Truong ideal 'crht'. With its help one can solve the decoding problem. For basics of the method @ref{Cooper philosophy}. SEE ALSO: sysNewton, sysBin EXAMPLE: example sysCRHT; shows an example " { int r=size(defset); ring @crht=(q,a),(Y(e..1),Z(1..e),X(r..1)),lp; ideal crht; int i,j; poly sum; int k; if ( size(#) > 0) { k = #[1]; } //---------------------- add check equations -------------------------- for (i=1; i<=r; i++) { sum=0; for (j=1; j<=e; j++) { sum=sum+Y(j)*Z(j)^defset[i]; } crht[i]=sum-X(i); } //--------------------- field equations on syndromes ------------------ for (i=1; i<=r; i++) { crht=crht,X(i)^(q^m)-X(i); } //------ restrictions on error-locations: n-th roots of unity ---------- for (i=1; i<=e; i++) { crht=crht,Z(i)^(n+1)-Z(i); } for (i=1; i<=e; i++) { crht=crht,Y(i)^(q-1)-1; } //--------- add Sala's additional conditions if necessary -------------- if ( k > 0 ) { for (i=1; i<=e; i++) { for (j=i+1; j<=e; j++) { crht=crht,Z(i)*Z(j)*p_poly(n,i,j); } } } export crht; return(@crht); } example { "EXAMPLE:"; echo=2; // binary cyclic [15,7,5] code with defining set (1,3) intvec v = option(get); list defset=1,3; // defining set int n=15; // length int e=2; // error-correcting capacity int q=2; // basefield size int m=4; // degree extension of the splitting field int sala=1; // indicator to add additional equations def A=sysCRHT(n,defset,e,q,m); setring A; A; // shows the ring we are working in print(crht); // the CRHT-ideal option(redSB); ideal red_crht=std(crht); // reduced Groebner basis print(red_crht); //============================ A=sysCRHT(n,defset,e,q,m,sala); setring A; print(crht); // CRHT-ideal with additional equations from Sala option(redSB); ideal red_crht=std(crht); // reduced Groebner basis print(red_crht); red_crht[5]; // general error-locator polynomial for this code option(set,v); } /////////////////////////////////////////////////////////////////////////////// proc sysCRHTMindist (int n, list defset, int w) "USAGE: sysCRHTMindist(n,defset,w); n,w are int, defset is list of int's @format - n length of the cyclic code, - defset is a list representing the defining set, - w is a candidate for the minimum distance @end format RETURN: the ring to work with the Sala's ideal for the minimum distance containing the ideal with name 'crht_md' THEORY: Based on 'defset' of the given cyclic code, the procedure constructs the corresponding Cooper-Reed-Heleseth-Truong ideal 'crht_md'. With its help one can find minimum distance of the code in the binary case. For basics of the method @ref{Cooper philosophy}. EXAMPLE: example sysCRHTMindist; shows an example " { int r=size(defset); ring @crht_md=2,Z(1..w),lp; ideal crht_md; int i,j; poly sum; //------------ add check equations -------------- for (i=1; i<=r; i++) { sum=0; for (j=1; j<=w; j++) { sum=sum+Z(j)^defset[i]; } crht_md[i]=sum; } //----------- locations are n-th roots of unity ------------ for (i=1; i<=w; i++) { crht_md=crht_md,Z(i)^n-1; } //------------ adding conditions on locations being different ------------ for (i=1; i<=w; i++) { for (j=i+1; j<=w; j++) { crht_md=crht_md,Z(i)*Z(j)*p_poly(n,i,j); } } export crht_md; return(@crht_md); } example { "EXAMPLE:"; echo=2; intvec v = option(get); // binary cyclic [15,7,5] code with defining set (1,3) list defset=1,3; // defining set int n=15; // length int d=5; // candidate for the minimum distance def A=sysCRHTMindist(n,defset,d); setring A; A; // shows the ring we are working in print(crht_md); // the Sala's ideal for mindist option(redSB); ideal red_crht_md=std(crht_md); print(red_crht_md); // reduced Groebner basis option(set,v); } /////////////////////////////////////////////////////////////////////////////// // slightly modified mod function static proc mod_ (int n, int m) { n=n mod m; if (n<=0){ return(n+m);} return(n); } /////////////////////////////////////////////////////////////////////////////// proc sysNewton (int n, list defset, int t, int q, int m, list #) "USAGE: sysNewton (n,defset,t,q,m,[tr]); n,t,q,m,tr int, defset is list int's @format - n is length, - defset is the defining set, - t is the number of errors, - q is basefield size, - m is degree extension of the splitting field, - if tr>0 it indicates that Newton identities in triangular form should be constructed @end format RETURN: the ring to work with the generalized Newton identities (in triangular form if applicable) containing the ideal with name 'newton' THEORY: Based on 'defset' of the given cyclic code, the procedure constructs the corresponding ideal 'newton' with the generalized Newton identities. With its help one can solve the decoding problem. For basics of the method @ref{Generalized Newton identities}. SEE ALSO: sysCRHT, sysBin EXAMPLE: example sysNewton; shows an example " { string s="ring @newton=("+string(q)+",a),("; int i,j; int flag; int tr; if (size(#)>0) { tr=#[1]; } for (i=n; i>=1; i--) { for (j=1; j<=size(defset); j++) { flag=1; if (i==defset[j]) { flag=0; break; } } if (flag) { s=s+"S("+string(i)+"),"; } } s=s+"sigma(1.."+string(t)+"),"; for (i=size(defset); i>=2; i--) { s=s+"S("+string(defset[i])+"),"; } s=s+"S("+string(defset[1])+")),lp;"; execute(s); ideal newton; poly sum; //------------ generate generalized Newton identities ---------- if (tr) { for (i=1; i<=t; i++) { sum=0; for (j=1; j<=i-1; j++) { sum=sum+sigma(j)*S(i-j); } newton=newton,S(i)+sum+number(i)*sigma(i); } } else { for (i=1; i<=t; i++) { sum=0; for (j=1; j<=t; j++) { sum=sum+sigma(j)*S(mod_(i-j,n)); } newton=newton,S(i)+sum; } } for (i=1; i<=n-t; i++) { sum=0; for (j=1; j<=t; j++) { sum=sum+sigma(j)*S(t+i-j); } newton=newton,S(t+i)+sum; } //----------- add field equations on sigma's -------------- for (i=1; i<=t; i++) { newton=newton,sigma(i)^(q^m)-sigma(i); } //----------- add conjugacy relations ------------------ for (i=1; i<=n; i++) { newton=newton,S(i)^q-S(mod_(q*i,n)); } newton=simplify(newton,2); export newton; return(@newton); } example { "EXAMPLE:"; echo = 2; // Newton identities for a binary 3-error-correcting cyclic code of //length 31 with defining set (1,5,7) int n=31; // length list defset=1,5,7; //defining set int t=3; // number of errors int q=2; // basefield size int m=5; // degree extension of the splitting field int tr=1; // indicator of triangular form of Newton identities def A=sysNewton(n,defset,t,q,m); setring A; A; // shows the ring we are working in print(newton); // generalized Newton identities //=============================== A=sysNewton(n,defset,t,q,m,tr); setring A; print(newton); // generalized Newton identities in triangular form } /////////////////////////////////////////////////////////////////////////////// // forms a list of special combinations needed for computation of Waring's //function static proc combinations_sum (int m, int n) { list result; list comb=combinations(m-1,n+m-1); int i,j,flag,count; list interm=comb; for (i=1; i<=size(comb); i++) { interm[i][1]=comb[i][1]-1; for (j=2; j<=m-1; j++) { interm[i][j]=comb[i][j]-comb[i][j-1]-1; } interm[i][m]=n+m-comb[i][m-1]-1; flag=1; count=2; while ((flag)&&(count<=m)) { if (interm[i][count] mod count != 0) {flag=0;} count++; } if (flag) { for (j=2; j<=m; j++) { interm[i][j]=interm[i][j] div j; } result[size(result)+1]=interm[i]; } } return(result); } /////////////////////////////////////////////////////////////////////////////// //if n=q^e*m, m and n are coprime, then return e static proc exp_count (int n, int q) { int flag=1; int result=0; while(flag) { if (n mod q != 0) {flag=0;} else {n=n div q; result++;} } return(result); } /////////////////////////////////////////////////////////////////////////////// proc sysBin (int v, list Q, int n, list #) "USAGE: sysBin (v,Q,n,[odd]); v,n,odd are int, Q is list of int's @format - v a number if errors, - Q is a defining set of the code, - n the length, - odd is an additional parameter: if set to 1, then the defining set is enlarged by odd elements, which are 2^(some power)*(some elment in the def.set) mod n @end format RETURN: the ring with the resulting system called 'bin' THEORY: Based on Q of the given cyclic code, the procedure constructs the corresponding ideal 'bin' with the use of the Waring function. With its help one can solve the decoding problem. For basics of the method @ref{Generalized Newton identities}. SEE ALSO: sysNewton, sysCRHT EXAMPLE: example sysBin; shows an example " { int odd; if (size(#)>0) { odd=#[1]; } //ring r=2,(sigma(1..v),S(1..n)),(lp(v),dp(n)); ring r=2,(S(1..n),sigma(1..v)),lp; list cyclot; ideal result; int i,j,k,s; list comb; poly sum_, mon; int count1, count2, upper, coef_, flag, gener; list Q_update; if (odd==1) { for (i=1; i<=n; i++) { cyclot[i]=0; } for (i=1; i<=size(Q); i++) { flag=1; gener=Q[i]; while(flag) { cyclot[gener]=1; gener=2*gener mod n; if (gener == Q[i]) {flag=0;} } } for (i=1; i<=n; i++) { if ((cyclot[i] == 1)&&(i mod 2 == 1)) {Q_update[size(Q_update)+1]=i;} } } else { Q_update=Q; } //---- form polynomials for the Bin system via Waring's function --------- for (i=1; i<=size(Q_update); i++) { comb=combinations_sum(v,Q_update[i]); sum_=0; for (k=1; k<=size(comb); k++) { upper=0; for (j=1; j<=v; j++) { upper=upper+comb[k][j]; } count1=0; for (j=2; j<=upper-1; j++) { count1=count1+exp_count(j,2); } count1=count1+exp_count(Q_update[i],2); count2=0; for (j=1; j<=v; j++) { for (s=2; s<=comb[k][j]; s++) { count2=count2+exp_count(s,2); } } if (count1count2) {coef_=0;} if (count1 == count2) {coef_=1;} mon=1; for (j=1; j<=v; j++) { mon=mon*sigma(j)^(comb[k][j]); } sum_=sum_+coef_*mon; } result=result,S(Q_update[i])-sum_; } ideal bin=simplify(result,2); export bin; return(r); } example { "EXAMPLE:"; echo = 2; // [31,16,7] quadratic residue code list l=1,5,7,9,19,25; // we do not need even synromes here def A=sysBin(3,l,31); setring A; print(bin); } /////////////////////////////////////////////////////////////////////////////// proc encode (matrix x, matrix g) "USAGE: encode (x, g); x a row vector (message), and g a generator matrix RETURN: corresponding codeword EXAMPLE: example encode; shows an example " { if (nrows(x)>1) {print("ERRORencode1!");} if (ncols(x)!=nrows(g)) {print("ERRORencode2!");} return(x*g); } example { "EXAMPLE:"; echo = 2; ring r=2,x,dp; matrix x[1][4]=1,0,1,0; matrix g[4][7]=1,0,0,0,0,1,1, 0,1,0,0,1,0,1, 0,0,1,0,1,1,1, 0,0,0,1,1,1,0; //encode x with the generator matrix g print(encode(x,g)); } /////////////////////////////////////////////////////////////////////////////// proc syndrome (matrix h, matrix c) "USAGE: syndrome (h, c); h a check matrix, c a row vector (codeword) RETURN: corresponding syndrome EXAMPLE: example syndrome; shows an example " { if (nrows(c)>1) {print("ERRORsyndrome1!");} if (ncols(c)!=ncols(h)) {print("ERRORsyndrome2!");} return(h*transpose(c)); } example { "EXAMPLE:"; echo = 2; ring r=2,x,dp; matrix x[1][4]=1,0,1,0; matrix g[4][7]=1,0,0,0,0,1,1, 0,1,0,0,1,0,1, 0,0,1,0,1,1,1, 0,0,0,1,1,1,0; //encode x with the generator matrix g matrix c=encode(x,g); // disturb c[1,3]=0; //compute syndrome //corresponding check matrix 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; print(syndrome(check,c)); c[1,3]=1; //now c is a codeword print(syndrome(check,c)); } /////////////////////////////////////////////////////////////////////////////// // (coordinatewise) star product of two vectors static proc star(matrix m, int i, int j) { matrix result[ncols(m)][1]; for (int k=1; k<=ncols(m); k++) { result[k,1]=m[i,k]*m[j,k]; } return(result); } /////////////////////////////////////////////////////////////////////////////// proc sysQE(matrix check, matrix y, int t, list #) "USAGE: sysQE(check,y,t,[fieldeq,formal]);check,y matrix;t,fieldeq,formal int @format - check is a parity check matrix of the code - y is a received word, - t the number of errors to be corrected, - if fieldeq=1, then field equations are added, - if formal=0, field equations on (known) syndrome variables are not added, in order to add them (note that the exponent should be equal to the number of elements in the INITIAL alphabet) one needs to set formal>0 for the exponent @end format RETURN: the ring to work with together with the resulting system called 'qe' THEORY: Based on 'check' of the given linear code, the procedure constructs the corresponding ideal that gives an opportunity to compute unknown syndrome of the received word y. After computing the unknown syndromes one is able to solve the decoding problem. For basics of the method @ref{Decoding method based on quadratic equations}. SEE ALSO: sysFL EXAMPLE: example sysQE; shows an example " { int fieldeq; int formal; if (size(#)>0) { fieldeq=#[1]; } if (size(#)>1) { formal=#[2]; } def br=basering; list rl=ringlist(br); int red=nrows(check); int n=ncols(check); int q=rl[1][1]; if (formal==0) { ring work=(q,a),(V(1..t),U(1..n)),dp; } else { ring work=(q,a),(V(1..t),U(1..n),s(1..red)),(dp(t),lp(n),dp(red)); } matrix check=imap(br,check); matrix y=imap(br,y); matrix h_full=genMDSMat(n,a); matrix h=submat(h_full,1..red,1..n); if (nrows(y)!=1) {print("ERROR1Pell");} if (ncols(y)!=n) {print("ERROR2Pell");} ideal result; list c; list a; list tmp,tmp2; int i,j,l,k; number sum,prod,sig; poly sum1,sum2,sum3; for (i=1; i<=n; i++) { c[i]=tmp; } matrix transf=inverse(transpose(h_full)); //------ expression matrix of check vectors w.r.t. the MDS basis ----------- for (i=1; i<=red ; i++) { a[i]=transpose(submat(check,i..i,1..n)); a[i]=transf*a[i]; } //----------- compute the structure constants ------------------------ matrix te[n][1]; for (i=1; i<=n; i++) { for (j=1; j<=t+1; j++) { if ((j0) { e=#[1]; } option(redSB); def br=basering; matrix h_full=genMDSMat(n,a); matrix z[1][ncols(h_full)]; //------------------ determine error-correction capacity ------------------- for (i=1; i<=ncodes; i++) { setring br; h=randomCheck(redun,n,1); "check matrix:"; print(h); if (e>0) { t=e; } else { tmp=mindist(h); dist=tmp; printf("d= %p",dist); t=(dist-1) div 2; } //------------- generate the template system ---------------------- def A=sysQE(h,z,t); setring A; matrix word,y,rec; ideal sys2,sys3; matrix h=imap(br,h); matrix g=dual_code(h); ideal sys=qe; print("The system is generated"); //------ modify the template according to every received word -------------- for (j=1; j<=ntrials; j++) { word=randomvector(n-redun,1); y=encode(transpose(word),g); print("Codeword:"); print(y); rec=errorRand(y,t,1); print("Received word:"); print(rec); sys2=add_synd(rec,h,redun,sys); option(redSB); sys3=std(sys2); print("The Groebenr basis of the QE system:"); print(sys3); } kill A; option(set,vopt); } } example { "EXAMPLE:"; echo = 2; int q=32; int n=25; int redun=n-11; int t=redun+1; ring r=(q,a),x,dp; // correct 2 errors in 2 random binary codes, 3 trials each decodeRandom(n,redun,2,3,2); } /////////////////////////////////////////////////////////////////////////////// proc decodeCode(matrix check, int ntrials, list #) "USAGE: decodeCode(check, ntrials, [e]); check matrix, ntrials,e int @format - check is a parity check matrix for the code, - ntrials is the number of received vectors per code to be corrected. - If e is given it sets the correction capacity explicitly. It should be used in case one expects some lower bound, otherwise the procedure tries to compute the real minimum distance to find out the error-correction capacity @end format RETURN: nothing; EXAMPLE: example decodeCode; shows an example " { intvec vopt = option(get); int n=ncols(check); int redun=nrows(check); int i,j; matrix h; int dist, t; ideal sys; int tmp; int e; if (size(#)>0) { e=#[1]; } option(redSB); def br=basering; matrix h_full=genMDSMat(n,a); matrix z[1][ncols(h_full)]; setring br; h=check; "check matrix:"; print(h); //------------------ determine error-correction capacity ------------------- if (e>0) { t=e; } else { tmp=mindist(h); dist=tmp; printf("d= %p",dist); t=(dist-1) div 2; } //------------- generate the template system ---------------------- def A=sysQE(h,z,t); setring A; matrix word,y,rec; ideal sys2,sys3; matrix h=imap(br,h); matrix g=dual_code(h); ideal sys=qe; print("The system is generated"); //--- modify the template according to every received word --------------- for (j=1; j<=ntrials; j++) { word=randomvector(n-redun,1); y=encode(transpose(word),g); print("Codeword:"); print(y); rec=errorRand(y,t,1); print("Received word:"); print(rec); sys2=add_synd(rec,h,redun,sys); option(redSB); sys3=std(sys2); print("Groebner basis of the QE system:"); print(sys3); } option(set,vopt); } example { "EXAMPLE:"; echo = 2; int q=32; int n=25; int redun=n-11; int t=redun+1; ring r=(q,a),x,dp; matrix check=randomCheck(redun,n,1); // correct 2 errors in using the code above, 3 trials decodeCode(check,3,2); } /////////////////////////////////////////////////////////////////////////////// // adding syndrome values to the template system static proc add_synd (matrix rec, matrix check, int redun, ideal sys) { ideal result=sys; matrix s[redun][1]=syndrome(check,rec); for (int i=1; i<=redun; i++) { result[i]=result[i]-s[i,1]; } return(result); } /////////////////////////////////////////////////////////////////////////////// // evaluate a polynomial at a given point static proc ev (poly f, matrix p) { if (ncols(p)>1) {ERROR("not a column vector");}; int m=size(p); poly temp=f; for (int i=1; i<=m; i++) { temp=subst(temp,var(i),p[i,1]); } return(number(temp)); } /////////////////////////////////////////////////////////////////////////////// // return index of an element in the ideal where it does not vanish at the //given point static proc find_index (ideal G, matrix p) { if (ncols(p)>1) {ERROR("not a column vector");}; int i=1; int n=size(G); while(i<=n) { if (ev(G[i],p)!=0) {return(i);} i++; } return(-1); } /////////////////////////////////////////////////////////////////////////////// // convert ideal to list static proc ideal2list (ideal id) { list l; for (int i=1; i<=size(id); i++) { l[i]=id[i]; } return(l); } /////////////////////////////////////////////////////////////////////////////// // convert list to ideal static proc list2ideal (list l) { ideal id; for (int i=1; i<=size(l); i++) { id[i]=l[i]; } return(id); } /////////////////////////////////////////////////////////////////////////////// // check whether given polynomial is divisible by some leading monomial of the //ideal static proc divisible (poly m, ideal G) { for (int i=1; i<=size(G); i++) { if (m/leadmonom(G[i])!=0) {return(1);} } return(0); } /////////////////////////////////////////////////////////////////////////////// proc vanishId (list points) "USAGE: vanishId (points); point is a list of matrices 'points' is a list of points for which the vanishing ideal is to be constructed RETURN: Vanishing ideal corresponding to the given set of points EXAMPLE: example vanishId; shows an example " { int m=size(points[1]); int n=size(points); ideal G=1; int i,k,j; list temp; poly h,cur; //------------- proceed according to Farr-Gao algorithm ---------------- for (k=1; k<=n; k++) { i=find_index(G,points[k]); cur=G[i]; for(j=i+1; j<=size(G); j++) { G[j]=G[j]-ev(G[j],points[k])/ev(G[i],points[k])*G[i]; } G=simplify(G,2); temp=ideal2list(G); temp=delete(temp,i); G=list2ideal(temp); for (j=1; j<=m; j++) { if (!divisible(var(j)*leadmonom(cur),G)) { attrib(G,"isSB",1); h=NF((var(j)-points[k][j,1])*cur,G); temp=ideal2list(G); temp=insert(temp,h); G=list2ideal(temp); G=sort(G)[1]; } } } attrib(G,"isSB",1); return(G); } example { "EXAMPLE:"; echo = 2; ring r=3,(x(1..3)),dp; //generate all 3-vectors over GF(3) list points=pointsGen(3,1); list points2=convPoints(points); //grasps the first 11 points list p=graspList(points2,1,11); print(p); //construct the vanishing ideal ideal id=vanishId(p); print(id); } /////////////////////////////////////////////////////////////////////////////// // construct the list of all vectors of length m with elements in p^e, where p // is theharacteristic proc pointsGen (int m, int e) { if (e>1) { list result; int count=1; int i,j; list l=ringlist(basering); int charac=l[1][1]; number a=par(1); list tmp; for (i=1; i<=charac^(e*m); i++) { result[i]=tmp; } if (m==1) { result[count][m]=0; count++; for (j=1; j<=charac^(e)-1; j++) { result[count][m]=a^j; count++; } return(result); } list prev=pointsGen(m-1,e); for (i=1; i<=size(prev); i++) { result[count]=prev[i]; result[count][m]=0; count++; for (j=1; j<=charac^(e)-1; j++) { result[count]=prev[i]; result[count][m]=a^j; count++; } } return(result); } if (e==1) { list result; int count=1; int i,j; list l=ringlist(basering); int charac=l[1][1]; list tmp; for (i=1; i<=charac^m; i++) { result[i]=tmp; } if (m==1) { for (j=0; j<=charac-1; j++) { result[count][m]=number(j); count++; } return(result); } list prev=pointsGen(m-1,e); for (i=1; i<=size(prev); i++) { for (j=0; j<=charac-1; j++) { result[count]=prev[i]; result[count][m]=number(j); count++; } } return(result); } } /////////////////////////////////////////////////////////////////////////////// // convert list to a column vector static proc list2vec (list l) { matrix m[size(l)][1]; for (int i=1; i<=size(l); i++) { m[i,1]=l[i]; } return(m); } /////////////////////////////////////////////////////////////////////////////// // convert all the point in the list with list2vec proc convPoints (list points) { for (int i=1; i<=size(points); i++) { points[i]=list2vec(points[i]); } return(points); } /////////////////////////////////////////////////////////////////////////////// // extracts elements from l in the range m..n proc graspList (list l, int m, int n) { list result; int count=1; for (int i=m; i<=n; i++) { result[count]=l[i]; count++; } return(result); } /////////////////////////////////////////////////////////////////////////////// // "characteristic" polynomial static proc xi_gen (matrix p, int e, int s) { poly prod=1; list rl=ringlist(basering); int charac=rl[1][1]; int l; for (l=1; l<=s; l++) { prod=prod*(1-(var(l)-p[l,1])^(charac^e-1)); } return(prod); } /////////////////////////////////////////////////////////////////////////////// // generating polynomials in Fitzgerald-Lax construction static proc gener_funcs (matrix check, list points, int e, ideal id, int s) { int n=ncols(check); if (n!=size(points)) {ERROR("Incompatible sizes of check and points");} ideal xi; int i,j; for (i=1; i<=n; i++) { xi[i]=xi_gen(points[i],e,s); } ideal result; int m=nrows(check); poly sum; for (i=1; i<=m; i++) { sum=0; for (j=1; j<=n; j++) { sum=sum+check[i,j]*xi[j]; } result[i]=NF(sum,id); } return(result); } /////////////////////////////////////////////////////////////////////////////// proc sysFL (matrix check, matrix y, int t, int e, int s) "USAGE: sysFL (check,y,t,e,s); check,y matrix, t,e,s int @format - check is a parity check matrix of the code, - y is a received word, - t the number of errors to correct, - e is the extension degree, - s is the dimension of the point for the vanishing ideal @end format RETURN: the system of Fitzgerald-Lax for the given decoding problem THEORY: Based on 'check' of the given linear code, the procedure constructs the corresponding ideal constructed with a generalization of Cooper's philosophy. For basics of the method @ref{Fitzgerald-Lax method}. SEE ALSO: sysQE EXAMPLE: example sysFL; shows an example " { list rl=ringlist(basering); int charac=rl[1][1]; int n=ncols(check); int m=nrows(check); list points=pointsGen(s,e); list points2=convPoints(points); list p=graspList(points2,1,n); ideal id=vanishId(p,e); ideal funcs=gener_funcs(check,p,e,id,s); ideal result; poly temp; int i,j,k; //--------------- add vanishing realtions --------------------- for (i=1; i<=t; i++) { for (j=1; j<=size(id); j++) { temp=id[j]; for (k=1; k<=s; k++) { temp=subst(temp,var(k),x_var(i,k,s)); } result=result,temp; } } //--------------- add field equations -------------------- for (i=1; i<=t; i++) { for (k=1; k<=s; k++) { result=result,x_var(i,k,s)^(charac^e)-x_var(i,k,s); } } for (i=1; i<=t; i++) { result=result,e(i)^(charac^e-1)-1; } result=simplify(result,8); //--------------- add check realtions -------------------- poly sum; matrix syn[m][1]=syndrome(check,y); for (i=1; i<=size(funcs); i++) { sum=0; for (j=1; j<=t; j++) { temp=funcs[i]; for (k=1; k<=s; k++) { temp=subst(temp,var(k),x_var(j,k,s)); } sum=sum+temp*e(j); } result=result,sum-syn[i,1]; } result=simplify(result,2); points=points2; export points; return(result); } example { "EXAMPLE:"; echo = 2; intvec vopt = option(get); list l=FLpreprocess(3,1,11,2,""); def r=l[1]; setring r; int s_work=l[2]; //the check matrix of [11,6,5] ternary code matrix h[5][11]=1,0,0,0,0,1,1,1,-1,-1,0, 0,1,0,0,0,1,1,-1,1,0,-1, 0,0,1,0,0,1,-1,1,0,1,-1, 0,0,0,1,0,1,-1,0,1,-1,1, 0,0,0,0,1,1,0,-1,-1,1,1; matrix g=dual_code(h); matrix x[1][6]; matrix y[1][11]=encode(x,g); //disturb with 2 errors matrix rec[1][11]=errorInsert(y,list(2,4),list(1,-1)); //the Fitzgerald-Lax system ideal sys=sysFL(h,rec,2,1,s_work); print(sys); option(redSB); ideal red_sys=std(sys); red_sys; // read the solutions from this redGB // the points are (0,0,1) and (0,1,0) with error values 1 and -1 resp. // use list points to find error positions; points; option(set,vopt); } /////////////////////////////////////////////////////////////////////////////// // preprocessing steps for the Fitzgerald-Lax scheme proc FLpreprocess (int p, int e, int n, int t, string minp) { ring r1=p,x,dp; int s=1; while(p^(s*e)=1; i--) { var_ord[count]=string("x("+string(i)+")"); count++; } for (i=t; i>=1; i--) { var_ord[count]=string("e("+string(i)+")"); count++; for (j=s; j>=1; j--) { var_ord[count]=string("x1("+string(s*(i-1)+j)+")"); count++; } } list rl; list tmp; if (e>1) { rl[1]=tmp; rl[1][1]=p; rl[1][2]=tmp; rl[1][2][1]=string("a"); rl[1][3]=tmp; rl[1][3][1]=tmp; rl[1][3][1][1]=string("lp"); rl[1][3][1][2]=1; rl[1][4]=ideal(0); } else { rl[1]=p; } rl[2]=var_ord; rl[3]=tmp; rl[3][1]=tmp; rl[3][1][1]=string("lp"); intvec v=1; for (i=1; i<=size(var_ord)-1; i++) { v=v,1; } rl[3][1][2]=v; rl[3][2]=tmp; rl[3][2][1]=string("C"); rl[3][2][2]=intvec(0); rl[4]=ideal(0); def r2=ring(rl); setring r2; list l=ringlist(r2); if (e>1) { execute(string("poly f="+minp)); ideal id=f; l[1][4]=id; } def r=ring(l); setring r; return(list(r,s)); } /////////////////////////////////////////////////////////////////////////////// // imitating two indeces static proc x_var (int i, int j, int s) { return(x1(s*(i-1)+j)); } /////////////////////////////////////////////////////////////////////////////// // random vector of length n with entries from p^e, p the characteristic static proc randomvector(int n, int e) { int i; matrix result[n][1]; for (i=1; i<=n; i++) { result[i,1]=asElement(random_prime_vector(e)); } return(result); } /////////////////////////////////////////////////////////////////////////////// // "convert" representation of an element from the field extension from vector //to an elelemnt static proc asElement(list l) { number s; int i; number w=1; if (size(l)>1) {w=par(1);} for (i=0; i<=size(l)-1; i++) { s=s+w^i*l[i+1]; } return(s); } /////////////////////////////////////////////////////////////////////////////// // random vector of length n with entries from p, p the characteristic static proc random_prime_vector (int n) { list rl=ringlist(basering); int i, charac; for (i=2; i<=rl[1][1]; i++) { if (rl[1][1] mod i ==0) { break; } } charac=i; list l; for (i=1; i<=n; i++) { l=l+list(random(0,charac-1)); } return(l); } /////////////////////////////////////////////////////////////////////////////// proc decodeRandomFL(int n, int redun, int p, int e, int t, int ncodes, int ntrials, string minpol) "USAGE: decodeRandomFL(redun,p,e,n,t,ncodes,ntrials,minpol); @format - n is length of codes generated, - redun = redundancy of codes generated, - p is the characteristic, - e is the extension degree, - t is the number of errors to correct, - ncodes is the number of random codes to be processed, - ntrials is the number of received vectors per code to be corrected, - minpol: due to some pecularities of SINGULAR one needs to provide minimal polynomial for the extension explicitly @end format RETURN: nothing EXAMPLE: example decodeRandomFL; shows an example " { intvec vopt = option(get); list l=FLpreprocess(p,e,n,t,minpol); def r=l[1]; int s_work=l[2]; export(s_work); setring r; int i,j; matrix h, g, word, y, rec; ideal sys, sys2, sys3; option(redSB); matrix z[1][n]; for (i=1; i<=ncodes; i++) { h=randomCheck(redun,n,e); g=dual_code(h); //---------------- generate the template system ----------------------- sys=sysFL(h,z,t,e,s_work); //------ modifying the template according to the received word --------- for (j=1; j<=ntrials; j++) { word=randomvector(n-redun,1); y=encode(transpose(word),g); print("Codeword:"); print(y); rec=errorRand(y,t,e); print("Received word"); print(rec); sys2=LF_add_synd(rec,h,sys); sys3=std(sys2); print("Groebner basis of the FL system:"); print(sys3); } } option(set,vopt); } example { "EXAMPLE:"; echo = 2; // correcting one error for one random binary code of length 25, // redundancy 14; 10 words are processed decodeRandomFL(25,14,2,1,1,1,10,""); } /////////////////////////////////////////////////////////////////////////////// // add syndrome values to the template system in FL static proc LF_add_synd (matrix rec, matrix check, ideal sys) { int redun=nrows(check); ideal result=sys; matrix s[redun][1]=syndrome(check,rec); for (int i=size(sys)-redun+1; i<=size(sys); i++) { result[i]=result[i]-s[i-size(sys)+redun,1]; } return(result); } /* ////////////// SOME RELATIVELY EASY EXAMPLES ////////////// /////////////////// THAT RUN AROUND ONE MINUTE //////////////// "EXAMPLE:"; echo = 2; int q=128; int n=120; int redun=n-30; ring r=(q,a),x,dp; decodeRandom(n,redun,1,1,6); int q=128; int n=120; int redun=n-20; ring r=(q,a),x,dp; decodeRandom(n,redun,1,1,9); int q=128; int n=120; int redun=n-10; ring r=(q,a),x,dp; decodeRandom(n,redun,1,1,19); int q=256; int n=150; int redun=n-10; ring r=(q,a),x,dp; decodeRandom(n,redun,1,1,22); ////////////// SOME HARD EXAMPLES ////////////////////// ////// THAT MAYBE WILL BE DOABLE LATER /////////////// 1.) These random instances are not doable in <=1000 sec. "EXAMPLE:"; echo = 2; int q=128; int n=120; int redun=n-40; ring r=(q,a),x,dp; decodeRandom(n,redun,1,1,6); redun=n-30; decodeRandom(n,redun,1,1,8); redun=n-20; decodeRandom(n,redun,1,1,12); redun=n-10; decodeRandom(n,redun,1,1,24); int q=256; int n=150; int redun=n-10; ring r=(q,a),x,dp; decodeRandom(n,redun,1,1,26); 2.) Generic decoding is hard! int q=32; int n=31; int redun=n-16; int t=3; ring r=(q,a),(V(1..n),U(n..1),s(redun..1)),(dp(n),lp(n),dp(redun)); matrix check[redun][n]= 1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1, 0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0, 0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1, 0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0, 0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0, 1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0, 1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1, 1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1, 0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1, 0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0, 0,1,0,1,0,0,1,0,0,1; matrix rec[1][n]; def A=sysQE(check,rec,t,1,2); setring A; print(qe); ideal red_qe=stdfglm(qe); */