# OK to post

# PART 1: EXERCISES FROM HILL'S BOOK:
# -----------------------

# - exercise 2.1 -
# (6,2,6) code:
# 000000
# 111111

# (3,8,1) code:
# 000
# 001
# 010
# 011
# 100
# 101
# 110
# 111

# (4,8,2) code (just added a parity check to previous code):
# 0000
# 0011
# 0101
# 0110
# 1001
# 1010
# 1100
# 1111

# A (5,3,4) code is not possible. We can assume that 00000 is a codeword WLOG,
# so there must then also exist a code with at least 4 1s, say 11110 WLOG. Every other
# code must also have at least 4 1s to remain distance 4 from the 0 codeword, but it is not
# possible to have 4 1s but still have distance 4 from 11110

# A (8,30,3) code is not possible. The sphere-packing bound for such a code is 30 * (1 + (8 choose 1)) = 30 * 9 = 270 < 256 = 2^8.
# It is clear that this bound is not satisfied, so such a code cannot exist.

# - exercise 2.2 -
# Partition the codewords in the given binary (n,M,d) code into 2 disjoint sets, those ending with 0 and those ending with 1.
# 1 of these sets must have at least M/2 elements, simply take this set and remove the last digit. This leaves the minimum distance
# unchanged since all codewords ended with the same digit, giving a (n-1, M', d) code where M'<=M/2.

# this also gives us that 2A(n-1,d) >= A(n,d)

# - exercise 2.3 -
# A(n+1,d+1)=A(n,d) if d is odd, so A(3,2)=A(2,1). A(2,1) is just the number of elements in (F_q)^2=q^2

# - exercise 2.4 -
# Given an element, v, from (F_2)^(n-1), there are 2 possible ways to extend it to an element of (F_2)^n,
# adding a 1 to the end or a 0 to the end (adding a digit anywhere else in the element is just equivalent to
# picking a different element of (F_2)^(n-1) and extending that). Adding both (v,1) and (v,0) for every element
# of (F_2)^(n-1) yields (F_2)^n. To generate E_n we instead only pick 1 possible extension. Only 1 of these possible extensions is in 
# E_n since either w(v) is odd, in which case add a 1 so w(v,1) is even or w(v) is even in which case add a 0
# so that w(v,0) is still even. This is equivalent to adding a parity check to each element of (F_2)^(n-1).

# E_n has 2^(n-1) members since |E_n|=|(F_2)^(n-1)| as shown earlier. Each element in E_n has minimum distance 2 to all
# other elements since adding a parity check to any (n,M,d) code yields a (n+1,M,d+1) code if d is odd.

# PART 3: CARDINALITY OF BINARY CODES:
# -----------------------

# all vectors in (F_q)^n
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;

#Nei(q,c): all the neighbors of the vector c in Fqn(q,n)
Nei:=proc(q,c) local n, i, a:
	n:=nops(c);
	{seq(seq([op(c[1..i-1]), a, op(c[i+1..n])], a=0..q-1), i=1..n)}
end;

#SP(q,c,t): the set of all vectors in Fqn(q,n) whose distance is <=t from c
SP:=proc(q,c,t) local S, s, i:
	S:={c};
	for i from 1 to t do:
		S:=S union {seq(op(Nei(q,s)), s in S)};
	od;

	return S;
end;

# greater random code
GRC:=proc(q,n,d) local S, A, v:
	A:=Fqn(q,n);
	S:={};

	while A <> {} do:
		v:=A[1];
		S:=S union {v};

		A:=A minus SP(q, v, d-1);
	od;

	return S;
end;

# sphere packing bound:
SPB:=proc(q,n,t) local i:
	trunc(q^n / add(binomial(n,i)*(q-1)^i, i=0..t));
end;

for n from 5 to 20 do:
    for d from 3 to 7 by 2 do:
        s:=nops(GRC(2,n,d));
        spb:=SPB(2,n,floor((d-1)/2));

        print(cat("n=", n, " d=", d, " - nops(GRC(2,n,d))=", s, " SPB(2,n,t)=", spb, " difference=", spb-s));
    od;
od;

# PART 4: CARDINALITY OF TERNARY CODES:
# -----------------------
for n from 5 to 10 do:
    for d from 3 to 7 by 2 do:
        s:=nops(GRC(3,n,d));
        spb:=SPB(3,n,floor((d-1)/2));

        print(cat("n=", n, " d=", d, " - nops(GRC(3,n,d))=", s, " SPB(3,n,t)=", spb, " difference=", spb-s));
    od;
od;

# PART 5: CADINALITY OF QUINARY CODES:
# -----------------------
for n from 5 to 7 do:
    for d from 3 to 7 by 2 do:
        s:=nops(GRC(5,n,d));
        spb:=SPB(5,n,floor((d-1)/2));

        print(cat("n=", n, " d=", d, " - nops(GRC(5,n,d))=", s, " SPB(5,n,t)=", spb, " difference=", spb-s));
    od;
od;

# PART 6: ParaBD() PROCEDURE:
# -----------------------

# returns the tuple (b,v,r,k,lambda) of a given block design, or FAIL if it is not a valid block design
ParaBD:=proc(BD, v) local b, T, r, rt, k, kt, lambda, lambdat, i, j, i1:
    b:=nops(BD);

    # generate incidence matrix:
    for i from 1 to v do:
        for j from 1 to b do:
            if member(i, BD[j]) then
                T[i,j]:=1;
            else
                T[i,j]:=0;
            fi;
        od;
    od;

    # check row sums are all equal
    r:=sum(T[1,i1], i1=1..b);
    for i from 1 to v do:
        rt:=sum(T[i,i1], i1=1..b); # r-temp

        if rt <> r then
            return FAIL;
        fi;
    od;

    # check column sums are all equal
    k:=sum(T[i1,1], i1=1..v);
    for i from 1 to b do:
        kt:=sum(T[i1,1], i1=1..v); # k-temp

        if kt <> k then
            return FAIL;
        fi;
    od;

    # check dot products of rows are all equal
    lambda:=sum(T[1,i1]*T[2,i1],i1=1..v);
    for i from 1 to v do:
        for j from i+1 to v do:
            lambdat:=sum(T[i,i1]*T[j,i1],i1=1..v); # lambda-temp

            if lambdat <> lambda then
                return FAIL;
            fi;
        od;
    od;

    return (b,v,r,k,lambda);
end;

BDfano:=proc():
	return {{1,2,4},{2,3,5},{3,4,6},{4,5,7},{5,6,1},{6,7,2},{7,1,3}};
end:

BDex212:=proc():
	return {
		{1,3,4,5,9},
		{2,4,5,6,10},
		{3,5,6,7,11},
		{1,4,6,7,8},
		{2,5,7,8,9},
		{3,6,8,9,10},
		{4,7,9,10,11},
		{1,5,8,10,11},
		{1,2,6,9,11},
		{1,2,3,7,10},
		{2,3,4,8,11}
	};
end:

# works!
ParaBD(BDfano(),7);
ParaBD(BDex212(),11);

# PART 7: PlotkinBound() PROCEDURE:
# -----------------------

# returns an upper bound for the size of a code given that n<2d
PlotkinBound:=proc(n,d):
    return 2 * floor(d / (2*d - n));
end;

for d from 2 to 20 do:
    for n from 2 to 2*d-1 do:
        pb:=PlotkinBound(n,d);
        spb:=SPB(2,n,floor((d-1)/2));

        if pb < spb then
            print(cat("Plotkin bound is smaller for n=",n,", d=", d));
        elif spb < pb then
            print(cat("Sphere packing bound is smaller for n=",n,", d=", d));
        else
            print(cat("Both bounds are equal for n=",n,", d=", d));
        fi;
    od;
od;

# it seems that the plotkin bound is almost always smaller than the sphere packing bound