#OK to post homework
#Ramesh Balaji,2/11/2024,hw6

#Alphabet {0,1,...q-1}, Fqn(q,n): {0,1,...,q-1}^n
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:

#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:

#MinD(C): The minimal (Hamming) distance of the code C
MinD:=proc(C) local i,j:
min( seq(seq(HD(C[i],C[j]),j=i+1..nops(C)), i=1..nops(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:

#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):

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

end:

# Part 1: 5.1, 5.2, 5.3, 5.5

# 5.1 Binary code (11, 24, 5) is not a linear code. Note that any linear
# code can be represented as a qary [n,k,d] code and every qary [n,k,d]
# code is a (n, q^k, d) code. Since q=2, we expect 24 = q^k but this is
# not possible as 24 is not a power of 2.
#
# 5.2 E_n's parameters are [n, n-1, 2]
# The minimum distance is 2 because the minimum Hamming distance between
# two even weight vectors must be at least 2, as there is only one vector with
# Hamming weight 0.
# Note that the total number of even weight vectors is equal to
# sum(binomial(n, k), k=even) = 2^(n-1). Then E_n is an (n, 2^(n-1), 2) code,
# so thus k=log_q(2^(n-1)) = n-1. A generating matrix for E_n in standard form
# is [I_(n-1)|1], where 1 is the 1 column vector of length n-1 (this ensures
# that the min HD between two vectors is always even).
#
#
# 5.3 We aim to prove that C is a lienar code. Then
# take arbitrary x, y in C and note xH^T = 0 and yH^T = 0.
# Note that 0 = 0 + 0 = xH^T + yH^T = (x+y)H^t. So x + y is in C.
# Take arbitrary c in GF(q) and x in C. Note that xH^T = 0.
# Note that c*0 = c*xH^T = (cx)H^T. So cx is in C. Thus C is a
# vector subspace of V(n,q) so it is a linear code.
#
# 5.4
# i. Let C' be the code C with a parity check added to every codeword. Now
# we aim to show that C' is linear code.
# Choose arbitrary x, y in C'. Then then there are 3 cases,
# either x and y have different parity, or they are both odd,
# or they are both even.
#
# Case 1: They both have different parity. Assume WLOG that x has
# an odd weight and y has an even weight. Then the parity check bit
# of x + y will be 1. Take x' to be x without the parity check bit
# and y' to be y without the parity check bit.
# By lemma 2.5 and 2.6, we have that w(x'+y')
# = d(x',y') = w(x') + w(y') - 2*w(x' intersect y')
# Since w(x') is odd and w(y') is even, w(x') + w(y') is
# odd. Note that 2*w(x' intersect y') is always even, so
# w(x') + w(y') - 2*w(x' intersect y') as an odd minus an even
# is always odd. Thus the parity of x' + y' is odd, so
# x' + y' with its parity bit added (which must be 1) will
# always be in C'.
#
# Case 2: They both have even parity. Then the parity bit
# must be zero as x and y's parity bits will both be zero.
# Take x' to be x without the parity check bit
# and y' to be y without the parity check bit.
# By lemma 2.5 and 2.6, we have that w(x'+y')
# = d(x',y') = w(x') + w(y') - 2*w(x intersect y)
# Since w(x'), w(y') and 2*w(x' intersect y') are all even,
# we know that w(x'+y') is even.
# Thus the parity of x' + y' is even, so
# x' + y' with its parity bit added (which must be 0) will
# always be in C'.
#
# Case 3: They both have odd parity. Then the parity bit
# must be zero as x and y's parity bits will both be zero.
# Take x' to be x without the parity check bit
# and y' to be y without the parity check bit.
# By lemma 2.5 and 2.6, we have that w(x'+y')
# = d(x',y') = w(x') + w(y') - 2*w(x intersect y)
# Since w(x'), w(y') are all odd, and 2*w(x' intersect y') is even
# we know that d(x',y') is odd, so
# we know that w(x'+y') is odd.
# Thus the parity of x' + y' is odd, so
# x' + y' with its parity bit added (which must be 1) will
# always be in C'.
#
# ii. The generator matrix is
# [ 1 0 0 0 1 1 1 0 ]
# [ 0 1 0 0 1 1 0 1 ]
# [ 0 0 1 0 1 0 1 1 ]
# [ 0 0 0 1 0 1 1 1 ]

# Check for ii:
r := {[1, 0, 0, 0, 1, 1, 1, 0],
      [0, 1, 0, 0, 1, 1, 0, 1],
      [0, 0, 1, 0, 1, 0, 1, 1],
      [0, 0, 0, 1, 0, 1, 1, 1]};

MinD(r);
MinD(LtoC(2, r));

r := 'r':

# Part 3

# return a random subset of k elements in the set st
randsubset := proc(st, k) local i, picker:
    picker := rand(1..nops(st));
    return {seq(st[picker()], i=1..k)};
end:


SearchLC := proc(q,n,k,K) local fqn, Mset, i, Mminw, ss:
    Mset := [];
    Mminw := [];
    fqn := Fqn(q,n);
    for i from 1 to K do:
        ss :=  randsubset(fqn, k);
        Mset := [op(Mset), ss];
        Mminw := [op(Mminw), MinW(q, ss)];
    od:

    return Mset[min[index](Mminw)];
end:

for q in [2, 3, 5] do
    for n from 4 to 10 do
        for k from 3 to 5 do
            print("q="||q||",n="||n||",k="||k);
            print(SearchLC(q, n, k, 1000));
        od:
    od:
od:

q := 'q':
n := 'n':
k := 'k':

# Part 4 (I tried to write this in PascalCase, for whatever reason)
# Honestly, it is not clear if my procedure yields the correct matrix.
# But I think it might?

ReductionR1 := proc(Mat, q, Row1, Row2) local MatOut, RowTemp:
    MatOut := Mat;

    # We could probably set Row1 = Row1 + Row2 and then Row2 = Row1 - Row2, and
    # then Row1 = Row1 - Row2, but I'm not sure if it would make a huge
    # difference in Maple.

    RowTemp := MatOut[Row1];

    MatOut[Row1] := MatOut[Row2];
    MatOut[Row2] := RowTemp;

    return MatOut;
end:

ReductionR2 := proc(Mat, q, Scalar, Row) local MatOut, ICol:
    MatOut := Mat;

    MatOut[Row] := [seq( (MatOut[Row][ICol] * Scalar) mod q , ICol=1..nops(MatOut[Row]) )];

    return MatOut;
end:

# We perform operation Scalar1*Row1 + Scalar2*Row2 -> Row2
ReductionR3 := proc(Mat, q, Scalar1, Row1, Scalar2, Row2) local MatOut, ICol:
    MatOut := Mat;

    MatOut[Row2] := [seq( (( (Scalar1*MatOut[Row1][ICol]) mod q)
                           + (Scalar2*MatOut[Row2][ICol] mod q) ) mod q,
                         ICol=1..nops(MatOut[Row1]) )];

    return MatOut;
end:

ReductionC1 := proc(Mat, q, Col1, Col2) local MatOut, NumRows, IRow, TempCell:
    MatOut := Mat;

    NumRows := nops(MatOut);

    for IRow from 1 to NumRows do:
        TempCell := MatOut[IRow][Col1];

        MatOut[IRow][Col1] := MatOut[IRow][Col2];
        MatOut[IRow][Col2] := TempCell;
    od:

    return MatOut;
end:

ReductionC2 := proc(Mat, q, Scalar, Col) local MatOut, NumRows, IRow, TempCell:
    MatOut := Mat;

    NumRows := nops(MatOut);

    for IRow from 1 to NumRows do:
        MatOut[IRow][Col] := (MatOut[IRow][Col] * Scalar) mod q;
    od:

    return MatOut;
end:

SGM := proc(Mat, q) local PivotRow, RowLen, MatOut, IRow, ICol:
    MatOut := Mat;

    PivotRow := 1;
    RowLen := nops(MatOut);

    # Downwards elimination
    for PivotRow from 1 to RowLen do:
        # If the location where the pivot should be is zero,
        # find a row where the location of the pivot in that row
        # is nonzero, and swap the two rows.
        if MatOut[PivotRow][PivotRow] = 0 then:
            for IRow from PivotRow + 1 to nops(MatOut) do:
                if MatOut[IRow][PivotRow] <> 0 then:
                    MatOut := ReductionR1(MatOut, q, PivotRow, IRow);
        print(Matrix(MatOut));
                    break;
                end:
            od:
        end:
        # If after all that we couldn't make the pivot value nonzero, then we
        # need to do some column swapping.
        if MatOut[PivotRow][PivotRow] = 0 then:
           for ICol from PivotRow to nops(MatOut[1]) do:
               if MatOut[PivotRow][ICol] <> 0 then:
                  print("Swapping Col1=" || PivotRow || " Col2=" || ICol);
                  MatOut := ReductionC1(MatOut, q, PivotRow, ICol);
                  break;
        print(Matrix(MatOut));
               end:
           end:
        end:
        # Make the pivot location to 1
        MatOut := ReductionR2(MatOut, q, (1/MatOut[PivotRow][PivotRow]) mod q, PivotRow);
        print(Matrix(MatOut));

        # Eliminate every nonzero entry in the column of the pivot for all rows
        # beyond the current pivot row.
        for IRow from PivotRow + 1 to RowLen do:
            if MatOut[IRow][PivotRow] <> 0 then:
                MatOut := ReductionR3(MatOut, q, (-1*MatOut[IRow][PivotRow]) mod q, PivotRow, 1, IRow);
        print(Matrix(MatOut));
            end:
        od:
        print(Matrix(MatOut));
    od:

    # Upwards elimination
    for PivotRow from RowLen to 1 by -1 do:
        # Eliminate every nonzero entry in the column of the pivot for all rows
        # above the current pivot row.
        for IRow from PivotRow - 1 to 1 by -1 do:
            if MatOut[IRow][PivotRow] <> 0 then:
                MatOut := ReductionR3(MatOut, q, (-1*MatOut[IRow][PivotRow]) mod q, PivotRow, 1, IRow);
        print(Matrix(MatOut));
            end:
        od:
        # print(Matrix(MatOut));
    end:


    return MatOut:
end: