LIB "decode.lib";
LIB "poly.lib";

// Coopers philosphy
// begin Cooper_decode

// binary cyclic [15,7,5] code with a defining set (1,3)
  
list defset=1,3; // defining set  
int n=15; // length
int e=2; // error-correcting capacity
int q=2; // basefield size
int m=4; // degree extension of the splitting field
int sala=1; // indicator to add additional equations
   
def A=CRHT(n,defset,e,q,m); 
setring A;
A; // shows the ring we are working in
print(crht); // the CRHT-ideal 
option(redSB);
ideal red_crht=std(crht);
// reduced Groebner basis
red_crht;
  
A=CRHT(n,defset,e,q,m,sala); 
setring A;  
print(crht); // the CRHT-ideal with additional equations from Sala
option(redSB);
ideal red_crht=std(crht);
// reduced Groebner basis
print(red_crht);
// general error-locator polynomial for this code
poly gen_err_loc_poly=red_crht[5];
" ";
gen_err_loc_poly;

// change the ring to work in the splitting field
list l=ringlist(basering);
l[1][4]=ideal(a4+a+1);
def B=ring(l);
setring B;
poly gen_err_loc_poly=imap(A,gen_err_loc_poly);

// now the code itself
int k=7;
int redun=n-k;
// its check matrix
matrix check[redun][n]=1,1,0,1,0,0,0,1,0,0,0,0,0,0,0,
0,1,1,0,1,0,0,0,1,0,0,0,0,0,0,
0,0,1,1,0,1,0,0,0,1,0,0,0,0,0,
0,0,0,1,1,0,1,0,0,0,1,0,0,0,0,
0,0,0,0,1,1,0,1,0,0,0,1,0,0,0,
0,0,0,0,0,1,1,0,1,0,0,0,1,0,0,
0,0,0,0,0,0,1,1,0,1,0,0,0,1,0,
0,0,0,0,0,0,0,1,1,0,1,0,0,0,1;

// generator matrix
matrix g=dual_code(check);
// some messege
matrix x[1][k]=1,1,0,1,0,0,1;
// encode it with the code
matrix y[1][n]=enc(x,g);
//disturb with 2 errors
matrix rec[1][n]=error(y,list(4,7),list(1,1));

// check rows for defining set (1,3)
matrix checkrow1[1][n];
matrix checkrow3[1][n];
int i;
number work=a;
for (i=0; i<=n-1; i++)
{
  checkrow1[1,i+1]=work^i;
}
work=a^3;
for (i=0; i<=n-1; i++)
{
  checkrow3[1,i+1]=work^i;
}
// compute syndromes
matrix s1mat=checkrow1*transpose(rec);
matrix s3mat=checkrow3*transpose(rec);
number s1=number(s1mat[1,1]);
number s3=number(s3mat[1,1]);

// evaluate general error-locator polynomial at these syndromes
poly specialized_gen=substitute(gen_err_loc_poly,X(1),s1,X(2),s3);
// find its roots
factorize(specialized_gen);
// check
a3;
a6;

// exercise: do the same for the message (1,0,0,0,0,0,1) and errors on positions 2 and 13

x=1,0,0,0,0,0,1;
y=enc(x,g);
//disturb with 2 errors
rec=error(y,list(2,13),list(1,1));

// compute syndromes
s1mat=checkrow1*transpose(rec);
s3mat=checkrow3*transpose(rec);
s1=number(s1mat[1,1]);
s3=number(s3mat[1,1]);

// evaluate general error-locator polynomial at these syndromes
specialized_gen=substitute(gen_err_loc_poly,X(1),s1,X(2),s3);
// find its roots
factorize(specialized_gen);
// check
a;
a12;

// end Cooper_decode

// begin Cooper_mindist

// binary cyclic [15,7] code with a defining set (1,3)
// proof that its mindist is 5 with the use of GB
  
list defset=1,3; // defining set  
int n=15; // length
int w=1; // candidate for the minimum distance
option(redSB);

while (1)
{      
  def A=CRHT_mindist_binary(n,defset,w); 
  setring A;
  A; // shows the ring we are working in
  print(crht_md); // the Sala's ideal for mindist
  ideal red_crht_md=std(crht_md);
  // reduced Groebner basis
  print(red_crht_md);
  if (red_crht_md[1]!=1)
  {
    break;
  }
  else
  {
    w++;
  }
}

"Minimum distance is: ",w;

// exercise
// can we leave out the polynomials p(15,Z_i,Z_j)=0, 1 <= i < j <= w ?

w=1; // candidate for the minimum distance
option(redSB);

while (1)
{      
  def A=CRHT_mindist_binary(n,defset,w); 
  setring A;  
  ideal crht_md_ex=ideal(crht_md[1..size(defset)+w]);
  ideal red_crht_md_ex=std(crht_md_ex);
  // reduced Groebner basis
  print(red_crht_md_ex);
  if (red_crht_md_ex[1]!=1)
  {
    break;
  }
  else
  {
    w++;
  }
}

"What we obtain is: ",w;

// end Cooper_mindist

// Newton identities
// begin Newton

// We want to decode 3 errors with binary cyclic [31,16,7] code with a defining set (1,5,7)

// Newton identities for a binary 3-error-correcting cyclic code of length 31
     
int n=31; // length
list defset=1,5,7; // definig set
int t=3; // number of errors
int q=2; // basefield size
int m=5; // degree extension of the splitting field
     
def A=Newton(n,defset,t,q,m);
setring A;
A; // shows the ring we are working in
print(newton); // generalized Newton identities

// change the ring to work in the splitting field
list l=ringlist(basering);
l[1][4]=ideal(a5+a2+1);
def B=ring(l);
setring B;
ideal newton=imap(A,newton);

// code description
int k=16;
int redun=n-k;
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,
0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,
0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,
0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,
0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,
0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,
1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,
1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,
1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,0,1,
0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,0,1,
0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,1,1,0,1,1,0,0,0,1,0,1,0,0,1,0,0,1,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,1,1,0,0,
0,1,0,1,0,0,1,0,0,1;

// generator matrix
matrix g=dual_code(check);
// some messege
matrix x[1][k]=0,1,1,1,0,1,1,1,0,0,1,1,1,0,0,1;
// encode it with the code
matrix y[1][n]=enc(x,g);
//disturb with 3 errors
matrix rec[1][n]=error(y,list(3,10,24),list(1,1,1));

// check rows for defining set (1,5,7)
matrix checkrow1[1][n];
matrix checkrow5[1][n];
matrix checkrow7[1][n];
int i;
number work=a;
for (i=0; i<=n-1; i++)
{
  checkrow1[1,i+1]=work^i;
}
work=a^5;
for (i=0; i<=n-1; i++)
{
  checkrow5[1,i+1]=work^i;
}
work=a^7;
for (i=0; i<=n-1; i++)
{
  checkrow7[1,i+1]=work^i;
}
// compute syndromes
matrix s1mat=checkrow1*transpose(rec);
matrix s5mat=checkrow5*transpose(rec);
matrix s7mat=checkrow7*transpose(rec);
number s1=number(s1mat[1,1]);
number s5=number(s5mat[1,1]);
number s7=number(s7mat[1,1]);

// substitute the known syndromes in the system
ideal specialize_newton;
for (i=1; i<=size(newton); i++)
{
  specialize_newton[i]=substitute(newton[i],S(1),s1,S(5),s5,S(7),s7);
}
option(redSB);
// find sigmas
ideal red_spec_newt=std(specialize_newton);
red_spec_newt;

// find the roots of erro-locator
ring solve=(2,a),Z,lp;minpoly=a5+a2+1;
poly error_loc=Z3+(a4+1)*Z2+(a^4+a^3+a^2)*Z+(a^3);
factorize(error_loc);
// check
a2;
a9;
a23;

// exercise
// try to solve for 2-error correcting [15,7,5] binary cyclic code with a defining set (1,3) formally
// and then use "closed formulas" to decode the message 1,1,0,1,0,0,1 from the first example

int n=15; // length
list defset=1,3; // defining set
int t=2; // number of errors
int q=2; // basefield size
int m=4; // degree extension of the splitting field
     
def A=Newton(n,defset,t,q,m);
setring A;
ideal red_newton=std(newton);
red_newton;
// obtained some closed formulas
// take e.g. 4 and 7
ideal closed_form=red_newton[4],red_newton[7];

// from the first example we know that S1=a^2, S3=a, so substitute!
closed_form[1]=substitute(closed_form[1],S(1),a2,S(3),a);
closed_form[2]=substitute(closed_form[2],S(1),a2,S(3),a);

// change the ring to work in the splitting field
list l=ringlist(basering);
l[1][4]=ideal(a4+a+1);
def B=ring(l);
setring B;
ideal closed_form=imap(A,closed_form);
// finally find concrete values for sigmas
closed_form=std(closed_form);
closed_form;

// find the roots of error-locator
ring solve=(2,a),Z,lp;minpoly=a4+a+1;
poly error_loc=Z2+(a2)*Z+(a^3+a);
factorize(error_loc);
// check
a3;
a6;

// end Newton

// Fitzgerald-Lax
// begin fitzgerald-lax

// we want to decode 2 errors with [11,6,5] ternary code

int n=11;
int k=6;
int redun=n-k;
int q=3; 
int t=2;

// some preparation
list l=FLpreprocess(q,1,n,t,"");
def r=l[1];
setring r;
// s in the definition of affine code
int s_work=l[2];
     
//the check matrix of [11,6,5] ternary code
matrix check[redun][n]=1,0,0,0,0,1,1,1,-1,-1,0,
                       0,1,0,0,0,1,1,-1,1,0,-1,
                       0,0,1,0,0,1,-1,1,0,1,-1,
                       0,0,0,1,0,1,-1,0,1,-1,1,
                       0,0,0,0,1,1,0,-1,-1,1,1;
matrix g=dual_code(check);
// some message
matrix x[1][k]=-1,0,1,-1,-1,1;
matrix y[1][n]=enc(x,g);
//disturb with 2 errors
matrix rec[1][n]=error(y,list(2,4),list(1,-1));

//the Fitzgerald-Lax system
ideal FL=FLsystem(check,rec,t,1,s_work);
print(FL);
option(redSB);
ideal red_FL=std(FL);
red_FL; // read the solutions from this redGB
// the points are (0,0,1) and (0,1,0) with error values 1 and -1 resp.
// use list points to find error positions;
points;

// end fitzgerald-lax
$