######################################################################
## Coding.txt Save this file as  Coding.txt  to use it,              #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read `Coding.txt` <Enter>                         #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger,                                      #
##  DoronZeil at gmail dot com                                       # 
######################################################################


print(`First Written: March 2023: tested for Maple 2020 `):
print(`Version  : March 2023:  `):
print():
print(`This is Coding.txt, A Maple package`):
print(`to study Coding theory following the wonderful book of Raymond Hill, "A first course in Coding Theory" `):
print(`and the wikipedia article on the Reed-Muller code`):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/Coding.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezra1();`):
print():
 

with(combinat):


ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(`  CheckISBN, CL, CodeISBN,  Coeff1, Coeff1G, DecodeHam2, DecodeRMp, EncodeRM1, Fano, IM, Fqn,  Ham2, Ham2pcm, HamPcm, HD,  GenCode,  MinD, MinWt, OneMore, OneStep, PCMrev, RaGM, RandS, Spa, SubCube, SubS, VecToTab`):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` Coding.txt: A Maple package for  `):

print(`The MAIN procedures are`):
print(` AllCodes, BestCode, Decode, DecodeHam, DecodeRM, Encode,  EncodeRM, Ham, GenCode, Fano, FanoCode, OneCode, OneCodeL, OneMoreL, PCM, SPb, ST `):

elif nargs=1 and args[1]=AllCodes then
print(`AllCodes(q,n,d): All the codes of maximum size of min. distance d in F_q^n (aka Fqn(q,n)). Try:`):
print(`AllCodes(2,5,3);`):



elif nargs=1 and args[1]=BestCode then
print(`BestCode(q,n,d,M): One q-ary code of length n, with min. distance d with  size M. Using brute force. Try:`):
print(`BestCode(2,5,3,4);`):

elif nargs=1 and args[1]=CheckISBN then
print(`CheckISBN(x,X)L checks the ISBN code`):


elif nargs=1 and args[1]=CL then
print(`CL(A,p): Given the generating matrix A of an [n,k] linear code over G(p), finds the list of coset leades. Try`):
print(`A:=OneCodeL(7,4,2,3); CL(A,2);`):

elif nargs=1 and args[1]=CodeISBN then
print(`The ISBN code for a decimal vector of length 9, X stands for 10.`):
print(`Try:`):
print(`CodeISBN([1,2,3,4,5,6,7,8,9],X);`):

elif nargs=1 and args[1]=Coeff1 then
print(`Coeff1(m,r,x,S): Given a word of length 2^m in a Reed-Muller code and a subset S of {1,..,m} with at most r members finds the coeff. of`):
print(`z[S] assuming that there are no errors. Try:`):
print(`m:=4; r:=2; f:=z[1]*z[3]+z[1]; L:=F1qn(2,m): x:=[seq(subs({seq(z[i1]=L[i][i1],i1=1..m)},f) mod 2,i=1..nops(L))]: Coeff1(m,r,x,{1,3});`):



elif nargs=1 and args[1]=Coeff1G then
print(`Coeff1G(m,r,x,S): Given a word of length 2^m in a Reed-Muller code and a subset S of {1,..,m} with at most r members finds the coeff. of`):
print(`z[S] assuming that there are at most 2^(m-r-1)-1 errors, using majority rule. Try:`):
print(`m:=4; r:=2; f:=z[1]*z[3]+z[1]; L:=F1qn(2,m): x:=[seq(subs({seq(z[i1]=L[i][i1],i1=1..m)},f) mod 2,i=1..nops(L))]: Coeff1G(m,r,x,{1,3});`):

elif nargs=1 and args[1]=Decode then
print(`Decode(A,y,p): decoding the message y with the linear code generated by A (in standard form) mod p. Try:`):
print(`x:=[1,1]; A:=[[1,0,1,1],[0,1,0,1]]: y:=Encode(A,x,2)+[0,1,0,0] mod 2;Decode(A,y,2);`):


elif nargs=1 and args[1]=DecodeHam then
print(`DecodeHam(r,p,y): decoding the message y with the linear code generated by the binary Happing code Ham(r,p).`):
print(`Try: r:=4; p:=2; A:=Ham(r,p): ra:=rand(0..p-1): x:=[seq(ra(),i=1..(p^r-1)/(p-1)-r)]:y:=Encode(A,x,p);  j:=rand(1..(p^r-1)/(p-1))(); b:=rand(1..p-1)():v:=y+b*[0$(j-1),1,0$(nops(y)-j)] mod p; y1:=DecodeHam(r,p,v); evalb(y=y1);`):


elif nargs=1 and args[1]=DecodeHam2 then
print(`DecodeHam2(r,y): decoding the message y with the linear code generated by the binary Happing code Ham(r,2).`):
print(`Try: r:=4; A:=Ham2(r): ra:=rand(0..1): x:=[seq(ra(),i=1..2^r-1-r)]:y:=Encode(A,x,2): j:=rand(1..2^r-1)(); v:=y+[0$(j-1),1,0$(2^r-j-1)] mod 2; y1:=DecodeHam2(r,v,2): evalb(y=y1);`):

elif nargs=1 and args[1]=DecodeRM then
print(`DecodeRM(m,r,x): Given a message of length 2^m in and Reed-Muller code, decodes it assuming up to 2^(m-r-1)-1 errors. `):
print(`To test that it can correct up to 2^(m-r-1)-1 errors. Try:`):
print(`ra:=rand(0..1): r:=2; m:=5; k:=add(binomial(m,i),i=0..r): x:=[seq(ra(),i=1..k)]: y:=EncodeRM(r,m,x); y:=y+[1$(2^(m-r-1)-1),0$(nops(y)-(2^(m-r-1)-1))]: evalb(x=DecodeRM(m,r,y)); `):
print(`To test that it cannot correct  2^(m-r-1) errors. Try:`):
print(`ra:=rand(0..1): r:=2; m:=5; k:=add(binomial(m,i),i=0..r): x:=[seq(ra(),i=1..k)]: y:=EncodeRM(r,m,x); y:=y+[1$2^(m-r-1),0$(nops(y)-2^(m-r-1))]: evalb(x=DecodeRM(m,r,y)); `):

elif nargs=1 and args[1]=DecodeRMp then
print(`DecodeRMp(m,r,x): Given a message of length 2^m in and Reed-Muller code, decodes it assuming no errors. Try:`):
print(`ra:=rand(0..1): r:=2; m:=5; k:=add(binomial(m,i),i=0..r): x:=[seq(ra(),i=1..k)]: y:=EncodeRM(r,m,x); evalb(x=DecodeRMp(m,r,y));`):

elif nargs=1 and args[1]=Encode then
print(`Encode(A,x,p): The coding of the vector x in the linear code generated by A mod p. Try:`):
print(`Encode([[1,0,0,0,1,0,1],[0,1,0,0,1,1,1],[0,0,1,0,1,1,0],[0,0,0,1,0,1,1]],[1,1,1,0],2);`):

elif nargs=1 and args[1]=EncodeRM then
print(`EncodeRM(r,m,x): Given a 0-1 vector of length k=Sum(binomial(m,i),i=0..r) outputs the vector of length 2^n, correpsonding to the value. Try:`):
print(`For the wikipedia exanmple, try:`):
print(`r:=2; m:=4; x:=[1,1,0,1,0,0,1,0,1,0,1]; EncodeRM(r,m,x);`):
print(``):
print(`Also try:`):
print(`ra:=rand(0..1): r:=2; m:=5; k:=add(binomial(m,i),i=0..r); x:=[seq(ra(),i=1..k)]; EncodeRM(r,m,x);`):

elif nargs=1 and args[1]=EncodeRM1 then
print(`EncodeRM1(m,S): Given a subset S of {1,...,m} outputs the list of values in F11n(2,m) given by the monomial Z[S]. Try:`):
print(`EncodeRM1(5,{1,3});`):

elif nargs=1 and args[1]=Fano then
print(`Fano(): The Fano plane `):


elif nargs=1 and args[1]=FanoCode then
print(`FanoCode(): The (7,16,3)  binary code based on the Fano plane`):

elif nargs=1 and args[1]=Fqn then
print(`Fqn(q,n): {0,1,..,q-1}^n. Try:`):
print(`Fqn(2,3);`):

elif nargs=1 and args[1]=Ham then
print(`Ham(r,p): The generating matrix In standard form) of the binary Hamming-code Ham(r,p). Try`):
print(`Ham(4,3);`):

elif nargs=1 and args[1]=HamPcm then
print(`HamPcm(r,p): The Parity Check matrix in standard form) of the binary Hamming-code Ham(r,p), with p prime. Try:`):
print(`HamPcm(4,2);`):


elif nargs=1 and args[1]=Ham2 then
print(`Ham2(r): A generating matrix, in standard form, of the binary Hamming-code Ham(r,2). Try:`):
print(`Ham2(3);`):

elif nargs=1 and args[1]=Ham2pcm then
print(`Ham2pcm(r): The Parity Check matrix in standard form) of the binary Hamming-code Ham(r,2). Try:`):
print(`Ham2pcm(4);`):

elif nargs=1 and args[1]=HD then
print(`HD(x,y): The Hamming distance between vectors x and y. Try:`):
print(`HD([1,0,2],[0,0,2]);`):


elif nargs=1 and args[1]=GenCode then
print(`GenCode(B,p): The code generated by the basis B over GP(p).`):
print(`Try:`):
print(`GenCode([[1,2,0],[2,1,1]],3);`):

elif nargs=1 and args[1]=IM then
print(`IM(S): The incidence matrix of the list of sets S. Try:`):
print(`IM(Fano());`):

elif nargs=1 and args[1]=MinD then
print(`MinD(C): Given a code C (a set of vectors of the same length) find the min. distance. Try:`):

elif nargs=1 and args[1]=MinWt then
print(`MinWt(C): Given a code C finds its minimal weight (that is not 0)`):

elif nargs=1 and args[1]=OneCode then
print(`OneCode(q,n,d): One q-ary code of length n, with min. distance d with many members. Using brute force. Try:`):
print(`OneCode(2,5,3);`):

elif nargs=1 and args[1]=OneCodeL then
print(`OneCodeL(n,k,p,d): Tries to find the generator matrix (in standard form) of ONE [n,k,d] code. May only give part of it. Try`):
print(`OneCodeL(7,4,2,3);`):

elif nargs=1 and args[1]=OneMore then
print(`OneMore(q,n,C,d): Given a subset of F_q^n, C, finds all possible extensions in F_q^n whose distance from all the members of C is at least d. Try:`):
print(`OneMore(2,5,{[0,0,0,0,0]},3);`):


elif nargs=1 and args[1]=OneMoreL then
print(`OneMoreL(A,n,k,p,d): Given a partial n by r, say matrix over GF(p) such that the min. wt. of the spanned subspace is d finds all the`):
print(`possible extensions that still has this property. Try:`):
print(`OneMoreL([[1,0,1]],7,4,2,3);`):

elif nargs=1 and args[1]=OneStep then
print(`OneStep(q,n,d,L): Given a list of lists of codes, q-ary codes with length n and min. distance d, extends it by one step, try`):
print(`OneStep(2,7,3,[[{[0,0,0,0,0,0,0]}]]);`):

elif nargs=1 and args[1]=PCM then
print(`PCM(A,p): Inputs char. p, and a generator matrix A`):
print(`of a code, outputs the Parity-Check Matrix`):
print(`Try: `):
print(`A:=OneCodeL(7,4,3,3); PCM(A,3);`):


elif nargs=1 and args[1]=PCMrev then
print(`PCMrev(A,p): Inputs char. p, and a parity-check matrix in standard form`):
print(`of a code, outputs the generating Matrix in standard form`):
print(`Try: `):
print(`A:=OneCodeL(7,4,3,3); M:=PCM(A,3); PCMrev(A,3);`):


elif nargs=1 and args[1]=RaGM then
print(`RaGM(n,k,p): A random generating matrix for a linear code lf length n and dimension k over GF(p)  (hence with p^n elements). Try:`):
print(`RaGM(7,4,2);`):

elif nargs=1 and args[1]=RandS then
print(`RandS(q,n,M): a random subset of Fqn(q,n) with M members. Try`):
print(`RandS(2,4,5);`):

elif nargs=1 and args[1]=Spa then
print(`Spa(L,n,p): inputs a list of vectors of length n, and finds the set of all linear combinations over F_p^n spanned by L. Try:`):
print(`Spa([[1,1,0],[0,1,1]],3,2);`);

elif nargs=1 and args[1]=SPb then
print(`SPb(q,n,t): The Sphere packing bound for q-ary codes of length n and correcting 2*t+1 errors. Try:`):
print(`SPb(2,5,2);`):

elif nargs=1 and args[1]=ST then
print(`ST(A,p): The Syndrom table of the Linear Code generated by the matrix A in characteristic p. Try:`):
print(`ST([[1,0,1,1],[0,1,0,1]],2);`):

elif nargs=1 and args[1]=SubCube then
print(`SubCube(m,S): The subcube of {0,1}^m where all the entries not in S are 0. Try:`):
print(`SubCube(4,{2,3});`):

elif nargs=1 and args[1]=SubCubeG then
print(`SubCubeG(m,S,v): The subcube of {0,1}^m where all the entries not in S are given by v. Try:`):
print(`SubCubeG(4,{2,3},[0,0]);`):

elif nargs=1 and args[1]=SubS then
print(`SubS(m,r): The list of subests of {1,...,m} with at most r elements, with decreasing cardinalities. Try:`):
print(`SubS(4,2);`):

elif nargs=1 and args[1]=VecToTab then
print(`VecToTab(v,m): Given a 0-1 vector v of length 2^m representing a code-word in the Reed-Muller code, translates it into a table format. Try: `):
print(`VecToTab([1,0,1,1],2);`):

else
 print(`There is no such thing as`, args):

fi:


end:



#HD(x,y): The Hamming distance between vectors x and y. Try:
#HD([1,0,2],[0,0,2]);

HD:=proc(x,y) local i,co:

co:=0:

for i from 1 to nops(x) do
 if x[i]<>y[i] then
   co:=co+1:
 fi:
od:
co:
end:

#Fqn(q,n): {0,1,..,q-1}^n. Try:
#Fqn(2,3);
Fqn:=proc(q,n) local gu,gu1,i:
option remember:

if n<0 then
 RETURN({}):
fi:

if n=0 then
 RETURN({[]}):
fi:

gu:=Fqn(q,n-1):

{seq(seq([op(gu1),i],gu1 in gu),i=0..q-1)}:

end:

#RandS(q,n,M): a random subset of Fqn(q,n) with M members. Try
#RandS(2,4,5);
RandS:=proc(q,n,M) local gu,mu,i,newguy:
mu:=Fqn(q,n):

gu:={}:
for i from 1 to M do
 newguy:=mu[rand(1..nops(mu))()]:
 gu:=gu union {newguy}:
 mu:=mu minus {newguy}:
od:
gu:
end:

#MinD(C): Given a code C (a set of vectors of the same length) find the min. distance. Try:
MinD:=proc(C) local i,j:
min(seq(seq(HD(C[i],C[j]),j=i+1..nops(C)),i=1..nops(C))):
end:

#SPb(q,n,t): The Sphere packing bound for q-ary codes of length n and correcting 2*t+1 errors. Try:
#SPb(2,5,2);
SPb:=proc(q,n,t) local i:trunc(q^n/add(binomial(n,i)*(q-1)^i,i=0..t)):end:


#OneMore(q,n,C,d): Given a subset of F_q^n, C, finds all possible extensions in F_q^n whose distance from all the members of C is at least d. Try:
#OneMore(2,5,{[0,0,0,0,0]},3);
OneMore:=proc(q,n,C,d) local gu,gu1,c,mu:
gu:=Fqn(q,n) minus C:

mu:={}:
for gu1 in gu do
 if min(seq(HD(gu1,c),c in C))>=d then
  mu:=mu union {C union {gu1}}:
 fi:
od:
mu:
end:



#AllCodes(q,n,d): All the codes of maximum size of min. distance d in F_q^n (aka Fqn(q,n)). Try:
#AllCodes(2,5,3);
AllCodes:=proc(q,n,d) local gu,lu,mu,gu1:
gu:={{[0$n]}}:

lu:={seq(op(OneMore(q,n,gu1,d)), gu1 in gu)}:

while lu<>{} do

mu:={seq(op(OneMore(q,n,gu1,d)), gu1 in lu)}:
gu:=lu:
lu:=mu:
od:
gu:
end:


#OneCode(q,n,d): One q-ary code of length n, with min. distance d with M members. Using Backtracking
OneCode:=proc(q,n,d) local C1,C2:
C1:={{[0$n]}}:
C2:=OneMore(q,n,C1[1],d):

if C2={} then
 RETURN(C1):
fi:

while C2<>{} do
C1:=C2[1]:
C2:=OneMore(q,n,C1,d):
od:
convert(C1,list):
end:


#Fano(): The Fano plane
Fano:=proc():
[{1,2,5},{1,3,7},{5,6,7},{1,4,6},{2,4,7},{3,4,5},{2,3,6}];
end:


#IM(S): Inputs a list of sets and outputs the 0-1 matrix that
#is M[i,j]=1 iff i belongs to S[j] 
IM:=proc(S) local T,i,j,V:
V:={seq(op(S[i]),i=1..nops(S))}:
for i from 1 to nops(S) do
 for j from 1 to nops(V) do
   if member(V[j],S[i]) then
    T[i,j]:=1:
   else
    T[i,j]:=0:
    fi:
  od:
od:
[seq([seq(T[i,j],j=1..nops(V))],i=1..nops(S))];
end:

#FanoCode(): The (7,16,3)  binary code based on the Fano plane
FanoCode:=proc() local A,B:
A:=IM(Fano()):
B:=subs({0=1,1=0},A):
{op(A),op(B),[0$7],[1$7]}:

end:

#The ISBN code for a decimal vector of length 9, X stands for 10.
#Try:
#CodeISBN([1,2,3,4,5,6,7,8,9]);
CodeISBN:=proc(x,X) local i,a:
if nops(x)<>9 then
 RETURN(FAIL):
fi:
a:=add(i*x[i],i=1..9) mod 11:
if a=10 then
 a:=X:
fi:
[op(x),a]:
end:


#CheckISBN(x,X)L checks the ISBN code
CheckISBN:=proc(x,X) local i,x1,y:

if nops(x)<>10 then
 RETURN(FAIL):
fi:

if x[10]=X then
 x1:=[op(1..9,x),10]:
 else
 x1:=x:
fi:

y:=add(i*x[i],i=1..10) mod 11:

if y<>0 then
 false:
else
 true:
fi:

end:


#GenCode(B,p): The code generated by the basis B over GP(p).
#Try:
#GenCode([[1,2,0],[2,1,1]],3);
GenCode:=proc(B,p) local B1, S, a, s; 
if nops(B) = 1 then 

   RETURN({seq(`mod`(a*B[1], p), a = 0 .. p-1)}) end if; B1 := 

   [op(1 .. nops(B)-1, B)]; S := GenCode(B1, p); {seq(seq(`mod`(\

  s+a*B[nops(B)], p), a = 0 .. p-1), s in S)} :
end:




#Spa(L,n,p): inputs a list of vectors, and finds the set of all linear combinations over F_p^n spanned by L
Spa:=proc(L,n,p) local S,L1,i,s,NG:
if L=[] then
 RETURN({[0$n]}):
fi:


L1:=[op(1..nops(L)-1,L)]:
S:=Spa(L1,n,p):
NG:=L[nops(L)]:
{seq(seq(s+i*NG mod p,i=0..p-1),s in S)}:
end:


#OneStep(q,n,d,L): Given a list of lists of codes, q-ary codes with length n and min. distance d, extends it by one step, try
#OneStep(2,7,3,[[{[0,0,0,0,0,0]}]]):
OneStep:=proc(q,n,d,L) local C,gu:

if L=[] then
 RETURN(FAIL):
fi:

if L[-1]=[] then
RETURN([op(1..nops(L)-2,L),[op(2..nops(L[-2]),L[-2])]]):
fi:

C:=L[-1][1]:


gu:=OneMore(q,n,C,d):


if gu<>{} then
 RETURN([op(L),convert(gu,list)]):
else
RETURN([op(1..nops(L)-1,L),[op(2..nops(L[-1]),L[-1])]]):
fi:
end:




#BestCode(q,n,d,M): One q-ary code of length n, with min. distance d with M members. Using Backtracking
BestCode:=proc(q,n,d,M) local L:
L:=[[{[0$n]}]]:



while nops(L[nops(L)][1])<M do
 L:=OneStep(q,n,d,L):
od:

L[nops(L)][1]:
end:


#MinWt(C): Given a code C finds its minimal weight
MinWt:=proc(C) local n,gu,w,C1:

if C={} then
 RETURN(FAIL):
fi:

gu:={seq(nops(w),w in C)}:
if nops(gu)<>1 then
 RETURN(FAIL):
fi:

n:=gu[1]:



C1:=convert(C,set) minus {[0$n]}:

min(seq(HD([0$n],w), w in C1))  :

end:



#RaGM(n,k,p): A random generating matrix for a linear code lf length n and dimension k over GF(p)  (hence with p^n elements). Try:
#RaGM(7,4,2);
RaGM:=proc(n,k,p) local ra,i,j:
ra:=rand(0..p-1):
[seq([0$(i-1),1,0$(k-i),seq(ra(),j=1..n-k)],i=1..k)]:
end:

#OneMoreL(A,n,k,p,d): Given a partial n by r, say matrix over GF(p) such that the min. wt. of the spanned subspace is d finds all the
#possible extensions that still has this property. Try:
#OneMoreL([[1,0]],5,3,2,3);
OneMoreL:=proc(A,n,k,p,d) local i,r,lu,lu1,gu,A1,A2:
r:=nops(A):

if r>=k then
 RETURN(FAIL):
fi:

A1:=[seq([0$(i-1),1,0$(k-i),op(A[i])],i=1..r)]:


lu:=Fqn(p,n-k):


gu:={}:
for lu1 in lu do
A2:=[op(A1),[0$r,1,0$(k-r-1),op(lu1)]]:

if MinWt(Spa(A2,n,p))=d then
gu:=gu union {[op(A),lu1]}:
fi:
od:
gu:


end:



#OneCodeL(n,k,p,d): Tries to find the generator matrix (in standard form) of ONE [n,k,d] code. May only give part of it. Try
#OneCode(7,4,2,3);
OneCodeL:=proc(n,k,p,d) local L,L1,L2,i:
L:={[[1$(n-k)]]}:
L1:=OneMoreL(L[1],n,k,p,d):



while L1<>{} and L1<>FAIL do
 L2:=OneMoreL(L1[1],n,k,p,d):
 L:=L1:
 L1:=L2:
od:
L:=L[1]:



[seq([0$(i-1),1,0$(k-i),op(L[i])],i=1..nops(L))]:
end:
 



PCMold:=proc(A,p) local i,M,M1,j,n,k:
option remember:
k:=nops(A):
n:=nops(A[1]):

if not (type(A,list) and nops(A)=k and nops(A[1])=n) then
 RETURN(FAIL):
fi:
for i from 1 to k do  
if [op(1..k,A[i])]<>[0$(i-1),1,0$(k-i)] then
 RETURN(FAIL):
fi:
od:
#[I_k|A]
#[-A^T,I_{n-k}]
M:=[]:
for i from 1 to n-k do

 M1:=[seq(-A[j][i+k] mod p,j=1..k),0$(i-1),1,0$(n-k-i)]:
M:=[op(M),M1]:
od:
M:
end:


#PCM(A,p):  Inputs char. p, and a generator matrix A
#of a code, outputs the Parity-Check Matrix
PCM:=proc(A,p) local i,n,k,A1:
option remember:
k:=nops(A):
n:=nops(A[1]):

if not (type(A,list) and nops(A)=k and nops(A[1])=n) then
 RETURN(FAIL):
fi:
for i from 1 to k do  
if [op(1..k,A[i])]<>[0$(i-1),1,0$(k-i)] then
 RETURN(FAIL):
fi:
od:

A1:=-Tran1([seq([op(k+1..nops(A[i]),A[i])],i=1..k)]) mod p:

[seq([op(A1[i]),0$(i-1),1,0$(n-k-i)],i=1..nops(A1))]:
end:


#PCMrev(A,p):  Inputs char. p, and a parity-check matrix in standard form outputs the corresponding
#geneator matrix in standard-form
PCMrev:=proc(A,p) local i,n,k,A1:
option remember:
n:=nops(A[1]):
k:=n-nops(A):


for i from 1 to n-k do  
if [op(k+1..n,A[i])]<>[0$(i-1),1,0$(n-k-i)] then
 RETURN(FAIL):
fi:
od:

A1:=-Tran1([seq([op(1..k,A[i])],i=1..n-k)]) mod p:

[seq([0$(i-1),1,0$(k-i),op(A1[i])],i=1..nops(A1))]:
end:




#CL(A,p): Given the generating matrix A of an [n,k] linear code over G(p), finds the list of coset leades. Try
#A:=OneCodeL(7,4,2,3); CL(A,2);
CL:=proc(A,p) local n,gu,lu,LED1,d1,gu1,lu1:
n:=nops(A[1]):

gu:=Fqn(p,n):


lu:=Spa(A,n,p):



LED1:=[[0$n]]:

gu:=gu minus lu:


while gu<>{} do
 d1:=MinWt(gu):


 for gu1 in gu while HD(gu1,[0$n])<>d1 do od:

 LED1:=[op(LED1),gu1]:
 gu:=gu minus {seq(lu1+gu1 mod p,lu1 in lu)}:
od:
LED1:
end:

#ST(A,p): The Syndrom table of the Linear Code generated by the matrix A in characteristic p. Try:
#ST([[1,0,1,1],[0,1,0,1]],2);
ST:=proc(A,p) local T,M,gu,y,lu1,lu,i1,j1:
option remember:
M:=PCM(A,p):
gu:=convert(CL(A,p),set):

lu:={}:

for y in gu do
 lu1:=[seq(add(M[i1][j1]*y[j1],j1=1..nops(y)) mod p,i1=1..nops(M))]:
 lu:=lu union {lu1}:
 T[lu1]:=y:
od:
lu,op(T):
end:


#Encode(A,x,p): The coding of the vector x in the linear code generated by A mod p. Try
#Encode([[1,0,0,0,1,0,1],[0,1,0,0,1,1,1],[0,0,1,0,1,1,0],[0,0,0,1,0,1,1]],[1,1,1,0],2);
Encode:=proc(A,x,p) local i:
add(A[i]*x[i],i=1..nops(x)) mod p:
end:

#Decode(A,y,p): decoding the message y with the linear code generated by A (in standard form) mod p. Try:
#x:=[1,1,1,0]; A:=[[1,0,0,0,1,0,1],[0,1,0,0,1,1,1],[0,0,1,0,1,1,0],[0,0,0,1,0,1,1]]: y:=Encode(A,x,2)+[0,1,0,0] mod 2;Decode(A,y,2);
Decode:=proc(A,y,p) local lu,M,e,i1,j1:
M:=PCM(A,p):
lu:=ST(A,p):
e:=[seq(add(M[i1][j1]*y[j1],j1=1..nops(y)) mod p,i1=1..nops(M))]:
if not member(e,lu[1]) then
 RETURN(FAIL):
fi:
y-lu[2][e] mod p:
end:



#Tran1(A): The transpose of the matrix A (given a list-of-lists)
Tran1:=proc(A) local i,j:
[seq([seq(A[i][j],i=1..nops(A))],j=1..nops(A[1]))]:
end:

#Ham2pcm(r): The Parity Check matrix in standard form) of the binary Hamming-code Ham(r,2). Try:
#Ham2pcm(4);
Ham2pcm:=proc(r) local gu,i:
gu:=Fqn(2,r) minus {[0$r]}:

gu:=gu minus {seq([0$(i-1),1,0$(r-i)],i=1..r)}:
gu:=convert(gu,list):

[op(gu),seq([0$(i-1),1,0$(r-i)],i=1..r)]:

end:


#Ham2(r): The generating matrix In standard form) of the binary Hamming-code Ham(r,2)
Ham2:=proc(r): PCMrev(Tran1(Ham2pcm(r)),2):end:




#DecodeHam2(r,y): decoding the message y with the linear code generated by the binary Happing code Ham(r,2).
#Try: r:=4; A:=Ham2(r): ra:=rand(0..1): x:=[seq(ra(),i=1..2^r-1-r)]:y:=Encode(A,x,2): v:=y+[1,0$(2^r-2)] mod 2; y1:=DecodeHam2(r,v,2): evalb(y=y1):

DecodeHam2:=proc(r,y) local M,e,i:
M:=Ham2pcm(r):
e:=add(M[i]*y[i],i=1..nops(y)) mod 2:

if e=[0$nops(e)] then
 RETURN(y):
fi:

for i from 1 to nops(M) while M[i]<>e do od:
if i=nops(M)+1 then
 RETURN(FAIL):
fi:


y+[0$(i-1),1,0$(nops(y)-i)] mod 2:

end:


#FNZ(v): Given a vector v finds the  first non-zero entry
FNZ:=proc(v) local i:
for i from 1 to nops(v) while v[i]=0 do od:
if i=nops(v)+1 then
 FAIL:
else
 v[i]:
fi:
v[i]:
end:





#HamPcm(r,p): The Parity Check matrix in standard form) of the binary Hamming-code Ham(r,p), with p prime. Try:
#HamPcm(4,2);
HamPcm:=proc(r,p) local gu,i,mu,gu1:
gu:=Fqn(p,r) minus {[0$r]}:

gu:=gu minus {seq([0$(i-1),1,0$(r-i)],i=1..r)}:
mu:={}:

for gu1 in gu do

if FNZ(gu1)=1 then
 mu:=mu union {gu1}:
fi:

od:

[op(mu),seq([0$(i-1),1,0$(r-i)],i=1..r)]:

end:

#Ham(r,p): The generating matrix In standard form) of the binary Hamming-code Ham(r,p)
Ham:=proc(r,p): PCMrev(Tran1(HamPcm(r,p)),p):end:



#IsMul(u,v,p): u/v mod p if it exists otherwise FAIL 
IsMul:=proc(u,v,p) local b:
if nops(u)<>nops(v) then
 RETURN(FAIL):
fi:

for b from 1 to p-1 do
if b*v-u mod p=[0$nops(v)] then
 RETURN(b):
fi:
od:
FAIL:
end:






#DecodeHam(r,p,y): decoding the message y with the linear code generated by the binary Happing code Ham(r,p).
#Try: r:=4; A:=Ham(r,2): ra:=rand(0..1): x:=[seq(ra(),i=1..2^r-1-r)]:y:=Encode(A,x,2): v:=y+[1,0$(2^r-2)] mod 2; y1:=DecodeHam(r,2,v): evalb(y=y1):

DecodeHam:=proc(r,p,y) local M,e,i,b:
M:=HamPcm(r,p):
e:=add(M[i]*y[i],i=1..nops(y)) mod p:

if e=[0$nops(e)] then
 RETURN(y):
fi:

for i from 1 to nops(M) while IsMul(e,M[i],p)=FAIL do od:
if i=nops(M)+1 then
 RETURN(FAIL):
fi:

b:=IsMul(e,M[i],p):

y-b*[0$(i-1),1,0$(nops(y)-i)] mod p:

end:



#SubS(m,r): The list of subests of {1,...,m} with at most r elements
SubS:=proc(m,r) local i,L:
L:=[seq(i,i=1..m)]:
[seq(op(choose(L,r-i)),i=0..r)]:
end:


#F1qn(q,n): {0,1,..,q-1}^n. as a list Try:
#F1qn(2,3);
F1qn:=proc(q,n) local gu,gu1,i:
option remember:

if n<0 then
 RETURN([]):
fi:

if n=0 then
 RETURN([[]]):
fi:

gu:=F1qn(q,n-1):

[seq(seq([i,op(gu1)],gu1 in gu),i=0..q-1)]:

end:

#EncodeRMold(r,m,x): Given a 0-1 vector of length k=Sum(binomial(m,i),i=0..r) outputs the vector of length 2^n, correpsonding to the value. Try:
#ra:=rand(0..1): r:=2; m:=5; k:=add(binomial(m,i),i=0..r): x:=[seq(ra(),i=1..k)]: RM(r,m,x);

EncodeRMold:=proc(r,m,x) local P,z,i,L1,L2,s,i1:
L1:=SubS(m,r):
if nops(L1)<>nops(x) then
 RETURN(FAIL):
fi:

P:=add(x[i]*mul(z[s], s in L1[i]),i=1..nops(L1)):


L2:=F1qn(2,m):

[seq(subs({seq(z[i1]=L2[i][i1],i1=1..m)},P) mod 2,i=1..nops(L2))]:
end:



#EncodeRM1(m,S): Given a subset S of {1,...,m} outputs the list of values in F11n(2,m) given by the monomial Z[S]. Try:
#EncodeRM1(5,{1,3}):
EncodeRM1:=proc(m,S) local P,s,L,i,i1,z:

P:=mul(z[s], s in S):
L:=F1qn(2,m):

[seq(subs({seq(z[i1]=L[i][i1],i1=1..m)},P) mod 2,i=1..nops(L))]:
end:

#EncodeRM(r,m,x): Given a 0-1 vector of length k=Sum(binomial(m,i),i=0..r) outputs the vector of length 2^n, correpsonding to the value. Try:
#ra:=rand(0..1): r:=2; m:=5; k:=add(binomial(m,i),i=0..r): x:=[seq(ra(),i=1..k)]: RM(r,m,x);
EncodeRM:=proc(r,m,x) local L1,i:
L1:=SubS(m,r):
if nops(L1)<>nops(x) then
 RETURN(FAIL):
fi:

add(x[i]*EncodeRM1(m,L1[i]),i=1..nops(L1)) mod 2:
end:



#SubCube(m,S): The subcube of {0,1}^m where all the entries not in S are 0. Try:
#SubCube(4,{2,3});
SubCube:=proc(m,S) local gu,i,S1:
option remember:
if not {seq(i, i in S)} minus {seq(i,i=1..m)}={} then
 RETURN(FAIL):
fi:

if S={} then
RETURN({[0$m]}):
fi:

if not member(m,S) then
 gu:=SubCube(m-1,S):
 RETURN([seq([op(gu[i]),0],i=1..nops(gu))]):
else
S1:=S minus {m}:
gu:=SubCube(m-1,S1) :
RETURN([seq([op(gu[i]),0],i=1..nops(gu)),seq([op(gu[i]),1],i=1..nops(gu))]):
fi:

end:

#VecToTab(v,m): Given a 0-1 vector v of length 2^m representing a code-word in the Reed-Muller code, translates it into a table format
#Try(VecToTab([1,0,1,1],2);
VecToTab:=proc(v,m) local T,gu,i:

if  not nops(v)=2^m then
 RETURN(FAIL):
fi:

gu:=F1qn(2,m):

for i from 1 to nops(gu) do
T[gu[i]]:=v[i]:
od:

op(T):

end:



#Coeff1(m,r,x,S): Given a word of length 2^m in a Reed-Muller code and a subset S of {1,..,m} with at most r members finds the coeff. of
#z[S] assuming that there are no errors. Try:
#m:=4; r:=2; f:=z[1]*z[3]+z[1]; L:=F1qn(2,m): x:=[seq(subs({seq(z[i1]=L[i][i1],i1=1..m)},f) mod 2,i=1..nops(L))]: Coeff1(m,r,x,[1,3]);
Coeff1:=proc(m,r,x,S) local T,L1,i,lu,lu1:
if not (nops(x)=2^m and r<=m and nops(S)<=r) then
 RETURN(FAIL):
fi:

L1:=F1qn(2,m):

for i from 1 to nops(L1) do
 T[L1[i]]:=x[i]:
od:

lu:=SubCube(m,S):

add(T[lu1],lu1 in lu) mod 2:

end:


#DecodeRMp(m,r,x): Given a message of length 2^m in and Reed-Muller code, decodes it assuming no errors. Try:
#ra:=rand(0..1): r:=2; m:=5; k:=add(binomial(m,i),i=0..r): x:=[seq(ra(),i=1..k)]: y:=EncodeRM(r,m,x); evalb(x=DecodeRMp(m,r,y)):
DecodeRMp:=proc(m,r,x) local lu,y,i,c,S,H:

if nops(x)<>2^m then
 RETURN(FAIL):
fi:
lu:=SubS(m,r):
y:=x:
H:=[]:
for i from 1 to nops(lu) do
S:=convert(lu[i],set):
c:=Coeff1(m,r,y,S):
H:=[op(H),c]:
if c=1 then
y:=y-EncodeRM1(m,S):
fi:
od:
H:

end:





#SubCubeG(m,S,v): The subcube of {0,1}^m where all the entries not in S are given by v. Try:
#SubCubeG(4,{2,3},[0,0]);
SubCubeG:=proc(m,S,v) local gu,i,S1:
option remember:
if not {seq(i, i in S)} minus {seq(i,i=1..m)}={} then
 RETURN(FAIL):
fi:

if not nops(S)+nops(v)=m then
 RETURN(FAIL):
fi:

if S={} then
RETURN({v}):
fi:

if v=[] then
 RETURN({op(Fqn(2,m))}):
fi:

if not member(m,S) then
 gu:=SubCubeG(m-1,S,[op(1..nops(v)-1,v)]):
 RETURN([seq([op(gu[i]),v[nops(v)]],i=1..nops(gu))]):
else
S1:=S minus {m}:
gu:=SubCubeG(m-1,S1,v) :
RETURN([seq([op(gu[i]),0],i=1..nops(gu)),seq([op(gu[i]),1],i=1..nops(gu))]):
fi:

end:



#Coeff1G(m,r,x,S): Given a word of length 2^m in a Reed-Muller code and a subset S of {1,..,m} with at most r members finds the coeff. of
#z[S] assuming that there are at most 2^(m-r-1)-1 errors, using majority rule. Try:
#m:=4; r:=2; f:=z[1]*z[3]+z[1]; L:=F1qn(2,m): x:=[seq(subs({seq(z[i1]=L[i][i1],i1=1..m)},f) mod 2,i=1..nops(L))]: Coeff1G(m,r,x,[1,3]);
Coeff1G:=proc(m,r,x,S) local T,L1,i,lu,lu1,L2,v,c0,c1,ka,katan:
if not (nops(x)=2^m and r<=m and nops(S)<=r) then
 RETURN(FAIL):
fi:

L1:=F1qn(2,m):

for i from 1 to nops(L1) do
 T[L1[i]]:=x[i]:
od:

L2:=convert(F1qn(2,m-nops(S)),set):

c0:=0:
c1:=0:

for v in L2 do



lu:=SubCubeG(m,S,v):

ka:=add(T[lu1],lu1 in lu) mod 2:


if ka=0 then
 c0:=c0+1:
else
 c1:=c1+1:
fi:
od:


katan:=min(c0,c1):

if katan>2^(m-r-1)-1 then
 RETURN(FAIL):
else
if c0>c1 then
  RETURN(0):
else
 RETURN(1):
fi:
fi:


end:



#DecodeRM(m,r,x): Given a message of length 2^m in and Reed-Muller code, decodes it assuming up to 2^(m-r-1)-1 errors. Try:
#ra:=rand(0..1): r:=2; m:=5; k:=add(binomial(m,i),i=0..r): x:=[seq(ra(),i=1..k)]: x:=x+[1,0$(2^m-1)]: y:=EncodeRM(r,m,x); evalb(x=DecodeRM(m,r,y)):
DecodeRM:=proc(m,r,x) local lu,y,i,c,S,H:

if nops(x)<>2^m then
 RETURN(FAIL):
fi:
lu:=SubS(m,r):
y:=x:
H:=[]:
for i from 1 to nops(lu) do
S:=convert(lu[i],set):
c:=Coeff1G(m,r,y,S):
if c=FAIL then
 RETURN(FAIL):
else
H:=[op(H),c]:
if c=1 then
y:=y-EncodeRM1(m,S):
fi:

fi:
od:
H:

end: