###########################################################
#
#             codes.mpl
#      basic error correcting code
#            procedures
#
#David Joyner, wdj@usna.edu, 6-98
#Reed-Muller generator matrix added 10-98
#dual_code_list added 1-1-99
#scalar_mult, projective_space, hamming_check 10-20-99
#weight_lexorder, isequivaelent, ..., automorphismgroup, 12-99
#dimension_code, minor modifications, 7-2001
#small bug in decode_hamming fixed, 10-2001
#seems to work in maple7 - however, use codes_v7.mpl instead
###########################################################

matrix_times_vector_modp:=proc(A::array,v::list,p::integer)
#returns A*v mod p, A a matrix, v a vector
#written as a list
local B,C;
C:=convert(map(B->modp(B,p),evalm(A&*v)),list);
RETURN(C);
end:

matrix_modp:=proc(A::array,p::integer)
#returns A mod p, A a matrix
local B,C;
C:=map(B->modp(B,p),evalm(A));
RETURN(evalm(C));
end:

cartesian_product:=proc(L1::list,L2::list)
#returns the cartesian product of 2 lists
 local i,j,L3;
 L3:=[];
  for i in L1 do
   for j in L2 do
    L3:=[op(L3),[i,j]];
   od;
  od;
 RETURN(L3);
end:

cartesian_power:=proc(L::list,n::integer)
#returns the n-th cartesian power
#of the list L
 local i,j,k,L0,L1;
 if n=0 then RETURN([]); fi;
 if n=1 then RETURN(L); fi;
 if n<0 then ERROR(`must have `,n,` >0`); fi;
 L0:=cartesian_product(L,L);
 if n=2 then RETURN(L0); fi;
 for k from 3 to n do
 L1:=[];
 for i in L0 do
   for j in L do
    L1:=[op(L1),[op(i),j]];
   od;
 od;
 L0:=L1;
 od;
 RETURN(L1);
end:

scalar_mult:=proc(L,p)
#returns scalar mult of L in F_p^k
#with 1st non-zero coeff =1
local i,j,c,k,i0,j0,L0;
k:=nops(L);
for i from 1 to k do
if L[i]<>0 then i0:=i; break; fi;
od;
j0:=L[i0]^(-1) mod p;
L0:=expand(j0*L) mod p;
RETURN(L0);
end:

which_scalar_mult:=proc(L1,L2,p)
#returns scalar c which multiplied by
# L1 is L2 in F_p^k: c*L1=L2
local i,j,c,k,i0,j0,L0,i1,i2,k2,k1;
k1:=nops(L1);
k2:=nops(L2):
if k1<>k2 then RETURN(`not multiples`); fi;
if scalar_mult(L1,p)<>scalar_mult(L2,p)
 then RETURN(`not multiples`); fi;
if scalar_mult(L1,p)=scalar_mult(L2,p)
 then 
  for i from 1 to k1 do
   if L1[i]<>0 then i1:=i; break; fi;
  od;
  for j from 1 to k1 do
   if L2[j]<>0 then i2:=j; break; fi;
  od;
 j0:=L1[i1]^(-1) mod p;
fi;
 RETURN(L2[i2]*j0 mod p);
end:

projective_space:=proc(k,p)
#contructs a list L of all
#vectors of length k modulo scalars in Fp
local i,j,L,v,v0,S,Fp,OO;
Fp:=[seq(i,i=0..p-1)];
S:=cartesian_power(Fp,k);
L:=[];
OO:=[seq(0,j=1..k)];
for v in S do
v0:=scalar_mult(v,p);
if not(member(v0,L)) and v0<>OO then
 L:=[op(L),v0];
fi;
od;
RETURN(L):
end:

golay12:=proc()
#returns a generator for the
#ternary (12,6) Golay code
local G,G0,I6;
I6:= array(identity, 1..6,1..6);
G0:=matrix([[0,1,1,1,1,1],[1,0,1,2,2,1],[1,1,0,1,2,2],
[1,2,1,0,1,2],[1,2,2,1,0,1],[1,1,2,2,1,0]]):
G:=concat(I6,G0);
RETURN(evalm(G));
end:

golay11:=proc()
#returns a generator for the
#perfect ternary (11,6) Golay code
local G,G0,I6;
I6:= array(identity, 1..6,1..6);
G0:=matrix([[1,1,1,1,1],[0,1,2,2,1],[1,0,1,2,2],
[2,1,0,1,2],[2,2,1,0,1],[1,2,2,1,0]]):
G:=concat(I6,G0);
RETURN(evalm(G));
end:

golay_check12:=proc()
#returns a parity check matrix for the
#ternary (12,6) Golay code
local G,G0,G1,H,I6,A;
I6:= array(identity, 1..6,1..6);
G0:=matrix([[0,1,1,1,1,1],[1,0,1,2,2,1],[1,1,0,1,2,2],
[1,2,1,0,1,2],[1,2,2,1,0,1],[1,1,2,2,1,0]]):
G1:=-transpose(G0):
H:=map(A->modp(A,3),concat(G1,I6));
RETURN(evalm(H));
end:

golay_check11:=proc()
#returns a parity check matrix for the
#perfect ternary (11,6) Golay code
local G,G0,G1,H,I5,A;
I5:= array(identity, 1..5,1..5);
G0:=matrix([[1,1,1,1,1],[0,1,2,2,1],[1,0,1,2,2],
[2,1,0,1,2],[2,2,1,0,1],[1,2,2,1,0]]):
G1:=-transpose(G0):
H:=map(A->modp(A,3),concat(G1,I5));
RETURN(evalm(H));
end:

golay23:=proc()
#returns a generator for the
#perfect binary (23,12) Golay code
local G,G0,I12;
I12:= array(identity, 1..12,1..12);
G0:=matrix([[0,1,1,1,1,1,1,1,1,1,1,1],[1,1,1,0,1,1,1,0,0,0,1,0],
[1,1,0,1,1,1,0,0,0,1,0,1],[1,0,1,1,1,0,0,0,1,0,1,1],
[1,1,1,1,0,0,0,1,0,1,1,0],[1,1,1,0,0,0,1,0,1,1,0,1],
[1,1,0,0,0,1,0,1,1,0,1,1],[1,0,0,0,1,0,1,1,0,1,1,1],
[1,0,0,1,0,1,1,0,1,1,1,0],[1,0,1,0,1,1,0,1,1,1,0,0],
[1,1,0,1,1,0,1,1,1,0,0,0],[1,0,1,1,0,1,1,1,0,0,0,1]]):
G0:=delcols(G0, 12..12);
G:=concat(I12,G0);
RETURN(evalm(G));
end:

golay24:=proc()
#returns a generator for the
#binary (24,12) Golay code
local G,G0,I12;
I12:= array(identity, 1..12,1..12);
G0:=matrix([[0,1,1,1,1,1,1,1,1,1,1,1],[1,1,1,0,1,1,1,0,0,0,1,0],
[1,1,0,1,1,1,0,0,0,1,0,1],[1,0,1,1,1,0,0,0,1,0,1,1],
[1,1,1,1,0,0,0,1,0,1,1,0],[1,1,1,0,0,0,1,0,1,1,0,1],
[1,1,0,0,0,1,0,1,1,0,1,1],[1,0,0,0,1,0,1,1,0,1,1,1],
[1,0,0,1,0,1,1,0,1,1,1,0],[1,0,1,0,1,1,0,1,1,1,0,0],
[1,1,0,1,1,0,1,1,1,0,0,0],[1,0,1,1,0,1,1,1,0,0,0,1]]):
G:=concat(I12,G0);
RETURN(evalm(G));
end:

golay_check23:=proc()
#returns a generator for the
#perfect binary (23,12) Golay code
local G0,G1,I11,A,H;
I11:= array(identity, 1..11,1..11);
G0:=matrix([[0,1,1,1,1,1,1,1,1,1,1,1],[1,1,1,0,1,1,1,0,0,0,1,0],
[1,1,0,1,1,1,0,0,0,1,0,1],[1,0,1,1,1,0,0,0,1,0,1,1],
[1,1,1,1,0,0,0,1,0,1,1,0],[1,1,1,0,0,0,1,0,1,1,0,1],
[1,1,0,0,0,1,0,1,1,0,1,1],[1,0,0,0,1,0,1,1,0,1,1,1],
[1,0,0,1,0,1,1,0,1,1,1,0],[1,0,1,0,1,1,0,1,1,1,0,0],
[1,1,0,1,1,0,1,1,1,0,0,0],[1,0,1,1,0,1,1,1,0,0,0,1]]):
G0:=delcols(G0, 12..12);
G1:=-transpose(G0):
H:=map(A->modp(A,2),concat(G1,I11));
RETURN(evalm(H));
end:

golay_check24:=proc()
#returns a generator for the
#binary (24,12) Golay code
local G0,G1,I12,A,H;
I12:= array(identity, 1..12,1..12);
G0:=matrix([[0,1,1,1,1,1,1,1,1,1,1,1],[1,1,1,0,1,1,1,0,0,0,1,0],
[1,1,0,1,1,1,0,0,0,1,0,1],[1,0,1,1,1,0,0,0,1,0,1,1],
[1,1,1,1,0,0,0,1,0,1,1,0],[1,1,1,0,0,0,1,0,1,1,0,1],
[1,1,0,0,0,1,0,1,1,0,1,1],[1,0,0,0,1,0,1,1,0,1,1,1],
[1,0,0,1,0,1,1,0,1,1,1,0],[1,0,1,0,1,1,0,1,1,1,0,0],
[1,1,0,1,1,0,1,1,1,0,0,0],[1,0,1,1,0,1,1,1,0,0,0,1]]):
G1:=-transpose(G0):
H:=map(A->modp(A,2),concat(G1,I12));
RETURN(evalm(H));
end:

### WARNING: semantics of type `string` have changed
### from V4 to V5
seg:=proc(s::string,m::integer,n::integer)
 #returns the substring of s
 #from positions m through n
# this can be combined with MAPLE's
# procedure SearchText
 local s0;
 if m>n then ERROR(`arg[2] must be less than arg[3]`);
 fi;
 s0:=substring(s,m..n);
 RETURN(s0);
end:

alphabet:=["A","B","C","D","E","F","G","H","I","J","K","L","M","N","O","P","Q","R","S","T","U","V","W","X","Y","Z","a","b","c","d","e","f","g","h","i","j","k","l","m","n","o","p","q","r","s","t","u","v","w","x","y","z"," 
",",",".","?","!"]:

### WARNING: semantics of type `string` have changed
### from V4 to V5
locate:=proc(s1::string,s2::string)
 #RETURNS list of positions where
 #s1 occurs in s2
 #case sensitive
 local i,l1,l2,L;
 L:=[]:
 l1:=length(s1);
 l2:=length(s2);
 for i from 1 to l2-l1+1 do
  if searchtext(s1,substring(s2,i..(i+l1-1)))<>0 then
    L:=[op(L),i];
  fi;
 od: RETURN(L);
end:

### WARNING: semantics of type `string` have changed
### from V4 to V5
Locate:=proc(s1::string,s2::string)
 #RETURNS list of positions where
 #s1 occurs in s2
 #case sensitive
 local i,l1,l2,L;
 L:=[]:
 l1:=length(s1);
 l2:=length(s2);
 for i from 1 to l2-l1+1 do
  if s1=substring(s2,i..(i+l1-1)) then
    L:=[op(L),i];
  fi;
 od: RETURN(L);
end:

statistics:=proc(mytext::string)
#returns the ranked list of most common
#characters and their percentages of occurrance
 local x,L;
 L:=[];
 for x in alphabet do
 L:=[op(L),[x,100*evalf(nops(Locate(x,mytext))/length(mytext))]]; 
 od:
 RETURN(sort(L,occur_order));
end:

check_matrix:=proc(g,n::integer,p::prime)
#enter gen poly of cyclic code and RETURN parity check matrix
local A,i,j,v,h;
if not(type(g, polynom(integer, x))) then 
 ERROR(`g should be a polynomial in x`);
fi;
if not(iscyclic(g,n,p)) then
 ERROR(`the polynomial doesn't generate a cyclic code`);
fi;
h:=Quo(x^n-1,g,x) mod p;
A:=generator_matrix(reciprocal(h,x),n,p);
RETURN(evalm(A));
end:

reciprocal:=proc(f,x)
#returns the reciprocal polynomial
local g;
 ### WARNING: degree(0,x) now returns -infinity
 g:=x^(degree(f))*subs(x=1/x,f);
 RETURN(expand(g));
end:

generator_matrix:=proc(g,n::integer,p::prime)
#enter gen poly of cyclic code and RETURN generator matrix
 local A,i,j,v;
 if not(iscyclic(g,n,p)) then
  ERROR(`the polynomial doesn't generate a cyclic code`);
 fi;
 ### WARNING: degree(0,x) now returns -infinity
 for i from 1 to n-degree(g) do
 v[i]:=[seq(0,j=1..(i-1)),op(poly2coeff_list(g)),
seq(0,j=(i+degree(g))..(n-1))];
 od;
 A:=matrix([seq(v[i],i=1..(n-degree(g)))]);
 RETURN(evalm(A));
end:

poly2coeff_list:=proc(f)
#RETURNS the list of coefficients of a polynomial f
local i,L;
L:=[];
### WARNING: degree(0,x) now returns -infinity
for i from 0 to degree(f) do
 L:=[op(L),coeff(f,x,i)];
od;
RETURN(L);
end: 

coeff_list2poly:=proc(L::list)
 #inverse of the above procedure
 local i,f;
 f:=add(L[i]*x^(i-1),i=1..(nops(L)));
 RETURN(f);
end:

iscyclic:=proc(g,n::integer,p::prime)
 #tests if g is the generating polynomial
 # of a cyclic code of length n over F_p
 RETURN(Divide(x^n-1,g) mod p);
end:

padic_coeff:=proc(p,k)
 #returns coeffs in p-adic expan of each 1<= n <= p^k
 local C,L,L0,i,j;
 L0:={[seq(0,i=1..k)]};
 L:=[seq(i-1,i=1..p)];
 C:=[seq(cartesian_power(L,k)[i],i=2..p^k)];
 RETURN(C);
end:

occur_order:=proc(a,b)
 #assumes a,b are pairs of the
 # form a=[c_a,f_a], c_a a character, 
 # f_a a real number,
 #returns true if f_a<f_b
  local n,i,j;
  if a[2]<b[2] then RETURN(`false`); fi;
  RETURN(`true`);
end:

weight_order:=proc(a,b)
 #assumes a,b are lists of equal length
 #returns true if wt(a)<wt(b)
 local n,i,j;
 n:=nops(a):
 for i from 1 to n do
 if weight_codeword(b)<weight_codeword(a) then RETURN(`false`); fi;
 od;
 RETURN(`true`);
end:

weight_lexorder:=proc(a,b)
 #assumes a,b are lists of equal length
 #returns true if wt(a)<wt(b) and lexorder(a)<lexorder(b)
 local n,i,j;
 n:=nops(a):
 for i from 1 to n do
 if weight_codeword(b)<weight_codeword(a) then RETURN(`false`); fi;
 if (weight_codeword(b)=weight_codeword(a) and
   lexorder(convert(a,string),convert(b,string)))
 then RETURN(`false`); fi;
od;
 RETURN(`true`);
end:

weight_codeword:=proc(w)
 #assumes w is a list
 #returns the weight of the word
 #(number of non-zero entries)
  local S,n,i;
 n:=nops(w):
 S:=[];
 for i from 1 to n do
  if w[i]<>0 then S:=[op(S),i]; fi;
 od;
  RETURN(nops(S));
end:

weight_enumerator_vector:=proc(G::array,p::integer)
 # returns the list of number of codewords of a
 # given weight of the code over F_p with generator
 # matrix G
 local L,g,A,i,w,C,cdim,c;
C:=code_list(G,p);
cdim:=coldim(G);
A[0]:=[[seq(0,i=1..cdim)]];
for i from 1 to cdim do
 A[i]:=[];
od;
for c in C do
 for i from 1 to cdim do
  if weight_codeword(c)=i then 
     A[i]:=[op(A[i]),c];
  fi;
 od;
od;
w:=[seq(nops(A[i]),i=0..cdim)];
RETURN(w);
end:

weight_enumeration_polynomial:=proc(G::array,p::integer,s,t)
 # returns the weight_enumeration_polynomial
 # of the code over F_p with generator
 # matrix G in the variables s,t
 local w,p0,p1,cdim,i;
cdim:=coldim(G);
w:=weight_enumerator_vector(G,p);
p1:=0;
for i from 0 to cdim do
 p1:=p1+w[i+1]*s^i*t^(cdim-i);
od;
p0:=(s,t)->p1;
RETURN(p0(s,t));
end:

min_distance:=proc(G::array,p::integer)
#returns minumum distance between two
#codewords in the code over F_p
#generated by the matrix G
local C,w,d;
C:=code_list(G,p);
C:=sort(C, weight_order);
d:=weight_codeword(C[2]);
#the 1st codeword is 0, so we take the second
#which has minimal non-zero weight
RETURN(d);
end:

##This proc seems to have a mistake.
#hamming:=proc(k)
#local G0,G1,C;
##supposedly returns the parity check matrix of a 
##binary Hamming code in GF(2^k)
# C:=padic_coeff(2,k);
# G0:=matrix(C):
# G1:=transpose(G0);
# RETURN(evalm(G1)):
#end:

hamming_check_binary:=proc(k)
 local G0,G1,G2,C,C0;
#returns the parity check matrix of the 
#binary Hamming code in GF(2)^k
#in Gauss eliminated form
 C0:=padic_coeff(2,k);
 C:=sort(C0,weight_order);
 G0:=matrix(C):
 G1:=transpose(G0);
 G2:=gausselim(G1);
 RETURN(evalm(G2)):
end:

hamming:=proc(k)
#use hamming_check instead
 local G0,G1,G2,C,C0;
#returns the parity check matrix of the 
#binary Hamming code in GF(2)^k
#in Gauss eliminated form
 C0:=padic_coeff(2,k);
 C:=sort(C0,weight_order);
 G0:=matrix(C):
 G1:=transpose(G0);
 G2:=gausselim(G1);
 RETURN(evalm(G2)):
end:

hamming_check:=proc(k,p)
 local G0,G1,G2,C,C0;
#returns the parity check matrix of the 
#p-ary Hamming code in GF(p)^k
#in Gauss eliminated form
 C0:=projective_space(k,p);
 C:=sort(C0,weight_order);
 G0:=matrix(C):
 G1:=transpose(G0);
 G2:=gausselim(G1);
 RETURN(evalm(G2)):
end:

hamming_check_nonstandard:=proc(k,p)
 local G0,G1,G2,C,C0;
#returns the parity check matrix of the 
#binary Hamming code in GF(2)^k
#in Gauss eliminated form
 C0:=projective_space(k,p);
 C:=sort(C0,weight_order);
 G0:=matrix(C):
 G1:=transpose(G0);
# G2:=gausselim(G1);
 RETURN(evalm(G1)):
end:

hamming_generator:=proc(k,p)
local H0,H1,G1,G2,Ik,C,C0;
 #returns the generator matrix of the 
 #binary Hamming code in GF(p)^k
 Ik:= array(identity, 1..((p^k-1)/(p-1)-k),1..((p^k-1)/(p-1)-k));
 H0:=hamming_check(k,p);
 H1:=delcols(H0,1..k);
 G1:=-transpose(H1):
 G2:=subs(p0=p,map(A->modp(A,p0),concat(G1,Ik)));
 RETURN(evalm(map(A->simplify(A),G2))):
end:

#hamming_generator2:=proc(k)
#local H0,H1,G1,G2,G3,Ik,C,C0;
# #returns the generator matrix of the 
# #binary Hamming code in GF(2^k)
# #in Gauss eliminated form
# Ik:= array(identity, 1..(2^k-k-1),1..(2^k-k-1));
# H0:=hamming(k);
# H1:=delcols(H0,1..k);
# G1:=-transpose(H1):
# G2:=subs(p0=2,map(A->modp(A,p0),concat(Ik,G1)));
# G3:=subs(p0=2,map(A->modp(A,p0),gausselim(G2)));
# RETURN(evalm(map(A->simplify(A),G3))):
#end:

RM_generator:=proc(n::integer)
local H,H1,i,I0,O0,G_RM;
H:=hamming(n,2);
O0 := matrix(n,1,[seq(0,i=1..n)]):
H1:=augment(H,O0):
I0 := matrix(1,2^n,[seq(1,i=1..2^n)]):
G_RM:=stackmatrix(I0,H1);
RETURN(evalm(G_RM));
end:

add_lists:=proc(L1::list,L2::list)
 #adds two lists of vectors
 #(represented as lists)
 local x,y,L;
 L:=[];
 for x in L1 do
  for y in L2 do
   L:=[op(L),x+y]:
  od;
 od;
 RETURN(L);
end:

add_many_lists:=proc(L0::list)
 #adds a list of lists of vectors
 #(represented as lists)
 local x,y,L,i,N;
 N:=nops(L0);
 if N=1 then RETURN(L0); fi;
 if N=2 then RETURN(add_lists(L0[1],L0[2])); fi;
 if N>2 then
  L:=[seq(L0[i],i=2..N)]; 
  RETURN(add_lists(L0[1],add_many_lists(L)));
 fi;
end:

code_list:=proc(G::array,p::integer)
 # returns the list of all code words
 # of the code over F_p with generator
 # matrix G
 local L,g,a,i,j,C,rdim;
 C:=[];
 rdim:=linalg[rowdim](G);
 for i from 1 to rdim do
  g[i]:=convert(linalg[row](G,i),list);
 od;
 for i from 1 to rdim do
  L[i]:=[seq(expand(modp(a[i]*g[i],p)),a[i]=0..(p-1))];
 od;
 C:=modp(add_many_lists([seq(L[i],i=1..rdim)]),p);
end:

dimension_code:=proc(G::array,p::integer)
 # returns the dim 
 # of the code over F_p with generator
 # matrix G
local r,dim;
Gaussjord(G,'r') mod p;dim:=r;
RETURN(dim);
end:

dual_code_list:=proc(G::array,p::integer)
  # returns the list of all code words
  # of the DUAL code over F_p of the code
  # with generator matrix G in STANDARD form
  local p0,A,L,g,a,i,j,C,dualdim,T,B,rdim,rdim1,G0,G1,cdim,Ik;
  C:=[];
  rdim:=rowdim(G);
  cdim:=coldim(G);
  Ik:= array(identity, 1..(cdim-rdim),1..(cdim-rdim));
  G0:=submatrix(G,1..rdim,(rdim+1)..cdim);
  T:=-transpose(G0);
#  G1:=subs(p0=p,map(A->modp(A,p0),concat(T,Ik)));
G1:=map(A->modp(A,p),concat(T,Ik));
  rdim1:=rowdim(G1);
#print(evalm(G1),rdim1);
  for i from 1 to rdim1 do
   g[i]:=convert(row(G1,i),list);
  od;
  for i from 1 to rdim1 do
   L[i]:=[seq(expand(modp(a[i]*g[i],p)),a[i]=0..(p-1))];
  od;
  C:=modp(add_many_lists([seq(L[i],i=1..rdim1)]),p);
end:

decode1:=proc(v::list,G::array,p::integer)
 #assumes G is the generator
 #matrix for a code in F_p
 # and that v has only 1 error
 local v_neg,C,x,syndrome_v;
 v_neg:=modp(expand(-v),p);
 C:=code_list(G,p);
 syndrome_v:=modp(add_lists([v_neg],C),p):
 for x in syndrome_v do
  if weight_codeword(x)<2 then
        print(`corrected codeword: `,modp(expand(x+v),p));
  fi;
 od;
end:

decode:=proc(v::list,G::array,p::integer,d::integer)
 #assumes G is the generator
 #matrix for a code in F_p
 # and that v has <= d/2 errors
 local v_neg,C,x,syndrome_v;
 v_neg:=modp(expand(-v),p);
 C:=code_list(G,p);
 syndrome_v:=modp(add_lists([v_neg],C),p):
 for x in syndrome_v do
  if weight_codeword(x)<=d/2 then
        print(`corrected codeword: `,modp(expand(x+v),p));
  fi;
 od;
end:

basis_vector:=proc(i::integer,n::integer)
 local j,v0,v1;
 #returns vector of length n with
 #1 in ith position, 0 elsewhere
 for j from 1 to n do
  if j=i then v0[j]:=1; 
 else v0[j]:=0;
 fi;
 od;
 v1:=[seq(v0[j],j=1..n)];
 RETURN(v1);
end:

decode_hamming_binary:=proc(v::list)
local k,i,colH,n,v0,G,H,e,syndrome_v,v_correct,error_index,colsH;
 #assumes v is a received word and that
 #the code is a binary hamming code
 #returns the decoded word
 n:=nops(v);
 k:=floor(0.001+evalf(log[2](n+1)));
 H:=hamming(k,2);
 for i from 1 to coldim(H) do
  colH[i]:=convert(col(H,i),list):
 od:
 colsH:=[seq(colH[i],i=1..coldim(H))];
 syndrome_v:=matrix_times_vector_modp(H,v,2);
 if member(syndrome_v,colsH) then
  member(syndrome_v,colsH,'error_index'):
  e:=basis_vector(error_index,n);
  v_correct:=modp(v+e,2);
  RETURN(v_correct);
 else RETURN(v);
 fi;
end:

decode_hamming:=proc(v::list,p)
local c,k,i,colH,n,N,v0,j,G,H,e,syndrome_v,v_correct,error_index,colsH;
 #assumes v is a received word and that
 #the code is a p-ary hamming code
 #returns the decoded word
 n:=nops(v);
 k:=floor(evalf(log[p]((p-1)*n+1))+0.001);
 H:=hamming_check(k,p);
#print(v,n,k,H);
 N:=linalg[coldim](H);
 if n<>N then ERROR(`wrong dims`); fi;
 for i from 1 to N do
  colH[i]:=scalar_mult(convert(col(H,i),list),p):
 od:
 colsH:=[seq(colH[j],j=1..N)];
 syndrome_v:=scalar_mult(matrix_times_vector_modp(H,v,p),p);
 if member(syndrome_v,colsH) then
  member(syndrome_v,colsH,'error_index'):
  
c:=which_scalar_mult(convert(linalg[col](H,error_index),list),matrix_times_vector_modp(H,v,p),p);
  e:=expand(c*basis_vector(error_index,n));
#print(error_index,e,syndrome_v);
  v_correct:=modp(v-e,p);
  RETURN(v_correct);
 else RETURN(v);
 fi;
end:

hamming_weight:=proc(x,y,r)
local sigma1,sigma22,sigma4,n,w;
 n:=2^r-1;
 sigma1:=x+y;
 sigma22:=x^2*y^2;
 sigma4:=x^4+y^4;
 w:=(n+1)^(-1)*(sigma1^n+n*sigma1^((n-1)/2)*(x-y)^((n+1)/2));
 RETURN(expand(w));
end:

isequivalent:=proc(G1,G2)
 #if G1,G2 are the gen mats of linear
 #*binary* codes then this returns true *and* an
 #associated col permutation if
 #the codes are equiv and false otherwise
 ######## SLOW!
 local N,i,j,G,C,C1,C2,S_n,x;
 print(`Uses group21_v5.mpl. Please be patient for output ...`);
 C1:=code_list(G1,2):
 C1:=sort(C1,weight_lexorder);
 N:=linalg[coldim](G1);
 if N<>linalg[coldim](G2) then RETURN(`false`); fi;
 S_n:=gp_elements(N,{[[1,2]],[[seq(i,i=1..N)]]});
 for x in S_n do
 G:=permcol_matrix(G2,x);
 C:=code_list(G,2):
 C:=sort(C,weight_lexorder);
 if evalb(C=C1) then RETURN([`true`,x]); fi;
 od;
 RETURN(`false`);
end:

isequivalent2:=proc(G1,G2)
 #if G1,G2 are the gen mats of linear
 #binary codes then this returns true *and* an
 #associated col permutation if
 #the codes are equiv and false otherwise
 ######## slow but about 5-8 times faster than the one above
 local N,i,j,G,C,G10,G20,S_n,x;
 print(`Uses group21_v5.mpl. Please be patient for output ...`);
 G10:=Gausselim(G1) mod 2;
 N:=linalg[coldim](G1);
 if N<>linalg[coldim](G2) then RETURN(`false`); fi;
 S_n:=gp_elements(N,{[[1,2]],[[seq(i,i=1..N)]]});
 for x in S_n do
  G:=permcol_matrix(G2,x);
  G20:=Gausselim(G) mod 2;
  if linalg[equal](G10,G20) then RETURN([`true`,x]); fi;
 od;
 RETURN(`false`);
end:

automorphismgroup:=proc(G)
 #if G is the gen mat of a linear
 #*binary* code then this returns the
 #list of col permutations for which
 #the codes are equal
 local N,i,j,C1,C2,G10,G20,S_n,x,L;
 print(`Uses group21_v5.mpl. Please be patient for output ...`);
 C1:=code_list(G,2):
 C1:=sort(C1,weight_order);
 N:=linalg[coldim](G);
 L:=[];
 S_n:=gp_elements(N,{[[1,2]],[[seq(i,i=1..N)]]});
 for x in S_n do
  G20:=permcol_matrix(G,x);
  C2:=code_list(G20,2):
  C2:=sort(C2,weight_order);
  if evalb(C2=C1)  then L:=[op(L),x];; fi;
 od;
 if nops(L)>0 then 
  RETURN(L); else RETURN(`false`);
 fi;
end:

print(`         Basic coding theory commands:`);
print(` matrix_times_vector_modp, matrix_modp, golay12, golay11, `);
print(` golay23, golay24, golay12_check, golay11_check, `);
print(` golay23_check, golay24_check, check_matrix, generator_matrix,`);
print(` iscyclic, weight_codeword, weight_enumerator_vector,`);
print(` hamming_generator, hamming_weight, code_list, dual_code_list,`);
print(` decode, decode_hamming, RM_generator, weight_order,`);
print(` weight_lexorder, isequivalent, isequivalent2, `);
print(` automorphismgroup ... loaded. Last update 10-26-2001.`);