Changeset fba86f8 in git
- Timestamp:
- Oct 8, 2008, 6:52:24 PM (15 years ago)
- Branches:
- (u'spielwiese', 'f6c3dc58b0df4bd712574325fe76d0626174ad97')
- Children:
- 7d56875cf7720719f5b6acd193768a4d5b57848f
- Parents:
- 02efaee061610b6730796d2e862e2ab63663c50a
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/decodegb.lib
r02efae rfba86f8 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: decodegb.lib,v 1. 5 2008-10-07 14:14:56bulygin Exp $";2 version="$Id: decodegb.lib,v 1.6 2008-10-08 16:52:24 bulygin Exp $"; 3 3 category="Coding theory"; 4 4 info=" 5 LIBRARY: decodedistGB.lib Generating and solving systems of polynomial equations for decoding and finding the minimum distance of linear codes6 AUTHOR S:Stanislav Bulygin, bulygin@mathematik.uni-kl.de5 LIBRARY: decodegb.lib Decoding and min distance of linear codes with GB 6 AUTHOR: Stanislav Bulygin, bulygin@mathematik.uni-kl.de 7 7 8 8 OVERVIEW: 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. 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, 15 cf. 16 @format 17 Stanislav Bulygin, Ruud Pellikaan: 'Decoding and finding the minimum 18 distance with Groebner bases: history and new insights', in 'Selected Topics 19 in Information and Coding Theory', World Scientific (2008) (to appear) (*). 20 @end format 21 For the vanishing ideal computation the algorithm of Farr and Gao is 22 implemented. 17 23 18 24 MAIN PROCEDURES: 19 sysCRHT( n,defset,e,q,m,#); generates the CRHT-ideal that followsCooper's philosophy20 sysCRHTMindist Binary(n,defset,w); generates theCRHT-ideal to find the minimum distance in the binary case21 sysNewton( n,defset,t,q,m,#); generates the ideal with the Generalized Newton identities22 sysBin( v,Q,n,#); generates Bin system as in the work of Augot et.al, cf. [*] for the reference23 encode(x,g); encodes given message with the given generator matrix24 syndrome(h, c);computes a syndrome w.r.t. the given check matrix25 sysQE( check,y,t,fieldeq,formal); generates the system of quadratic equations as in [*]26 errorInsert( y,pos,val);inserts errors in a word27 errorRand(y,num,e); 28 randomCheck(m,n,e); 29 genMDSMat(n,a); generates an MDS (actually an RS) matrix30 mindist(check); computes the minimum distance of a code31 decode(rec); decoding of a word using the system of quadratic equations32 solveForRandom(redun,ncodes,ntrials,#);a procedure for manipulation with random codes33 solveForCode(check,ntrials,#); a procedure for manipulation with the given codes34 vanishId(points); computes the vanishing ideal for the given set of points35 sysFL( check,y,t,e,s);generates the Fitzgerald-Lax system36 solveForRandomFL(n,redun,p,e,t,ncodes,ntrials,minpol); a procedure formanipulation with random codes via Fitzgerald-Lax25 sysCRHT(..); generates the CRHT-ideal as in Cooper's philosophy 26 sysCRHTMindist(..); CRHT-ideal to find the minimum distance in the binary case 27 sysNewton(..); generates the ideal with the generalized Newton identities 28 sysBin(..); generates Bin system using Waring function 29 encode(x,g); encodes given message x with the given generator matrix g 30 syndrome(h,c); computes a syndrome w.r.t. the given check matrix 31 sysQE(..); generates the system of quadratic equations for decoding 32 errorInsert(..); inserts errors in a word 33 errorRand(y,num,e); inserts random errors in a word 34 randomCheck(m,n,e); generates a random check matrix 35 genMDSMat(n,a); generates an MDS (actually an RS) matrix 36 mindist(check); computes the minimum distance of a code 37 decode(rec); decoding of a word using the system of quadratic equations 38 decodeRandom(..); a procedure for manipulation with random codes 39 decodeCode(..); a procedure for manipulation with the given code 40 vanishId(points); computes the vanishing ideal for the given set of points 41 sysFL(..); generates the Fitzgerald-Lax system 42 decodeRandomFL(..); manipulation with random codes via Fitzgerald-Lax 37 43 38 44 … … 76 82 77 83 /////////////////////////////////////////////////////////////////////////////// 78 // the po olynomial for Sala's restrictions84 // the polynomial for Sala's restrictions 79 85 static proc p_poly(int n, int a, int b) 80 86 { … … 89 95 /////////////////////////////////////////////////////////////////////////////// 90 96 91 proc 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 95 RETURN: the ring to work with the CRHT-ideal (with Sala's additions), the ideal itself is exported with the name 'crht' 96 EXAMPLE: example sysCRHT; shows an example 97 proc sysCRHT (int n, list defset, int e, int q, int m, list #) 98 "USAGE: sysCRHT(n,defset,e,q,m,[k]); n,e,q,m,k int, defset list of int's 99 @format 100 - n length of the cyclic code, 101 - defset is a list representing the defining set, 102 - e the error-correcting capacity, 103 - q field size 104 - m degree extension of the splitting field, 105 - if k>0 additional equations representing the fact that every two 106 error positions are either different or at least one of them is zero 107 @end format 108 RETURN: 109 @format 110 the ring to work with the CRHT-ideal (with Sala's additions), 111 containig an ideal with name 'crht' 112 @end format 113 EXAMPLE: example sysCRHT; shows an example 97 114 " 98 115 { … … 102 119 int i,j; 103 120 poly sum; 121 int k; 122 if ( size(#) > 0) 123 { 124 k = #[1]; 125 } 104 126 105 //------------ add check equations-----------------127 //---------------------- add check equations -------------------------- 106 128 for (i=1; i<=r; i++) 107 129 { … … 114 136 } 115 137 116 //------------ field equations on syndromes-----------138 //--------------------- field equations on syndromes ------------------ 117 139 for (i=1; i<=r; i++) 118 140 { … … 120 142 } 121 143 122 //------ ------ restrictions on error-locations: n-th roots of unity --------------144 //------ restrictions on error-locations: n-th roots of unity ---------- 123 145 for (i=1; i<=e; i++) 124 146 { … … 131 153 } 132 154 133 //------------ add Sala's additional conditions if necessary ----------------- 134 if (#) 155 //--------- add Sala's additional conditions if necessary -------------- 156 if ( k > 0 ) 157 135 158 { 136 159 for (i=1; i<=e; i++) … … 144 167 export crht; 145 168 return(@crht); 146 } example147 { 148 169 } 170 example 171 { "EXAMPLE:"; echo=2; 149 172 // binary cyclic [15,7,5] code with defining set (1,3) 150 151 list defset=1,3; // defining set 152 153 int n=15; // length154 int e=2; // error-correcting capacity155 int q=2; // basefield size156 int m=4; // degree extension of the splitting field157 int sala=1; // indicator to add additional equations173 intvec v = option(get); 174 175 list defset=1,3; // defining set 176 int n=15; // length 177 int e=2; // error-correcting capacity 178 int q=2; // basefield size 179 int m=4; // degree extension of the splitting field 180 int sala=1; // indicator to add additional equations 158 181 159 182 def A=sysCRHT(n,defset,e,q,m); 160 183 setring A; 161 A; // shows the ring we are working in162 print(crht); // the CRHT-ideal184 A; // shows the ring we are working in 185 print(crht); // the CRHT-ideal 163 186 option(redSB); 164 ideal red_crht=std(crht); 165 // reduced Groebner basis 187 ideal red_crht=std(crht); // reduced Groebner basis 166 188 print(red_crht); 167 189 … … 169 191 A=sysCRHT(n,defset,e,q,m,sala); 170 192 setring A; 171 print(crht); // theCRHT-ideal with additional equations from Sala193 print(crht); // CRHT-ideal with additional equations from Sala 172 194 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 183 proc 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 186 RETURN: the ring to work with the Sala's ideal for the minimum distance, the ideal itself is exported with the name 'crht_md' 187 EXAMPLE: example sysCRHTMindistBinary; shows an example 195 ideal red_crht=std(crht); // reduced Groebner basis 196 print(red_crht); 197 red_crht[5]; // general error-locator polynomial for this code 198 option(set,v); 199 } 200 201 /////////////////////////////////////////////////////////////////////////////// 202 203 204 proc sysCRHTMindist (int n, list defset, int w) 205 "USAGE: sysCRHTMindist(n,defset,w); n,w are int, defset is list of int's 206 @format 207 - n length of the cyclic code, 208 - defset is a list representing the defining set, 209 - w is a candidate for the minimum distance 210 @end format 211 RETURN: 212 @format 213 the ring to work with the Sala's ideal for the minimum distance 214 containing the ideal with name 'crht_md' 215 @end format 216 EXAMPLE: example sysCRHTMindist; shows an example 188 217 " 189 218 { … … 223 252 export crht_md; 224 253 return(@crht_md); 225 } example 254 } 255 example 226 256 { 227 257 "EXAMPLE:"; echo=2; 258 intvec v = option(get); 228 259 // binary cyclic [15,7,5] code with defining set (1,3) 229 260 230 list defset=1,3; // defining set 261 list defset=1,3; // defining set 262 int n=15; // length 263 int d=5; // candidate for the minimum distance 264 265 def A=sysCRHTMindist(n,defset,d); 266 setring A; 267 A; // shows the ring we are working in 268 print(crht_md); // the Sala's ideal for mindist 269 option(redSB); 270 ideal red_crht_md=std(crht_md); 271 print(red_crht_md); // reduced Groebner basis 231 272 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); 273 option(set,v); 243 274 } 244 275 … … 253 284 /////////////////////////////////////////////////////////////////////////////// 254 285 255 proc 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 259 RETURN: the ring to work with the generalized Newton identities (in triangular form if applicable), 260 the ideal itself is exported with the name 'newton' 286 proc sysNewton (int n, list defset, int t, int q, int m, list #) 287 "USAGE: sysNewton (n,defset,t,q,m,[tr]); n,t,q,m,tr int, defset list int's 288 @format 289 - n is length, 290 - defset is the defining set, 291 - t is the number of errors, 292 - q is basefield size, 293 - m is degree extension of the splitting field, 294 - if tr>0 it indicates that Newton identities in triangular 295 form should be constructed 296 @end format 297 RETURN: 298 @format 299 the ring to work with the generalized Newton identities (in triangular 300 form if applicable) containing the ideal with name 'newton' 301 @end format 261 302 EXAMPLE: example sysNewton; shows an example 262 303 " … … 265 306 int i,j; 266 307 int flag; 308 int tr; 309 310 if (size(#)>0) 311 { 312 tr=#[1]; 313 } 314 267 315 for (i=n; i>=1; i--) 268 316 { … … 295 343 296 344 //------------ generate generalized Newton identities ---------- 297 if ( #)345 if (tr) 298 346 { 299 347 for (i=1; i<=t; i++) … … 342 390 export newton; 343 391 return(@newton); 344 } example 392 } 393 example 345 394 { 346 395 "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 396 // Newton identities for a binary 3-error-correcting cyclic code of 397 //length 31 with defining set (1,5,7) 398 399 int n=31; // length 350 400 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 401 int t=3; // number of errors 402 int q=2; // basefield size 403 int m=5; // degree extension of the splitting field 404 int tr=1; // indicator of triangular form of Newton identities 405 355 406 356 407 def A=sysNewton(n,defset,t,q,m); 357 408 setring A; 358 A; // shows the ring we are working in359 print(newton); // generalized Newton identities409 A; // shows the ring we are working in 410 print(newton); // generalized Newton identities 360 411 361 412 //=============================== 362 A=sysNewton(n,defset,t,q,m,tr iangular);413 A=sysNewton(n,defset,t,q,m,tr); 363 414 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 415 print(newton); // generalized Newton identities in triangular form 416 } 417 418 /////////////////////////////////////////////////////////////////////////////// 419 // forms a list of special combinations needed for computation of Waring's 420 //function 369 421 static proc combinations_sum (int m, int n) 370 422 { … … 417 469 418 470 419 proc 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 422 RETURN: keeps the ring with the resulting system, which ideal is called 'bin' 471 proc sysBin (int v, list Q, int n, list #) 472 "USAGE: sysBin (v,Q,n,[odd]); v,n,odd are int, Q is list of int's 473 @format 474 - v a number if errors, 475 - Q is a generating set of the code, 476 - n the length, 477 - odd is an additional parameter: if 478 set to 1, then the generating set is enlarged by odd elements, 479 which are 2^(some power)*(some elment in the gen.set) mod n 480 @end format 481 THEORY: 482 @format 483 The system using Waring's function due to Augot et.al. 484 is built as in [*]. 485 @end format 486 RETURN: the ring with the resulting system called 'bin' 423 487 EXAMPLE: example sysBin; shows an example 424 488 " 425 489 { 490 int odd; 491 if (size(#)>0) 492 { 493 odd=#[1]; 494 } 495 426 496 //ring r=2,(sigma(1..v),S(1..n)),(lp(v),dp(n)); 427 497 ring r=2,(S(1..n),sigma(1..v)),lp; … … 433 503 int count1, count2, upper, coef_, flag, gener; 434 504 list Q_update; 435 if ( #==1)505 if (odd==1) 436 506 { 437 507 for (i=1; i<=n; i++) … … 460 530 } 461 531 462 //---- ---------- form polynomials for the Bin system via Waring's function --------------532 //---- form polynomials for the Bin system via Waring's function --------- 463 533 for (i=1; i<=size(Q_update); i++) 464 534 { … … 583 653 /////////////////////////////////////////////////////////////////////////////// 584 654 585 proc 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 655 proc sysQE(matrix check, matrix y, int t, list #) 656 "USAGE: sysQE(check,y,t,[fieldeq,formal]);check,y matrix;t,fieldeq,formal int 657 @format 658 - check is the check matrix of the code 659 - y is a received word, 660 - t the number of errors to be corrected, 661 - if fieldeq=1, then field equations are added, 662 - if formal=0, field equations on (known) syndrome variables 663 are not added, in order to add them (note that the exponent should 664 be as a number of elements in the INITIAL alphabet) one 590 665 needs to set formal>0 for the exponent 591 RETURN: the ring to work with together with the resulting system 666 @end format 667 THEORY: The system of quadratic equations is built as in [*]. 668 RETURN: the ring to work with together with the resulting system called 'qe' 592 669 EXAMPLE: example sysQE; shows an example 593 670 " 594 671 { 672 int fieldeq; 673 int formal; 674 if (size(#)>0) 675 { 676 fieldeq=#[1]; 677 } 678 if (size(#)>1) 679 { 680 formal=#[2]; 681 } 682 595 683 def br=basering; 596 684 list rl=ringlist(br); … … 632 720 matrix transf=inverse(transpose(h_full)); 633 721 634 //------ ----- expression matrix of check vectors w.r.t. the MDS basis --------------722 //------ expression matrix of check vectors w.r.t. the MDS basis ----------- 635 723 tim=rtimer; 636 724 for (i=1; i<=red ; i++) … … 706 794 } 707 795 708 //----- ---- form the quadratic equations according to the theory ---------------796 //----- form the quadratic equations according to the theory ----------- 709 797 for (i=1; i<=n; i++) 710 798 { … … 731 819 ideal qe=result; 732 820 export qe; 733 return(work); 734 //exportto(Top,h_full); 821 return(work); 735 822 } example 736 823 { 737 824 "EXAMPLE:"; echo = 2; 825 intvec v = option(get); 826 738 827 //correct 2 errors in [7,3] 8-ary code RS code 739 828 int t=2; int q=8; int n=7; int redun=4; … … 744 833 matrix x[1][3]=0,0,1,0; 745 834 matrix y[1][7]=encode(x,g); 835 746 836 //disturb with 2 errors 747 matrix rec[1][7]=errorInsert(y,list(2,4),list(1,a)); 837 matrix rec[1][7]=errorInsert(y,list(2,4),list(1,a)); 838 748 839 //generate the system 749 def A=sysQE(h,rec,t ,0,0);840 def A=sysQE(h,rec,t); 750 841 setring A; 751 842 print(qe); 843 752 844 //let us decode 753 845 option(redSB); 754 846 ideal sys_qe=std(qe); 755 847 print(sys_qe); 848 849 option(set,v); 756 850 } 757 851 … … 759 853 760 854 proc 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 855 "USAGE: errorInsert(y,pos,val); y is matrix, pos,val list of int's 856 @format 857 - y is a (code) word, 858 - pos = positions where errors occured, 859 - val = their corresponding values 860 @end format 762 861 RETURN: corresponding received word 763 EXAMPLE: example error ; shows an example862 EXAMPLE: example errorInsert; shows an example 764 863 " 765 864 { … … 793 892 794 893 proc 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)) 894 "USAGE: errorRand(y, num, e); y matrix, num,e int 895 @format 896 - y is a (code) word, 897 - num is the number of errors, 898 - e is an extension degree (if one wants values to be from GF(p^e)) 899 @end format 797 900 RETURN: corresponding received word 798 901 EXAMPLE: example errorRand; shows an example … … 854 957 /////////////////////////////////////////////////////////////////////////////// 855 958 856 proc 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)) 959 proc randomCheck(int m, int n, int e) 960 "USAGE: randomCheck(m, n, e); m,n,e are int 961 @format 962 - m x n are dimensions of the matrix, 963 - e is an extension degree (if one wants values to be from GF(p^e)) 964 @end format 859 965 RETURN: random check matrix 860 966 EXAMPLE: example randomCheck; shows an example … … 867 973 for (i=1; i<=m; i++) 868 974 { 869 temp=randomvector(n-m,e ,#);975 temp=randomvector(n-m,e); 870 976 for (j=1; j<=n-m; j++) 871 977 { … … 893 999 894 1000 proc 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. 896 NOTE: 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. 1001 "USAGE: genMDSMat(n, a); n is int, a is number 1002 @format 1003 - n x n are dimensions of the MDS matrix, 1004 - a is a primitive element of the field. 1005 @end format 1006 NOTE: 1007 @format 1008 An MDS matrix is constructed in the following way. We take a to be a 1009 generator of the multiplicative group of the field. Then we construct 1010 the Vandermonde matrix with this a. 1011 @end format 898 1012 ASSUME: extension field should already be defined 899 1013 RETURN: a matrix with the MDS property … … 926 1040 927 1041 proc mindist (matrix check) 928 "USAGE: mindist (check, q); check is a check matrix, q = field size 929 RETURN: minimum distance of the code together with the time needed for its calculation 1042 "USAGE: mindist (check, q); check matrix, q int 1043 @format 1044 - check is a check matrix, 1045 - q is the field size 1046 @end format 1047 RETURN: 1048 @format 1049 minimum distance of the code together with the time needed for its 1050 calculation 1051 @end format 930 1052 EXAMPLE: example mindist; shows an example 931 1053 " 932 1054 { 1055 intvec vopt = option(get); 1056 933 1057 int n=ncols(check); int redun=nrows(check); int t=redun+1; 934 1058 … … 947 1071 matrix z[1][n]; 948 1072 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 -------------------- 1073 def A=sysQE(check,z,count); 1074 1075 //proceed with solving the system w.r.t zero vector until some solutions 1076 //are found 952 1077 while (flag) 953 1078 { 954 A=sysQE(check,z,count ,0,0);1079 A=sysQE(check,z,count); 955 1080 setring A; 956 1081 ideal temp=qe; … … 979 1104 } 980 1105 list result=list(count,timsolve); 1106 1107 option(set,vopt); 981 1108 return(result); 982 } example 1109 } 1110 example 983 1111 { 984 1112 "EXAMPLE:"; echo = 2; … … 999 1127 1000 1128 proc 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 1003 ASSUME: Errors in rec should be correctable, otherwise the output is unpredictable 1129 "USAGE: decode(check, rec, t); check, rec matrix, t int 1130 @format 1131 - check is the check matrix of the code, 1132 - rec is a received word, 1133 - t is an upper bound for the number of errors one wants to correct 1134 @end format 1135 ASSUME: 1136 @format 1137 Errors in rec should be correctable, otherwise the output is 1138 unpredictable 1139 @end format 1004 1140 RETURN: a codeword that is closest to rec 1005 1141 EXAMPLE: example decode; shows an example 1006 1142 " 1007 { 1143 { 1144 intvec vopt = option(get); 1145 1008 1146 def br=basering; 1009 1147 int n=ncols(check); 1010 1148 1011 1149 int count=1; 1012 def A=sysQE(check,rec,count ,0,0);1150 def A=sysQE(check,rec,count); 1013 1151 while(1) 1014 1152 { 1015 A=sysQE(check,rec,count ,0,0);1153 A=sysQE(check,rec,count); 1016 1154 setring A; 1017 1155 matrix h_full=genMDSMat(n,a); … … 1043 1181 setring br; 1044 1182 matrix real_syn=imap(A,real_syn); 1183 1184 option(set,vopt); 1045 1185 return(rec-transpose(real_syn)); 1046 } example 1186 } 1187 example 1047 1188 { 1048 1189 "EXAMPLE:"; echo = 2; … … 1074 1215 1075 1216 1076 proc 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 1217 proc decodeRandom(int n, int redun, int ncodes, int ntrials, list #) 1218 "USAGE: decodeRandom(redun,q,ncodes,ntrials,[e]); all parameters int 1219 @format 1220 - redun is a redundabcy of a (random) code, 1221 - q is the field size, 1222 - ncodes is the number of random codes to be processed, 1223 - ntrials is the number of received vectors per code to be corrected 1224 - If e is given it sets the correction capacity explicitly. It 1225 should be used in case one expects some lower bound, 1226 otherwise the procedure tries to compute the real minimum distance 1227 to find out the error-correction capacity 1228 @end format 1082 1229 RETURN: nothing; 1083 EXAMPLE: example solveForRandom; shows an example1230 EXAMPLE: example decodeRandom; shows an example 1084 1231 " 1085 1232 { 1233 intvec vopt = option(get); 1234 1086 1235 int i,j; 1087 1236 matrix h; … … 1089 1238 ideal sys; 1090 1239 list tmp; 1240 int e; 1241 if (size(#)>0) 1242 { 1243 e=#[1]; 1244 } 1091 1245 1092 1246 option(redSB); … … 1102 1256 "check matrix:"; 1103 1257 print(h); 1104 if ( #>0)1105 { 1106 t= #;1258 if (e>0) 1259 { 1260 t=e; 1107 1261 } else { 1108 1262 tim=rtimer; … … 1118 1272 1119 1273 //------------- generate the template system ---------------------- 1120 def A=sysQE(h,z,t ,0,0);1274 def A=sysQE(h,z,t); 1121 1275 setring A; 1122 1276 matrix word,y,rec; … … 1128 1282 timdec3=timdec3+rtimer-tim3; 1129 1283 1130 //------ ------- modify the template according to every received word -------------------1284 //------ modify the template according to every received word -------------- 1131 1285 for (j=1; j<=ntrials; j++) 1132 1286 { … … 1142 1296 timdec2=timdec2+rtimer-tim2; 1143 1297 kill A; 1298 option(set,vopt); 1144 1299 } 1145 1300 printf("Time for mindist: %p", timdist); … … 1148 1303 printf("Time for GB in decoding: %p", timdec); 1149 1304 printf("Time for sysQE in decoding: %p", timdec3); 1150 } example 1305 } 1306 example 1151 1307 { 1152 1308 "EXAMPLE:"; echo = 2; … … 1155 1311 1156 1312 // correct 2 errors in 5 random binary codes, 50 trials each 1157 solveForRandom(n,redun,5,50,2); 1158 } 1159 1160 /////////////////////////////////////////////////////////////////////////////// 1161 1162 1163 proc 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 1313 decodeRandom(n,redun,5,50,2); 1314 } 1315 1316 /////////////////////////////////////////////////////////////////////////////// 1317 1318 1319 proc decodeCode(matrix check, int ntrials, list #) 1320 "USAGE: decodeCode(check, ntrials, [e]); check matrix, ntrials,e int 1321 @format 1322 - check is a check matrix for the code, 1323 - ntrials is the number of received vectors per code to be 1324 corrected. 1325 - If e is given it sets the correction capacity explicitly. It 1326 should be used in case one expects some lower bound, 1327 otherwise the procedure tries to compute the real minimum distance 1328 to find out the error-correction capacity 1329 @end format 1168 1330 RETURN: nothing; 1169 EXAMPLE: example solveForCode; shows an example1331 EXAMPLE: example decodeCode; shows an example 1170 1332 " 1171 1333 { 1334 intvec vopt = option(get); 1335 1172 1336 int n=ncols(check); 1173 1337 int redun=nrows(check); … … 1177 1341 ideal sys; 1178 1342 list tmp; 1343 int e; 1344 if (size(#)>0) 1345 { 1346 e=#[1]; 1347 } 1179 1348 1180 1349 option(redSB); … … 1188 1357 1189 1358 //------------------ determine error-correction capacity ------------------- 1190 if ( #>0)1191 { 1192 t= #;1359 if (e>0) 1360 { 1361 t=e; 1193 1362 } else { 1194 1363 tim=rtimer; … … 1204 1373 1205 1374 //------------- generate the template system ---------------------- 1206 def A=sysQE(h,z,t ,0,0);1375 def A=sysQE(h,z,t); 1207 1376 setring A; 1208 1377 matrix word,y,rec; … … 1214 1383 timdec3=timdec3+rtimer-tim3; 1215 1384 1216 //--- ---------- modify the template according to every received word -------------------1385 //--- modify the template according to every received word --------------- 1217 1386 for (j=1; j<=ntrials; j++) 1218 1387 { … … 1233 1402 printf("Time for GB in decoding: %p", timdec); 1234 1403 printf("Time for sysQE in decoding: %p", timdec3); 1235 } example 1404 1405 option(set,vopt); 1406 } 1407 example 1236 1408 { 1237 1409 "EXAMPLE:"; echo = 2; … … 1241 1413 1242 1414 // correct 2 errors in using the code above, 50 trials 1243 solveForCode(check,50,2);1415 decodeCode(check,50,2); 1244 1416 } 1245 1417 … … 1273 1445 } 1274 1446 1275 ////////////////////////////////////////////////////////////////////////////////////// 1276 // return index of an element in the ideal where it does not vanish at the given point 1447 /////////////////////////////////////////////////////////////////////////////// 1448 // return index of an element in the ideal where it does not vanish at the 1449 //given point 1277 1450 static proc find_index (ideal G, matrix p) 1278 1451 { … … 1312 1485 } 1313 1486 1314 //////////////////////////////////////////////////////////////////////////////////// 1315 // checl whether given polynomial is divisible by some leading monomial of the ideal 1487 /////////////////////////////////////////////////////////////////////////////// 1488 // check whether given polynomial is divisible by some leading monomial of the 1489 //ideal 1316 1490 static proc divisible (poly m, ideal G) 1317 1491 { … … 1326 1500 1327 1501 proc vanishId (list points) 1328 "USAGE: vanishId (points,e); points is a list of points, for which the vanishing ideal is to be constructed. 1502 "USAGE: vanishId (points); point is a list of matrices 1503 @format 1504 - points is a list of points for which the vanishing ideal is to be 1505 constructed 1506 @end format 1329 1507 RETURN: Vanishing ideal corresponding to the given set of points 1330 1508 EXAMPLE: example vanishId; shows an example … … 1386 1564 } 1387 1565 1388 ////////////////////////////////////////////////////////////////////////////////////////////////// 1389 // construct the list of all vectors of length m with elements in p^e, where p is characteristics 1566 /////////////////////////////////////////////////////////////////////////////// 1567 // construct the list of all vectors of length m with elements in p^e, where p 1568 //is characteristics 1390 1569 proc pointsGen (int m, int e) 1391 1570 { … … 1548 1727 1549 1728 proc 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 1729 "USAGE: sysFL (check,y,t,e,s); check,y matrix, t,e,s int 1730 @format 1731 - check is a check matrix of the code, 1732 - y is a received word, 1733 - t the number of errors to correct, 1734 - e is the extension degree, 1735 - s is the dimension of the point for the vanishing ideal 1736 @end format 1552 1737 RETURN: the system of Fitzgerald-Lax for the given decoding problem 1553 1738 EXAMPLE: example sysFL; shows an example … … 1620 1805 export points; 1621 1806 return(result); 1622 } example 1807 } 1808 example 1623 1809 { 1624 1810 "EXAMPLE:"; echo = 2; 1811 intvec vopt = option(get); 1625 1812 1626 1813 list l=FLpreprocess(3,1,11,2,""); … … 1646 1833 option(redSB); 1647 1834 ideal red_sys=std(sys); 1648 red_sys; // read the solutions from this redGB 1835 red_sys; 1836 // read the solutions from this redGB 1649 1837 // the points are (0,0,1) and (0,1,0) with error values 1 and -1 resp. 1650 1838 // use list points to find error positions; 1651 points; 1839 points; 1840 1841 option(set,vopt); 1652 1842 } 1653 1843 … … 1691 1881 rl[1][2][1]=string("a"); 1692 1882 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)+")"); 1883 rl[1][3][1]=tmp; 1695 1884 rl[1][3][1][1]=string("lp"); 1696 1885 rl[1][3][1][2]=1; … … 1704 1893 rl[3]=tmp; 1705 1894 rl[3][1]=tmp; 1706 //rl[3][1][1]=string("dp("+string((t-1)*(s+1)+s)+"),lp("+string(s+1)+")");1707 1895 rl[3][1][1]=string("lp"); 1708 1896 intvec v=1; … … 1743 1931 /////////////////////////////////////////////////////////////////////////////// 1744 1932 // random vector of length n with entries from p^e, p the characteristic 1745 static proc randomvector(int n, int e , int #)1933 static proc randomvector(int n, int e) 1746 1934 { 1747 1935 int i; … … 1749 1937 for (i=1; i<=n; i++) 1750 1938 { 1751 result[i,1]=asElement(random_prime_vector(e ,#));1939 result[i,1]=asElement(random_prime_vector(e)); 1752 1940 } 1753 1941 return(result); 1754 1942 } 1755 1943 1756 //////////////////////////////////////////////////////////////////////////////////////////// 1757 // "convert" representation of an element from the field extension from vector to an elelemnt 1944 /////////////////////////////////////////////////////////////////////////////// 1945 // "convert" representation of an element from the field extension from vector 1946 //to an elelemnt 1758 1947 static proc asElement(list l) 1759 1948 { … … 1771 1960 /////////////////////////////////////////////////////////////////////////////// 1772 1961 // random vector of length n with entries from p, p the characteristic 1773 static 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 } 1962 static proc random_prime_vector (int n) 1963 { 1964 list rl=ringlist(basering); 1965 int i, charac; 1966 for (i=2; i<=rl[1][1]; i++) 1967 { 1968 if (rl[1][1] mod i ==0) 1969 { 1970 break; 1971 } 1972 } 1973 charac=i; 1974 1782 1975 list l; 1783 int i;1976 1784 1977 for (i=1; i<=n; i++) 1785 1978 { … … 1791 1984 /////////////////////////////////////////////////////////////////////////////// 1792 1985 1793 proc 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 1986 proc decodeRandomFL(int n, int redun, int p, int e, int t, int ncodes, int ntrials, string minpol) 1987 "USAGE: decodeRandomFL(redun,p,e,n,t,ncodes,ntrials,minpol); 1988 @format 1989 - n is length of codes generated, redun = redundancy of codes 1990 generated, 1991 - p is characteristics, 1992 - e is the extension degree, 1993 - t is the number of errors to correct, 1994 - ncodes is the number of random codes to be processed, 1995 - ntrials is the number of received vectors per code to be corrected, 1996 - minpol: due to some pecularities of SINGULAR one needs to provide 1997 minimal polynomial for the extension explicitly 1998 @end format 1799 1999 RETURN: nothing 1800 EXAMPLE: example solveForRandomFL; shows an example 1801 { 2000 EXAMPLE: example decodeRandomFL; shows an example 2001 " 2002 { 2003 intvec vopt = option(get); 2004 1802 2005 list l=FLpreprocess(p,e,n,t,minpol); 1803 2006 … … 1818 2021 for (i=1; i<=ncodes; i++) 1819 2022 { 1820 h=randomCheck(redun,n,e ,1);2023 h=randomCheck(redun,n,e); 1821 2024 g=dual_code(h); 1822 2025 tim2=rtimer; 1823 2026 tim3=rtimer; 1824 2027 1825 //------------ generate the template system ------------------2028 //---------------- generate the template system ----------------------- 1826 2029 sys=sysFL(h,z,t,e,s_work); 1827 2030 timdec3=timdec3+rtimer-tim3; 1828 2031 1829 //------ ------- modifying the template according to the received word ---------------2032 //------ modifying the template according to the received word --------- 1830 2033 for (j=1; j<=ntrials; j++) 1831 2034 { … … 1844 2047 printf("Time for GB in decoding: %p", timdec); 1845 2048 printf("Time for generating Fitzgerald-Lax system during decoding: %p", timdec3); 1846 } example 2049 2050 option(set,vopt); 2051 } 2052 example 1847 2053 { 1848 2054 "EXAMPLE:"; echo = 2; 1849 2055 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,""); 2056 // correcting one error for one random binary code of length 25, 2057 //redundancy 14; 300 words are processed 2058 decodeRandomFL(25,14,2,1,1,1,300,""); 1852 2059 } 1853 2060 … … 1874 2081 int q=128; int n=120; int redun=n-30; 1875 2082 ring r=(q,a),x,dp; 1876 solveForRandom(n,redun,1,1,6);2083 decodeRandom(n,redun,1,1,6); 1877 2084 1878 2085 int q=128; int n=120; int redun=n-20; 1879 2086 ring r=(q,a),x,dp; 1880 solveForRandom(n,redun,1,1,9);2087 decodeRandom(n,redun,1,1,9); 1881 2088 1882 2089 int q=128; int n=120; int redun=n-10; 1883 2090 ring r=(q,a),x,dp; 1884 solveForRandom(n,redun,1,1,19);2091 decodeRandom(n,redun,1,1,19); 1885 2092 1886 2093 int q=256; int n=150; int redun=n-10; 1887 2094 ring r=(q,a),x,dp; 1888 solveForRandom(n,redun,1,1,22);2095 decodeRandom(n,redun,1,1,22); 1889 2096 1890 2097 ////////////// SOME HARD EXAMPLES ////////////////////// … … 1896 2103 int q=128; int n=120; int redun=n-40; 1897 2104 ring r=(q,a),x,dp; 1898 solveForRandom(n,redun,1,1,6);2105 decodeRandom(n,redun,1,1,6); 1899 2106 1900 2107 redun=n-30; 1901 solveForRandom(n,redun,1,1,8);2108 decodeRandom(n,redun,1,1,8); 1902 2109 1903 2110 redun=n-20; 1904 solveForRandom(n,redun,1,1,12);2111 decodeRandom(n,redun,1,1,12); 1905 2112 1906 2113 redun=n-10; 1907 solveForRandom(n,redun,1,1,24);2114 decodeRandom(n,redun,1,1,24); 1908 2115 1909 2116 int q=256; int n=150; int redun=n-10; 1910 2117 ring r=(q,a),x,dp; 1911 solveForRandom(n,redun,1,1,26);2118 decodeRandom(n,redun,1,1,26); 1912 2119 1913 2120
Note: See TracChangeset
for help on using the changeset viewer.