# Joseph Koutsoutis, 05-01-2024, proj2

with(NumberTheory):
with(combinat):

Help:=proc(): 
 print(` MacWm(f,n,k,q,z), GetDualGenerator(n,q,x,g), HD(u,v), CtoL(n,x,g), WeightEnumeratorBCH(n,d,q,t) `):
 print(` GetEnumerators(k,K,B), AveAndMoms(f,x,N), EncodeBCH(d,q,x), HillBCHMatrix(n,d,q) `):
 print(` HillBCHPCM(n,d,q), HillWeightEnumeratorBCH(n,d,q,t) `):
end:

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:


GetDualGenerator := proc(n,q,x,g) local i,h,r:
 h := Quo(x^n-1,g,x) mod q:
 r := degree(h,x):
 h := h * coeff(h,x,0)^(-1) mod q:
 h := add(coeff(h,x,i)* x^(r-i), i=0..r):
end:


# This gets weight enumerators for the more standard definition of BCH codes (https://en.wikipedia.org/wiki/BCH_code).
# This currently only works for GF(p) where p is prime, but this will eventually be usable for prime powers.
WeightEnumeratorBCH := proc(n,d,q,t) local g, alpha,i,used_dual,x_len,M,zero,ptr,go_back,word,T,x,dist,enumerator,val:
 alpha := PrimitiveRoot(q):


 g := 1: # g is our generator polynomial
 for i from 1 to d-1 do:
  g := g * (x - alpha&^i) mod q:
 od:
 
 x_len := n+1-d:
 if x_len <= n - x_len then:
  M := CtoL(n,x,g):
  used_dual := false:
 else:
  g := GetDualGenerator(n,q,x,g):
  M := CtoL(n,x,g):
  x_len := n - x_len:
  used_dual := true:
 fi:
 
 x := [seq(0,i=1..x_len)]:
 zero := [seq(0,1..n)]:
 word := zero:
 ptr := x_len:
 go_back := false:
 T[0] := 1:
 while ptr <> 0 do:
  if not go_back then:
    x[ptr] += 1:
    word := (word + M[ptr]) mod q:
    if x[ptr] = q then:
     x[ptr] := 0:
     ptr -= 1:
     go_back := true:  #gc():
    else:
     dist := HD(word,zero):
     if not type(T[dist], integer) then T[dist] := 0: fi:
     T[dist] += 1:
    fi:
  else:
   x[ptr] += 1:
   word := (word + M[ptr]) mod q:
   if x[ptr] = q then:
    x[ptr] := 0:
    ptr -= 1:
   else:
    go_back := false:
    ptr := x_len:
    dist := HD(word,zero):
    if not type(T[dist], integer) then T[dist] := 0: fi:
    T[dist] += 1:
   fi:
  fi:
 od:
 
 enumerator := 0:
 for i,val in eval(T) do:
  enumerator += val * t^i:
 od:

 if used_dual then: 
  return(simplify(expand(MacWm(enumerator,n,x_len,q,t)))):
 else 
  return(enumerator): 
 fi:
end:

EncodeBCH := proc(d,q,x) local S,eqs,i,j,n,y:
 n := nops(x) + d - 1:
 eqs := {seq(add(i^j * x[i], i=1..nops(x)) + add(i^j * y[i], i=(nops(x)+1)..n), j=0..(d-2))}:
 S := msolve(eqs, q):
 subs(S, [op(x), seq(y[i], i=(nops(x)+1)..n)]):
end:

HillBCHMatrix := proc(n,d,q) local x_len,x,i,M:
 x_len := n+1-d:
 M := [seq([], i=1..x_len)]:
 for i from 1 to x_len do:
  x := [0$(i-1), 1, 0$(x_len-i)]:
  M[i] := EncodeBCH(d,q,x):
 od:
 M:
end:

HillBCHPCM := proc(n,d,q) local M,i,j:
 [seq([seq(j^i, j=1..n)], i=0..d-2)]:
end:

HillWeightEnumeratorBCH := proc(n,d,q,t) local i,used_dual,x_len,M,zero,ptr,go_back,word,T,x,dist,enumerator,val:
 
 x_len := n+1-d:
 if x_len <= n - x_len then:
  M := HillBCHMatrix(n,d,q):
  used_dual := false:
 else:
  M := HillBCHPCM(n,d,q):
  x_len := n - x_len:
  used_dual := true:
 fi:
 
 x := [seq(0,i=1..x_len)]:
 zero := [seq(0,1..n)]:
 word := zero:
 ptr := x_len:
 go_back := false:
 T[0] := 1:
 while ptr <> 0 do:
  if not go_back then:
    x[ptr] += 1:
    word := (word + M[ptr]) mod q:
    if x[ptr] = q then:
     x[ptr] := 0:
     ptr -= 1:
     go_back := true:  #gc():
    else:
     dist := HD(word,zero):
     if not type(T[dist], integer) then T[dist] := 0: fi:
     T[dist] += 1:
    fi:
  else:
   x[ptr] += 1:
   word := (word + M[ptr]) mod q:
   if x[ptr] = q then:
    x[ptr] := 0:
    ptr -= 1:
   else:
    go_back := false:
    ptr := x_len:
    dist := HD(word,zero):
    if not type(T[dist], integer) then T[dist] := 0: fi:
    T[dist] += 1:
   fi:
  fi:
 od:
 
 enumerator := 0:
 for i,val in eval(T) do:
  enumerator += val * t^i:
 od:

 if used_dual then: 
  return(simplify(expand(MacWm(enumerator,n,x_len,q,t)))):
 else 
  return(enumerator): 
 fi:
end:


#MacWm(f,n,k,q,z): The MacWilliams transform takes the weight-enumerator of a linear [n,k]-code over GF(q) and outputs
#the weight-enumerator of the dual code (that is [n,n-k])
MacWm:=proc(f,n,k,q,z)
subs(z=(1-z)/(1+(q-1)*z),f)*(1+(q-1)*z)^n/q^k:
end:


CtoL:=proc(n,x,g) local r,L,i:
r:=degree(g,x):
L:=[seq(coeff(g,x,i),i=0..r)]:
[seq([0$i,op(L), 0$(n-r-1-i)],i=0..n-r-1)]:
end:

#AveAndMoms(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list whose first entry is the average 
#(alias expectation, alias mean)
#followed by the variance, the third moment (about the mean) and
#so on, until the N-th moment (about the mean).
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#AveAndMoms(((1+x)/2)^100,x,4);

AveAndMoms:=proc(f,x,N) local mu,F,memu1,gu,i:
mu:=simplify(subs(x=1,f)):

if mu=0 then
print(f, `is neither a prob. generating function nor can it be made so`):
RETURN(FAIL):
fi:

F:=f/mu:


memu1:=simplify(subs(x=1,diff(F,x))):

gu:=[memu1]:

F:=F/x^memu1:

F:=x*diff(F,x):
for i from 2 to N do
 F:=x*diff(F,x):
 gu:=[op(gu),simplify(subs(x=1,F))]:
od:

gu:

end:


GetEnumerators := proc(k,K,B,print_to_file,use_hill) local p,d,n,enumerator,avg_and_var:
 p := k:
 while p <= K do:
  for d from 3 to p-1 by 2 do:
   for n from d to p-1 do:
    if p^(n-d+1) <= B or p^(d-1) <= B then: 
     if use_hill then:
      enumerator := HillWeightEnumeratorBCH(n,d,p,t):
     else:
      enumerator := WeightEnumeratorBCH(n,d,p,t):
     fi:
     avg_and_var := AveAndMoms(enumerator,t,2): 
     
     if print_to_file then:
      fprintf(`database.txt`, `d=%d, n=%d, p=%d, mean=%e, var=%e, stdev=%e, enumerator=%a\n`,  
                       d, n, p, avg_and_var[1], avg_and_var[2], sqrt(avg_and_var[2]), enumerator):
      fflush(`database.txt`):
     else:
      printf(`d=%d, n=%d, p=%d, mean=%e, var=%e, stdev=%e, enumerator=%a\n`,  
                       d, n, p, avg_and_var[1], avg_and_var[2], sqrt(avg_and_var[2]), enumerator):
     fi:
    fi:
   od:
  od:
  p := nextprime(p):
 od:
 return():
end: