############################################################################
## SMCboole.txt: Save this file as SMCboole.txt. To use it, stay in the    #
## same directory, get into Maple (by typing: maple <Enter> )              #
## and then type:  read `SMCboole.txt` : <Enter>                           #
## Thene follow the instructions given there                                #
##                                                                         #
## Written by Blair Seidler and Doron Zeilberger, Rutgers University       #
############################################################################ 

#Created: Nov. 2022
#This version:  Apr. 2023
#SMCboole.txt: A Maple package to implement Blair Seidler's and Doron Zeilberger's
#Symbolic Moment Calculus applied to Boolean functions
#Please report bugs to DoronZeil@gmail.com



print(`Created: Nov. 2022`):
print(`This version: Apr. 3, 2023`):
print(`This is SMCboole.txt, a Maple package to implement`):
print(`Blair Seidler's Doron Zeilberger's Symbolic Moment Calculus for Boolean functions/`):
print(`as it is described in theirarticle:`):
print(` Symbolic Moment Calculus III: Boolean Functions `):
print(``):
print(`Written by Blair Seidler and Doron Zeilberger`):
lprint(``):
print(`Please report bugs to DoronZeil at gmail dot com`):
lprint(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/smcIII.html`):

print(`-----------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);`):
 print(``):
print(`-----------------------------`):



print(`-----------------------------`):
 print(`For a list of the SUPPORTING procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);`):
 print(``):
print(`-----------------------------`):


print(`-----------------------------`):
 print(`For a list of the  NUMERIC procedures type ezraN();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);`):
 print(``):
print(`-----------------------------`):


print(`-----------------------------`):
 print(`For a list of the  SIMULATION procedures type ezraSi();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);`):
 print(``):
print(`-----------------------------`):


with(combinat):



ezraSi:=proc(): 
if args=NULL then
print(` The  Simulation procedures are: `):
else
 ezra(args):
fi:
end:
 
 
ezraN:=proc(): 
if args=NULL then
print(` The  Numeric procedures are:  MomLn `):
else
 ezra(args):
fi:
end:
 
ezra1:=proc(): 
if args=NULL then
print(` The supporting procedures are:  Conf, Conf1, ExpaS, ExpkC, Imags, Kids1, Kids11, LtoM, Mats,  MOMrk1, NuVerts, NuVerts1, Reps, Trans1, Vecsk, Wt `):
else
 ezra(args):
fi:
end:

ezra:=proc()
if args=NULL then
print(`The main procedures are: Aek1, Aek2, AekM, Cor, Cov, MOMrk, Moms,  SMOMrk, VARr `):


elif nargs=1 and args[1]=Aek1 then
print(`Aek1(n,k): The  formula, due to Thotsaporn Thanatipanonda, for the mean of the r.v. "number of k-dim subcubes contained in f", for a random Boolean function with n variables. n is symbolic, but k ins numeric (but arbitrary) `):
print(`It should give the same output as Moms([k],n);  (after they are both simplified).Try: `):
print(`seq(simplify(Aek1(n,i)-Moms([i],n)),i=1..4);`):

elif nargs=1 and args[1]=Aek2 then
print(`Aek2(n,k): The  formula, due to Thotsaporn Thanatipanonda, for the second moment of the r.v. "number of k-dim subcubes contained in f", for a random Boolean function with n variables. n is symbolic, but k ins numeric (but arbitrary) `):
print(`It should give the same output as Moms([k,k],n);  (after they are both simplified).Try: `):
print(`seq(simplify(Aek2(n,i)-Moms([i,i],n)),i=1..4);`):

elif nargs=1 and args[1]=AekM then
print(`AekM(n,k1,k2): The  formula for the mixed moment of the r.v.'s "number of k1-dim subcubes contained in f" and "number of k2-dim subcubes contained in f", for a random Boolean function with n variables. n is symbolic, but k1<k2 are numeric (but arbitrary) `):
print(`It should give the same output as Moms([k1,k2],n);  (after they are both simplified).Try: `):
print(`seq(simplify(AekM(n,1,i)-Moms([1,i],n)),i=2..4);`):

elif nargs=1 and args[1]=Conf then
print(`Conf(A,r): Inputs a list of positive integers, of size k, say, outputs ALL lists of length k subsets of {1, ...,r}   [S1,...,Sk] such that Si has cardinality A[i] and their union is {1,..,r}. Try: `):
print(`Conf([2,3,2],5);`):

elif nargs=1 and args[1]=Conf1 then
print(`Conf1(A,r): Inputs a list of positive integers, of size k, say, outputs ALL lists of length k subsets of {1, ...,r}   [S1,...,Sk] such that Si has cardinality A[i].`):
print(`Try: Conf1([2,3,2],5);`):


elif nargs=1 and args[1]=Cor then
print(`Cor(d1,d2,n): The Correlation betweed X[d1] and X[d2], where X[r] is the r.v. defined on the sample-space of all Boolean functions (with the uniform distribution. Try:`):
print(`Cor(2,3,n);`):

elif nargs=1 and args[1]=Cov then
print(`Cov(d1,d2,n): The Covariance betweed X[d1] and X[d2], where X[r] is the r.v. defined on the sample-space of all Boolean functions (with the uniform distribution. Try:`):
print(`Cov(2,3,n);`):


elif nargs=1 and args[1]=ExpaS then
print(`ExpaS(F): Given a subcube where 0 indicates a wild-card, finds all the vertices beloning to it. Try:`):
print(`ExpaS([1,-1,0,1,0]);`):

elif nargs=1 and args[1]=ExpkC then
print(`ExpkC(v): Given a vector v of {0,1,-1} finds the set of {1,-1} vectors where 0 can be either 1 or -1. Try:`):
print(`ExpkC([1,0,-1,0]);`):

elif nargs=1 and args[1]=Imags then
print(`Imags(M): all the images of the matrix M under permuting the rows and the columns. try: `):
print(`Imags([[1,0,1],[1,1,1],[1,0,1]]);`):


elif nargs=1 and args[1]=Kids1 then
print(`Kids1(M): Inputs a 0-1 matrix M and outputs the set of all {0,1,-1} matrices M1 such that abs(M1[i][j])=M[i][j], in other words where 1 is re;aced by {1,-1}. Try:`):
print(`Kids1([[1,0],[0,1]]);`):

elif nargs=1 and args[1]=Kids11 then
print(`Kids11(M,i,j): Given a matrix M (represented as list of lists) with entries {0,-1,1} and a location [i,j] with M[i][j]=1 outputs the set of`):
print(`two matrices where the [i,j] entry is either 1 or of -1.  Try:`):
print(`Kids11([[0,1,0],[1,1,0]],2,1)`):

elif nargs=1 and args[1]=LtoM then
print(`LtoM(L,r): given a list of subsets L of {1,...,r} of size k, that cover {1,..,r}, outputs the 1-0 matrix M (given as a list of lists) such that  M[i][j]=0 iff j belongs to L[i]`):
print(`Try: `):
print(`LtoM([{1,2},{1,3},{2,3}],3);`):

elif nargs=1 and args[1]=Mats then
print(`Mats(n,L):  Inputs a positive ineteger n and a list L of non-negative integers of length k, say, outputs all k by n {-1,0,1} matrices `):
print(`where there are exactly L[1] 0's in the first row, exactly L[2] 0's in the second row, ... exactly L[k] 0's in the k-th row. Try:`):
print(`Mats(3,[1,1,1]);`):

elif nargs=1 and args[1]=MomLn then
print(`MomLn(L,n): For Numeric n, The L-mixed moment, i.e. E[X[L[1]]*X[L[2]]*...*X[L[k]]), where X[r]:=number of sub-cubes of dimension r.  for Boolean functions on n variables. VERY SLOW. For CHECKING PURPOSES ONLY. Try:`):
print(`MomLn([1,1,1],3);`):

elif nargs=1 and args[1]=MOMrk1 then
print(`MOMrk(r,k,n): The k-th moment of X[r] (NOT about the mean). Try:`):
print(`MOMrk(1,3,n);`):


elif nargs=1 and args[1]=MOMrk then
print(`MOMrk(r,k,n): The k-th moment of X[r] ABOUT THE MEAN. Try:`):
print(`MOMrk(1,3,n);`):

elif nargs=1 and args[1]=Moms then
print(`Moms(A,n): Inputs a list A of positive integers of length k, say, and a SYMBOL n, outputs an expression  as a polynomial in n and 2^n,`):
print(`for E[X[A[1]]X[A[2]]...X[A[k]]] where X[r] is the r.v. number of r-dim subcubes of a Boolean function f in n variables. Try`):
print(`Moms([1,1],n);`):

elif nargs=1 and args[1]=NuVerts then
print(`NuVerts(M,SP): inputs a matrix M of {0,1,-1} of size k by r, say, and a set-partition SP of {1,..,k} outputs the total number of vertices in that configuration. Try`):
print(`NuVerts([[1,0,-1],[1,1,-1]],{{1},{2}});`):

elif nargs=1 and args[1]=NuVerts1 then
print(`NuVerts1(M): inputs a matrix M in {0,1,-1}, outputs the number of vetrices in the union of subcubes given by its rows. Try:`):
print(`NuVerts1([[0,1],[1,0]]);`):

elif nargs=1 and args[1]=Reps then
print(`Reps(A,r): Given a list of positive intgegers A, of length k, say,  and a positive integer r outputs a set of pairs [M,multiplicitly], where M is  `):
print(` a k by by r 0-1 matrices with A[i] 0's in row i, and multipilicity is the number of its images under permuting rows and columns that still has the same A. Try:`):
print(`Reps([2,3],3);`):


elif nargs=1 and args[1]=SMOMrk then
print(`SMOMrk(r,k,n):  SCALED  k-th moment of X[r] about the mean. Try:`):
print(`SMOMrk(1,3,n);`):

elif nargs=1 and args[1]=Trans1 then
print(`Trans1(M,pi1,pi2); permutes the rows of M according to pi1 and the columns according to pi2. Try:`):
print(`Trans1([[a11,a12],[a21,a22]],[2,1],[2,1]);`):

elif nargs=1 and args[1]=VARr then
print(`VARr(n,r): The expression, in n, for the variance of the r.v. "number of r-dim subcubes of a random Boolean function with n variables. Try:`):
print(`VARr(n,4): `):

elif nargs=1 and args[1]=Vecsk then
print(`Vecsk(k):{-1,0,1}^k. Try: `):
print(`Vecsk(3);`):

elif nargs=1 and args[1]=Wt then
print(`Wt(M): Given a matrix M represening subcubes of the nops(M[1]) unit cubes, finds the weight 1/2^(Number of Vertices in their Union). Try`):
print(`Wt([[0,-1,-1],[-1,0,-1],[0,-1,-1]]);`):

else

print(`there is no help for`, args):

fi:

end:



ez:=proc(): 
print(` Cov(d1,d2,n) `):
end:


####NUMERIC


#Vecsk(k):{-1,0,1}^k
Vecsk:=proc(k) local gu,gu1:
option remember:

if k=0 then
RETURN({[]}):
fi:
gu:=Vecsk(k-1):

{seq([op(gu1),-1],gu1 in gu),seq([op(gu1),0],gu1 in gu),seq([op(gu1),1],gu1 in gu)}:
end:

#Mats(n,L):  Inputs a positive ineteger n and a list L of non-negative integers of length k, say, outputs all k by n {-1,0,1} matrices 
#where there are exactly L[1] 0's in the first row, exactly L[2] 0's in the second row, ... exactly L[k] 0's in the k-th row. Try:
#Mats(3,[1,1,1]);
Mats:=proc(n,L) local k, gu,lu,L1,gu1,gu11,i,lu1:
option remember:

k:=nops(L):

if min(op(L))<0 then
 RETURN({}):
fi:

if n=0 then
if L=[0$k] then
 RETURN({[[]$k]}):
else
 RETURN({}):
fi:
fi:
 
lu:=Vecsk(k):

gu:={}:

for lu1 in lu do


L1:=[]:
for i from 1 to k do
 if lu1[i]=0 then
  L1:=[op(L1),L[i]-1]:
 else
 L1:=[op(L1),L[i]]:
 fi:
od:

gu1:=Mats(n-1,L1):

for gu11 in gu1 do
gu:=gu union {[seq([op(gu11[i]),lu1[i]],i=1..k)]}:
od:
od:
gu:

end:

#ExpaS(F): Given a subcuble where 0 indicates a wild-card, finds all the vertices beloning to it. Try
#ExpS([1,-1,0,1,0]);
ExpaS:=proc(F) local gu,gu1,F1,n:
option remember:
n:=nops(F):
if n=0 then
 RETURN({[]}):
fi:

F1:=[op(1..nops(F)-1,F)]:

gu:=ExpaS(F1):

if F[n]=0 then
 RETURN({seq([op(gu1),-1],gu1 in gu), seq([op(gu1),1],gu1 in gu)}):
else
 RETURN({seq([op(gu1),F[n]],gu1 in gu)}):
fi:

end:

#Wt(M): Given a matrix M represening subcubes of the nops(M[1]) unit cubes, finds the weight 1/2^(Number of Vertices in their Union). Try
#Wt([[0,-1,-1],[-1,0,-1],[0,-1,-1]]);
Wt:=proc(M) local i:

1/2^nops({seq(op(ExpaS(M[i])),i=1..nops(M))}):

end:





#MomLn(L,n): The L-mixed moment, i.e. E[X[L[1]]*X[L[2]]*...*X[L[k]]), where X[r]:=number of sub-cubes of dimension r.  for Boolean functions on n variables
#Try: MomLn([1,1,1],3);
MomLn:=proc(L,n) local gu,gu1:
gu:=Mats(n,L):
add(Wt(gu1),gu1 in gu):
end:

#Aek1(n,k): Aek's formula for mean
Aek1:=proc(n,k) local i:
simplify(expand(binomial(n,k)*2^(n-k)/2^(2^k))):
end:


#Aek2(n,k): Aek's formula
Aek2:=proc(n,k) local i:
normal(
simplify(
add(
n!/i!/(k-i)!^2/(n-2*k+i)!*2^(n-i)/2^(2^(k+1))*(2^(2^i)-1),i=0..k)+
(binomial(n,k)*2^(n-k))^2/2^(2^(k+1)))
):
end:

#AekM(n,k1,k2): Modification to Aek's formula for mixed moments
AekM:=proc(n,k1,k2) local i:
if(k1>=k2) then error("k1 must be smaller than k2"): fi:
normal(
simplify(
add(
n!/i!/(k1-i)!/(k2-i)!/(n-k1-k2+i)!*2^(n-i)/2^(2^(k1)+2^(k2))*(2^(2^i)-1),i=0..k1)+
(binomial(n,k1)*2^(n-k1))*(binomial(n,k2)*2^(n-k2))/2^(2^(k1)+2^(k2)))
):
end:

##END NUMERIC

#SYMBOLIC
with(combinat):

#Conf1(A,r): Inputs a list of positive integers, of size k, say, outputs ALL lists of length k subsets of {1, ...,r}   [S1,...,Sk] such that Si has cardinality A[i].
#Try: Conf1([2,3,2],5);
Conf1:=proc(A,r) local k,A1,gu,mu,gu1,i1:

k:=nops(A):
if k=0 then
 RETURN({[]}):
fi:

A1:=[op(1..k-1,A)]:
gu:=Conf1(A1,r):
mu:=choose(r,A[k]):

{seq(seq([op(gu1),mu[i1]],i1=1..nops(mu)),gu1 in gu)}:
end:




#Conf(A,r): Inputs a list of positive integers, of size k, say, outputs ALL lists of length k subsets of {1, ...,r}   [S1,...,Sk] such that Si has cardinality A[i] and their union is {1,..,r}.
#Try: Conf([2,3,2],5);
Conf:=proc(A,r) local mu,mu1,gu,i:
mu:=Conf1(A,r) :

gu:={}:

for mu1 in mu do
if {seq(op(mu1[i]),i=1..nops(mu1))}={seq(i,i=1..r)} then
 gu:=gu union {mu1}:
fi:
od:

gu:

end:


#LtoM(L,r): given a list of subsets L of {1,...,r} of size k, that cover {1,..,r}, outputs the 1-0 matrix M (given as a list of lists) such that  M[i][j]=0 iff j belongs to L[i]
#Try:
#LtoM([{1,2},{1,3},{2,3}],3);
LtoM:=proc(L,r) local k,i,j,T,s:
k:=nops(L):

 if {seq(op(L[i]),i=1..k)}<>{seq(i,i=1..r)} then
  RETURN(FAIL):
fi:

for i from 1 to k do
for j from 1 to r do
 T[i,j]:=1:
od:
od:

for i from 1 to k do
 for s in L[i] do
  T[i,s]:=0:
 od:
od:

[seq([seq(T[i,j],j=1..r)],i=1..k)]:

end:





#Trans1(M,pi1,pi2); permutes the rows of M according to pi1 and the columns according to pi2
Trans1:=proc(M,pi1,pi2) local M1,i,j:
 M1:=[seq(M[pi1[i]],i=1..nops(pi1))]:
[seq([seq(M1[i][pi2[j]],j=1..nops(pi2))],i=1..nops(M1))]:
end:

#Imags(M): all the images of the matrix M under permuting the rows and the columns
#Try:
#Imags([[1,0,1],[1,1,1],[1,0,1]]);
Imags:=proc(M) local k,r,gu1,gu2,pi1,pi2:
k:=nops(M):
r:=nops(M[1]):
gu1:=permute(k):
gu2:=permute(r):

{seq(seq(Trans1(M,pi1,pi2),pi1 in gu1),pi2 in gu2)}:

end:

#Reps(A,r): Given a list of positive intgegers A and a positive integer r outputs a set of representatives of all nops(A) by r 0-1 matrices with A[i] 0's in row i. Try:
#Reps([2,3],3);
Reps:=proc(A,r) local mu,mu1,gu,F:
mu:=Conf(A,r):
mu:={seq(LtoM(mu1,r),mu1 in mu)}:

gu:={}:
while mu<>{} do
mu1:=mu[1]:
F:=Imags(mu1) intersect mu:
gu:=gu union {[mu1,nops(F)]}:
mu:=mu minus F:
od:
gu:

end:

#Kids11(M,i,j): Given a matrix M (represented as list of lists) with entries {0,-1,1} and a location [i,j] with M[i][j]=1 outputs the set of
#two matrices where the [i,j] entry is either 1 or of -1. Do
#Kids11([[0,1,0],[1,1,0]],2,1)
Kids11:=proc(M,i,j):

if M[i][j]=0 then
 RETURN(FAIL):
fi:

{M,[op(1..i-1,M),[op(1..j-1,M[i]),-1,op(j+1..nops(M[i]),M[i])],op(i+1..nops(M),M)]}:

end:

#Kids1old(M): Inputs a 0-1 matrix M and outputs the set of all {0,1,-1} matrices M1 such that abs(M1[i][j])=M[i][j], in other words where 1 is re;aced by {1,-1}. Try:
#Kids1old([[1,0],[0,1]]);

Kids1old:=proc(M) local gu,i,j,gu1:
gu:={M}:

for i from 1 to nops(M) do
 for j from 1 to nops(M[i]) do
  if M[i][j]=1 then
   gu:={seq(op(Kids11(gu1,i,j)),gu1 in gu)}:
  fi:
 od:
od:
gu:
end:


#ExpkC(v): Given a vector v of {0,1,-1} finds the set of {1,-1} vectors where 0 can be either 1 or -1. Try:
#ExpkC([1,0,-1,0]);
ExpkC:=proc(v) local gu,k,gu1:
option remember:
if v=[] then
 RETURN({[]}):
fi:

k:=nops(v):
 gu:=ExpkC([op(1..k-1,v)]):
if v[k]<>0 then
 RETURN({seq([op(gu1),v[k]],gu1 in gu)}):
else
RETURN({seq([op(gu1),1],gu1 in gu),seq([op(gu1),-1],gu1 in gu)}):
fi:

end:


#NuVerts1(M): inputs a matrix M in {0,1,-1}, outputs the number of vetrices in the union of subcubes given by its rows. Try:
#NuVerts1([[0,1],[1,0]]);
NuVerts1:=proc(M) local i:
nops({seq(op(ExpkC(M[i])),i=1..nops(M))}):
end:



#NuVerts(M,SP): inputs a matrix M of {0,1,-1} of size k by r, say and a set partition SP of {1,..,k} outputs the total number of vertices in that configuration. Try
#NuVerts([[1,0,-1],[1,1,-1]],{{1},{2}});
NuVerts:=proc(M,SP) local M1,i,k,lu,S,j:

k:=nops(M):

if {seq(op(SP[i]),i=1..nops(SP))}<>{seq(i,i=1..k)} then
 RETURN(FAIL):
fi:


lu:=0:

for i from 1 to nops(SP) do
 S:=SP[i]:
 M1:=[seq(M[S[j]],j=1..nops(S))]:
 lu:=lu + NuVerts1(M1):
od:

lu:

end:

#MomsOld(A,n): Inputs a list A of positive integers of length k, say, and a symbol n, outputs an expression for E[X[A[1]]X[A[2]]...X[A[k]]] where X[r] is the r.v. number of r-dim subcubes of a Boolean function f in n variables. Try
#MomsOld([1,1],n);
MomsOld:=proc(A,n) local gu,k,r,lu,mu,mu1,gu1,i,ku,lu1,gu11,Gu,fu,fu1:
option remember:
k:=nops(A):

 mu:=setpartition({seq(i,i=1..k)}):


Gu:=0:
for r from  max(op(A)) to convert(A,`+`) do
gu:=0:
 lu:=Reps(A,r):


  for i from 1 to nops(lu) do
  lu1:=lu[i][1]:
  ku:=lu[i][2]:

    gu1:=0:

fu:=Kids1old(lu1):

for fu1 in fu do

  for mu1 in mu do


gu11:=expand(binomial(2^(n-r),nops(mu1))*nops(mu1)!/2^NuVerts(fu1,mu1)):
   gu1:=gu1+gu11:

 od:
od:
gu:=gu+ku*gu1:
od:

Gu:=Gu+gu*expand(binomial(n,r))
od:  
factor(Gu):
 

end:



#Kids1(M): Inputs a 0-1 matrix M and outputs the set of all {0,1,-1} matrices M1 such that abs(M1[i][j])=M[i][j], in other words where 1 is re;aced by {1,-1}. Try:
#Kids1([[1,0],[0,1]]);

Kids1:=proc(M) local gu,i,j,gu1:
gu:={M}:

for i from 2 to nops(M) do
 for j from 1 to nops(M[i]) do
  if M[i][j]=1 then
   gu:={seq(op(Kids11(gu1,i,j)),gu1 in gu)}:
  fi:
 od:
od:
gu:
end:


#Moms(A,n): Inputs a list A of positive integers of length k, say, and a symbol n, outputs an expression for E[X[A[1]]X[A[2]]...X[A[k]]] where X[r] is the r.v. number of r-dim subcubes of a Boolean function f in n variables. Try
#Moms([1,1],n);
Moms:=proc(A,n) local gu,k,r,lu,mu,mu1,gu1,i,ku,lu1,gu11,Gu,fu,fu1:
option remember:
k:=nops(A):

 mu:=setpartition({seq(i,i=1..k)}):


Gu:=0:
for r from  max(op(A)) to convert(A,`+`) do
gu:=0:
 lu:=Reps(A,r):


  for i from 1 to nops(lu) do
  lu1:=lu[i][1]:
  ku:=lu[i][2]:

    gu1:=0:


fu:=Kids1(lu1):

for fu1 in fu do

  for mu1 in mu do


gu11:=expand(binomial(2^(n-r),nops(mu1))*nops(mu1)!*2^(r-A[1])/2^NuVerts(fu1,mu1)):
   gu1:=gu1+gu11:

 od:
od:
gu:=gu+ku*gu1:
od:

Gu:=Gu+gu*expand(binomial(n,r))
od:  
factor(Gu):
 

end:


#Cov(d1,d2,n): The Covariance betweed X[d1] and X[d2]
Cov:=proc(d1,d2,n) local mu1,mu2,mu12:
option remember:
if d1=d2 then
  mu12:=Aek2(n,d1):
else
  mu12:=AekM(n,d1,d2):
fi:
mu1:=Aek1(n,d1):
mu2:=Aek1(n,d2):
#mu1:=Moms([d1],n):
#mu2:=Moms([d2],n):
mu12-mu1*mu2:
#normal(mu12-mu1*mu2):
end:


#Cor(d1,d2,n): The Correlation betweed X[d1] and X[d2]
Cor:=proc(d1,d2,n) local gu,mu:
option remember:
gu:=Cov(d1,d2,n):
mu:=sqrt(Cov(d1,d1,n)*Cov(d2,d2,n)):
#simplify(gu/mu,symbolic):
simplify(gu/mu):
end:

#VARr(n,r): The expression, in n, for the variance of the r.v. "number of r-dim subcubes of a random Boolean function with n variables. Try:
#VARr(n,4)
VARr:=proc(n,r) 
factor(simplify(expand(normal(Aek2(n,r)-(2^(n-r)*binomial(n,r)/2^(2^r))^2)))):
end:


#MOMrk1(r,k,n): The k-th moment of X[r]. Try:
#MOMrk1(1,3,n);
MOMrk1:=proc(r,k,n) option remember:
Moms([r$k],n):
end:

#MOMrk(r,k,n): The k-th moment about the mean of X[r]. Try:
#MOMrk(1,4,n);
MOMrk:=proc(r,k,n) local x,mu,f,i: 
mu:=2^(n-r)*binomial(n,r)/2^(2^r):
f:=expand((x-mu)^k):

factor(normal(
add(coeff(f,x,i)*MOMrk1(r,i,n),i=1..degree(f,x))+coeff(f,x,0)
)
):
end:



#SMOMrk(r,k,n): The SCALED k-th moment about the mean of X[r]. Try:
SMOMrk:=proc(r,k,n):
factor(simplify(MOMrk(r,k,n)/VARr(n,r)^(k/2))): 
end:



#END SYMBOLIC