######################################################################
## proj1.txt: Save this file as  proj1.txt to use it,                #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
#then type: read `proj1.txt`: <Enter>                                #
## Then follow the instructions given there                          #
##                                                                   #
## Written by  Himanshu Chandra, Nkhalo Malawo, Lucy Martinez, and   #
## Doron Zeilberger						     #
######################################################################

#VERSION OF April 30 2024

with(combinat):
with(DocumentTools):



print(`First Written: April 2024: tested for Maple 20 `):
print():
print(`This is proj1.txt, a Maple package `):
print(`accompanying Himanshu Chandra, Nkhalo Malawo, and Lucy Martinez's  article: `):
print(`NAME OF ARTICLE`):

print(`-------------------------------`):


print(`Please report all bugs to: lm1154@scarletmail.rutgers.edu  .`):
print():
print(`-------------------------------`):
print(`For general help, and a list of the MAIN functions,`):
print(` type "Help();". For specific help type "Help(procedure_name);" `):

print(`-------------------------------`):
print(`For a list of the supporting functions type:  Help1();`):
print(`For help with a specific procedure, type "Help(procedure_name);"`):
print():
print(`-------------------------------`):



Help1:=proc()
if args=NULL then

print(` The supporting procedures are: HD, LtoC, CtoL, CtoDual, AllFactors`):

else
 Help(args):
fi:
end:

Help:=proc()

if args=NULL then

print(`proj1.txt: A Maple package for Weight-Enumerators of all Cyclic Linear Codes $[n,k]$ over $GF(q)$ up to a Given Size`):
print(`The main procedures are: DualCodes, WtEnumerator, MacWm, AveAndMoms, Di,`):
print(`Rqnk, Span1, WtEnu, Stat, CreateTable, CreateTable2`):


elif nargs=1 and args[1]=HD then
print(`HD(u,v): The Hamming distance between two words (of the same length). Try:`):
print(`HD([0, 1, 1], [1, 0, 1]); `):


elif nargs=1 and args[1]=LtoC then
print(`LtoC(q,M): Inputs a list of basis vectors for our linear code over GF(q)`):
print(`and outputs all the codewords (the actual subset of GF(q)^n)`):
print(`with q^k elements. Try:`):
print(`LtoC(2, [[1, 1, 0], [0, 1, 1]]); `):


elif nargs=1 and args[1]=CtoL then
print(`CtoL(n,x,g): Inputs a generator polynomial g(x) over R_n (mod x^n-1) and`):
print(`outputs the generator matrix for the cyclic code <g(x)> over R^n (mod x^n-1). Try:`):
print(`CtoL(3,x,1+x); `):


elif nargs=1 and args[1]=CtoDual then
print(`CtoDual(n,q,x,g): Inputs a positive integer n, a prime q, a variable x, `):
print(`and a polynomial g (a generator of a cyclic code) mod q, and`):
print(`outputs the generator matrix for the dual code of CtoL(n,x,g). Try:`):
print(`CtoDual(3,2,x,1+x); `):


elif nargs=1 and args[1]=AllFactors then
print(`AllFactors(n,q,x): Inputs positive integer n, a prime q, and outputs`):
print(`all the factors of x^n-1 mod q. Try: `):
print(`AllFactors(2,2,x); `):


elif nargs=1 and args[1]=DualCodes then
print(`DualCodes(n,q,a): Generates all cyclic linear [n,n-a] dual codes over GF(q) of length n`):
print(`and outputs the Dual Codes followed by their corresponding generating polynomial g(x) of the Code C. Try:`):
print(`CLC(3,2,1); `):


elif nargs=1 and args[1]=WtEnumerator then
print(`WtEnumerator(q,C,z): The weight-enumerator, in z, of the linear code C`):
print(` over GF(q). Try:`):
print(`WtEnumerator(2,{[0, 0, 0], [1, 1, 1]} ,z)`):



elif nargs=1 and args[1]=MacWm then
print(`MacWm(n,q,a,Cdual,z): Inputs a dual code C over GF(q) of length n,`):
print(`computes its weight enumerator, and outputs the weight enumerator polynomial`):
print(`for the code C using the MacWilliams identity. Try:`):
print(`MacWm(3,2,1, {[0, 0, 0], [1, 1, 1]},z); `):


elif nargs=1 and args[1]=AveAndMoms then
print(`AveAndMoms(f,z,N): Given a probability generating function`):
print(`f(z) (a polynomial, rational function and possibly beyond)`):
print(`returns a list whose first entry is the average, followed by the variance,`):
print(`the third moment, so on, until the N-th moment (about the mean). Try:`):
print(`AveAndMoms(((1+z)/2)^100,z,2); `):


elif nargs=1 and args[1]=Di then
print(`Di(a,b) outputs true with prob. a/b. Try: `):
print(`Di(1,3); `):


elif nargs=1 and args[1]=Rqnk then
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 nargs=1 and args[1]=Span1 then
print(`Span1(M,n,q): inputs a base M for a vector subspace over GF(q)^n, finds all the vectors. Try:`):
print(`M:=Rqnk(2,10,3); Span1(M,10,2); `):


elif nargs=1 and args[1]=WtEnu then
print(`WtEnu(M,q,z): The  weight-enumerator according to the variable z, of the subspace of GF(q)^n`):
print(`generated by M (where n:=nops(M[1]). Try:`):
print(`M:=Rqnk(2,10,3);  WtEnu(M,2,z); `):


elif nargs=1 and args[1]=Stat then
print(`Stat(n,q,a,z): Inputs a positive integer n, a prime q, small a, and a variable z,`):
print(`generates all cyclic linear [n,n-a] dual codes over GF(q) of length n,`):
print(`DualCodes(n,q,a), and outputs a list of lists with the following information`):
print(`the generating polynomial (g(x)) of the code C, weight enumerator of the code C,`):
print(`the average weight for each code C, and the standard deviation for each code C. Try:`):
print(`Stat(7,2,3,z); `):


elif nargs=1 and args[1]=CreateTable then
print(`CreateTable(n,q,a,z): Inputs a positive integer n, a prime q, small a and variable z,`):
print(`runs Stat(i, q, a, z) for 1<=i<=n, outputs a table with 4 columns where`):
print(`the first column outputs the triple [n,q,k] corresponding to the [n,k]-code C`):
print(`(or [n,n-a]) over GF(q), the second column outputs the generating polynomial g(x)`):
print(`of the code C, the third column outputs the weight enumerator of the code C in z,`):
print(`and the fourth column outputs a pair {[Avg Wt, Std]} where Avg Wt is the average`):
print(`weight of the [n,q,k]-code and Std is the standard deviation of the weight. Try:`):
print(`CreateTable(30,2,3,z); `):


elif nargs=1 and args[1]=CreateTable2 then
print(`CreateTable2(n,q,a,z): Inputs a positive integer n>3, a prime q, small a and variable z,`):
print(`runs Stat(i, q, a, z) for 1<=i<=n, outputs a table with 5 columns where`):
print(`the first fourth columns are as in CreateTable(n,q,a,z) and`):
print(`the fifth column outputs a pair {[ZAvg Wt, ZStd]} where ZAvg Wt is the average`):
print(`weight of the [n,q,k]-code and Std is the standard deviation of the weight`):
print(`using the CalabiWilf.txt file from Zeilberger. Try:`):
print(`CreateTable2(10,2,3,z); `):


else
 print(`There is no such thing as`, args):
fi:

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:



# LtoC(q,M): Inputs a list of basis vectors for our linear code over GF(q)
# and 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:



# CtoL(n,x,g): Inputs a generator polynomial g(x) over R_n (mod x^n-1)
# and outputs the generator matrix for the cyclic code <g(x)> over R^n (mod x^n-1)
CtoL:=proc(n,x,g) local r,L,i:
r:=degree(g,x):
L:=[seq(coeff(g,x,i),i=0..r)]:
#n-(i+r+1): the first 0$i has i elements, then r+1, then we need n- those

#this will be a n-r by n matrix
[seq( [0$i,op(L),0$(n-r-1-i)], i=0..n-r-1)]:
end:



# CtoDual(n,q,x,g): Inputs a positive integer n, a prime q, a variable x, 
# and a polynomial g (a generator of a cyclic code) mod q
# and outputs the generator matrix for the dual code of CtoL(n,x,g)
CtoDual:=proc(n,q,x,g) local i,h,k:
h:=Quo(x^n-1,g,x) mod q:
k:=degree(h,x):
h:=add(coeff(h,x,i)* x^(k-i), i=0..k):
CtoL(n,x,h):
end:



# AllFactors(n,q,x): Inputs positive integer n, a prime q, and outputs
# all the factors of x^n-1 mod q except for 1 and x^n-1
AllFactors:=proc(n,q,x) local R,S,s,i,j:
R:=Sqrfree(x^n-1) mod q:
S:=map(a->Berlekamp(a[1], x), R[2]) mod q:
S:=map(a->convert(a, list), S):
S:=[seq([seq(op(S[i]), j in 1..R[2][i][2])], i in 1..nops(S))]:
S:=ListTools:-Flatten(S):
S:={op(powerset(S))} minus {[]}:
S:={seq(expand(convert(s,`*`)) mod q,s in S)}:
S minus {expand(x^n - 1) mod q}:
end:

 

#DualCodes(n,q,a): Generates all cyclic linear [n,n-a] dual codes over GF(q) of length n
# and outputs the Dual Codes followed by their corresponding
# generating polynomial g(x) of the Code C
DualCodes:=proc(n,q,a) local polys,C,g,temp,i,j,h,S,M: 
if a=n then 
 return {[0$n]}:
fi:

# need to add when a=0? This means g(x)=1 which should output nothing
# Code is everything in Rn which means the dual code is nothing


polys:=AllFactors(n,q,x):
polys:=polys union {1}:

S:={}:
for h in polys do
 g:=Quo(x^n-1,h,x) mod q:
 if degree(g)=a then
  S:=S union {g}:
 fi:
od:

C:={}:
for g in S do
 #for every g(x) with deg(g)=a, compute the generating matrix of the dual code
 M:={op(CtoDual(n,q,x,g))}:
 #includes the generating poly g(x) in the output
 C:= C union {[LtoC(q,M),g]}:
od:
C:
end:



#WtEnumerator(q,C,t): The weight-enumerator, in t, of the linear code C over GF(q)
WtEnumerator:=proc(q,C,t) local v,zero,n,C1:
 n:=nops(C[1]):
 zero:=[0$n]:
 add(t^HD(zero,v), v in C):
end:



#MacWm(n,q,a,Cdual,z): Inputs a Dual Code C over GF(q) of length n,  
# computes its weight enumerator, and outputs the weight enumerator polynomial
# for the code C using the MacWilliams identity
MacWm:=proc(n,q,a,Cdual,z) local k,WCd:
k:=n-a:
#computes the weight enumerator for the dual code of C
WCd:=WtEnumerator(q,Cdual,z):
expand(simplify(1/(q^(n-k))*(1+(q-1)*z)^n*subs(z=(1-z)/(1+(q-1)*z),WCd))):
end:



#AveAndMoms(f,z,N): Given a probability generating function
#f(z) (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+z)/2)^100,z,4);
AveAndMoms:=proc(f,z,N) local mu,F,memu1,gu,i:
mu:=simplify(subs(z=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(z=1,diff(F,z))):
gu:=[memu1]:

F:=F/z^memu1:
F:=z*diff(F,z):
for i from 2 to N do
 F:=z*diff(F,z):
 gu:=[op(gu),simplify(subs(z=1,F))]:
od:
gu:
end:



#Di(a,b)  outputs true with prob. a/b. Try: Di(1,3);
Di:=proc(a,b) 
 evalb(rand(1..b)()<=a):
end:


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



#Span1(M,n,q): inputs a base M for a vector subspace over GF(q)^n, finds all the vectors. Try:
#M:=Rqnk(2,10,3); Span1(M,10,2);
Span1:=proc(M,n,q) local k,gu,v,gu1,i,M1:
k:=nops(M):

if k=0 then
 RETURN({[0$n]}):
fi:
v:=M[k]:
M1:=[op(1..k-1,M)]:
gu:=Span1(M1,n,q):
{seq(seq(gu1+i*v mod q, i=0..q-1),gu1 in gu)}:
end:



#Wt1(v): the number of non-zero entries in the vector v
Wt1:=proc(v) local i,co:
co:=0:
for i from 1 to nops(v) do
 if v[i]<>0 then
   co:=co+1:
 fi:
od:
co:
end:



#WtEnu(M,q,z): The  weight-enumerator according to the variable z, of the subspace of GF(q)^n generated by M (where n:=nops(M[1]). Try:
#M:=Rqnk(2,10,3); WtEnu(M,2,z);
WtEnu:=proc(M,q,z) local n,gu,v:
n:=nops(M[1]):
gu:=Span1(M,n,q) minus {[0$n]}:
1+add(z^Wt1(v),v in gu):
end:



#Stat(n,q,a,z): inputs a positive integer n, a prime q, small a, and a variable z,
# generates all cyclic linear [n,n-a] dual codes over GF(q) of length n,
# DualCodes(n,q,a), and outputs a list of lists with the following information
# the generating polynomial (g(x)) of the code C, weight enumerator of the code C,
# the average weight for each code C, and the standard deviation for each code C
# For example, Stat(7,2,3,z); returns
# [[x^3 + x + 1, z^7 + 7*z^4 + 7*z^3 + 1, 7/2, sqrt(7)/2], [x^3 + x^2 + 1, z^7 + 7*z^4 + 7*z^3 + 1, 7/2,sqrt(7)/2]]
Stat:=proc(n,q,a,z) local Cduals,i,c,S,L,sd:
 S:=[]:
 Cduals:=DualCodes(n,q,a):
 #adds the generating polynomials g(x) for each [n,k]-code over GF(q)
 for i from 1 to nops(Cduals) do
  S:=[op(S),[Cduals[i,2]]]:
 od: 

 #adds the Weight Enumerator Polynomial to the list S for each Code
 #second element in the list of lists S
 for i from 1 to nops(S) do
  S[i]:=[op(S[i]),MacWm(n,q,a,Cduals[i,1],z)]:
 od:
 
 #computes the average and the variance, then uses variance to compute 
 #standard deviation 
 for i from 1 to nops(S) do
  L:=AveAndMoms(S[i,2],z,2):
  sd:=sqrt(L[2]):
  S[i]:=[op(S[i]),op([L[1],sd])]:
 od:
S:
end:



#CreateTable(n,q,a,z): Inputs a positive integer n, a prime q, small a and variable z,
# runs Stat(i, q, a, z) for 1<=i<=n, outputs a table with 4 columns where
# the first column outputs the triple [n,q,k] corresponding to the [n,k]-code C
# (or [n,n-a]) over GF(q), the second column outputs the generating polynomial g(x)
# of the code C, the third column outputs the weight enumerator of the code C in z,
# and the fourth column outputs a pair {[Avg Wt, Std]} where Avg Wt is the average
# weight of the [n,q,k]-code and Std is the standard deviation of the weight
CreateTable:=proc(n,q,a,z) local T,M,C,Our,O,Z,i,j,P,xml,str,L,f:
 T:=[["{[n,q,k]}","Generator Polynomial","Weight Numerator","{Avg Wt,Std}"]]: 
 for i from 1 to n do
  #Stat(n,q,a,z)
  Our:=Stat(i, q, a, z):
  P:=[]:
  L:=[i,q,i-a]:
  for j from 1 to nops(Our) do
   O:=Our[j]:
   P:=[op(P), [{L},O[1],O[2],{[O[3],O[4]]}]]:
  od:
  T := [op(T), op(P)]: 
 od: 
#Tabulate(T, width = 50):
xml:=subsindets(Tabulate((T),'output'=':-XML'),"pagebreak"="none",()->"pagebreak"="row"):
str:=ContentToString(Layout:-Worksheet(Layout:-Group(xml))):

## If you want the table alone in a new GUI tab
Worksheet:-Display(str):

## If you want it embedded in this same worksheet
#InsertContent(str):
end:



#CreateTable2(n,q,a,z): Inputs a positive integer n>3, a prime q, small a and variable z,
# runs Stat(i, q, a, z) for 1<=i<=n, outputs a table with 5 columns where
# the first column outputs the triple [n,q,k] corresponding to the [n,k]-code C
# (or [n,n-a]) over GF(q), the second column outputs the generating polynomial g(x)
# of the code C, the third column outputs the weight enumerator of the code C in z,
# the fourth column outputs a pair {[Avg Wt, Std]} where Avg Wt is the average
# weight of the [n,q,k]-code and Std is the standard deviation of the weight,
# the fifth column outputs a pair {[ZAvg Wt, ZStd]} where ZAvg Wt is the average
# weight of the [n,q,k]-code and Std is the standard deviation of the weight
# using the CalabiWilf.txt file from Zeilberger
CreateTable2:=proc(n,q,a,z) local T,M,C,Our,O,Z,i,j,P,xml,str,L,f:
 T:=[["{[n,q,k]}","Generator Polynomial","Weight Numerator","{AvgWt,Std}","{ZAvgWt,ZStd}"]]: 
 for i from 4 to n do
  #Rqnk(q,n,k) where k=n-a
  M:=Rqnk(q, i, i-a):
  #Span1(M,n,q) where M=Rqnk(q,n,k)
  C:=Span1(M, i, q):

  
  #WtEnu(M,q,x) where M=Rqnk(q,n,k)
  f:=WtEnu(M,q,z):
  Z:=AveAndMoms(f,z,2):
  
  #Stat(n,q,a,z)
  Our:=Stat(i, q, a, z):
  P:=[]:
  L:=[i,q,i-a]:
  for j from 1 to nops(Our) do
   O:=Our[j]:
   P:=[op(P), [{L},O[1],O[2],{[O[3],O[4]]}]]:
   #P:=[op(P), [{L},O[1],O[2],{[O[3],O[4]]},{[Z[1],Z[2]]}]]:
  od:
  T := [op(T), op(P)]: 
 od: 
#Tabulate(T, width = 50):
xml:=subsindets(Tabulate((T),'output'=':-XML'),"pagebreak"="none",()->"pagebreak"="row"):
str:=ContentToString(Layout:-Worksheet(Layout:-Group(xml))):

## If you want the table alone in a new GUI tab
Worksheet:-Display(str):

## If you want it embedded in this same worksheet
#InsertContent(str):
end: