/*  Reed-Solomon Codes  and encoding
*/

//  An example field.  
q := 9;
n := q-1;		// Length of the code.
k := 4; 		// Dimension of the code.
Fq<a> := GF(q);        //Creates the field.

// Let C be the Reed-Solomon code of dimension k and length n=k-1
// constructed in the standard form identified in Sec. 5 
// of the course notes.
// 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.

// Careful Magma starts iterations on the right.

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]];

//  Check the the product is zero.
print "The generator matrix", G; 
print "The check matrix", H; 
print "Their product", G*H;


// ********************************************
// For encoding it is useful to think of RS codes as cyclic codes,
// generated by an ideal in the polynomial ring over Fq.
P<x>  := PolynomialRing(Fq);

gx := &*[x-a^i : i in [1..n-k] ] ;  // Generator poly for C.

gvec := Vec_n ! [Coefficient(gx,i) :  i in [0..n-1] ];

print "\n";
print "The generator polynomial";  gx;
print "The generator polynomial lives in"; Parent(gx);
print "The corresponding vector"; gvec;
print "Check that this is in the nullspace of H"; 
gvec*H;

// Here we check that a multiple of gx is also in C.
gxmult :=  x^(2)*gx;
gxmultvec:= Vec_n ! [Coefficient(gxmult,i) : i in [0..n-1] ];
print "\n";
print "A multiple of the generator polynomial"; gxmult;
print "The corresponding vector"; gxmultvec;
print "Check that this is in the nullspace of H"; 
gxmultvec*H;

// **********************************
//  Systematic encoding.

// Get a random message vector.
m := Random(Vec_k);

// Here is non-systematic encoding.
mx := P ! [m[i] : i in [1..k] ];
cx := gx*mx;
c := Vec_n ! [Coefficient(cx,i) : i in [0..n-1] ];
print "\n";
print "The message vector" ; m;
print "The message polynomial "; mx;
print "The code poynomial" ; cx;
print "The code vector";  c;
print "Check its product with the check matrix";
c*H;

// Here is systematic encoding.

mxp := mx*x^(n-k);
qx, rx := Quotrem(mxp, gx);
cx := mxp- rx;
c := Vec_n ! [ Coefficient(cx,i) : i in [0.. n-1] ];
print "\n";
print "The message vector" ; m;
print "The message polynomial, times x^(n-k) "; mxp;
print "The remander  poynomial" ; rx;
print "The code poynomial" ; cx;
print "The code vector";  c;
print "Check its product with the check matrix";
c*H;