Help:=proc(): print(` MacWm(f,n,k,q,z), LtoC(q,M), HD(u,v), CtoL(n,x,g), AllFactors(n,q,x), WtEnumerator(q,M,t),  AveAndMoms(f,x,N), Alpha(f,x,N),`):end:


with(combinat):

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

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


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

#WtEnumerator(q,M,t): The weight-enumerator, in t, of the linear code generated by the basis M over GF(q)
WtEnumerator := proc(q,M,t) local v,C,zero,n:
 n := nops(M[1]):
 zero := [0$n]:
 C := LtoC(q,M):
 add(t^HD(zero,v), v in C):
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:


#AllFactors(n,q,x): all the factors of x^n-1 mod q (q is a prime)
AllFactors:=proc(n,q,x) local S,s:
S:=Factor(x^n-1) mod q:
S:=powerset({op(S)}) minus {{},{op(S)}}:

{seq(expand(convert(s,`*`)) mod q,s in S)}:
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:


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