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