# OK to post

# PART 1: VALENTINE AGES:
# -----------------------

# Alex Varjabedian: 19
# Aurora Hiveley: 23
# Gloria Liu: 21
# Himanshu Chandra: 22
# Joseph Koutsoutis: 20
# Kaylee Weatherspoon: 22
# Lucy Martinez: 25
# Natalya Ter-Saakov: 24
# Pablo Blanco: NOT CONSISTENT
# Ramesh Balaji: 19
# Ryan Badi: 20
# Shaurya Baranwal: 19

# PART 2: antiPCM() PROCEDURE:
# -----------------------
antiPCM:=proc(q,H) local n,k,M,i,j:
    n:=nops(H[1]);
    k:=n-nops(H);

    for i from 1 to n-k do:
        for j from 1 to k do:
            M[j, k+i]:=-H[i][j] mod q;
        od;
    od;

    for i from 1 to k do:
        for j from 1 to k do:
            if i=j then
                M[i,j]:=1;
            else
                M[i,j]:=0;
            fi;
        od;
    od;

    return [seq([seq(M[i,j], j=1..n)], i=1..k)];
end;

# PART 3: DecodeLT() PROCEDURE:
# -----------------------
DecodeLT:=proc(q,M) local T,ST,H,n,i,y,k:
    n:=nops(M[1]);
    k:=nops(M);

    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:
    
    ST:=SynT(q,M)[1];
    H:=PCM(q,M);

    T:=table();
    for y in Fqn(q,n) do:
        T[y]:=y-ST[Syn(q,H,y)] mod q;
    od;

    return op(T);
end;

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

# dot product
DP:=proc(q,u,v) local i,n:
    n:=nops(v);
    return add(u[i]*v[i],i=1..n) mod q;
end;

#start code from Aurora Hiveley
# finds a vector of smallest weight among A (collection of given vectors)
minW := proc(A) local n,v,w,minw,minv: 
n:= nops(A[1]):
minw := HD(A[1],[0$n]):
minv := A[1]: # initialize min weight vectors

for v in A do 
    w:= HD(v,[0$n]):
    if w < minw then
        minw := w :
        minv := v :
    fi:
od:

minv;
end:

#SAah(q,n,M): inputs a basis M of a linear [n,nops(M),d] code outputs Slepian's Standard Array
#as a matrix of vectors containing all the vectors in GF(q)^n (alas Fqn(q,n)) such that
#the first row is an ordering of the members of the actual code (LtoC(q,M)) with
#[0$n] the first entry and the first columns are the coset reprenatives

SAah:=proc(q,n,M) local SL,C,A,a,r1,r,j:
# copied from class
 C:=LtoC(q,M):
 C:=C minus {[0$n]}:

# added
 r1 := [[0$n],op(C)]:
 SL := [r1]:
 A:=Fqn(q,n) minus {op(r1)}: # changed from class

 while A<>{} do
    # find coset representative of min weight
    a := minW(A) :

    # build next row
    r := []:
    for j from 1 to nops(r1) do
        r := [op(r), a + r1[j] mod q]:
    od:
    # add new row to array
    SL := [op(SL), r];
    
    # print(SL);
    # update available vectors
    A := A minus {op(SL[-1])}:

 od:

SL;
end:

# Syn(q,H,y): the syndrome of the received vector y
Syn:=proc(q,H,y) local i:
    [seq(DP(q,y,H[i]),i=1..nops(H))]:
end:

# SynT(q,M): inputs q and a basis (in standard form) M of some linear code over GF(q)
# outputs the syndrome table mapping each possible syndrome to its corresponding 
# coset leader in its standard array
SynT:=proc(q,M) local n,T,A,H,S,i,s:
    if SF(q,M)<>M then
        return FAIL;
    fi;

    n:=nops(M[1]);
    A:=SAah(q,n,M);
    H:=PCM(q,M);
    S:={};

    for i from 1 to nops(A) do:
        s:=Syn(q,H,A[i][1]);
        S:=S union {s};
        T[s]:=A[i][1];
    od;

    return op(T),S;
end;

# 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,i,M1,c:
    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);
    return {seq(seq(c + i*M[k] mod q, i=0..q-1), c in C)};
end;

# hamming distance
HD:=proc(u,v) local i, c:
	# TODO: type checking
	
	if nops(u) <> nops(v) then
		return FAIL;
	fi;

	c:=0;
	for i from 1 to nops(u) do
		if u[i] <> v[i] then
			c:=c+1;
		fi;
	od;

	return c;
end;

# MinW(q,M): The minimal weight of the linear code generated by M over GF(q)
MinW:=proc(q,M) local n,c,C:
    n:=nops(M[1]);
    C:=LtoC(q,M);

    return min(seq(HD(c, [0$n]), c in C minus {[0$n]}));
end;

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

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

# GLC1: tries to add a new number 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;

    return 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 give 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;

    return M;
end;

# PCM(q,M) inputs a non-empty basis describing some linear code (assumed in standard form)
# over GF(q)^n outputs the parity check matrix
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:

    return [seq([seq(H[i,j],j=1..n)],i=1..n-k)]:
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
SF:=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;