#LtoC(q,M): inputs a list of basis vectors for our linear code over GF(q)
#outputs all the codewords (the actual subset of GF(q)^n with q^k elements
LtoC:=proc(q,M) local n,k,C,c,i,M1:
option remember:
k:=nops(M):
n:=nops(M[1]):
if k=1 then
 RETURN({seq(i*M[1] mod q,i=0..q-1) }):
fi:
M1:=M[1..k-1]:
C:=LtoC(q,M1):
{seq(seq(c+i*M[k] mod q,i=0..q-1),c in C)}:
end:


#HD(u,v): The Hamming distance between two words (of the same length)
HD:=proc(u,v) local i,co:
co:=0:
for i from 1 to nops(u) do
  if u[i]<>v[i] then
      co:=co+1:
  fi:
od:
co:
end:

#GLC1(q,M,d): tries to add a new member to the current basis M out of the other vectors
#that still has min. weight d
GLC1:=proc(q,M,d) local n,A,c,M1:
n:=nops(M[1]):
A:=Fqn(q,n) minus LtoC(q,M):
for c in A do
 M1:=[op(M),c]:
 if MinW(q,M1)=d then
   RETURN(M1):
 fi:
od:
FAIL:
end:


#GLC(q,n,d): inputs q and n and d (for min. distance) greedily and randomly tries to
#get as large a linear code given by a basis as possible)
GLC:=proc(q,n,d) local M,M1:
M:=[[1$d,0$(n-d)]]:
M1:=GLC1(q,M,d):
while M1<>FAIL do
 M:=M1:
 M1:=GLC1(q,M,d):
od:
M;
end:


# gets a row from a matrix
GetRow:=proc(M,r) local v,i,n:
    v:=[];
    n:=nops(M[1]);

    for i from 1 to n do:
        v:=[op(v),M[r,i]];
    od;

    return v;
end:

# sets a row in a matrix
SetRow:=proc(M,r,v) local i,M1:
    M1:=copy(M);
    
    for i from 1 to nops(v) do:
        M1[r,i]:=v[i];
    od;

    return M1;
end:

# gets a column from a matrix
GetCol:=proc(M,c) local v,i,n:
    v:=[];
    n:=nops(M);
    
    for i from 1 to n do:
        v:=[op(v),M[i,c]];
    od;

    return v;
end:

# sets a column in a matrix
SetCol:=proc(M,c,v) local i,M1:
    M1:=copy(M);
    
    for i from 1 to nops(v) do:
        M1[i,c]:=v[i];
    od;

    return M1;
end:


# returns matrix M in standard form
SFde:=proc(q,M) local k,n,i,j,S,rj,cj,ri:
    k:=nops(M);
    n:=nops(M[1]);
    S:=copy(M);

    for i from 1 to k do: # algorithm is iterated from 1 to k
        # if S_ii = 0, then we need to perform a swap:
        if S[i,i] = 0 then 
            for j from i+1 to k while S[j,i] = 0 do od; # look for available row

            if j<=k then # swap rows
                rj:=GetRow(S,j);
                S:=SetRow(S,j, GetRow(S,i));
                S:=SetRow(S,i,rj);

            else # look to swap columns

                for j from i+1 to n while S[i,j] = 0 do od; # look for available col

                # swap cols
                cj:=GetCol(S,j);
                S:=SetCol(S,j, GetCol(S,i));
                S:=SetCol(S,i,cj);
            fi;
        fi;

        # scale row to have leading entry 1
        ri:=GetRow(S,i);
        ri:=(ri*(ri[i]&^(-1) mod q)) mod q;
        S:=SetRow(S,i,ri);

        for j from 1 to k do:
            if j <> i then
                rj:=GetRow(S,j);
                rj:=(rj - (rj[i] * ri mod q)) mod q;
                S:=SetRow(S,j,rj);
            fi;
        od;
    od;

    return S;
end:


#PCM(q,M): inputs the mod. q and a non-empty basis (a list of lists)
#n=nops(M[1]) M is k by n matrix and H=PCM(q,M) is an n-k by n matrix
#describing the basis of some linear code over GF(q)^n outputs of dimension k=nops(M)
#the (n-k) by n
PCM:=proc(q,M) local k,n,i,j,H:
option remember:
k:=nops(M):
n:=nops(M[1]):

if [seq([op(1..k,M[i])],i=1..k)]<>[seq([0$(i-1),1,0$(k-i)],i=1..k)] then
 print(`Not standard form , please use SFde(q,M) first `):
 RETURN(FAIL):
fi:

for i from 1 to n-k do
 for j from 1 to k do
  H[i,j]:=-M[j][i+k]  mod q:
 od:
 for j from k+1 to n do
  if j-k=i then
   H[i,j]:=1:
  else
   H[i,j]:=0:
  fi:
od:
od:

[seq([seq(H[i,j],j=1..n)],i=1..n-k)]:
end:
#____________________________________


#HW8#1

WtEnumerator:=proc(q,M,t) local C,v:
C:=LtoC(q,M):
add(t^HD([0$nops(v)], v), v in C):
end:


#HW8#2

#Relation between WtEnumerator(q,M,t), WtEnumerator(q,H,t)

M:=GLC(2,5,3);
M:=SFde(3, M);
H:=PCM(3, M);


for i from 1 to 10 do:
print(WtEnumerator(i,H,3), WtEnumerator(i,M,3));
od:

with(NumberTheory):

for i from 1 to 10 do:
print(PrimeFactors(WtEnumerator(i,H,3)), PrimeFactors(WtEnumerator(i,M,3)));
od:


for i from 1 to 10 do:
print(Totient(WtEnumerator(i,H,3)), Totient(WtEnumerator(i,M,3)));
od:


for i from 1 to 10 do:
print(evalf(WtEnumerator(i,H,3)/WtEnumerator(i,M,3)));
od:



#Observations:
#When one is even, so is the other, and when one is odd, so is the other.
#I haven't noticed any association in terms of prime factors. Some pairs share a factor, some dont
#Nothing interesting from Totient
#Does not appear to be an interesting pattern in the ratios of the two quantities.