# OK to post

# PART 1: Hampcm() PROCEDURE:
# -----------------------
Hampcm:=proc(r,q) local A,a,i,j,H,skip,ncols:
	A:=Fqn(q,r) minus {[0$r], seq([0$i, 1, 0$(r-i-1)], i=0..r-1)};

	ncols:=0;
	for a in A do:
		skip:=false;

		# check if a has leading nonzero entry 1:
		for i from 1 to r do:
			if a[i] <> 0 then
				skip:=evalb(a[i] <> 1);
				break;
			fi;
		od;

		if skip = false then
			ncols:=ncols+1;
			for i from 1 to r do:
				H[i, ncols]:=a[i];
			od;
		fi;
	od;

	for i from 1 to r do:
		for j from 1 to r do:
			if i = j then
				H[i, j+ncols]:=1;
			else
				H[i, j+ncols]:=0;
			fi;
		od;
	od;

	ncols:=ncols+r;
	return [seq([seq(H[i,j], j=1..ncols)], i=1..r)];
end;

# PART 2: Ham() PROCEDURE:
# -----------------------
Ham:=proc(r,q):
	return antiPCM(q,Hampcm(r,q));
end;

# PART 3: DecodeHamming() PROCEDURE:
# -----------------------
DecodeHamming:=proc(q,r,v) local H,syn,b,i,j,colequals,u:
	H:=Hampcm(r,q);
	syn:=Syn(q,H,v);

	# deal with all 0 vector, since it wont work 
	if v = [0$r] then
		return v;
	fi;

	for i from 1 to nops(syn) while syn[i] = 0 do od;

	b:=syn[i];
	syn:=syn * (b&^(-1) mod q) mod q;

	for j from 1 to nops(H[1]) do:
		colequals:=true;

		for i from 1 to  nops(H) do:
			if H[i][j] <> syn[i] then
				colequals:=false;
				break;
			fi;
		od;

		if colequals = true then
			break;
		fi;
	od;

	u:=v;
	u[j]:=u[j]-b mod q;
	return u;
end;

############################################################## OLD CODE ##############################################################

Fqn:=proc(q,n) local S,a,v:
option remember:
if n=0 then
 RETURN({[]}):
fi:

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

{seq(seq([op(v),a],a=0..q-1), v in S)}:

end:

antiPCM:=proc(q,H) local A,B,M,n,rows,k,co,tempCol:
 rows:=nops(H): # this is n-k
 n:=nops(H[1]): # number of H columns
 k:= n-rows: 

 B:=extractMCols(1,k,H):                         # this is -A^T, we wish to return [I_k | A].
 tempCol := extractMCols(1,1,B):                 # take the first column, negate it
 A:= [[seq(op(tempCol[k]), k=1..nops(tempCol))]]:  # make the first column into a row
 
 for co from 2 to k do:
   tempCol := extractMCols(co,co,B): # take the (co)-th column 
   A:= [op(A), [seq(op(tempCol[k]), k=1..nops(tempCol))]]: # add the co-th column as a row
 od:

 A:= -A mod q: # we took the transpose, now we negate.

 co:=1:
 M:=A:
 for co from 1 to k do: #we will append the identity matrix row by row
  M[co] := [0$(co-1), 1, 0$(k-co),op(A[co])]:
 od:

 M:
end:

#DP(q,u,v):
DP:=proc(q,u,v) local i,n:
n:=nops(u):
add(u[i]*v[i],i=1..n) mod q:
end:
 
#Syn(q,H,y): The syndrom of the transmitted vector y if the PCM is H
#(a row vector of length n-k) 
Syn:=proc(q,H,y) local i:
[seq(DP(q,y,H[i]),i=1..nops(H))]:
end: