# OK to post homework
# Joseph Koutsoutis, 02-04-2024, Assignment 5

with(ListTools):

#1 done

#Exercise 2.1:
#Example of a (6,2,6)-code: {[0,0,0,0,0,0],[1,1,1,1,1,1]}
#Example of a (3,8,1)-code: (F_2)^3
#Example of a (4,8,2)-code: This can be obtained by adding an overall
#                           parity check to the elements of (F_2)^3
#Example of a (5,3,4)-code: By corollary 2.8, a (5,3,4)-code exists iff
#                           a (4,3,3)-code exists, so we will inspect these codes instead.
#                           Assume for contradiction such a code exists. By lemma 2.3, we
#                           may assume that one of the elements is [0,0,0,0]. Hence this code
#                           contains no elements with one or two 1s. Since this code can only
#                           contain one element with four 1s, at least one of the remaining
#                           two elements must have three 1s. Suppose x,y have three 1s and z has four.
#                           By lemma 2.6: d(x,z)  = 3 + 4 - 6 = 1 
#                                         d(x,y) <= 3 + 3 - 4 = 2.
#                           This is a contradiction, so a (5,3,4)-code does not exist.
#Example of a (8,30,3)-code: By the sphere backing bound: M <= 2^8 / (1 + 8) < 29.
#                            Therefore no (8,30,3)-code exists.

#Exercise 2.2:
#Suppose that there exists a binary (n,M,d)-code. Fix an arbitrary position of the code.
#Since this is a binary code, either 1s or 0s must account for at least half of the entries in
#this position. Without loss of generality assume that there are at least as many 0s as 1s in this
#position. Take all of the members of the (n,M,d)-code with 0s in this fixed position and remove their
#zeros from this fixed position. Let's say that this new code is a (n-1,M',d')-code. Since all of these
#members had a 0 in the removed position and we started with a (n,M,d)-code, this implies that d' >= d.
#Since we took the entries with 0s, this implies that M' >= M/2. Hence we have constructed the desired code.

#Let M = A_2(n,d). By the result we just showed, we may construct a (n-1,ceil(A_2(n,d)/2),d)-code. 
#Therefore A_2(n-1,d) >= A_2(n,d) / 2, implying A_2(n,d) <= 2A_2(n-1,d).

#Exercise 2.3: 
#Suppose q>=2. Suppose we are given a (3,A_q(3,2),2)-code. Observe that if we choose some arbitrary position
#and remove it from all members of our code, we are left with a (2,A_q(3,2),1)-code. There can be at most
#q^2 entries here, implying A_q(3,2) <= q^2.

#Next we show that a (3,q^2,2)-code exists. Start with the (2,q^2,1)-code (F_q)^2. Next take each member and
#add a third entry that is equal to the sum of the first two entries mod q. Next take any two members x,y of 
#this (3,q^2,d)-code that we just constructed. Assume wlog x,y are elements such that their first two entries
#are only off only in the first position. Let's say that the first two entries of x are x_1,x_2 and the first
#two entries of y are y_1,y_2. We find that: x_1 + x_2 = y_1 + y_2 (mod q) iff x_1 = y_1 (mod q).
#Since x_1 != y_1, d(x,y) = 2, and we have constructed a (3,q^2,2)-code.

#Exercise 2.4:
#We want to show that E_n is the code obtained by adding an overall parity check to (F_2)^(n-1).
#Observe that for any x in E_n, the first n-1 positions of x must match one of the elements found in (F_2)^(n-1).
#These first n-1 entries force the nth position of x to have the same parity as the parity of sum of the first n-1 entries,
#so x is found in the code obtained by adding an overall parity check to (F_2)^(n-1).
#The other direction holds since adding a parity check means that the overall parity is even since even + even = even
#and odd + odd = even.
#Hence E_n is (n,2^(n-1),d)-code. Now suppose x,y are members of E_n. Assume wlog that the distance between these two
#elements in the first n-1 entries is 1. Since this is a binary code, the parities of x,y must differ, so d(x,y) = 2.
#Therefore E_n is a (n,2^(n-1),2)-code.

#2 done

#3 

#These 5 functions were written in lectures 4+5:
SPB:=proc(q,n,t) local i:
trunc(q^n/add(binomial(n,i)*(q-1)^i,i=0..t)):
end:

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:

Nei:=proc(q,c) local n,i,a: 
n:=nops(c):
{seq(seq([op(1..i-1,c),a,op(i+1..n,c)], a=0..q-1) , i=1..n)}:
end:

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:
S:
end:

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:
S:
end:

#RunGrcAndCompareToSPB(q,nUpperBound) takes as input integers q and nUpperBound and prints the sphere-packing bound,
#the result of GRC(q,n,d), and the difference between these two values for n in [5,...,nUpperBound] and d in [3,5,7]
RunGrcAndCompareToSPB:=proc(q,nUpperBound)
 local d,n,bound,grc_result:
 for d in [3,5,7] do:
  for n from 5 to nUpperBound do:
   bound := SPB(q,n,(d-1)/2):
   grc_result := nops(GRC(q,n,d)):
   print(sprintf(`q=%d, d=%d, n=%d, SPB=%d, GRC=%d, SPB-GRC=%d`, q, d, n, bound, grc_result, bound-grc_result)):
  od:
 od:
end:

#Here is the result of running RunGrcAndCompareToSPB(2,20) (the third column is the distance from the sphere-packing bound):

# n,d    GRC(q,n,d)    Diff
# 5,3             4       1
# 6,3             8       1
# 7,3            16       0
# 8,3            16      12
# 9,3            32      19
#10,3            64      29
#11,3           128      42
#12,3           256      59
#13,3           512      73
#14,3          1024      68
#15,3          2048       0
#16,3          2048    1807
#17,3          4096    3185
#18,3          8192    5605
#19,3         16384    9830
#20,3         32768   17164
# 5,5             2       0
# 6,5             2       0
# 7,5             2       2
# 8,5             4       2
# 9,5             4       7
#10,5             8      10
#11,5            16      14
#12,5            16      35
#13,5            32      57
#14,5            64      90
#15,5           128     142
#16,5           256     222
#17,5           512     339
#18,5           512    1012
#19,5          1024    1720
#20,5          2048    2921
# 5,7             1       0
# 6,7             1       0
# 7,7             2       0
# 8,7             2       0
# 9,7             2       1
#10,7             2       3
#11,7             4       4
#12,7             4       9
#13,7             8      13
#14,7            16      18
#15,7            32      24
#16,7            32      62
#17,7            64      93
#18,7           128     137
#19,7           256     195
#20,7           512     264

#4

#Here is the result of running RunGrcAndCompareToSPB(3,10) (the third column is the distance from the sphere-packing bound):

# n,d    GRC(q,n,d)    Diff
# 5,3             9      13
# 6,3            24      32
# 7,3            72      73
# 8,3           198     187
# 9,3           519     516
#10,3          1390    1421
# 5,5             3       1
# 6,5             3       6
# 7,5             8      14
# 8,5            18      32
# 9,5            39      81
#10,5            83     210
# 5,7             1       0
# 6,7             1       2
# 7,7             3       2
# 8,7             3       8
# 9,7             3      20
#10,7             9      41

#5

#Here is the result of running RunGrcAndCompareToSPB(5,7) (the third column is the distance from the sphere-packing bound):

# n,d    GRC(q,n,d)    Diff
# 5,3            74      74
# 6,3           265     360
# 7,3          1113    1580
# 5,5             5      12
# 6,5            13      45
# 7,5            39     175
# 5,7             1       2
# 6,7             1       9
# 7,7             5      24

#6

#this function was written in lecture 5:
CV:=proc(S,n) local v,i:
v:=[]:
for i from 1 to n do
 if member(i,S) then
  v:=[op(v),1]:
 else
  v:=[op(v),0]:
 fi:
od:
v:
end:


#ParaBD(BD,v) takes as input a block design BD and v varieties and outputs FAIL if BD
#is not a block design and (b,v,r,k,lambda) otherwise.
#I will assume that v>=2.
ParaBD:=proc(BD,v)
 local incidence,s,i,j,b,r,k,lambda:
 b := nops(BD):
 incidence := Transpose([seq(CV(s,v), s in BD)]):

 r := sum(incidence[1][j], j=1..b):
 for i from 2 to v do:
  if sum(incidence[i][j], j=1..b) <> r then RETURN(FAIL): fi: 
 od:

 i := 'i': #I don't know why but the sum in the next line doesn't work if I don't set this.
 k := sum(incidence[i][1], i=1..v):
 for j from 2 to b do:
  if sum(incidence[i][j], i=1..v) <> k then RETURN(FAIL): fi:
 od:

 lambda := DotProduct(incidence[1], incidence[2]):
 for i from 1 to v do:
  for j from i+1 to v do:
   if DotProduct(incidence[i], incidence[j]) <> lambda then RETURN(FAIL): fi:
  od:
 od:

 (b,v,r,k,lambda):
end:

#7

PlotkinBound:=proc(n,d)
 2 * trunc(d / (2*d - n)):
end:

#CompareSphereAndPlotkin(d_lower,d_upper) prints whether the sphere-packing or Plotkin
#bound is larger for d_lower <= d <= d_upper and 2 <= n < 2d
CompareSphereAndPlotkin:=proc(d_lower,d_upper) local n,d,s,p:
 for d from d_lower to d_upper do:
  for n from 2 to 2*d-1 do:
   s := SPB(2,n,trunc((d-1)/2)):
   p := PlotkinBound(n,d):
   if s > p then print(sprintf(`S, d=%d, n=%d`, d, n)):
   elif p > s then print(sprintf(`P, d=%d, n=%d`, d, n)):
   else print(sprintf(`T, d=%d, n=%d`, d, n)) fi:
  od:
 od:
end:

#After running CompareSphereAndPlotkin(2,20), the plotkin Bound was only higher when n=5 and d=3.
#The actual values were 6 for the Plotkin and 5 for the sphere-packing.

#8 I didn't attempt this challenge.