//Makes the F_n+1 matrices of the Barrolleta/Villanueva paper


//first makes the fn+1 matrices from simplex paper, Lemma 5  


//value for  n
WD:=WeightDistribution;
PG:=PermutationGroup;
Dim:=Dimension;

//input value for n 
d:=n;

f:= GaloisField(2);
vv := VectorSpace(f, n);
gln:=GL(n,f);
E<z>:=ext<f|n>;
z:=PrimitiveElement(E);
m<x>:=MinimalPolynomial(z);m;
vve:={x:x in vv};
zd:=z^n;
L:=(2^n-1); "Length of simplex code=",L;
v:=L;
fn:=Floor(L/n)-1;
"fn for simplex code=",fn;
Fn:=Floor(2^n/(n+1)) -1;
"Fn for RM code=",Fn;
//Z:=Integers();L:=Z!L;

SS:=SymmetricGroup(L);
//Fq:=f;Fqs:={x:x in Fq|x ne 0};
HC:=HammingCode(GF(2),n);C:=Dual(HC);
G:=GeneratorMatrix(C); gpts:=[Transpose(G)[i]:i in [1..L]];
vn:=VectorSpace(f,L);
g:=SS;
 npts:=Seqset(gpts) diff Seqset(Basis(vv));
 fpts:=Basis(vv) cat Setseq(npts);
ns:=&cat[Eltseq(fpts[i]):i in [1..#fpts]];
G:=Transpose(Matrix(L,n,ns));
//standard form
C1:=LinearCode(G);gpts:=fpts;
"weight distribution of the simplex code C1 of length 2^n -1=",L;
 WD(C1);
 
 


zz:=z^(1+n);
T:={z^i:i in [0..n]};
T0:=[Eltseq(z^i,f):i in [0..n]];
Trs:=[T];TR:=[T0];
for r:=1 to Fn do
Tr:={x*zz^r:x in T};
Trs:=Append(Trs,Tr);TR:=Append(TR,[Eltseq(x,f):x in Tr]);
end for;
ETR:=[[vv!TR[i][x]:x in [1..#TR[i]]]:i in [1..#TR]];
////



vdp:=VectorSpace(f,d+1);

TRn:=[];
for i:=1 to Fn+1 do
b:=TR[i];
nb:=[[f!1] cat b[j]:j in [1..d+1]];
TRn:=Append(TRn,nb);
end for;


MTRn:=[];
for i:=1 to Fn+1 do
b:=TRn[i];
nb:=[vdp!b[1]] cat [vdp!b[1]+vdp!b[j]:j in [2..d+1]];
MTRn:=Append(MTRn,nb);
end for;


KK:=KMatrixSpace(GF(2),n+1,n+1);
Tmats:=[KK!TRn[i]:i in [1..Fn+1]];

"Tmats have first column 1s";
Tinv:=[KK!MTRn[i]:i in [1..Fn+1]];

"Tinv is the sequence mats*";
TRpds:=[Tinv[i]^-1:i in [1..Fn+1]];

"TRpds is the F_n-PD set of F_n+1 matrices";

vdpe:={x: x in vdp};
sp:=[x:x in vdpe|x[1] eq 1];
ze:=vdp!0;
ba:=Basis(vdp);
pts:=[ba[1]] cat [ba[1]+ba[i]:i in [2..#ba]];
rpts:=Seqset(sp) diff Seqset(pts);
pts:=pts cat Setseq(rpts);
NG:=Transpose(Matrix(#pts,n+1,pts));
"NG is a generating matrix for RM(1,n), info set e_1, e_1+e_i";

NC:=LinearCode(NG);
"NC is RM(1,",n,") of length", Length(NC);
"weight distribution of  RM(1,n) of length 2^n ";
 WD(NC);


GLN:=GL(n+1,f);
SSS:=SymmetricGroup(2^(d));


//ian:=matspds;
ian:=TRpds;
pian:=[];
cpian:=[];icpian:=[];
for l:=1 to #ian do
ma:=ian[l];
aa:=[];
cian:=[];
ician:=[];
for i:=1 to #pts do
x:=pts[i]*ma;
for j:=1 to #pts do
if pts[j] in sub<vdp|x> then
for y in f do
if pts[j] eq y*x then
cian:=Append(cian,y);
ician:=Append(ician,y^-1);
end if;
end for;

aa:=Append(aa,j);
end if;
end for;end for;
aap:=SSS!aa;
pian:=Append(pian,aap);
cpian:=Append(cpian,cian);
icpian:=Append(icpian,ician);
end for;

 paut:=PG(NC);
C:=NC;v:=2^d;
 t:=Fn;au:=pian;
 t;
CG:=GeneratorMatrix(C);
mons:=[MonomialMatrix(cpian[i],pian[i]):i in [1..#ian]];
for i:=1 to #mons do
for j:=1 to n+1 do
if not (NG[j]*mons[i] in C) then
j,i;
end if;
end for;end for;

diags:=[DiagonalMatrix(icpian[i]):i in [1..#ian]];
pigs:=[diags[i]*mons[i]:i in [1..#ian]];

t:=Fn;
vd:=VectorSpace(f,2^n-n-1);
A:=&cat[[-CG[i][j]:j in [n+2..2^n]]:i in [1..n+1]]; 
MA:=Matrix(n+1,2^n-n-1,A);
AT:=Transpose(MA);
SH:=&cat[Eltseq(AT[i]) cat Eltseq(Basis(vd)[i]):i in [1..2^n-1-n]];
H:=Matrix(2^n-1-n,2^n,SH);
HT:=Transpose(H);

"5 tests decoding random sent words that have  Fn=",Fn," errors";
for i:=1 to 5 do
e:=Random(C);
e1:=e;

sl:={1..2^n};
ts:=[];
for i:=1 to t do

ts:=Append(ts,Random(sl diff Seqset(ts)));
end for;
//ts;
#Seqset(ts) eq t, "that there are", t ,"places for errors";
for i:=1 to t do
e[ts[i]]:=e[ts[i]]+1;
//Random(Fq)*e[ts[i]]+Random(Fq);
end for;
Weight(e1+e),"=no. of errors";

"after up to t errors, received vector y";
for i:=1 to #mons do
w:=Weight(e*mons[i]*HT);
i,"th member a of PD set;",w,"=weight of y^a";
if w le t then
j:=i;break;
end if;
end for;
y:=e*mons[j];
sy:=vdp![y[i]:i in [1..n+1]];
my:=sy*CG;
r:=my*mons[j]^-1;

r eq e1,"that the corrected received word is the sent word";
"--------";
end for;

omons:=mons;

/**********************************/
" The general  method for RM(r,n) using (n,M,3) code";
/**********************/
BP3:=function(u);
vu:=VectorSpace(GF(2),u);
VV:=[x:x in vu];
D3:=[VV[1]];
for i:=2 to #VV do
flag:=true;
for j:=1 to #D3 do
dist:=Distance(VV[i],D3[j]);
if (dist lt 3) then 
  flag:=false;
end if;
end for;
if flag eq true then
D3:=Append(D3,VV[i]);
end if;
end for;


return(D3);
end function;
/************************/
CT:=BP3(n);
"(n,M,3)=","(",n,#CT,3,")";
SCT:=[Eltseq(CT[i]):i in [1..#CT]];
aa:=KK!Id(GLN);
nmats:=[];
for i:=1 to #CT do
cc:=vdp!([1] cat SCT[i]);
ba:=aa; ba[1]:=cc;
nmats:=Append(nmats,ba);
end for;


//now make these into permutations on the RM code
ian:=nmats;

pian:=[];
cpian:=[];icpian:=[];
for l:=1 to #ian do
ma:=ian[l];
aa:=[];
cian:=[];
ician:=[];
for i:=1 to #pts do
x:=pts[i]*ma;
for j:=1 to #pts do
if pts[j] in sub<vdp|x> then
for y in f do
if pts[j] eq y*x then
cian:=Append(cian,y);
ician:=Append(ician,y^-1);
end if;
end for;

aa:=Append(aa,j);
end if;
end for;end for;
aap:=SSS!aa;
pian:=Append(pian,aap);
cpian:=Append(cpian,cian);
icpian:=Append(icpian,ician);
end for;

 t:=#nmats -1;
 au:=pian;
 t;
//CG:=GeneratorMatrix(C);
nmons:=[MonomialMatrix(cpian[i],pian[i]):i in [1..#ian]];
for i:=1 to #nmons do
for j:=1 to n+1 do
if not (NG[j]*mons[i] in C) then
j,i;
end if;
end for;end for;

diags:=[DiagonalMatrix(icpian[i]):i in [1..#ian]];
pigs:=[diags[i]*nmons[i]:i in [1..#ian]];

"5 tests decoding random sent words that have  #nmats-1=",t," errors";
for i:=1 to 5 do
e:=Random(C);
e1:=e;

sl:={1..2^n};
ts:=[];
for i:=1 to t do

ts:=Append(ts,Random(sl diff Seqset(ts)));
end for;
//ts;
#Seqset(ts) eq t, "that there are", t ,"places for errors";
for i:=1 to t do
e[ts[i]]:=e[ts[i]]+1;
//Random(Fq)*e[ts[i]]+Random(Fq);
end for;
Weight(e1+e),"=no. of errors";

"after up to t errors, received vector y";
for i:=1 to #nmons do
w:=Weight(e*nmons[i]*HT);
i,"th member a of PD set;",w,"=weight of y^a";
if w le t then
j:=i;break;
end if;
end for;
y:=e*nmons[j];
sy:=vdp![y[i]:i in [1..n+1]];
my:=sy*CG;
r:=my*nmons[j]^-1;

r eq e1,"that the corrected received word is the sent word";
"--------";
end for;
