// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// Decoding Reed-Solomon Codes with the
// Berlekamp-Massey algorithm.
// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
// Let C be the Reed-Solomon code of dimension k and length n=k-1
// constructed in the standard form identified in Sec. 65
// of the course notes.
// An example field.
/*
q := 9;
n := q-1; // Length of the code.
k := 4; // Dimension of the code.
Fq := GF(q); //Creates the field.
*/
q := 16;
n := q-1; // Length of the code.
k := 9; // Dimension of the code.
Fq := GF(q); //Creates the field.
// The following generate the generator matrix, G, and the check
// matrix, H, for C.
Vec_n := VectorSpace(Fq, n); // n dim vec space.
Vec_k := VectorSpace(Fq, k); // k dim vec space.
MatH := KMatrixSpace(Fq, n, n-k); // Space for check matrix.
MatG := KMatrixSpace(Fq, k, n); // Space for gen matrix.
G := MatG ! [ a^((i)*(j-1)) : j in [1..n], i in [1..k]];
H := MatH ! [ a^((i-1)*(j-1)) : j in [1..n-k], i in [1..n]];
//G := MatG ! [ a^((i-1)*(j-1)) : j in [1..n], i in [1..k]];
//H := MatH ! [ a^((i-1)*(j)) : j in [1..n-k], i in [1..n]];
G*H;
// ********************************************
// Creation of a message vector.
// Systematic encoding to a code vector.
// ********************************************
P := PolynomialRing(Fq);
gx := &*[x-a^(i-1) : i in [1..n-k] ] ; // Generator poly for C.
m := Random(Vec_k);
mx := P ! [m[i] : i in [1..k] ];
mxp := mx*x^(n-k);
qx, rx := Quotrem(mxp, gx);
cx := mxp- rx;
c := Vec_n ! [ Coefficient(cx,i) : i in [0.. n-1] ];
// ********************************************
// Corruption by an error vector of bounded weight.
// ********************************************
// The maximum numbers of errors that can be corrected is (d-1)/2.
// The code below gives some weight less than d/2.
// You might want to alter to give exactly some prescribed weight.
// Or you could simulate a channel by corrupting a symbol with
// a given probability.
t := Floor((n-k)/2);
locs := Random(Subsets({1..n},t)); // A random set of t locations.
e := Vec_n !0;
for l in locs do
e[l] := Random(Fq); // Random values at the given locations, possibly 0.
end for;
// The vector recieved.
v := c+e;
// The syndrome vector.
syn := v*H;
// ************************************************
// Some functions introduced in Sec. 5 of the notes.
// ********************************************
// We only use the function Syn() in the decoding algorithm
Syn := function(syn, fx)
print("here"); print Degree(fx), Degree(syn);
if Degree(fx) ge Degree(syn) then
print "Polynomial is of larger degree than the known syndromes";
return 0;
end if;
fvec := Parent(syn) ! [ Coefficient(fx,i) : i in [0.. n-k-1] ];
return InnerProduct(syn, fvec);
end function;
Spandiscfail := function(syn, fx)
d := Fq ! 0; i := -1;
while d eq 0 and i lt Degree(syn)-Degree(fx)-1 do
i:= i+1;
d := Syn(syn, x^i*fx);
end while;
if d eq 0 then
print "Span appears to be infinite";
return -1, 0, -1;
else
return i, d, i + Degree(fx);
end if;
end function;
Disc := func< syn, fx, shft | Syn(syn, x^shft*fx) >;
// ************************************************
// The Berlekamp-Massey algorithm
// ************************************************
// The Update procedure:
Update := procedure(~info, m , syn)
fx := info[1,1];
c := Degree(fx) - 1;
r := m - c -1;
mu := Syn(syn, x^r*fx);
print "m is", m, "f is ", fx;
print "shift of g is ", c, "shift of f is ", r, "discrepancy of f is", mu;
if (mu eq 0) then
Q := RMatrixSpace(Parent(fx),2,2) ! [1,0, 0,1];
elif (r le c ) then
Q := RMatrixSpace(Parent(fx),2,2) ! [1, -mu*x^(c-r), 0,1];
else
Q := RMatrixSpace(Parent(fx),2,2) ! [x^(r-c), -mu, mu^(-1), 0];
end if;
print "Update matrix is ", Q;
info := Q*info;
print "New info is", info;
end procedure;
// Initialization:
MatPoly := RMatrixSpace(P,2,1);
info := MatPoly ![1,0];
// Run the algorithm:
for m in [0..n-k-1] do
Update(~info, m, syn);
end for;
// **********************************************
// Find the error locations
// **********************************************
// Two functions to find the locations and the values
// (here we use the formula from Sec. 8 of the notes.
FindLocations := function(fx, vec_n)
elocs := vec_n !0;
for i in [0..Degree(vec_n)-1 ] do
if (Evaluate(fx,a ^i) eq 0) then
elocs[i+1] := 1;
end if;
end for;
return elocs;
end function;
FindVals := function(fx, gx, erlocs)
evals := erlocs;
fpr := Derivative(fx);
// Degree(evals) is just the dimension of the space it lives in.
for i in [0..Degree(evals) -1 ] do
if (evals[i+1] eq (Fq ! 1)) then
evals[i+1] := (Evaluate(fpr,a^i)*Evaluate(gx,a^i))^(-1) ;
end if;
end for;
return evals;
end function;
// Find the locations and values and check the results.
e_locs := FindLocations(info[1,1], Vec_n);
print "The locations of the errors are the nonzero entries in \n", e_locs;
e_out := FindVals(info[1,1], info[2,1], e_locs);
print "The error vector as determined by the algorithm is \n ", e_out;
print "The original error vector was \n ", e;