# Joseph Koutsoutis, 05-01-2024, proj2

print(`This is proj2.txt, A Maple package`):
print(`Version: May 2024  `):
print():
print(`----------------------------------------------------`):

print():
print(`To see the major procedures in this package, type "Help()"`):

print():
print(`For  a list of the SUPPORTING functions type "Help1();" `):
print(`For explanation of each procedure, type "WhatIs(procedure_name)"`):


Help:=proc(): 
 print(` HillWeightEnumeratorBCH(n,d,q,t), HillWeightEnumeratorBCHFast1(n,d,q,t), HillWeightEnumeratorBCHFast2(n,d,q,t) `):
 print(` WeightEnumeratorRandLinearCode(n,d,q), GetEnumerators(k,K,B,print_to_file,use_fast,filename), ProveTheorem06BaseCase() `):
 print(` GetRandLinearCodeMeansAndVars(k,K,B,T,print_to_file,filename), TestConjecture08(k,K,num_moments) `)
end:

Help1:=proc():
 print(` HD(u,v), GetEnumeratorFromM(M,x_len,n,q), R1op(M,i1,i2), C1op(M,j1,j2,k) `):
 print(` Step1(M1,i,g,k,n), Step2(M1, q, g, i), Step3(M1, q, i, k) SF(q,M) `):
 print(` PCM(q,M), EncodeHillBCH(d,q,x), HillBCHMatrix(n,d,q), MacWm(f,n,k,q,z) `):
 print(` AveAndMoms(f,x,N), Alpha(f,x,N), Di(a,b), Rqnk(q,n,k) `):
end:

WhatIs := proc(procedure):
 if procedure = HD then:
  print(`HD(u,v): takes as input 2 vectors u,v of the same length and outputs the Hamming distance`):

 elif procedure = GetEnumeratorFromM then:
  print(`GetEnumeratorFromM(M,x_len,n,q): takes as input a x_len x n matrix M with entries in GF(q) and computes`):
  print(`the weight enumerator for the code generated by this matrix`):

 elif procedure = R1op then:
  print(`R1op(M,i1,i2): swaps M[i1] and M[i2]`):

 elif procedure = C1op then:
  print(`C1op(M,j1,j2,k): swaps column j1 of M with column j2.`):

 elif procedure = Step1 then:
  print(`This corresponds to step 1 in an algorithm included in`):
  print(`Hill's nice textbook "A first course in coding theory"`):
  print(`which puts a generating matrix into standard form`):

 elif procedure = Step2  then:
  print(`This corresponds to step 2 in an algorithm included in`):
  print(`Hill's nice textbook "A first course in coding theory"`):
  print(`which puts a generating matrix into standard form`):

 elif procedure = Step3 then:
  print(`This corresponds to step 3 in an algorithm included in`):
  print(`Hill's nice textbook "A first course in coding theory"`):
  print(`which puts a generating matrix into standard form`):

 elif procedure = SF then:
  print(`SF(q,M): takes as input a matrix M with entries in GF(q) and converts it to standard form`)
:
 elif procedure = PCM then:
  print(`PCM(q,M): takes as input as matrix M with entries in GF(q)`):
  print(`and outputs the corresponding parity check matrix`):

 elif procedure = EncodeHillBCH then:
  print(`#EncodeHillBCH(d,q,x): takes as input a prime p, a pos. integer d,`):
  print(`and a vector x entries in GF(q) such that nops(x) + d - 1 <= q-1.`):
  print(`We will return a vector of length k+d-1 by extending x so that it solves the following equations:`):
  print(`seq(add(i^j * x[i], i=1..nops(x)) + add(i^j * y[i], i=(nops(x)+1)..n), j=0..(d-2))`):

 elif procedure = HillBCHMatrix then:
  print(`HillBCHMatrix(n,d,q): create a matrix of size n+1-d x n`):
  print(`with entries in GF(q) (q prime) that is a generating matrix for the`):
  print(`class of BCH codes described in Hill's nice textbook`):
  print(`"A First Course in Coding Theory"`):

 elif procedure = MacWm then:
  print(`MacWm(f,n,k,q,z): The MacWilliams transform takes the weight-enumerator of a linear [n,k]-code over GF(q) and outputs`):
  print(`the weight-enumerator of the dual code (that is [n,n-k])`):

 elif procedure = AveAndMoms then:
  print(`Credit: This method is from Dr. Z's nice Maple package Stat.txt found here:`):
  print(`https://sites.math.rutgers.edu/~zeilberg/EM24/Stat.txt`):
  print(`AveAndMoms(f,x,N): Given a probability generating function`):
  print(`f(x) (a polynomial, rational function and possibly beyond)`):
  print(`returns a list whose first entry is the average `):
  print(`(alias expectation, alias mean)`):
  print(`followed by the variance, the third moment (about the mean) and`):
  print(`so on, until the N-th moment (about the mean).`):
  print(`If f(1) is not 1, than it first normalizes it by dividing`):
  print(`by f(1) (if it is not 0) .`):
  print(`For example, try:`):
  print(`AveAndMoms(((1+x)/2)^100,x,4);`):

 elif procedure = Alpha then:
  print(`Credit: This method is from Dr. Z's nice Maple package Stat.txt found here:`):
  print(`https://sites.math.rutgers.edu/~zeilberg/EM24/Stat.txt`):
  print(`Alpha(f,x,N): Given a probability generating function`):
  print(`f(x) (a polynomial, rational function and possibly beyond)`):
  print(`returns a list, of length N, whose `):
  print(`(i) First entry is the average`):
  print(`(ii): Second entry is the variance`):
  print(`for i=3...N, the i-th entry is the so-called alpha-coefficients`):
  print(`that is the i-th moment about the mean divided by the`):
  print(`variance to the power i/2 (For example, i=3 is the skewness`):
  print(`and i=4 is the Kurtosis)`):
  print(`If f(1) is not 1, than it first normalizes it by dividing`):
  print(`by f(1) (if it is not 0) .`):
  print(`For example, try:`):
  print(`Alpha(((1+x)/2)^100,x,4);`):

 elif procedure = Di then:
  print(`Credit: This method is from Dr. Z's nice Maple package and article found here:`):
  print(`https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/calabi.html`):
  print(`Di(a,b)  outputs true with prob. a/b. Try: Di(1,3);`):

 elif procedure = Rqnk then:
  print(`Credit: This method is from Dr. Z's nice Maple package and article found here:`):
  print(`https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/calabi.html`):
  print(`Rqnk(q,n,k): a (uniformly at) random n row-echelon matrices over GF(q)^n representing k-dim subspaces. Try`):
  print(`Rqnk(2,10,5);`):

 elif procedure = HillWeightEnumeratorBCH then:
  print(`HillWeightEnumeratorBCH(n,d,q,t): returns a weight enumerator as a polynomial in t`):
  print(`for the [n,n-d+1,d] code over GF(q). We require d <= n <= q-1, where q is prime.`):

 elif procedure = HillWeightEnumeratorFast1 then:
  print(`HillWeightEnumeratorBCHFast1(n,d,q,t): returns a weight enumerator as a polynomial in t`):
  print(`for the [n,n-d+1,d] code over GF(q). This algorithm is faster than HillWeightEnumeratorBCH(n,d,q,t)`):
  print(`and corresponds to algorithm 1 presented in the corresponding article`):

 elif procedure = HillWeightEnumeratorFast2 then:
  print(`HillWeightEnumeratorBCHFast2(n,d,q,t): returns a weight enumerator as a polynomial in t`):
  print(`for the [n,n-d+1,d] code over GF(q). This algorithm is faster than HillWeightEnumeratorBCH(n,d,q,t)`):
  print(`and corresponds to algorithm 2 presented in the corresponding article`):

 elif procedure = WeightEnumeratorRandLinearCode then:
  print(`WeightEnumeratorRandLinearCode(n,d,q): Finds the weight enumerator of a random`):
  print(`[n, n-d+1] code over GF(q). Here d does NOT mean distance.`):

 elif procedure = GetEnumerators then:
  print(`GetEnumerators(k,K,B,print_to_file,use_fast,filename): This prints out all weight enumerators of Hill's BCH codes`):
  print(`over all primes p in the range k <= p <= K with at most B codewords (if B=-1, there is no limit).`):
  print(`Note Hill specified that d should be odd, so we will print the weight enumerators such that 3 <= d <= n <= p-1.`):
  print(`If print_to_file is set to true, we will print to the file specified by argument filename, which defaults to`):
  print(`"database.txt" if no name is provided.`):

 elif procedure = ProveTheorem06BaseCase then:
  print(`ProveTheorem06BaseCase(): This does the computation to prove`):
  print(`the base case of theorem 0.6 in the article accompanying this file`):

 elif procedure = GetRandLinearCodeMeansAndVars then:
  print(`GetRandLinearCodeMeansAndVars(k,K,B,T,print_to_file,filename): `):
  print(`For each prime p in the range k<=p<=K and odd d where 3 <= d <= n and length n where n <= p-1, `):
  print(`we generate T random linear codes and get the average of the means and the average of the variances`):
  print(`as long as this code or its dual has at most B words.`):
  print(`If print_to_file is set to true, we will print to the file specified by argument filename, which defaults to`):
  print(`"database2.txt" if no name is provided.`):

 elif procedure = TestConjecture08 then:
  print(`TestConjecture08(): For each prime k <= p <= K, we compute the Hill BCH weight enumerators for all`):
  print(`3 <= n <= p-1 and 1 <= d <= n+1-num_moments. We then check to see if for each fixed p and n,`):
  print(`the first num_moments standardized moments of the weight distributions`):
  print(`are the same for all d in the specified range.`):

 else:
  print(`procedure name not recognized`):
 fi:
end:


#HD(u,v): takes as input 2 vectors u,v of the same length
#and outputs the Hamming distance
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:


#GetEnumeratorFromM(M,x_len,n,q): takes as input a x_len x n matrix M with entries in GF(q) and computes
#the weight enumerator for the code generated by this matrix
GetEnumeratorFromM := proc(M,x_len,n,q) local x,zero,word,ptr,T,enumerator,dist,i,go_back,val:
 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: 
    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:
 
 return(enumerator):
end:


#R1op(M,i1,i2): swaps M[i1] and M[i2]
R1op := proc(M,i1,i2)
 local u,v:
 u := M[i1]: v := M[i2]:
 M[i1] := v: M[i2] := u:
end:

#C1op(M,j1,j2,k): swaps column j1 of M with column j2.
C1op := proc(M,j1,j2,k)
 local v, i:
 for i from 1 to k do:
  v := M[i,j1]:
  M[i,j1] := M[i,j2]:
  M[i,j2] := v:
 od:
end:


#This corresponds to step 1 in an algorithm included in
#Hill's nice textbook "A first course in coding theory"
#which puts a generating matrix into standard form
Step1 := proc(M1,i,g,k,n) local j,g1:
   j := i+1:
   while j <= k do:
    g1 := M1[j,i]:
    if g1 <> 0 then:
     R1op(M1,i,j):
     RETURN(g1):
    fi:
    j += 1:
   od:

   j := i+1:
   while j <= n do:
    g1 := M1[i,j]:
    if g1 <> 0 then:
     C1op(M1,i,j,k):
     RETURN(g1):
    fi:
    j += 1:
   od:
end:


#This corresponds to step 2 in an algorithm included in
#Hill's nice textbook "A first course in coding theory"
#which puts a generating matrix into standard form
Step2 := proc(M1, q, g, i) local g_inv:
 g_inv := g&^(-1) mod q:
 M1[i] := (M1[i] * g_inv) mod q:
end:


#This corresponds to step 3 in an algorithm included in
#Hill's nice textbook "A first course in coding theory"
#which puts a generating matrix into standard form
Step3 := proc(M1, q, i, k) local j:
 for j from 1 to k do:
  if j <> i then:
   M1[j] := (M1[j] - M1[j,i] * M1[i]) mod q:
  fi:
 od:
end:


#SF(q,M): takes as input a matrix M with entries in GF(q)
#and converts it to standard form
SF := proc(q,M) 
 local k,i,j,n,M1,g,g1,g_inv:
 k := nops(M):
 n := nops(M[1]):
 M1 := M:
 M1 := Array(M1): 
 
 for i from 1 to k do:
  g := M1[i,i]:
  if g = 0 then:
   g := Step1(M1,i,g,k,n):
  fi:
  Step2(M1,q,g,i):
  Step3(M1,q,i,k):
 od:

 convert(M1,list,nested=true):
end: 


#PCM(q,M): takes as input as matrix M with entries in GF(q)
#and outputs the corresponding parity check matrix
PCM:=proc(q,M) local k,n,i,j,H:
k:=nops(M):
n:=nops(M[1]):
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:

[seq([seq(H[i,j],j=1..n)],i=1..n-k)]:
end:


#EncodeHillBCH(d,q,x): takes as input a prime p, a pos. integer d,
#and a vector x entries in GF(q) such that nops(x) + d - 1 <= q-1.
#We will return a vector of length k+d-1 by extending x so that it solves the following equations:
#seq(add(i^j * x[i], i=1..nops(x)) + add(i^j * y[i], i=(nops(x)+1)..n), j=0..(d-2)):
EncodeHillBCH := 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(n,d,q): create a matrix of size n+1-d x n
#with entries in GF(q) (q prime) that is a generating matrix for the
#class of BCH codes described in Hill's nice textbook
#"A First Course in Coding Theory"
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] := EncodeHillBCH(d,q,x):
 od:
 M:
end:


#HillBCHPCM(n,d,q): create a matrix of size d-1 x n
#with entries in GF(q) (q prime) that is a parity check matrix for the
#class of BCH codes described in Hill's nice textbook
#"A First Course in Coding Theory"
HillBCHPCM := proc(n,d,q) local M,i,j:
 [seq([seq(j^i, j=1..n)], i=0..d-2)]:
end:


#HillWeightEnumeratorBCH(n,d,q,t): returns a weight enumerator as a polynomial in t
#for the [n,n-d+1,d] code over GF(q). We require d <= n <= q-1, where q is prime.
HillWeightEnumeratorBCH := proc(n,d,q,t) local i,used_dual,x_len,M,enumerator:
 
 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:
 
 enumerator := GetEnumeratorFromM(M,x_len,n,q):

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


#HillWeightEnumeratorBCHFast1(n,d,q,t): returns a weight enumerator as a polynomial in t
#for the [n,n-d+1,d] code over GF(q). This algorithm is faster than HillWeightEnumeratorBCH(n,d,q,t)
#and corresponds to algorithm 1 presented in the corresponding article
HillWeightEnumeratorBCHFast1 := proc(n,d,q,t) local counts,i,k,enumerator,first:
 counts := Array([seq(q^i - 1, i=1..n-d+1)]):
 k := d + 1:
 enumerator := 1:
 while k <= n+1 do:
  first := counts[k-d]:
  
  for i from 0 to n-k+1 do:
   counts[k-d+i] -= first * binomial(k-1+i, k-1):
  od:

  enumerator += t^(k-1) * first * binomial(n, k-1):

  k += 1:
 od:
 enumerator:
end:



#HillWeightEnumeratorBCHFast2(n,d,q,t): returns a weight enumerator as a polynomial in t
#for the [n,n-d+1,d] code over GF(q). This algorithm is faster than HillWeightEnumeratorBCH(n,d,q,t)
#and corresponds to algorithm 2 presented in the corresponding article
HillWeightEnumeratorBCHFast2 := proc(n,d,q,t) local counts,i,k,enumerator,first:
 counts := Array([seq(q^i - 1, i=1..n-d+1)]):
 k := d + 1:
 enumerator := 1:
 while k <= n+1 do:
  first := counts[k-d]:
  
  for i from 1 to n-k+1 do:
    first := first * (k+i-1) / i:
   counts[k-d+i] -= first:
  od:

  enumerator += t^(k-1) * first:

  k += 1:
 od:
 enumerator:
end:


#WeightEnumeratorRandLinearCode(n,d,q): Finds the weight enumerator of a random
#[n, n-d+1] code over GF(q). Here d does NOT mean distance.
WeightEnumeratorRandLinearCode := proc(n,d,q) local M,k:
 k := n-(d-1):
 M := Rqnk(q,n,k):
 
 if k <= d-1 then:
  GetEnumeratorFromM(M,k,n,q):
 else:
  expand(simplify(MacWm(GetEnumeratorFromM(PCM(q,SF(q,M)),d-1,n,q), n, d-1, q, t))):
 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:


#Credit: This method is from Dr. Z's nice Maple package Stat.txt found here:
#https://sites.math.rutgers.edu/~zeilberg/EM24/Stat.txt
#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:


#Credit: This method is from Dr. Z's nice Maple package Stat.txt found here:
#https://sites.math.rutgers.edu/~zeilberg/EM24/Stat.txt
#Alpha(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list, of length N, whose 
#(i) First entry is the average
#(ii): Second entry is the variance
#for i=3...N, the i-th entry is the so-called alpha-coefficients
#that is the i-th moment about the mean divided by the
#variance to the power i/2 (For example, i=3 is the skewness
#and i=4 is the Kurtosis)
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#Alpha(((1+x)/2)^100,x,4);
Alpha:=proc(f,x,N) local gu,i:
 gu:=AveAndMoms(f,x,N):
 if gu=FAIL then
  RETURN(gu):
 fi:

if gu[2]=0 then
print(`The variance is 0`):
  RETURN(FAIL):
fi:

[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..N)]:

end:


#Credit: This method is from Dr. Z's nice Maple package and article found here:
#https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/calabi.html
#Di(a,b)  outputs true with prob. a/b. Try: Di(1,3);
Di:=proc(a,b) evalb(rand(1..b)()<=a):end:


#Credit: This method is from Dr. Z's nice Maple package and article found here:
#https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/calabi.html
#Rqnk(q,n,k): a (uniformly at) random n row-echelon matrices over GF(q)^n representing k-dim subspaces. Try
#Rqnk(2,10,5);
Rqnk:=proc(q,n,k) local ra1,gu1,gu2,gu,lu,i1,i:

if k>n  then
 RETURN({}):
fi:

ra1:=rand(0..q-1):

if k=1 then
 
if n=1 then
  RETURN([[1]]):
  else
   if Di(q-1,q^n-1) then
     RETURN([[1,0$(n-1)]]):
   else
    gu:=Rqnk(q,n-1,1):
   RETURN([[ra1(),op(gu[1])]]):
 fi:

fi:
fi:

if Di(q^k-1,q^n-1) then
 
gu1:=Rqnk(q,n-1,k-1):

gu:=RETURN([[1,0$(n-1)],seq([0,op(gu1[i1])],i1=1..nops(gu1))]):

else
gu2:=Rqnk(q,n-1,k):
lu:=[seq(ra1(),i=1..k)]:

RETURN([seq([lu[i1],op(gu2[i1])],i1=1..k)]):
fi:

end:


#ProveTheorem06BaseCase(): This does the computation to prove
#the base case of theorem 0.6 in the article accompanying this file
ProveTheorem06BaseCase := proc() local f,a0,ad:
 print(`the following code:`):
 print(` a0 := 1: `):
 print(` ad := n * (q-1): `):
 print(` f :=  a0 + ad * t^d + (q^2 - (a0 + ad))*t^(d+1): `):
 print(` f := subs(d=n-1, f): `):
 print(` simplify(AveAndMoms(f,t,2)) `):
 print(` outputs: `):
 a0 := 1:
 ad := n * (q-1):
 f :=  a0 + ad * t^d + (q^2 - (a0 + ad))*t^(d+1):
 f := subs(d=n-1, f):
 print(simplify(AveAndMoms(f,t,2))):
end:


#GetEnumerators(k,K,B,print_to_file,use_fast,filename): This prints out all weight enumerators of Hill's BCH codes
#over all primes p in the range k <= p <= K with at most B codewords (if B=-1, there is no limit).
#Note Hill specified that d should be odd, so we will print the weight enumerators such that 3 <= d <= n <= p-1.
#If print_to_file is set to true, we will print to the file specified by argument filename, which defaults to
#`database.txt` if no name is provided.
GetEnumerators := proc(k,K,B,print_to_file,use_fast,filename:=`database.txt`) 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 or B=-1 then: 
     if use_fast then:
      enumerator := HillWeightEnumeratorBCHFast2(n,d,p,t):
     else:
      enumerator := HillWeightEnumeratorBCH(n,d,p,t):
     fi:
     avg_and_var := AveAndMoms(enumerator,t,2): 
     
     if print_to_file then:
      fprintf(filename, `d=%d, n=%d, p=%d, mean=%e, var=%e, enumerator=%a\n`,  
                       d, n, p, avg_and_var[1], avg_and_var[2], enumerator):
      fflush(filename):
     else:
      printf(`d=%d, n=%d, p=%d, mean=%e, var=%e, enumerator=%a\n`,  
                       d, n, p, avg_and_var[1], avg_and_var[2], enumerator):
     fi:
    fi:
   od:
  od:
  p := nextprime(p):
 od:
 if print_to_file then FileTools[Text]:-Close(filename): fi:
 return():
end:


#GetRandLinearCodeMeansAndVars(k,K,B,T,print_to_file,filename): 
#For each prime p in the range k<=p<=K and odd d where 3 <= d <= n and length n where n <= p-1, 
#we generate T random linear codes and get the average of the means and the average of the variances
#as long as this code or its dual has at most B words.
#If print_to_file is set to true, we will print to the file specified by argument filename, which defaults to
#`database2.txt` if no name is provided.
GetRandLinearCodeMeansAndVars := proc(k,K,B,T,print_to_file,filename:=`database2.txt`) local p,d,n,enumerator,avg_and_var,i:
 p := nextprime(k-1):
 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 or B=-1 then:
      avg_and_var := [0,0]:
      for i from 1 to T do:
       enumerator := WeightEnumeratorRandLinearCode(n,d,p):
       avg_and_var += AveAndMoms(enumerator, t, 2):
      od:
     avg_and_var /= T: 
     
     if print_to_file then:
      fprintf(filename, `d=%d, n=%d, p=%d, mean=%e, var=%e\n`,  
                       d, n, p, avg_and_var[1], avg_and_var[2]):
      fflush(filename):
     else:
      printf(`d=%d, n=%d, p=%d, mean=%e, var=%e\n`,  
                       d, n, p, avg_and_var[1], avg_and_var[2]):
     fi:
    fi:
   od:
  od:
  p := nextprime(p):
 od:
 if print_to_file then FileTools[Text]:-Close(filename): fi:
 return():
end:


#TestConjecture08(): For each prime k <= p <= K, we compute the Hill BCH weight enumerators for all
#3 <= n <= p-1 and 1 <= d <= n+1-num_moments. We then check to see if for each fixed p and n,
#the first num_moments standardized moments of the weight distributions
#are the same for all d in the specified range.
TestConjecture08 := proc(k,K,num_moments) local p,n,enumerator,avg_and_var_def,avg_and_var,d:
 p := nextprime(k-1):
  while p <= K do:
   for n from 3 to p-1 do:
    if n+1-num_moments > 0 then:
     enumerator := HillWeightEnumeratorBCHFast2(n,1,p,t):
     avg_and_var_def := Alpha(enumerator,t,num_moments): fi:
     for d from 2 to n+1-num_moments do:
        enumerator := HillWeightEnumeratorBCHFast2(n,d,p,t):
        avg_and_var := Alpha(enumerator,t,num_moments):    

        if avg_and_var_def <> avg_and_var then: return(FAIL): fi:
    od:
   od:
   p := nextprime(p):
  od:
 return(SUCCESS):
end: