1############################################################################# 2## 3#A decoders.gi GUAVA library Reinald Baart 4#A &Jasper Cramwinckel 5#A &Erik Roijackers 6## &David Joyner 7## This file contains functions for decoding codes 8## 9## Decodeword (unrestricted) modified 11-3-2004 10## Decodeword (linear) created 10-2004 11## Hamming, permutation, and BCH decoders moved into this file 10-2004 12## GeneralizedReedSolomonDecoderGao added 11-2004 13## GeneralizedReedSolomonDecoder added 11-2004 14## bug fix in GRS decoder 11-29-2004 15## added GeneralizedReedSolomonListDecoder and GRS<functions>, 12-2004 16## added PermutationDecodeNC, CyclicDecoder 5-2005 17## revisions to Decodeword, 11-2005 (fixed bug spotted by Cayanne McFarlane) 18## 1-2006: added bit flip decoder 19## 7-2007: Fixed several bugs in CyclicDecoder found by 20## Punarbasu Purkayastha <ppurka@umd.edu> 21## 2009-3: Bugfix for Decodeword (J. Fields and wdj) 22## 23 24 25############################################################################# 26## 27#F Decode( <C>, <v> ) . . . . . . . . . decodes the vector(s) v from <C> 28## 29## v can be a "codeword" or a list of "codeword"s 30## 31 32InstallMethod(Decode, "method for unrestricted code, codeword", true, 33 [IsCode, IsCodeword], 0, 34function(C, v) 35 local c; 36 if v in C then 37 return Codeword(v, C); 38 fi; 39 c := Codeword(v, C); 40 if HasSpecialDecoder(C) then 41 return InformationWord(C, SpecialDecoder(C)(C, c) ); 42 elif IsCyclicCode(C) or IsLinearCode(C) then 43 return Decode(C, c); 44 else 45 Error("no decoder present"); 46 fi; 47end); 48 49InstallOtherMethod(Decode, 50 "method for unrestricted code, list of codewords", 51 true, [IsCode, IsList], 0, 52function(C, L) 53 local i; 54 return List(L, i->Decode(C, i)); 55end); 56 57InstallMethod(Decode, "method for linear code, codeword", true, 58 [IsLinearCode, IsCodeword], 0, 59function(C, v) 60 local ok, c, S, syn, index, corr, Gt, i, x, F; 61 if v in C then 62 return InformationWord(C,v); 63 fi; 64 c := Codeword(v, C); 65 if HasSpecialDecoder(C) then 66 return InformationWord(C, SpecialDecoder(C)(C, c) ); 67 fi; 68 F := LeftActingDomain(C); 69 S := SyndromeTable(C); 70 syn := Syndrome(C, c); 71 ok := false; 72 for index in [1..Length(S)] do 73 if IsBound(S[index]) and S[index][2] = syn then 74 ok := true; 75 break; 76 fi; 77 od; 78 if not ok then # this should never happen! 79 Error("In Decodeword: index not found"); 80 fi; 81#This is a hack. The subtraction operation for codewords is causing an error 82#and rather than trying to understand the method selection process, I'm brute- 83#forcing things... 84 85# corr := VectorCodeword(c - S[index][1]); # correct codeword 86 corr := VectorCodeword(c) - VectorCodeword(S[index][1]); 87 x := SolutionMat(GeneratorMat(C), corr); 88 return Codeword(x,LeftActingDomain(C)); # correct "message" 89end); 90 91 92 93############################################################################# 94## 95#F Decodeword( <C>, <v> ) . . . . . . . decodes the vector(s) v from <C> 96## 97## v can be a "codeword" or a list of "codeword"s 98## 99 100#nearest neighbor 101InstallMethod(Decodeword, "method for unrestricted code, codeword", true, 102 [IsCode, IsCodeword], 0, 103function(C, v) 104 local r, c, Cwords, NearbyWords, dist; 105 if v in C then 106 return v; 107 fi; 108 r := Codeword(v, C); 109 if r in C then return r; fi; 110 if HasSpecialDecoder(C) then 111 return SpecialDecoder(C)(C, r); 112 elif IsCyclicCode(C) or IsLinearCode(C) then 113 return Decodeword(C, r); 114 else 115 dist:=Int((MinimumDistance(C)-1)/2); 116 Cwords:=Elements(C); 117 NearbyWords:=[]; 118 for c in Cwords do 119 if WeightCodeword(r-c) <= dist then 120 NearbyWords:=Concatenation(NearbyWords,[r]); 121 fi; 122 od; 123 return NearbyWords; 124 fi; 125end); 126 127 128InstallOtherMethod(Decodeword, 129 "method for unrestricted code, list of codewords", 130 true, [IsCode, IsList], 0, 131function(C, L) 132 local i; 133 return List(L, i->Decodeword(C, i)); 134end); 135 136#syndrome decoding 137InstallMethod(Decodeword, "method for linear code, codeword", true, 138 [IsLinearCode, IsCodeword], 0, 139function(C, v) 140 local ok, c, c0, S, syn, index, corr, Gt, i, x, F; 141 if v in C then 142 return v; 143 fi; 144 c := Codeword(v, C); 145 if HasSpecialDecoder(C) then 146 c0:=SpecialDecoder(C)(C, c); 147 if Size(c0) = Dimension(C) then 148 return Codeword(c0*C); 149 fi; 150 if Size(c0) = Size(c) then 151 return Codeword(c0); 152 fi; 153 return c0; 154 fi; 155 F := LeftActingDomain(C); 156 S := SyndromeTable(C); 157 syn := Syndrome(C, c); 158 ok := false; 159 for index in [1..Length(S)] do 160 if IsBound(S[index]) and S[index][2] = syn then 161 ok := true; 162 break; 163 fi; 164 od; 165 if not ok then # this should never happen! 166 Error("In Decodeword: index not found"); 167 fi; 168 corr := VectorCodeword(c - S[index][1]); # correct codeword 169 return Codeword(corr,F); 170end); 171 172 173################################################################## 174## 175## PermutationDecode( <C>, <v> ) - decodes the vector v to <C> 176## 177## Uses Leon's AutomorphismGroup in the binary case, 178## PermutationAutomorphismGroup in the non-binary case, 179## performs permutation decoding when possible. 180## Returns original vector and prints "fail" when not possible. 181## 182 183InstallMethod(PermutationDecode, "attribute method for linear codes", true, 184 [IsLinearCode,IsList], 0, 185function(C,v) 186 #C is a linear [n,k,d] code, 187 #v is a vector with <=(d-1)/2 errors 188local M,G,perm,F,P,t,p0,p,pc,i,k,pv,c,H,d,G0,H0; 189 G:=GeneratorMat(C); 190 H:=CheckMat(C); 191 d:=MinimumDistance(C); 192 k:=Dimension(C); 193 G0 := List(G, ShallowCopy); 194 perm:=PutStandardForm(G0); 195 H0 := List(H, ShallowCopy); 196 PutStandardForm(H0,false); 197 F:=LeftActingDomain(C); 198 if F=GF(2) then 199 P:=AutomorphismGroup(C); 200 else 201 P:=PermutationAutomorphismGroup(C); 202 fi; 203 t:=Int((d-1)/2); 204 p0:=0; 205 if WeightCodeword(H0*v)=0 then return(v); fi; 206 for p in P do 207 pv:=Permuted(v,p); 208 if WeightCodeword(Codeword(H0*pv))<=t then 209 p0:=p; 210 break; 211 fi; 212 od; 213 if p0=0 then Print("fail \n"); return(v); fi; 214 pc:=TransposedMat(G0)*List([1..k],i->pv[i]); 215 c:=Permuted(pc,p0^(-1)); 216 return Codeword(Permuted(c,perm)); 217end); 218 219 220################################################################## 221## 222## 223## PermutationDecodeNC( <C>, <v>, <P> ) - decodes the 224## vector v to <C> 225## 226## Performs permutation decoding when possible. 227## Returns original vector and prints "fail" when not possible. 228## NC version does Not Compute the aut gp so one must 229## input the permutation automorphism group 230## <P> subset SymmGp(n), n=wordlength(<C>) 231## 232#### wdj,2-7-2005 233 234InstallMethod(PermutationDecodeNC, "attribute method for linear codes", true, 235 [IsLinearCode,IsCodeword,IsGroup], 0, 236function(C,v,P) 237 #C is a linear [n,k,d] code, 238 #v is a vector with <=(d-1)/2 errors 239local M,G,G0, perm,F,t,p0,p,pc,i,k,pv,c,H,H0, d; 240 G:=GeneratorMat(C); 241 H:=CheckMat(C); 242 d:=MinimumDistance(C); 243 k:=Dimension(C); 244 G0 := List(G, ShallowCopy); 245 perm:=PutStandardForm(G0); 246 H0 := List(H, ShallowCopy); 247 PutStandardForm(H0,false); 248 F:=LeftActingDomain(C); 249 t:=Int((d-1)/2); 250 p0:=0; 251 if WeightCodeword(H0*v)=0 then return(v); fi; 252 for p in P do 253 pv:=Permuted(v,p); 254 if WeightCodeword(Codeword(H0*pv))<=t then 255 p0:=p; 256# Print(p0," = p0 \n",pv,"\n",H*pv,"\n"); 257 break; 258 fi; 259 od; 260 if p0=0 then Print("fail \n"); return(v); fi; 261 pc:=TransposedMat(G0)*List([1..k],i->pv[i]); 262 c:=Permuted(pc,p0^(-1)); 263 return Codeword(Permuted(c,perm)); 264end); 265 266############################################################################# 267## 268#F CyclicDecoder( <C>, <w> ) . . . . . . . . . . . . decodes cyclic codes 269## 270InstallMethod(CyclicDecoder, "method for cyclic code, codeword", true, 271 [IsCode, IsCodeword], 0, 272function(C,w) 273local d, g, wpol, s, ds, cpol, cc, c, i, m, e, x, n, ccc, r; 274 if not(IsCyclicCode(C)) then 275 Error("\n\n Code must be cyclic"); 276 fi; 277 if Codeword(w) in C then return Codeword(w); fi; ## bug fix 7-6-2007 278 n:=WordLength(C); 279 d:=MinimumDistance(C); 280 g:=GeneratorPol(C); 281 x:=IndeterminateOfUnivariateRationalFunction(g); 282 wpol:=PolyCodeword(w); 283 for i in [0..(n-1)] do 284 s:=x^i*wpol mod g; 285 ds:=DegreeOfLaurentPolynomial(s); 286 if ds<Int((d-1)/2) then 287 m:=i; 288 e:=x^(n-m)*s mod (x^n-1); 289 cpol:=wpol-e; 290 cc:=CoefficientsOfUnivariatePolynomial(cpol); 291 r:=Length(cc); 292 ccc:=Concatenation(cc,List([1..(n-r)],k-> Zero(LeftActingDomain(C)) )); 293 c:=Codeword(ccc); 294 #return InformationWord( C, c ); 295 return c; 296 fi; 297 od; 298 return "fail"; 299end); 300 301############################################################################# 302## 303#F BCHDecoder( <C>, <r> ) . . . . . . . . . . . . . . . . decodes BCH codes 304## 305InstallMethod(BCHDecoder, "method for code, codeword", true, 306 [IsCode, IsCodeword], 0, 307function (C, r) 308 local F, q, n, m, ExtF, x, a, t, ri_1, ri, rnew, si_1, si, snew, 309 ti_1, ti, qi, sigma, i, cc, cl, mp, ErrorLocator, zero, 310 Syndromes, null, pol, ExtSize, ErrorEvaluator, Fp; 311 F := LeftActingDomain(C); 312 q := Size(F); 313 n := WordLength(C); 314 m := OrderMod(q,n); 315 t := QuoInt(DesignedDistance(C) - 1, 2); 316 ExtF := GF(q^m); 317 x := Indeterminate(ExtF); 318 a := PrimitiveUnityRoot(q,n); 319 zero := Zero(ExtF); 320 r := PolyCodeword(Codeword(r, n, F)); 321 if Value(GeneratorPol(C), a) <> zero then 322 return Decode(C, r); ##LR - inf loop !!! 323 fi; 324 # Calculate syndrome: this simple line is faster than using minimal pols. 325 Syndromes := List([1..2*QuoInt(DesignedDistance(C) - 1,2)], 326 i->Value(r, a^i)); 327 if Maximum(Syndromes) = Zero(F) then # no errors 328 return Codeword(r mod (x^n-1), C); 329 fi; 330 # Use Euclidean algorithm: 331 ri_1 := x^(2*t); 332 ri := LaurentPolynomialByCoefficients( 333 ElementsFamily(FamilyObj(ExtF)), 334 Syndromes, 0); 335 rnew := ShallowCopy(ri); 336 si_1 := x^0; si := 0*x; snew := 0*x; 337 ti_1 := 0*x; ti := x^0; sigma := x^0; 338 while Degree(rnew) >= t do 339 rnew := (ri_1 mod ri); 340 qi := (ri_1 - rnew) / ri; 341 snew := si_1 - qi*si; 342 sigma := ti_1 - qi*ti; 343 ri_1 := ri; ri := rnew; 344 si_1 := si; si := snew; 345 ti_1 := ti; ti := sigma; 346 od; 347 # Chien search for the zeros of error locator polynomial: 348 ErrorLocator := []; 349 null := a^0; 350 ExtSize := q^m-1; 351 for i in [0..ExtSize-1] do 352 if Value(sigma, null) = zero then 353 AddSet(ErrorLocator, (ExtSize-i) mod n); 354 fi; 355 null := null * a; 356 od; 357 # And decode: 358 if Length(ErrorLocator) = 0 then 359 Error("not decodable"); 360 fi; 361 x := Indeterminate(F); 362 if q = 2 then # error locator is not necessary 363 pol := Sum(List([1..Length(ErrorLocator)], i->x^ErrorLocator[i])); 364 return Codeword((r - pol) mod (x^n-1), C); 365 else 366 pol := Derivative(sigma); 367 Fp := One(F)*(x^n-1); 368 #Print(ErrorLocator, "\n"); 369 ErrorEvaluator := List(ErrorLocator,i-> 370 Value(rnew,a^-i)/Value(pol, a^-i)); 371 pol := Sum(List([1..Length(ErrorLocator)], i-> 372 -ErrorEvaluator[i]*x^ErrorLocator[i])); 373 return Codeword((r - pol) mod (x^n-1) , C); 374 fi; 375end); 376 377############################################################################# 378## 379#F HammingDecoder( <C>, <r> ) . . . . . . . . . . . . decodes Hamming codes 380## 381## Generator matrix must have all unit columns 382## 383InstallMethod(HammingDecoder, "method for code, codeword", true, 384 [IsCode, IsCodeword], 0, 385function(C, r) 386 local H, S,p, F, fac, e,z,x,ind, i,Sf; 387 S := VectorCodeword(Syndrome(C,r)); 388 r := ShallowCopy(VectorCodeword(r)); 389 F := LeftActingDomain(C); 390 p := PositionProperty(S, s->s<>Zero(F)); 391 if p <> fail then 392 z := Z(Characteristic(S[p]))^0; 393 if z = S[p] then 394 fac := One(F); 395 else 396 fac := S[p]/z; 397 fi; 398 Sf := S/fac; 399 H := CheckMat(C); 400 ind := [1..WordLength(C)]; 401 for i in [1..Redundancy(C)] do 402 ind := Filtered(ind, j-> H[i][j] = Sf[i]); 403 od; 404 e := ind[1]; 405 r[e] := r[e]-fac; # correct error 406 fi; 407 #x := SolutionMat(GeneratorMat(C), r); 408 #return Codeword(x); 409 return Codeword(r, C); 410end); 411 412############################################################################# 413## 414#F GeneralizedReedSolomonDecoderGao( <C>, <v> ) . . decodes 415## generalized Reed-Solomon codes 416## using S. Gao's method 417## 418InstallMethod(GeneralizedReedSolomonDecoderGao,"method for code, codeword", true, 419 [IsCode, IsCodeword], 0, 420function(C,vec) 421 local vars,a,b,n,i,g0,g1,geea,u,q,v,r,g,f,F,R,x,k,GcdExtEuclideanUntilBound; 422 if C!.name<>" generalized Reed-Solomon code" then 423 Error("\N This method only applies to GRS codes.\n"); 424 fi; 425 426GcdExtEuclideanUntilBound:=function(F,f,g,d,R) 427## S Gao's version 428## R is a poly ring 429local vars,u,v,r,q,i,num,x; 430 vars:=IndeterminatesOfPolynomialRing(R); 431 x:=vars[1]; # define local x 432 u:=List([1..3],i->Zero(F)); ## Need u[3] as a temporary variable 433 v:=List([1..3],i->Zero(F)); 434 r:=List([1..3],i->Zero(F)); 435 q:=Zero(F); 436 u[1]:=One(F); u[2]:=Zero(F); v[1]:=Zero(F); v[2]:=One(F); 437 r[2]:=f; r[1]:=g; ### applied with r2=f=g1, r1=g=g0 438 while DegreeIndeterminate(r[2],x)> d-1 do ### assume Degree(0) < 0. 439 q := EuclideanQuotient(R,r[1],r[2]); 440 r[3]:=EuclideanRemainder(R,r[1],r[2]); 441 u[3]:=u[1] - q*u[2]; 442 v[3]:=v[1] - q*v[2]; 443 r[1]:=r[2]; r[2] :=r[3]; 444 u[1]:=u[2]; u[2] :=u[3]; 445 v[1]:=v[2]; v[2] :=v[3]; 446 od; 447 return([r[2],u[2],v[2]]); 448end; 449 450 a:=C!.points; 451 R:=C!.ring; 452 k:=C!.degree; 453 F:=CoefficientsRing(R); 454 b:=VectorCodeword(vec); 455 vars:=IndeterminatesOfPolynomialRing(R); 456 x:=vars[1]; # define local x 457 n:=Length(a); 458 if Size(Set(a)) < n then 459 Error("`Points in 1st vector must be distinct.`\n\n"); 460 fi; 461 g0:=One(F)*Product([1..n],i->x-a[i]); 462 g1:=InterpolatedPolynomial(F,a,b); 463 geea:=GcdExtEuclideanUntilBound(F,g1,g0,(n+k)/2,R); 464 u:=geea[2]; v:=geea[3]; g:=geea[1]; 465 if v=Zero(F) then return("fail"); fi; 466 if v=One(F) then 467 q:=g; 468 r:=Zero(F); 469 else 470 q:=EuclideanQuotient(R,g,v); 471 r:=EuclideanRemainder(R,g,v); 472 fi; 473 if ((r=Zero(F) or LeadingCoefficient(r)=Zero(F)) and (Degree(q) < k)) then 474 f:=q; 475 else 476 f:="fail"; ## this does not occur if num errors < (mindist - 1)/2 477 fi; 478 if f="fail" then return(f); else 479 return Codeword(List(a,i->Value(f,i)),C); 480 fi; 481end); 482 483########################################################## 484# 485# Input: Pts=[x1,..,xn], l = highest power, 486# L=[h_1,...,h_ell] is list of powers 487# r=[r1,...,rn] is received vector 488# Output: Computes matrix described in Algor. 12.1.1 in [JH] 489# 490GRSLocatorMat:=function(r,Pts,L) 491 local a,j,b,ell,DiagonalPower,add_col_mat,add_row_mat,block_matrix; 492 493 add_col_mat:=function(M,N) ## "AddColumnsToMatrix" 494 #N is a matrix with same rowdim as M 495 #the fcn adjoins N to the end of M 496 local i,j,S,col,NT; 497 col:=MutableTransposedMat(M); #preserves M 498 NT:=MutableTransposedMat(N); #preserves N 499 for j in [1..DimensionsMat(N)[2]] do 500 Add(col,NT[j]); 501 od; 502 return MutableTransposedMat(col); 503 end; 504 505 add_row_mat:=function(M,N) ## "AddRowsToMatrix" 506 #N is a matrix with same coldim as M 507 #the fcn adjoins N to the bottom of M 508 local i,j,S,row; 509 row:=ShallowCopy(M);#to preserve M; 510 for j in [1..DimensionsMat(N)[1]] do 511 Add(row,N[j]); 512 od; 513 return row; 514 end; 515 516 block_matrix:=function(L) ## slightly simpler syntax IMHO than "BlockMatrix" 517 # L is an array of matrices of the form 518 # [[M1,...,Ma],[N1,...,Na],...,[P1,...,Pa]] 519 # returns the associated block matrix 520 local A,B,i,j,m,n; 521 n:=Length(L[1]); 522 m:=Length(L); 523 A:=[]; 524 if n=1 then 525 if m=1 then return L[1][1]; fi; 526 A:=L[1][1]; 527 for i in [2..m] do 528 A:=add_row_mat(A,L[i][1]); 529 od; 530 return A; 531 fi; 532 for j in [1..m] do 533 A[j]:=L[j][1]; 534 od; 535 for j in [1..m] do 536 for i in [2..n] do 537 A[j]:=add_col_mat(A[j],L[j][i]); 538 od; 539 od; 540 B:=A[1]; 541 for j in [2..m] do 542 B:= add_row_mat(B,A[j]); 543 od; 544 return B; 545 end; 546 547 DiagonalPower:=function(r,j) 548 local R,n,i; 549 n:=Length(r); 550 R:=DiagonalMat(List([1..n],i->r[i]^j)); 551 return R; 552 end; 553 554 ell:=Length(L); 555 a:=List([1..ell],j->DiagonalPower(r,(j-1))*VandermondeMat(Pts,L[j])); 556 b:=List([1..ell],j->[1,j,a[j]]); 557 return (block_matrix([a])); 558end; 559 560############################################################## 561# 562# 563# Input: Pts=[x1,..,xn], 564# L=[h_1,...,h_ell] is list of powers 565# r=[r1,...,rn] is received vector 566# 567# Compute kernel of matrix in alg 12.1.1 in [JH]. 568# Choose a basis vector in kernel. 569# Output: the lists of coefficients of the polynomial Q(x,y) in alg 12.1.1. 570# 571GRSErrorLocatorCoeffs:=function(r,Pts,L) 572 local a,j,b,vec,e,QC,i,lengths,ker,ell; 573 ell:=Length(L); 574 e:=GRSLocatorMat(r,Pts,L); 575 ker:=TriangulizedNullspaceMat(TransposedMat(e)); 576 if ker=[] then Print("Decoding fails.\n"); return []; fi; 577 vec:=ker[Length(ker)]; 578 QC:=[]; 579 lengths:=List([1..ell],i->Sum(List([1..i],j->1+L[j]))); 580 QC[1]:=List([1..lengths[1]],j->vec[j]); 581 for i in [2..ell] do 582 QC[i]:=List([(lengths[i-1]+1)..lengths[i]],j->vec[j]); 583 od; 584 return QC; 585end; 586 587 588####################################################### 589# 590# Input: List L of coefficients ell, R = F[x] 591# Pts=[x1,..,xn], 592# Output: list of polynomials Qi as in Algor. 12.1.1 in [JH] 593# 594GRSErrorLocatorPolynomials:=function(r,Pts,L,R) 595 local q,p,i,ell; 596 ell:=Length(L)+1; ## ?? Length(L) instead ?? 597 q:=GRSErrorLocatorCoeffs(r,Pts,L); 598 if q=[] then Print("Decoding fails.\n"); return []; fi; 599 p:=[]; 600 for i in [1..Length(q)] do 601 p:=Concatenation(p,[CoefficientToPolynomial(q[i],R)]); 602 od; 603 return p; 604end; 605 606########################################################## 607# 608# Input: List L of coefficients ell, R = F[x] 609# Pts=[x1,..,xn], 610# Output: interpolating polynomial Q as in Algor. 12.1.1 in [JH] 611# 612GRSInterpolatingPolynomial:=function(r,Pts,L,R) 613 local poly,i,Ry,F,y,Q,ell; 614 ell:=Length(L)+1; ## ?? Length(L) instead ?? ##### not used ??? 615Q:=GRSErrorLocatorPolynomials(r,Pts,L,R); 616 if Q=[] then Print("Decoding fails.\n"); return 0; fi; 617 F:=CoefficientsRing(R); 618 y:=IndeterminatesOfPolynomialRing(R)[2]; 619# Ry:=PolynomialRing(F,[y]); 620# poly:=CoefficientToPolynomial(Q,Ry); 621 poly:=Sum(List([1..Length(Q)],i->Q[i]*y^(i-1))); 622 return poly; 623end; 624 625############################################################################# 626## 627#F GeneralizedReedSolomonDecoder( <C>, <v> ) . . decodes 628## generalized Reed-Solomon codes 629## using the interpolation algorithm 630## 631InstallMethod(GeneralizedReedSolomonDecoder,"method for code, codeword", true, 632 [IsCode, IsCodeword], 0, 633function(C,vec) 634local v,R,k,P,z,F,f,s,t,L,n,Qpolys,vars,x,c,y; 635 636 v:=VectorCodeword(vec); 637 R:=C!.ring; 638 P:=C!.points; 639 k:=C!.degree; 640 F:=CoefficientsRing(R); 641 vars:=IndeterminatesOfPolynomialRing(R); 642 x:=vars[1]; 643#y:=vars[2]; 644 n:=Length(v); 645 t:=Int((n-k)/2); 646 L:=[n-1-t,n-t-k]; 647 Qpolys:=GRSErrorLocatorPolynomials(v,P,L,R); 648 f:=-Qpolys[2]/Qpolys[1]; 649 c:=List(P,s->Value(f,[x],[s])); 650 return Codeword(c,n,F); 651end); 652 653 654############################################################################# 655## 656#F GeneralizedReedSolomonListDecoder( <C>, <v> , <ell> ) . . ell-list decodes 657## generalized Reed-Solomon codes 658## using M. Sudan's algorithm 659## 660## 661## Input: v is a received vector (a GUAVA codeword) 662## C is a GRS code 663## ell>0 is the length of the decoded list (should be at least 664## 2 to beat GeneralizedReedSolomonDecoder 665## or Decoder with the special method of interpolation decoding) 666## 667InstallMethod(GeneralizedReedSolomonListDecoder,"method for code, codeword, integer", 668 true, [IsCode, IsCodeword, IsInt], 0, 669function(C,v,ell) 670local f,h,g,x,R,R2,L,F,t,i,c,Pts,k,n,tau,Q,divisorsf,div, 671 CodewordList,p,vars,y,degy, divisorsdeg1; 672 R:=C!.ring; 673 F:=CoefficientsRing(R); 674 vars:=IndeterminatesOfPolynomialRing(R); 675 x:=vars[1]; 676 Pts:=C!.points; 677 n:=Length(Pts); 678 k:=C!.degree; 679 tau:=Int((n-k)/2); 680 L:=List([0..ell],i->n-tau-1-i*(k-1)); 681 y:=X(F,vars);; 682 R2:=PolynomialRing(F,[x,y]); 683 vars:=IndeterminatesOfPolynomialRing(R2); 684 Q:=GRSInterpolatingPolynomial(v,Pts,L,R2); 685 divisorsf:=DivisorsMultivariatePolynomial(Q,R2); 686 divisorsdeg1:=[]; 687 CodewordList:=[]; 688 for div in divisorsf do 689 degy:=DegreeIndeterminate(div,y); 690 if degy=1 then ######### div=h*y+g 691 g:=Value(div,vars,[x,Zero(F)]); 692 h:=Derivative(div,y); 693 if DegreeIndeterminate(h,x)=0 then 694 f:= -h^(-1)*g*y^0; 695 divisorsdeg1:=Concatenation(divisorsdeg1,[f]); 696 if g=Zero(F)*x then 697 c:=List(Pts,p->Zero(F)); 698 else 699 c:=List(Pts,p->Value(f,[x,y],[p,Zero(F)])); 700 fi; 701 CodewordList:=Concatenation(CodewordList,[Codeword(c,C)]); 702 fi; 703 fi; 704 od; 705 return CodewordList; 706end); 707 708############################################################################# 709## 710#F NearestNeighborGRSDecodewords( <C>, <v> , <dist> ) . . . finds all 711## generalized Reed-Solomon codewords 712## within distance <dist> from v 713## *and* the associated polynomial, 714## using "brute force" 715## 716## Input: v is a received vector (a GUAVA codeword) 717## C is a GRS code 718## dist>0 is the distance from v to search 719## Output: a list of pairs [c,f(x)], where wt(c-v)<dist 720## and c = (f(x1),...,f(xn)) 721## 722## ****** very slow***** 723## 724InstallMethod(NearestNeighborGRSDecodewords,"method for code, codeword, integer", 725 true, [IsCode, IsCodeword, IsInt], 0, 726function(C,r,dist) 727# "brute force" decoder 728local k,F,Pts,v,p,x,f,NearbyWords,c,a; 729 k:=C!.degree; 730 Pts:=C!.points; 731 F:=LeftActingDomain(C); 732 NearbyWords:=[]; 733 for v in F^k do 734 a := Codeword(v); 735 f:=PolyCodeword(a); 736 x:=IndeterminateOfLaurentPolynomial(f); 737 c:=Codeword(List(Pts,p->Value(f,[x],[p]))); 738 if WeightCodeword(r-c) <= dist then 739 NearbyWords:=Concatenation(NearbyWords,[[c,f]]); 740 fi; 741od; 742return NearbyWords; 743end); 744 745############################################################################# 746## 747#F NearestNeighborDecodewords( <C>, <v> , <dist> ) . . . finds all 748## codewords in a linear code C 749## within distance <dist> from v 750## using "brute force" 751## 752## Input: v is a received vector (a GUAVA codeword) 753## C is a linear code 754## dist>0 is the distance from v to search 755## Output: a list of c in C, where wt(c-v)<dist 756## 757## ****** very slow***** 758## 759InstallMethod(NearestNeighborDecodewords,"method for code, codeword, integer", 760 true, [IsCode, IsCodeword, IsInt], 0, 761function(C,r,dist) 762# "brute force" decoder 763local k,F,Pts,v,p,x,f,NearbyWords,c,a; 764 F:=LeftActingDomain(C); 765 NearbyWords:=[]; 766 for v in C do 767 c := Codeword(v); 768 if WeightCodeword(r-c) <= dist then 769 NearbyWords:=Concatenation(NearbyWords,[c]); 770 fi; 771od; 772return NearbyWords; 773end); 774 775 776############################################################################# 777## 778#F BitFlipDecoder( <C>, <v> ) . . . decodes *binary* LDPC codes using bit-flipping 779## or Gallager hard-decision 780## ** often fails to work if C is not LDPC ** 781## 782## Input: v is a received vector (a GUAVA codeword) 783## C is a binary LDPC code 784## Output: a c in C, where wt(c-v)<mindis(C) 785## 786## ****** fast ***** 787## 788checksetJ:=function(H,w) 789 local i,I,n,k; 790 I:=[]; 791 n:=Length(H[1]); 792 for i in [1..n] do 793 if H[w][i]<>0*Z(2) then 794 I:=Concatenation(I,[i]); 795 fi; 796 od; 797 return I; 798end; 799 800checksetI:=function(H,u) 801 return checksetJ(TransposedMat(H),u); 802end; 803 804BFcounter:=function(I,s) 805 local S,i; 806 S:=Codeword(List(I, i -> List(s)[i])); 807 return WeightCodeword(S); 808end; 809 810bit_flip:=function(H,v) 811 local i,j,s,q,n,k,rho,gamma,I_j,tH; 812 tH:=TransposedMat(H); 813 q:=ShallowCopy(v); 814 s:=H*q; 815 n:=Length(H[1]); 816 k:=n-Length(H); 817 rho:=Length(Support(Codeword(H[1]))); 818 #gamma:= Length(Support(Codeword(TransposedMat(H)[1]))); 819 for j in [1..n] do 820 gamma:= Length(Support(Codeword(tH[j]))); 821 I_j:=checksetI(H,j); 822 if BFcounter(I_j,s)>gamma/2 then 823 q[j]:=q[j]+Z(2); 824 break; 825 fi; 826 od; 827 return Codeword(q); 828end; 829 830InstallMethod(BitFlipDecoder,"method for code, codeword, integer", 831 true, [IsCode, IsCodeword], 0, 832function(C,v) 833 local H,r,qnew,q; 834 H:=CheckMat(C); 835 r:=ShallowCopy(v); 836 if Length(Support(Codeword(H*v)))=0 then 837 return v; 838 fi; 839 q:=bit_flip(H,r); 840 if Length(Support(Codeword(H*q)))=0 then 841 return Codeword(q); 842 fi; 843 while Length(Support(Codeword(H*q)))<Length(Support(Codeword(H*r))) do 844 #while q<>r do 845 qnew:=q; 846 r:=q; 847 q:=bit_flip(H,qnew); 848 Print(" ",Length(Support(Codeword(H*q)))," ",Length(Support(Codeword(H*r))),"\n"); 849 od; 850 return Codeword(r); 851end); 852 853