Changeset c34952 in git
- Timestamp:
- Aug 11, 2008, 10:31:56 AM (16 years ago)
- Branches:
- (u'spielwiese', '17f1d200f27c5bd38f5dfc6e8a0879242279d1d8')
- Children:
- 1e5b760d17d4a47feae5716b7f581fd22e98960b
- Parents:
- 8e5a056876d28989ee075ac8260ef617777a486b
- File:
-
- 1 edited
-
Singular/LIB/decodedistGB.lib (modified) (96 diffs)
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/decodedistGB.lib
r8e5a056 rc34952 1 1 /////////////////////////////////////////////////////////////////////////////// 2 version="$Id: decodedistGB.lib,v 1. 1 2008-07-15 09:39:28 SingularExp $";2 version="$Id: decodedistGB.lib,v 1.2 2008-08-11 08:31:56 bulygin Exp $"; 3 3 category="Coding theory"; 4 4 info=" 5 5 LIBRARY: decodedistGB.lib Generating and solving systems of polynomial equations for decoding and finding the minimum distance of linear codes 6 6 AUTHORS: Stanislav Bulygin, bulygin@mathematik.uni-kl.de, 7 Ruud Pellikaan, g.r.pellikaan@tue.nl 7 Ruud Pellikaan, g.r.pellikaan@tue.nl 8 8 9 9 OVERVIEW: 10 10 In this library we generate several systems used for decoding cyclic codes. And finding the minimum distance 11 11 Namely, we work with the Cooper's philosophy and generalized Newton identities. 12 The original method of quadratic equations is worked out here as well. 12 The original method of quadratic equations is worked out here as well. 13 13 We also (for comparison) enable to work with the system of Fitzgerald-Lax. 14 14 We provide also some auxiliary functions for further manipulations and decoding. 15 For an overview of the methods mentioned above, cf. Stanislav Bulygin, Ruud Pellikaan: 'Decoding and finding the minimum distance with 15 For an overview of the methods mentioned above, cf. Stanislav Bulygin, Ruud Pellikaan: 'Decoding and finding the minimum distance with 16 16 Groebner bases: history and new insights', in 'Selected Topics in Information and Coding Theory', World Scientific (2008) (to appear). 17 17 18 MAIN PROCEDURES: 18 MAIN PROCEDURES: 19 19 CRHT(n,defset,e,q,m,#); generates the CRHT-ideal that follows Cooper's philosophy, Sala's extentions are available 20 20 CRHT_mindist_binary(n,defset,w); generates the ideal from Cooper's philosophy to find the minimum distance in the binary case … … 29 29 MDSmat(n,a); generates an MDS (actually an RS) matrix 30 30 mindist(check); computes the minimum distance of the code via solving systems of quadratic equations 31 decode(rec); decoding of a word using the systems of quadratic equations 31 decode(rec); decoding of a word using the systems of quadratic equations 32 32 solve_for_random(redun,ncodes,ntrials,#); a procedure for manipulation with random codes 33 33 solve_for_code(check,ntrials,#); a procedure for manipulation with a given codes 34 vanish_id(points); computes the vanishing ideal for the given set of points 34 vanish_id(points); computes the vanishing ideal for the given set of points. The algorithm of Farr and Gao is implemented 35 35 FLsystem(check,y,t,e,s); generates the Fitzgerald-Lax system 36 36 FL_solve_for_random(n,redun,p,e,t,ncodes,ntrials,minpol); a procedure for manipulation with random codes via Fitzgerald-Lax 37 38 39 KEYWORDS: Cyclic code; Linear code; Decoding; 37 38 39 KEYWORDS: Cyclic code; Linear code; Decoding; 40 40 Minimum distance; Groebner bases 41 41 "; … … 89 89 } 90 90 } 91 return(result); 91 return(result); 92 92 } 93 93 … … 110 110 " 111 111 { 112 int r=size(defset); 112 int r=size(defset); 113 113 ring @crht=(q,a),(Y(e..1),Z(1..e),X(r..1)),lp; 114 114 ideal crht; 115 115 int i,j; 116 116 poly sum; 117 117 118 118 // check equations 119 119 for (i=1; i<=r; i++) … … 126 126 crht[i]=sum-X(i); 127 127 } 128 128 129 129 // restrictions on syndromes 130 130 for (i=1; i<=r; i++) … … 132 132 crht=crht,X(i)^(q^m)-X(i); 133 133 } 134 134 135 135 // n-th roots of unity 136 136 for (i=1; i<=e; i++) … … 138 138 crht=crht,Z(i)^(n+1)-Z(i); 139 139 } 140 140 141 141 for (i=1; i<=e; i++) 142 142 { 143 143 crht=crht,Y(i)^(q-1)-1; 144 } 145 144 } 145 146 146 if (#) 147 147 { … … 153 153 } 154 154 } 155 } 156 export crht; 157 return(@crht); 158 } 159 example 155 } 156 export crht; 157 return(@crht); 158 } example 160 159 { 161 160 "EXAMPLE:"; echo=2; 162 161 // binary cyclic [15,7,5] code with defining set (1,3) 163 162 164 163 list defset=1,3; // defining set 165 164 166 165 int n=15; // length 167 166 int e=2; // error-correcting capacity … … 169 168 int m=4; // degree extension of the splitting field 170 169 int sala=1; // indicator to add additional equations 171 172 def A=CRHT(n,defset,e,q,m); 170 171 def A=CRHT(n,defset,e,q,m); 173 172 setring A; 174 173 A; // shows the ring we are working in 175 print(crht); // the CRHT-ideal 174 print(crht); // the CRHT-ideal 176 175 option(redSB); 177 176 ideal red_crht=std(crht); 178 177 // reduced Groebner basis 179 178 print(red_crht); 180 179 181 180 //============================ 182 A=CRHT(n,defset,e,q,m,sala); 183 setring A; 181 A=CRHT(n,defset,e,q,m,sala); 182 setring A; 184 183 print(crht); // the CRHT-ideal with additional equations from Sala 185 184 option(redSB); … … 199 198 " 200 199 { 201 int r=size(defset); 200 int r=size(defset); 202 201 ring @crht_md=2,Z(1..w),lp; 203 202 ideal crht_md; 204 203 int i,j; 205 204 poly sum; 206 205 207 206 // check equations 208 207 for (i=1; i<=r; i++) … … 214 213 } 215 214 crht_md[i]=sum; 216 } 217 218 215 } 216 217 219 218 // n-th roots of unity 220 219 for (i=1; i<=w; i++) 221 220 { 222 221 crht_md=crht_md,Z(i)^n-1; 223 } 224 225 222 } 223 224 226 225 for (i=1; i<=w; i++) 227 226 { … … 231 230 } 232 231 } 233 234 export crht_md; 235 return(@crht_md); 236 } 237 example 232 233 export crht_md; 234 return(@crht_md); 235 } example 238 236 { 239 237 "EXAMPLE:"; echo=2; 240 238 // binary cyclic [15,7,5] code with defining set (1,3) 241 239 242 240 list defset=1,3; // defining set 243 241 244 242 int n=15; // length 245 int d=5; // candidate for the minimum distance 246 247 def A=CRHT_mindist_binary(n,defset,d); 243 int d=5; // candidate for the minimum distance 244 245 def A=CRHT_mindist_binary(n,defset,d); 248 246 setring A; 249 247 A; // shows the ring we are working in … … 252 250 ideal red_crht_md=std(crht_md); 253 251 // reduced Groebner basis 254 print(red_crht_md); 252 print(red_crht_md); 255 253 } 256 254 … … 265 263 t is the number of errors, q is basefield size, m is degree extension of the splitting field 266 264 if triangular>0 it indicates that Newton identities in triangular form should be constructed 267 RETURN: a ring to work with the generalized Newton identities (in triangular form if applicable), 265 RETURN: a ring to work with the generalized Newton identities (in triangular form if applicable), 268 266 the ideal itself is exported with the name 'newton' 269 267 EXAMPLE: example Newton; shows an example 270 268 " 271 { 269 { 272 270 string s="ring @newton=("+string(q)+",a),("; 273 271 int i,j; … … 281 279 { 282 280 flag=0; 283 break; 281 break; 284 282 } 285 283 } … … 294 292 s=s+"S("+string(defset[i])+"),"; 295 293 } 296 s=s+"S("+string(defset[1])+")),lp;"; 297 294 s=s+"S("+string(defset[1])+")),lp;"; 295 298 296 execute(s); 299 300 ideal newton; 297 298 ideal newton; 301 299 poly sum; 302 303 300 301 304 302 // generate generalized Newton identities 305 303 if (#) … … 315 313 } 316 314 } else 317 { 315 { 318 316 for (i=1; i<=t; i++) 319 317 { … … 334 332 } 335 333 newton=newton,S(t+i)+sum; 336 } 337 334 } 335 338 336 // field equations on sigma's 339 337 for (i=1; i<=t; i++) … … 341 339 newton=newton,sigma(i)^(q^m)-sigma(i); 342 340 } 343 341 344 342 // conjugacy relations 345 343 for (i=1; i<=n; i++) … … 349 347 newton=simplify(newton,2); 350 348 export newton; 351 return(@newton); 352 } 353 example 349 return(@newton); 350 } example 354 351 { 355 352 "EXAMPLE:"; echo = 2; 356 353 // Newton identities for a binary 3-error-correcting cyclic code of length 31 with defining set (1,5,7) 357 354 358 355 int n=31; // length 359 356 list defset=1,5,7; //defining set … … 361 358 int q=2; // basefield size 362 359 int m=5; // degree extension of the splitting field 363 int triangular=1; // indicator of triangular form of Newton identities 364 360 int triangular=1; // indicator of triangular form of Newton identities 361 365 362 def A=Newton(n,defset,t,q,m); 366 363 setring A; 367 364 A; // shows the ring we are working in 368 365 print(newton); // generalized Newton identities 369 366 370 367 //=============================== 371 368 A=Newton(n,defset,t,q,m,triangular); 372 setring A; 369 setring A; 373 370 print(newton); // generalized Newton identities in triangular form 374 371 } … … 395 392 count++; 396 393 } 397 if (flag) 394 if (flag) 398 395 { 399 396 for (j=2; j<=m; j++) … … 401 398 interm[i][j]=interm[i][j] div j; 402 399 } 403 result[size(result)+1]=interm[i]; 400 result[size(result)+1]=interm[i]; 404 401 } 405 402 } … … 421 418 422 419 static proc Bin (int v, list Q, int n, int#) 423 "USAGE: Bin (v, Q, n, #); v a number if errors, Q is a generating set of the code, n the length, # is additional parameter is 420 "USAGE: Bin (v, Q, n, #); v a number if errors, Q is a generating set of the code, n the length, # is additional parameter is 424 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 425 422 RETURN: keeps a ring with the resulting system, which ideal is called 'bin' … … 462 459 Q_update=Q; 463 460 } 464 461 465 462 for (i=1; i<=size(Q_update); i++) 466 463 { … … 475 472 } 476 473 count1=0; 477 for (j=2; j<=upper-1; j++) 474 for (j=2; j<=upper-1; j++) 478 475 { 479 476 count1=count1+exp_count(j,2); … … 498 495 sum_=sum_+coef_*mon; 499 496 } 500 result=result,S(Q_update[i])-sum_; 497 result=result,S(Q_update[i])-sum_; 501 498 } 502 499 ideal bin=simplify(result,2); 503 500 export bin; 504 501 return(r); 505 } 506 example 502 } example 507 503 { 508 504 "EXAMPLE:"; echo = 2; … … 524 520 if (ncols(x)!=nrows(g)) {print("ERRORenc2!");} 525 521 return(x*g); 526 } 527 example 522 } example 528 523 { 529 524 "EXAMPLE:"; echo = 2; … … 535 530 0,0,0,1,1,1,0; 536 531 //encode x with the generator matrix g 537 print(enc(x,g)); 532 print(enc(x,g)); 538 533 } 539 534 … … 546 541 if (nrows(c)>1) {print("ERRORsyndrome1!");} 547 542 if (ncols(c)!=ncols(h)) {print("ERRORsyndrome2!");} 548 return(h*transpose(c)); 549 } 550 example 543 return(h*transpose(c)); 544 } example 551 545 { 552 546 "EXAMPLE:"; echo = 2; … … 566 560 print(syndrome(check,c)); 567 561 c[1,3]=1; 568 //now c is a codeword 562 //now c is a codeword 569 563 print(syndrome(check,c)); 570 564 } … … 581 575 582 576 proc QE(matrix check, matrix y, int t, int fieldeq, int formal) 583 "USAGE: QE(check, y, t, fieldeq, formal); check is the check matrix of the code 584 y is a received word, t the number of errors to be corrected, 585 if fieldeq=1, then field equations are added; if formal=0, fields equations on (known) syndrome variables 577 "USAGE: QE(check, y, t, fieldeq, formal); check is the check matrix of the code 578 y is a received word, t the number of errors to be corrected, 579 if fieldeq=1, then field equations are added; if formal=0, fields equations on (known) syndrome variables 586 580 are not added, in order to add them (note that the exponent should be as a number of elements in the INITIAL alphabet) one 587 581 needs to set formal>0 for the exponent … … 590 584 " 591 585 { 592 def br=basering; 586 def br=basering; 593 587 list rl=ringlist(br); 594 588 595 589 int red=nrows(check); 596 590 int n=ncols(check); 597 591 int q=rl[1][1]; 598 592 599 593 if (formal==0) 600 { 594 { 601 595 ring work=(q,a),(V(1..t),U(1..n)),dp; 602 596 } else … … 604 598 ring work=(q,a),(V(1..t),U(1..n),s(1..red)),(dp(t),lp(n),dp(red)); 605 599 } 606 600 607 601 matrix check=imap(br,check); 608 602 matrix y=imap(br,y); 609 603 610 604 matrix h_full=MDSmat(n,a); 611 605 matrix h=submat(h_full,1..red,1..n); 612 606 if (nrows(y)!=1) {print("ERROR1Pell");} 613 607 if (ncols(y)!=n) {print("ERROR2Pell");} 614 608 615 609 ideal result; 616 610 617 611 list c; 618 612 list a; … … 625 619 c[i]=tmp; 626 620 } 627 621 628 622 int tim=rtimer; 629 623 matrix transf=inverse(transpose(h_full)); 630 624 631 625 tim=rtimer; 632 626 for (i=1; i<=red ; i++) … … 635 629 a[i]=transf*a[i]; 636 630 } 637 631 638 632 tim=rtimer; 639 matrix te[n][1]; 633 matrix te[n][1]; 640 634 for (i=1; i<=n; i++) 641 635 { … … 643 637 { 644 638 if ((j<i)&&(i<=t+1)) {c[i][j]=c[j][i];} 645 else 639 else 646 640 { 647 if (i+j<=n+1) 641 if (i+j<=n+1) 648 642 { 649 643 c[i][j]=te; 650 c[i][j][i+j-1,1]=1; 644 c[i][j][i+j-1,1]=1; 651 645 } 652 646 else … … 654 648 c[i][j]=star(h_full,i,j); 655 649 c[i][j]=transf*c[i][j]; 656 } 650 } 657 651 } 658 652 } 659 653 } 660 661 662 tim=rtimer; 654 655 656 tim=rtimer; 663 657 if (formal==0) 664 658 { … … 699 693 result=result,V(j)^q-V(j); 700 694 } 701 } 695 } 702 696 for (i=1; i<=n; i++) 703 697 { … … 718 712 } 719 713 result=result,sum1-sum3; 720 } 721 714 } 715 722 716 result=simplify(result,2); 723 717 724 718 ideal qe=result; 725 719 export qe; 726 720 return(work); 727 //exportto(Top,h_full); 728 } 729 example 721 //exportto(Top,h_full); 722 } example 730 723 { 731 724 "EXAMPLE:"; echo = 2; … … 747 740 option(redSB); 748 741 ideal sys_qe=std(qe); 749 print(sys_qe); 742 print(sys_qe); 750 743 } 751 744 … … 763 756 } 764 757 return(result); 765 } 766 example 758 } example 767 759 { 768 760 "EXAMPLE:"; echo = 2; … … 776 768 matrix y[1][7]=enc(x,g); 777 769 print(y); 778 770 779 771 //disturb with 2 errors 780 772 matrix rec[1][7]=error(y,list(2,4),list(1,a)); … … 806 798 } 807 799 if (flag) {pos[i]=temp;break;} 808 } 809 } 810 800 } 801 } 802 811 803 for (i=1; i<=num; i++) 812 804 { 813 805 flag=1; 814 806 while(flag) 815 { 816 tempnum=randomvector(1,e); 817 if (tempnum!=0) {flag=0;} 818 } 819 val[i]=tempnum; 820 } 821 807 { 808 tempnum=randomvector(1,e); 809 if (tempnum!=0) {flag=0;} 810 } 811 val[i]=tempnum; 812 } 813 822 814 for (i=1; i<=size(pos); i++) 823 815 { … … 825 817 } 826 818 return(result); 827 } 828 example 819 } example 829 820 { 830 821 "EXAMPLE:"; echo = 2; … … 837 828 matrix x[1][3]=0,0,1,0; 838 829 matrix y[1][7]=enc(x,g); 839 830 840 831 //disturb with 2 random errors 841 832 matrix rec[1][7]=error_rand(y,2,3); 842 833 print(rec); 843 print(rec-y); 834 print(rec-y); 844 835 } 845 836 … … 857 848 for (i=1; i<=m; i++) 858 849 { 859 temp=randomvector(n-m,e,#); 850 temp=randomvector(n-m,e,#); 860 851 for (j=1; j<=n-m; j++) 861 852 { 862 853 rand[i,j]=temp[j,1]; 863 } 854 } 864 855 } 865 856 result=concat(rand,unitmat(m)); 866 857 return(result); 867 } 868 example 869 { 870 "EXAMPLE:"; echo = 2; 858 } example 859 { 860 "EXAMPLE:"; echo = 2; 871 861 int redun=5; int n=15; 872 862 ring r=2,x,dp; 873 863 874 864 //generate random check matrix for a [15,5] binary code 875 865 matrix h=random_check(redun,n,1); 876 866 print(h); 877 867 878 868 //corresponding generator matrix 879 869 matrix g=dual_code(h); 880 print(g); 870 print(g); 881 871 } 882 872 … … 885 875 An MDS matrix is constructed in the following way. We take a to be a generator of the multiplicative group of the field. 886 876 Then we construct the Vandermonde matrix with this a. 887 ASSUME: extension field should already be defined 877 ASSUME: extension field should already be defined 888 878 RETURN: a matrix with the MDS property 889 EXAMPLE: example random_check; shows an example 879 EXAMPLE: example random_check; shows an example 890 880 " 891 881 { … … 897 887 { 898 888 result[j+1,i+1]=(a^i)^j; 899 } 889 } 900 890 } 901 891 return(result); 902 } 903 example 892 } example 904 893 { 905 894 "EXAMPLE:"; echo = 2; 906 895 int q=16; int n=15; 907 896 ring r=(q,a),x,dp; 908 897 909 898 //generate an MDS (Vandermonde) matrix 910 899 matrix h_full=MDSmat(n,a); 911 print(h_full); 900 print(h_full); 912 901 } 913 902 … … 916 905 "USAGE: mindist (check, q); check is a check matrix, q = field size 917 906 RETURN: minimum distance of the code together with the time needed for its calculation 918 EXAMPLE: example mindist; shows an example 907 EXAMPLE: example mindist; shows an example 919 908 " 920 909 { 921 910 int n=ncols(check); int redun=nrows(check); int t=redun+1; 922 923 def br=basering; 924 list rl=ringlist(br); 911 912 def br=basering; 913 list rl=ringlist(br); 925 914 int q=rl[1][1]; 926 927 ring work=(q,a),(V(1..t),U(1..n)),dp; 928 matrix check=imap(br,check); 929 915 916 ring work=(q,a),(V(1..t),U(1..n)),dp; 917 matrix check=imap(br,check); 918 930 919 ideal temp; 931 920 int count=1; … … 933 922 int flag2; 934 923 int i, tim, timsolve; 935 matrix z[1][n]; 924 matrix z[1][n]; 936 925 option(redSB); 937 926 def A=QE(check,z,count,0,0); … … 943 932 tim=rtimer; 944 933 temp=std(temp); 945 timsolve=timsolve+rtimer-tim; 934 timsolve=timsolve+rtimer-tim; 946 935 flag2=1; 947 936 setring work; 948 temp=imap(A,temp); 937 temp=imap(A,temp); 949 938 for (i=1; i<=n; i++) 950 939 { 951 if 952 (temp[i]!=U(n-i+1)) 940 if 941 (temp[i]!=U(n-i+1)) 953 942 { 954 943 flag2=0; 955 944 } 956 945 } 957 if (!flag2) 946 if (!flag2) 958 947 { 959 948 flag=0; 960 949 } 961 else 950 else 962 951 { 963 952 count++; 964 953 } 965 954 } 966 list result=list(count,timsolve); 967 return(result); 968 } 969 example 955 list result=list(count,timsolve); 956 return(result); 957 } example 970 958 { 971 959 "EXAMPLE:"; echo = 2; 972 //determine a minimum distance for a [7,3] binary code 973 int q=8; int n=7; int redun=4; int t=redun+1; 974 ring r=(q,a),x,dp; 975 960 //determine a minimum distance for a [7,3] binary code 961 int q=8; int n=7; int redun=4; int t=redun+1; 962 ring r=(q,a),x,dp; 963 976 964 //generate random check matrix 977 965 matrix h=random_check(redun,n,1); 978 966 print(h); 979 list l=mindist(h); 980 print(l[1]); 981 //time for the comutation in secs 982 print(l[2]); 967 list l=mindist(h); 968 print(l[1]); 969 //time for the comutation in secs 970 print(l[2]); 983 971 } 984 972 … … 988 976 ASSUME: Errors in rec should be correctable, otherwise the output is unpredictable 989 977 RETURN: a codeword that is closest to rec 990 EXAMPLE: example decode; shows an example 978 EXAMPLE: example decode; shows an example 991 979 " 992 { 980 { 993 981 def br=basering; 994 int n=ncols(check); 995 996 int count=1; 982 int n=ncols(check); 983 984 int count=1; 997 985 def A=QE(check,rec,count,0,0); 998 986 while(1) 999 987 { 1000 988 A=QE(check,rec,count,0,0); 1001 setring A; 989 setring A; 1002 990 matrix h_full=MDSmat(n,a); 1003 991 matrix rec=imap(br,rec); 1004 992 option(redSB); 1005 ideal qe_red=std(qe); 993 ideal qe_red=std(qe); 1006 994 if (qe_red[1]!=1) 1007 995 { 1008 996 break; 1009 997 } 1010 else 998 else 1011 999 { 1012 1000 count++; … … 1014 1002 setring br; 1015 1003 } 1016 1017 setring A; 1018 1004 1005 setring A; 1006 1019 1007 //obtain a codeword 1020 1008 //this works only if our code is indeed can correct these errors … … 1023 1011 { 1024 1012 syn[i,1]=-qe_red[n-i+1]+lead(qe_red[n-i+1]); 1025 } 1026 1027 matrix real_syn=inverse(h_full)*syn; 1013 } 1014 1015 matrix real_syn=inverse(h_full)*syn; 1028 1016 setring br; 1029 1017 matrix real_syn=imap(A,real_syn); 1030 return(rec-transpose(real_syn)); 1031 } 1032 example 1033 { 1034 "EXAMPLE:"; echo = 2; 1035 //correct 1 error in [15,7] binary code 1018 return(rec-transpose(real_syn)); 1019 } example 1020 { 1021 "EXAMPLE:"; echo = 2; 1022 //correct 1 error in [15,7] binary code 1036 1023 int t=1; int q=16; int n=15; int redun=10; 1037 ring r=(q,a),x,dp; 1038 1024 ring r=(q,a),x,dp; 1025 1039 1026 //generate random check matrix 1040 matrix h=random_check(redun,n,1); 1027 matrix h=random_check(redun,n,1); 1041 1028 matrix g=dual_code(h); 1042 1029 matrix x[1][n-redun]=0,0,1,0,1,0,1; 1043 1030 matrix y[1][n]=enc(x,g); 1044 1031 print(y); 1045 1032 1046 1033 // find out the minimum distance of the code 1047 1034 list l=mindist(h); 1048 1035 1049 1036 //disturb with errors 1050 1037 "Correct ",(l[1]-1) div 2," errors"; 1051 1038 matrix rec[1][n]=error_rand(y,(l[1]-1) div 2,1); 1052 print(rec); 1053 1039 print(rec); 1040 1054 1041 //let us decode 1055 1042 matrix dec_word=decode(h,rec); 1056 print(dec_word); 1043 print(dec_word); 1057 1044 } 1058 1045 … … 1064 1051 if # is given it sets the correction capacity explicitly. It should be used in case one expexts some lower bound, 1065 1052 otherwise the procedure tries to compute the real minimum distance to find out the error-correction capacity 1066 RETURN: nothing; 1067 EXAMPLE: example solve_for_random; shows an example 1053 RETURN: nothing; 1054 EXAMPLE: example solve_for_random; shows an example 1068 1055 " 1069 1056 { … … 1076 1063 option(redSB); 1077 1064 def br=basering; 1078 matrix h_full=MDSmat(n,a); 1065 matrix h_full=MDSmat(n,a); 1079 1066 matrix z[1][ncols(h_full)]; 1080 1067 int n=ncols(h_full); … … 1094 1081 dist=tmp[1]; 1095 1082 printf("d= %p",dist); 1096 t=(dist-1) div 2; 1097 } 1098 tim2=rtimer; 1083 t=(dist-1) div 2; 1084 } 1085 tim2=rtimer; 1099 1086 tim3=rtimer; 1100 1087 def A=QE(h,z,t,0,0); … … 1103 1090 ideal sys2,sys3; 1104 1091 matrix h=imap(br,h); 1105 matrix g=dual_code(h); 1092 matrix g=dual_code(h); 1106 1093 ideal sys=qe; 1107 1094 print("The system is generated"); 1108 timdec3=timdec3+rtimer-tim3; 1095 timdec3=timdec3+rtimer-tim3; 1109 1096 for (j=1; j<=ntrials; j++) 1110 1097 { 1111 word=randomvector(n-redun,1); 1098 word=randomvector(n-redun,1); 1112 1099 y=enc(transpose(word),g); 1113 rec=error_rand(y,t,1); 1114 sys2=add_synd(rec,h,redun,sys); 1115 1100 rec=error_rand(y,t,1); 1101 sys2=add_synd(rec,h,redun,sys); 1102 1116 1103 tim=rtimer; 1117 sys3=std(sys2); 1104 sys3=std(sys2); 1118 1105 sys3; 1119 timdec=timdec+rtimer-tim; 1120 } 1121 timdec2=timdec2+rtimer-tim2; 1122 } 1106 timdec=timdec+rtimer-tim; 1107 } 1108 timdec2=timdec2+rtimer-tim2; 1109 } 1123 1110 printf("Time for mindist: %p", timdist); 1124 1111 printf("Time for GB in mindist: %p", timdist); 1125 1112 printf("Time for decoding: %p", timdec2); 1126 1113 printf("Time for GB in decoding: %p", timdec); 1127 printf("Time for QE in decoding: %p", timdec3); 1128 } 1129 example 1114 printf("Time for QE in decoding: %p", timdec3); 1115 } example 1130 1116 { 1131 1117 "EXAMPLE:"; echo = 2; 1132 1118 int q=32; int n=25; int redun=n-11; int t=redun+1; 1133 1119 ring r=(q,a),x,dp; 1134 1120 1135 1121 // correct 2 errors in 5 random binary codes, 50 trials each 1136 solve_for_random(n,redun,5,50,2); 1122 solve_for_random(n,redun,5,50,2); 1137 1123 } 1138 1124 1139 1125 1140 1126 proc solve_for_code(matrix check, int ntrials, int #) 1141 "USAGE: solve_for_code(check, ntrials); 1127 "USAGE: solve_for_code(check, ntrials); 1142 1128 check is a check matrix for the code, ntrials = number of received vectors per code to be corrected 1143 1129 if # is given it sets the correction cpacity explicitly. It should be used in case one expexts some lower bound 1144 1130 otherwise the procedure tries to compute the real minimum distance to find out the error-correction capacity 1145 RETURN: nothing; 1146 EXAMPLE: example solve_for_code; shows an example 1131 RETURN: nothing; 1132 EXAMPLE: example solve_for_code; shows an example 1147 1133 " 1148 1134 { … … 1157 1143 option(redSB); 1158 1144 def br=basering; 1159 matrix h_full=MDSmat(n,a); 1145 matrix h_full=MDSmat(n,a); 1160 1146 matrix z[1][ncols(h_full)]; 1161 1147 int n=ncols(h_full); … … 1173 1159 dist=tmp[1]; 1174 1160 printf("d= %p",dist); 1175 t=(dist-1) div 2; 1176 } 1177 tim2=rtimer; 1161 t=(dist-1) div 2; 1162 } 1163 tim2=rtimer; 1178 1164 tim3=rtimer; 1179 1165 def A=QE(h,z,t,0,0); … … 1182 1168 ideal sys2,sys3; 1183 1169 matrix h=imap(br,h); 1184 matrix g=dual_code(h); 1170 matrix g=dual_code(h); 1185 1171 ideal sys=qe; 1186 1172 print("The system is generated"); 1187 timdec3=timdec3+rtimer-tim3; 1173 timdec3=timdec3+rtimer-tim3; 1188 1174 for (j=1; j<=ntrials; j++) 1189 1175 { 1190 word=randomvector(n-redun,1); 1176 word=randomvector(n-redun,1); 1191 1177 y=enc(transpose(word),g); 1192 rec=error_rand(y,t,1); 1193 sys2=add_synd(rec,h,redun,sys); 1194 1178 rec=error_rand(y,t,1); 1179 sys2=add_synd(rec,h,redun,sys); 1180 1195 1181 tim=rtimer; 1196 sys3=std(sys2); 1182 sys3=std(sys2); 1197 1183 sys3; 1198 timdec=timdec+rtimer-tim; 1199 } 1200 timdec2=timdec2+rtimer-tim2; 1201 1184 timdec=timdec+rtimer-tim; 1185 } 1186 timdec2=timdec2+rtimer-tim2; 1187 1202 1188 printf("Time for mindist: %p", timdist); 1203 1189 printf("Time for GB in mindist: %p", timdist); 1204 1190 printf("Time for decoding: %p", timdec2); 1205 1191 printf("Time for GB in decoding: %p", timdec); 1206 printf("Time for QE in decoding: %p", timdec3); 1207 } 1208 example 1192 printf("Time for QE in decoding: %p", timdec3); 1193 } example 1209 1194 { 1210 1195 "EXAMPLE:"; echo = 2; … … 1212 1197 ring r=(q,a),x,dp; 1213 1198 matrix check=random_check(redun,n,1); 1214 1215 // correct 2 errors in using the code above, 50 trials 1216 solve_for_code(check,50,2); 1199 1200 // correct 2 errors in using the code above, 50 trials 1201 solve_for_code(check,50,2); 1217 1202 } 1218 1203 … … 1227 1212 return(result); 1228 1213 } 1229 1214 1230 1215 1231 1216 static proc add_synd (matrix rec, matrix check, int redun, ideal sys) … … 1234 1219 matrix s[redun][1]=syndrome(check,rec); 1235 1220 for (int i=1; i<=redun; i++) 1236 1237 { 1238 result[i]=result[i]-s[i,1]; 1221 1222 { 1223 result[i]=result[i]-s[i,1]; 1239 1224 } 1240 1225 return(result); … … 1290 1275 for (int i=1; i<=size(G); i++) 1291 1276 { 1292 if (m/leadmonom(G[i])!=0) {return(1);} 1277 if (m/leadmonom(G[i])!=0) {return(1);} 1293 1278 } 1294 1279 return(0); … … 1298 1283 "USAGE: vanish_id (points,e); points is a list of point that define, where polynomials from the vanishing ideal will vanish, 1299 1284 RETURN: Vanishing ideal corresponding to the given set of points 1300 EXAMPLE: example vanish_id; shows an example 1285 EXAMPLE: example vanish_id; shows an example 1301 1286 " 1302 1287 { 1303 1288 int m=size(points[1]); 1304 1289 int n=size(points); 1305 1290 1306 1291 ideal G=1; 1307 1292 int i,k,j; … … 1309 1294 poly h,cur; 1310 1295 for (k=1; k<=n; k++) 1311 { 1296 { 1312 1297 i=find_index(G,points[k]); 1313 cur=G[i]; 1298 cur=G[i]; 1314 1299 for(j=i+1; j<=size(G); j++) 1315 1300 { … … 1319 1304 temp=ideal2list(G); 1320 1305 temp=delete(temp,i); 1321 G=list2ideal(temp); 1306 G=list2ideal(temp); 1322 1307 for (j=1; j<=m; j++) 1323 1308 { … … 1325 1310 { 1326 1311 attrib(G,"isSB",1); 1327 h=NF((x(j)-points[k][j,1])*cur,G); 1312 h=NF((x(j)-points[k][j,1])*cur,G); 1328 1313 temp=ideal2list(G); 1329 1314 temp=insert(temp,h); 1330 1315 G=list2ideal(temp); 1331 G=sort(G)[1]; 1316 G=sort(G)[1]; 1332 1317 } 1333 1318 } … … 1335 1320 attrib(G,"isSB",1); 1336 1321 return(G); 1337 } 1338 example 1322 } example 1339 1323 { 1340 1324 "EXAMPLE:"; echo = 2; 1341 1325 ring r=3,(x(1..3)),dp; 1342 1326 1343 1327 //generate all 3-vectors over GF(3) 1344 1328 list points=points_gen(3,1); 1345 1329 1346 1330 list points2=conv_points(points); 1347 1331 1348 1332 //grasps the first 11 points 1349 1333 list p=grasp_list(points2,1,11); 1350 1334 print(p); 1351 1335 1352 1336 //construct the vanishing ideal 1353 1337 ideal id=vanish_id(p); … … 1375 1359 count++; 1376 1360 for (j=1; j<=charac^(e)-1; j++) 1377 { 1361 { 1378 1362 result[count][m]=a^j; 1379 1363 count++; … … 1381 1365 return(result); 1382 1366 } 1383 list prev=points_gen(m-1,e); 1367 list prev=points_gen(m-1,e); 1384 1368 for (i=1; i<=size(prev); i++) 1385 1369 { … … 1396 1380 return(result); 1397 1381 } 1398 1382 1399 1383 if (e==1) 1400 1384 { … … 1403 1387 int i,j; 1404 1388 list l=ringlist(basering); 1405 int charac=l[1][1]; 1389 int charac=l[1][1]; 1406 1390 list tmp; 1407 1391 for (i=1; i<=charac^m; i++) … … 1412 1396 { 1413 1397 for (j=0; j<=charac-1; j++) 1414 { 1398 { 1415 1399 result[count][m]=number(j); 1416 1400 count++; … … 1418 1402 return(result); 1419 1403 } 1420 list prev=points_gen(m-1,e); 1404 list prev=points_gen(m-1,e); 1421 1405 for (i=1; i<=size(prev); i++) 1422 1406 { … … 1430 1414 return(result); 1431 1415 } 1432 1416 1433 1417 } 1434 1418 … … 1467 1451 { 1468 1452 poly prod=1; 1469 list rl=ringlist(basering); 1453 list rl=ringlist(basering); 1470 1454 int charac=rl[1][1]; 1471 1455 int l; … … 1509 1493 " 1510 1494 { 1511 list rl=ringlist(basering); 1512 int charac=rl[1][1]; 1495 list rl=ringlist(basering); 1496 int charac=rl[1][1]; 1513 1497 int n=ncols(check); 1514 int m=nrows(check); 1498 int m=nrows(check); 1515 1499 list points=points_gen(s,e); 1516 1500 list points2=conv_points(points); 1517 list p=grasp_list(points2,1,n); 1518 ideal id=vanish_id(p,e); 1519 ideal funcs=gener_funcs(check,p,e,id,s); 1520 1501 list p=grasp_list(points2,1,n); 1502 ideal id=vanish_id(p,e); 1503 ideal funcs=gener_funcs(check,p,e,id,s); 1504 1521 1505 ideal result; 1522 1506 poly temp; 1523 1507 int i,j,k; 1524 1525 //vanishing realtions 1508 1509 //vanishing realtions 1526 1510 for (i=1; i<=t; i++) 1527 { 1511 { 1528 1512 for (j=1; j<=size(id); j++) 1529 1513 { 1530 temp=id[j]; 1514 temp=id[j]; 1531 1515 for (k=1; k<=s; k++) 1532 1516 { 1533 1517 temp=subst(temp,x(k),x_var(i,k,s)); 1534 } 1518 } 1535 1519 result=result,temp; 1536 1520 } 1537 } 1538 1521 } 1522 1539 1523 //field equations 1540 1524 for (i=1; i<=t; i++) … … 1549 1533 result=result,e(i)^(charac^e-1)-1; 1550 1534 } 1551 1535 1552 1536 result=simplify(result,8); 1553 1537 1554 1538 //check realtions 1555 1539 poly sum; 1556 matrix syn[m][1]=syndrome(check,y); 1540 matrix syn[m][1]=syndrome(check,y); 1557 1541 for (i=1; i<=size(funcs); i++) 1558 { 1542 { 1559 1543 sum=0; 1560 1544 for (j=1; j<=t; j++) … … 1565 1549 temp=subst(temp,x(k),x_var(j,k,s)); 1566 1550 } 1567 sum=sum+temp*e(j); 1551 sum=sum+temp*e(j); 1568 1552 } 1569 1553 result=result,sum-syn[i,1]; 1570 1554 } 1571 1555 1572 1556 result=simplify(result,2); 1573 1557 1574 1558 points=points2; 1575 export points; 1559 export points; 1576 1560 return(result); 1577 } 1578 example 1561 } example 1579 1562 { 1580 1563 "EXAMPLE:"; echo = 2; 1581 1564 1582 1565 list l=FLpreprocess(3,1,11,2,""); 1583 1566 def r=l[1]; 1584 1567 setring r; 1585 int s_work=l[2]; 1586 1568 int s_work=l[2]; 1569 1587 1570 //the check matrix of [11,6,5] ternary code 1588 1571 matrix h[5][11]=1,0,0,0,0,1,1,1,-1,-1,0, … … 1597 1580 matrix rec[1][11]=error(y,list(2,4),list(1,-1)); 1598 1581 1599 //the Fitzgerald-Lax system 1582 //the Fitzgerald-Lax system 1600 1583 ideal sys=FLsystem(h,rec,2,1,s_work); 1601 1584 print(sys); … … 1605 1588 // the points are (0,0,1) and (0,1,0) with error values 1 and -1 resp. 1606 1589 // use list points to find error positions; 1607 points; 1590 points; 1608 1591 } 1609 1592 … … 1615 1598 { 1616 1599 s++; 1617 } 1600 } 1618 1601 list var_ord; 1619 1602 int i,j; … … 1633 1616 count++; 1634 1617 } 1635 } 1636 1618 } 1619 1637 1620 list rl; 1638 1621 list tmp; 1639 1622 1640 1623 if (e>1) 1641 1624 { … … 1653 1636 rl[1]=p; 1654 1637 } 1655 1638 1656 1639 rl[2]=var_ord; 1657 1640 1658 1641 rl[3]=tmp; 1659 rl[3][1]=tmp; 1642 rl[3][1]=tmp; 1660 1643 //rl[3][1][1]=string("dp("+string((t-1)*(s+1)+s)+"),lp("+string(s+1)+")"); 1661 1644 rl[3][1][1]=string("lp"); … … 1665 1648 v=v,1; 1666 1649 } 1667 rl[3][1][2]=v; 1650 rl[3][1][2]=v; 1668 1651 rl[3][2]=tmp; 1669 rl[3][2][1]=string("C"); 1652 rl[3][2][1]=string("C"); 1670 1653 rl[3][2][2]=intvec(0); 1671 1672 rl[4]=ideal(0); 1673 1654 1655 rl[4]=ideal(0); 1656 1674 1657 def r2=ring(rl); 1675 1658 setring r2; 1676 list l=ringlist(r2); 1659 list l=ringlist(r2); 1677 1660 if (e>1) 1678 1661 { 1679 execute(string("poly f="+minp)); 1680 ideal id=f; 1662 execute(string("poly f="+minp)); 1663 ideal id=f; 1681 1664 l[1][4]=id; 1682 1665 } 1683 1666 1684 1667 def r=ring(l); 1685 setring r; 1686 1668 setring r; 1669 1687 1670 return(list(r,s)); 1688 1671 } 1689 1672 1690 1673 static proc x_var (int i, int j, int s) 1691 { 1674 { 1692 1675 return(x1(s*(i-1)+j)); 1693 1676 } 1694 1677 1695 1678 static proc randomvector(int n, int e, int #) 1696 { 1679 { 1697 1680 int i; 1698 1681 matrix result[n][1]; … … 1701 1684 result[i,1]=asElement(random_prime_vector(e,#)); 1702 1685 } 1703 return(result); 1686 return(result); 1704 1687 } 1705 1688 … … 1709 1692 int i; 1710 1693 number w=1; 1711 if (size(l)>1) {w=par(1);} 1694 if (size(l)>1) {w=par(1);} 1712 1695 for (i=0; i<=size(l)-1; i++) 1713 1696 { … … 1721 1704 if (#==1) 1722 1705 { 1723 list rl=ringlist(basering); 1706 list rl=ringlist(basering); 1724 1707 int charac=rl[1][1]; 1725 1708 } else { … … 1736 1719 1737 1720 proc FL_solve_for_random(int n, int redun, int p, int e, int t, int ncodes, int ntrials, string minpol) 1738 "USAGE: FL_solve_for_random(redun,p,e,n,t,ncodes,ntrials,minpol); n = length of codes generated, redun = redundancy of codes generated, 1739 p = characteristics, e is the extension degree, 1721 "USAGE: FL_solve_for_random(redun,p,e,n,t,ncodes,ntrials,minpol); n = length of codes generated, redun = redundancy of codes generated, 1722 p = characteristics, e is the extension degree, 1740 1723 q = number of errors to correct, ncodes = number of random codes to be processed 1741 1724 ntrials = number of received vectors per code to be corrected … … 1744 1727 EXAMPLE: example FLsystem; shows an example 1745 1728 { 1746 list l=FLpreprocess(p,e,n,t,minpol); 1747 1729 list l=FLpreprocess(p,e,n,t,minpol); 1730 1748 1731 def r=l[1]; 1749 int s_work=l[2]; 1732 int s_work=l[2]; 1750 1733 export(s_work); 1751 1734 setring r; 1752 1735 1753 1736 int i,j; 1754 1737 matrix h, g, word, y, rec; … … 1756 1739 int dist, tim, tim2, tim3, timdist, timdec, timdist2, timdec2, timdec3; 1757 1740 ideal sys, sys2, sys3; 1758 list tmp; 1759 1760 option(redSB); 1761 matrix z[1][n]; 1762 1741 list tmp; 1742 1743 option(redSB); 1744 matrix z[1][n]; 1745 1763 1746 for (i=1; i<=ncodes; i++) 1764 1747 { 1765 h=random_check(redun,n,e,1); 1748 h=random_check(redun,n,e,1); 1766 1749 g=dual_code(h); 1767 tim2=rtimer; 1768 tim3=rtimer; 1769 sys=FLsystem(h,z,t,e,s_work); 1750 tim2=rtimer; 1751 tim3=rtimer; 1752 sys=FLsystem(h,z,t,e,s_work); 1770 1753 timdec3=timdec3+rtimer-tim3; 1771 1754 1772 1755 for (j=1; j<=ntrials; j++) 1773 1756 { … … 1775 1758 y=enc(transpose(word),g); 1776 1759 rec=error_rand(y,t,e); 1777 sys2=LF_add_synd(rec,h,sys); 1760 sys2=LF_add_synd(rec,h,sys); 1778 1761 tim=rtimer; 1779 1762 sys3=std(sys2); 1780 1763 timdec=timdec+rtimer-tim; 1781 } 1764 } 1782 1765 timdec2=timdec2+rtimer-tim2; 1783 } 1784 1766 } 1767 1785 1768 printf("Time for decoding: %p", timdec2); 1786 1769 printf("Time for GB in decoding: %p", timdec); 1787 printf("Time for generating Fitzgerald-Lax system during decoding: %p", timdec3); 1788 } 1789 example 1770 printf("Time for generating Fitzgerald-Lax system during decoding: %p", timdec3); 1771 } example 1790 1772 { 1791 1773 "EXAMPLE:"; echo = 2; 1792 1774 1793 1775 // decoding for one random binary code of length 25, redundancy 14; 300 words are processed 1794 1776 FL_solve_for_random(25,14,2,1,1,1,300,""); … … 1806 1788 return(result); 1807 1789 } 1790 1791 1792 /* 1793 ////////////// SOME HARD EXAMPLES ////////////////////// 1794 ////// THAT MAYBE WILL BE DOABLE LATER /////////////// 1795 1796 1.) These random instances are not doable in <=1000 sec. 1797 1798 "EXAMPLE:"; echo = 2; 1799 int q=128; int n=120; int redun=n-40; 1800 ring r=(q,a),x,dp; 1801 solve_for_random(n,redun,1,1,6); 1802 1803 redun=n-30; 1804 solve_for_random(n,redun,1,1,8); 1805 1806 redun=n-20; 1807 solve_for_random(n,redun,1,1,12); 1808 1809 redun=n-10; 1810 solve_for_random(n,redun,1,1,24); 1811 1812 int q=256; int n=150; int redun=n-10; 1813 ring r=(q,a),x,dp; 1814 solve_for_random(n,redun,1,1,26); 1815 1816 1817 2.) Generic decoding is hard! 1818 1819 int q=32; int n=31; int redun=n-16; int t=3; 1820 ring r=(q,a),(V(1..n),U(n..1),s(redun..1)),(dp(n),lp(n),dp(redun)); 1821 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, 1822 0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1, 1823 0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0, 1824 0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0, 1825 0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0, 1826 0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0, 1827 0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1, 1828 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1, 1829 0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0, 1830 0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0, 1831 0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0, 1832 1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0, 1833 0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0, 1834 1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1, 1835 1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0, 1836 0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1, 1837 0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 1838 0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0, 1839 0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1, 1840 0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0, 1841 0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0, 1842 0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0, 1843 0,1,0,1,0,0,1,0,0,1; 1844 matrix rec[1][n]; 1845 1846 def A=QE(check,rec,t,1,2); 1847 setring A; 1848 print(qe); 1849 ideal red_qe=stdfglm(qe); 1850 1851 */
Note: See TracChangeset
for help on using the changeset viewer.
