# restart; read `c:/Users/Asus/Google Drive
# /Aek/MomentReview/Boolean.txt`;
# read `c:/Users/Asus/Google Drive/Guess.txt`;

with(combinat,powerset):
with(combinat,choose):
with(combinat,stirling1):
with(combinat,stirling2):

#print(`Welcome to Boolean, a Maple package`):
#print(`written by Thotsaporn "Aek" Thanatipanonda`):
#print(`First Version: November 13, 2013 `):
#print(`Second Version: April 15, 2018`):
#print(`Third Version: March 7, 2020`):

Help := proc()
print(`Part 1: Boolean Function`);
print(``);
print(`Section 1: Direct Computation`);
print(`The functions are`):
print(`Boolean(f,k,n),`);
print(`Insert1(S,p), ZO(n),`);
print(`Mo(n,k,r), MoMean(n,k,r),`);
print(`GenBoo(n,k,q), TestMo(n,k),`);
print(`SubGenBoo(n,m,k,q), `);
print(``);
print(`GFInd(n,k,q), `);
print(`TestPlot(n,k), TestGFInd(k,N),`);
print(`MoGFk1(r,n), MoMeanGFk1(r,n),`);
print(`MoGFk2(r,n), MoMeanGFk2(r,n),`);
print(`TestGen(r),`);
print(``);

print(`Section 2: Overlapping Stage Approach`);
print(`The functions are`):
print(`OSA(n,k,r), Overlap(n,k,r), `);
print(`Squares(n,k), Inserts(A,B,P),`);
print(`InsertOne(A,B,P)`);
print(``);
print(`TestOverlap(n,k,r,s)`);
print(``);

print(`Section 3: Normality for k=0`);
print(`The functions are`):
print(`Section 3.1:  `);
print(`generating functions`);
print(`Genk0(n,q), FallMok0(r,n),`);
print(`ToMo(r,n), InverMo(r,n),`);
print(`GenMeank0(n,q), FallMoMean(r,n),`);
print(`ToMoMean(r,n), ForFallMM(r,n),`);
print(``);
print(`Section 3.2: Moment Calculations`);
print(`Dim0For(r,n), Dim0MeanFor(r,n),`);
print(`InSum1(U,L), OutSum1(r,t),`);
print(`MyDim0(r,n), MyDim0Mean1(r,n), `);
print(`InSum2(r,t,u), OutSum2(r,t),`);
print(`MyDim0Mean2(r,n),`);
print(``);
print(`NewId(t,r)`);
print(``);

print(`Section 4: Formulas for all k`);
print(`The functions are`):
print(`ForMo1(n,k), ForMo2(n,k),`);
print(`ForVar(n,k), ForMo3k1(n)`);
print(``);
print(`Section 4.1: Helper functions`);
print(`Binom(n,k), MulNo(n,S),`);
print(`Fall(a,k),`);
print(``);

print(`Section 5: Work with independence`);
print(`MyH`);
print(``);
print(``);
print(``);
print(``);
print(``);





print(`Part 2: Boolean Board`);
print(``);
print(`Section 1: Empirical Counting`);
print(`Functions:`);
print(`ListBB(m,n), BBoard(s,m,n), `);
print(`CountD(T,m,n), `);
print(``);
print(`Section 2: Counting the rth-moment`);
print(`Functions:`);
print(`MoDoBo(r,m,n), MoMeDoBo(r,m,n),`);
print(`FactMeDoBo(r,m,n), GenBB(m,n,q),`);
print(`ForMoBB(m,r),`);
print(``);

print(`Section 3: Formula for moments 1-by-n case`);
print(`Functions:`);
print(`CenGenBB1(n,q), RecCenGenBB1(n,q),`);
print(`MyTaylor(r,n), LeadTermBB1(r),`);
print(`SecLeadTermBB1(r),`);
print(``);
print(`Section 4: Case 2-by-n board`);
print(`Functions:`);
print(`GenA(n,q), GenB(n,q),`);
print(``);
end:





####################
# Part 1 
# Boolean Function
####################

#################################
# Section 1: Direct Computation
#################################

# Input: Boolean S and the size k 
# that I want to count.
# Output: number of squares of size k of f.
# Try: Boolean({[0,0,0],[0,1,0],[1,1,1]},1,3);
Boolean := proc(f,k,n) option remember;
local i,c,S,SS,C,p,count;

count := 0;
C := choose(n,k);
SS := {};

for c in C do
	p := sort([op({seq(i,i=1..n)} minus {op(c)})]);
	SS := SS union Insert1(ZO(k),p);
od:

for S in SS do
	if S subset f then
		count := count+1;
	fi:
od:
count;
end:


# Try: Insert1({[0],[1]},[1]);
Insert1 := proc(S,p) option remember;
local s,a,T1,T2;

if p = [] then
	return({S});
fi:

a := p[1];
T1 := {seq([op(s[1..a-1]),
0,op(s[a..nops(s)])], s in S)};
T2 := {seq([op(s[1..a-1]),
1,op(s[a..nops(s)])], s in S)};

Insert1(T1,p[2..nops(p)]) union
Insert1(T2,p[2..nops(p)]);
end:


# Input: length n
# Output: all 0,1 of length n.
# Try: ZO(2);
ZO := proc(n) option remember;
local s;

if n = 0 then
	return({[]});
fi:

{seq(op([[op(s),0],[op(s),1]])
, s in ZO(n-1))};
end:


#####################################
# Average on sample space of length n
# of Boolean^r by brute force.
# Try:
# [seq(Mo(n,0,1),n=1..4)]; 
# seq(Mo(n,1,2),n=1..5);
Mo := proc(n,k,r) option remember;
local f, P;

P := powerset(ZO(n));
add(Boolean(f,k,n)^r, f in P)/nops(P);
end:

# Average on sample space of length n
# of (Boolean-E[x])^r by brute force.
# Try: seq(print(MoMean(n,1,1)),n=1..5);
MoMean := proc(n,k,r) option remember;
local f, P;

P := powerset(ZO(n));
add((Boolean(f,k,n)-Mo(n,k,1))^r,
f in P)/nops(P);
end:



#############################################
# Input: the size k of cube, 
# statement of length n and variable q
# Output: Generating function in q
# Try: 
# n:=4: A:=GenBoo(n,0,q);
# B:=add( binomial(2^n,i)*q^i,i=0..2^n)/2^(2^n):
# A-B;
# [seq(GenBoo(n,1,q),n=1..4)];
GenBoo := proc(n,k,q) option remember;
local f,P,A;

P := powerset(ZO(n));
A := add( q^Boolean(f,k,n), f in P)/nops(P);
#sort(nops(P)*A);
end:


# Test generating functions with moments.
# Input: the size k of cube 
# and statement of length n 
# Try: TestMo(1,0);
TestMo := proc(n,k) option remember;
local i,L,R,f,r;

R := [seq(Mo(n,k,r),r=0..10)];
L := [];
f := GenBoo(n,k,q);

for i from 0 to 10 do
	L := [op(L),subs(q=1,f)];
	f := q*diff(f,q);
od:

[R,L-R];
end:

# Difficult to find formula
# Input: m statements of length n, 
# the size k of cube, and variable q
# Output: Generating function in q
# Try: [seq(SubGenBoo(2,m,1,q),m=0..4)];
SubGenBoo := proc(n,m,k,q) option remember;
local f,P,A;

P := choose(ZO(n),m);
A := add( q^Boolean(f,k,n), f in P)/nops(P);
sort(nops(P)*A);
end:


###################################################
# Generatig function by 
# assume independence
# Try: GFInd(2,1,q);
GFInd := proc(n,k,q) option remember;
local m,p,M,A;

p := binomial(n,k)*(2^k)!/2^(k*2^k)/2^((n-k)*(2^k-1));
M := m -> binomial(m,2^k);

A:= (1/2)^(2^n)
*add( binomial(2^n,m)*(p*q+1-p)^M(m)  
,m=0..2^n);

sort(simplify(A));
end:

# Try: TestPlot(4,1);
TestPlot := proc(n,k) option remember;
local i,q,A,B,S,T;

A := GenBoo(n,k,q);
S := [seq(coeff(A,q,i),i=0..degree(A,q))];
B := GFInd(n,k,q);
T := [seq(coeff(B,q,i),i=0..degree(A,q)+3)];

with(Statistics):
Statistics:-ColumnGraph([S,T]);
end:

# Try: TestGFInd(1,6);
TestGFInd := proc(k,N) option remember;
local i,n,q,A,S;

for n from N-2 to N do 
	A := GFInd(n,k,q);
	S[n] := [seq(coeff(A,q,i),i=0..130)];
od:

with(Statistics):
Statistics:-ColumnGraph([seq(S[n],n=N-2..N)]);
end:


# Try: [seq(MoGFk1(r,n),r=0..3)];
MoGFk1 := proc(r,n) option remember;
local i,m,p,q,F,Mo;

F:=1/2^(2^n)*sum( binomial(2^n,m)
*(p*q+1-p)^Binom(m,2) ,m=0..2^n);

for i from 1 to r do
	F := q*diff(F,q);
od:	

Mo := eval(subs(q=1,F));
simplify(subs(p=n/2^n,Mo));
end:


# Try: [seq(MoMeanGFk1(r,n),r=0..3)];
MoMeanGFk1 := proc(r,n) option remember;
local mu1,A;

mu1 := MoGFk1(1,n);
A:=add(binomial(r,i)*MoGFk1(i,n)*(-mu1)^(r-i)  ,i=0..r);
sort(factor(simplify(A)));
end:


# Try: [seq(MoGFk2(r,n),r=0..3)];
MoGFk2 := proc(r,n) option remember;
local i,m,p,q,F,Mo;

F:=1/2^(2^n)*sum( binomial(2^n,m)
*(p*q+1-p)^Binom(m,4) ,m=0..2^n);

for i from 1 to r do
	F := q*diff(F,q);
od:	

Mo := eval(subs(q=1,F));
simplify(subs(p=Binom(n,2)*3!/2^(3*n),Mo));
end:


# Try: [seq(MoMeanGFk2(r,n),r=0..3)];
MoMeanGFk2 := proc(r,n) option remember;
local mu1,A;

mu1 := MoGFk2(1,n);
A:=add(binomial(r,i)*MoGFk2(i,n)*(-mu1)^(r-i)  ,i=0..r);
sort(factor(simplify(A)));
end:


# Test higher moment abou the mean
# of boolean function with k=1.
# This time in generating function form
# Try: TestGen(2);
TestGen := proc(r) option remember;
local n,q,i,m,p,mu,G,A;

mu := n*(2^n-1)/8;

G :=1/2^(2^n)*sum( binomial(2^n,m)
*(p*q+1-p)^Binom(m,2) ,m=0..2^n)/q^mu;

for i from 1 to r do
	G := q*diff(G,q);
od:	

A:= eval(subs(q=1,G));
factor(simplify(subs(p=n/2^n,A)));
end:



########################################
# Section 2: Overlapping Stage Approach
######################################## 


# Output: r_th moment of boolean function.
# Try: OSA(3,1,1);
# [seq(OSA(n,0,2),n=1..10)];
# seq([seq(OSA(n,k,2),k=0..n)],n=1..5);
# seq([seq(OSA(n,k,3),k=0..n)],n=1..4);
# seq([seq(OSA(n,k,4),k=0..n)],n=1..4);
OSA := proc(n,k,r)
local a,s,t,TT;

TT := [0$(r*2^k)];

for a in Overlap(n,k,r) do
	s := nops(`union`(seq(t, t in a)));  
	TT[s] := TT[s]+1;
od:

add( 2^(2^n-s)*TT[s],s=1..nops(TT))/2^(2^n);
end:


# Try: Overlap(2,1,1);
Overlap := proc(n,k,r) option remember;
local a,s,t;

if r =1 then
	return([seq({a}, a in Squares(n,k))]);
fi:

[seq(seq(s union t
,s in Overlap(n,k,1)),t in Overlap(n,k,r-1))];
end:


# Output: all the squares of length k 
# in word of length n.
# Try: Squares(2,1);
Squares := proc(n,k)
local a,PP,P;

PP := choose(n,k); 
{seq(seq(Inserts(a,ZO(k),P)
, P in PP),a in ZO(n-k))};
end:


# Output: Insert each memeber of B 
# at position P of A.
# Try: Inserts([10,20,30],{[1,2],[7,6]},[4,5]);
Inserts := proc(A,B,P)
local b,T;

T := sort(P,`<`);
{seq( InsertOne(A,b,T),b in B)};
end:

# Output: Insert each element of B 
# at position P of A.
# Try: InsertOne([10,20,30],[1,2],[4,5]);
InsertOne := proc(A,B,P)
local i,ret;

ret := A;

for i from 1 to nops(P) do
	ret := [op(1..P[i]-1,ret),B[i],
		   op(P[i]..nops(ret),ret)]; 
od:

ret;
end:


# Testing Formula
# See some examples.
# Input: size of the overlap.
# Try: TestOverlap(2,1,3,4);
TestOverlap := proc(n,k,r,s)
local a,t,ret;

ret := [[]$r];

for a in Overlap(n,k,r) do 
	if s = nops(`union`(seq(t, t in a))) then 
		ret[nops(a)] := [op(ret[nops(a)]),a];
	fi:
od:

print(ret);
[seq(nops(ret[i]),i=1..r)];
end:







########################################
# Section 3: Proving normality for k=0
########################################

# Section 3.1:
# Generating Functions


# Generating functions of this case k=0
# Try: Genk0(5,q);
Genk0 := proc(n,q) option remember;
local i;

# expand(((1+q)/2)^(2^n));
add( binomial(2^n,i)*q^i/2^(2^n),i=0..2^n);
end:


# Binomial moment of this case k=0
# Try: FallMok0(2,n);
FallMok0 := proc(r,n) option remember;
local i;

mul(2^n-i,i=0..r-1)/2^r;
end:

# Convert the falling moment to straight moment 
# Try: seq(ToMo(r,n),r=1..4);
ToMo := proc(r,n) option remember;
local i,A;

A := add( stirling2(r,i)*FallMok0(i,n)  ,i=0..r);
A := factor(A);
[A,Dim0For(r,n)-A];
end:


# Testing the Inversion formula 
# Convert the straight moment to binomial moment
# Try: InverMo(3,n);
InverMo := proc(r,n) option remember;
local L,R;

L := BiMok0(r,n)*r!;
R := add(stirling1(r,i)*Dim0For(i,n),i=0..r);

[L, expand(simplify(L-R))];
end:

#########################################################3
# Generating functions about 
# the mean of this case k=0.
# Try: GenMeank0(5,q);
GenMeank0 := proc(n,q) option remember;
local i,A,B;

A:=expand(((q^(-1/2)+q^(1/2))/2)^(2^n));
#B:=add( binomial(2^n,i)*q^(i-2^(n-1))/2^(2^n),i=0..2^n);
#[A,B];
end:

# Try: [seq(FallMoMean(r,n),r=0..5)];
FallMoMean := proc(r,n) option remember;
local A;

A := GenMeank0(n,q);

if r=0 then
	subs(q=1,A);
else	
	subs(q=1,diff(A,q$r));
fi:
end:

# Convert the falling moment to straight moment 
# Try: seq(ToMoMean(r,n),r=0..4);
ToMoMean := proc(r,n) option remember;
local i,A;

A := add( stirling2(r,i)*FallMoMean(i,n),i=0..r);
A := factor(A);
[A,Dim0MeanFor(r,n)-A];
end:



#########################################
# Section 3.2: Moment Calculations

# Output: Formula for k=0 cube with 
# numeric r moment  and symbolic n
# Try: seq(Dim0For(r,n),r=1..6);
Dim0For := proc(r,n) option remember;
local i;

factor(add( stirling2(r,i)
*Binom(2^n,i)*i!/2^i, i=0..r));
end:


# Input: numeric r moment and symbolic n
# Output: Moment about the mean for k=0
# Try: seq(Dim0MeanFor(r,n),r=1..6);
Dim0MeanFor := proc(r,n) option remember;
local i,A;

A := add(  (-1)^(r-i)*binomial(r,i)*Dim0For(i,n)
*Dim0For(1,n)^(r-i),i=0..r);

factor(A);
end:
##################################################


# inner sum of the first way
# [seq(InSum1(30,L),L=0..30)];
InSum1 := proc(U,L) option remember;
local j;

add( stirling2(U,j)*stirling1(j,L)/2^j  ,j=L..U);
end:


# Coeff of 2^(n*(r-t)) in E[(X-mu)^r]
# outer sum of the second way
# [seq(OutSum1(30,t),t=0..30)];
# [seq(OutSum1(30,t)-OutSum2(30,t),t=0..30)];
OutSum1 := proc(r,t) option remember;
local u;

add( (-1/2)^(r-i)*binomial(r,i)*InSum1(i,i-t), i=t..r);
end:


# Output: List of coeff of 2^(n*k) in E[X^r]
# Try: 
# [seq(simplify(MyDim0(r,n)-Dim0For(r,n)),r=0..10)];
MyDim0 := proc(r,n) option remember;
local k;
add( InSum1(r,k)*2^(n*k),k=0..r);
end:


# Formula for Dim0MeanFor
# Try:
# [seq(simplify(MyDim0Mean1(r,n)-Dim0MeanFor(r,n)),r=0..20)];
MyDim0Mean1 := proc(r,n) option remember;
local t;
add( OutSum1(r,t)*2^(n*(r-t)),t=0..r);
end:


# inner sum of the second way
# [seq(InSum2(30,20,j),j=0..21)];
InSum2 := proc(r,t,j) option remember;
local i;

add( (-1/2)^(r-i)*binomial(r,i)*stirling2(i,j)
*stirling1(j,i-t), i=max(t,j)..j+t);

end:


# Coeff of 2^(n*(r-t)) in E[(X-mu)^r]
# outer sum of the second way
# [seq(OutSum2(30,t),t=0..30)];
OutSum2 := proc(r,t) option remember;
local u;

add( InSum2(r,t,j)/2^j, j=0..r);
end:


# Try:
# [seq(simplify(MyDim0Mean2(r,n)-Dim0MeanFor(r,n)),r=0..20)];
MyDim0Mean2 := proc(r,n) option remember;
local i,t,u;
add( OutSum2(r,t)*2^(n*(r-t)),t=0..r);
end:


########################################################

# Input: Numeric t and r where 0<=t<=r.
# Output: coeff. of 2^(n*(r-t)) of E[X^r]
# Try: r:=37:[seq(NewId(t,r),t=0..r)];

NewId := proc(t,r) option remember;
local i,j;

sum( (-1/2)^(r-i)*binomial(r,i)*
sum(1/2^j*stirling2(i,j)*stirling1(j,i-t)
,j=i-t..i),i=t..r)   ;

end:



################################
# Section 4: Formulas for all k
################################

# Formula for the first moment
# Try: seq(ForMo1(n,1),n=1..4);
# seq(Mo(n,1,1),n=1..4);
ForMo1 := proc(n,k)
simplify(Binom(n,k)*2^(n-k)/2^(2^k));
end:


# Formula for the second moment
# Try: seq(ForMo2(n,1),n=1..4);
# seq(Mo(n,1,2),n=1..4);
# ForMo2(n,1);
ForMo2 := proc(n,k) option remember;
local i;

simplify(
add(MulNo(n,[i,k-i,k-i,n-2*k+i])
*2^(n-i)*(2^(2^i)-1)/2^(2^(k+1)),i=0..k)
+(Binom(n,k)*2^(n-k))^2/2^(2^(k+1))
);

end:


# Input: symbolic n and numeric k
# Try: ForVar(n,1);
# factor(ForMo2(n,1)-ForMo1(n,1)^2);
# seq(simplify(ForMo2(n,k)
# -ForMo1(n,k)^2-ForVar(n,k)),k=0..5);
ForVar := proc(n,k) option remember;
local i;

simplify(
add(MulNo(n,[i,k-i,k-i,n-2*k+i])
*2^(n-i)*(2^(2^i)-1)/2^(2^(k+1)),i=0..k)
);

end:

# Third moment and k=1
# Try: seq(ForMo3k1(n),n=1..6);
# [seq(Mo(n,1,3),n=1..4)]; 
ForMo3k1 := proc(n)
local i,a,b,c,S;

S := [0$6];
a := 2^(n-1)*n;
b := 2^n*n*(n-1)*3;
c := 2^n*n*(n-1)*(n-2)*4;

# left, #
S[2] := a;
S[3] := b;
S[4] := c+3*a*(a-1);
S[5] := b*(a-3*n+2);
S[6] := a^3-S[2]-S[3]-S[4]-S[5];

simplify(add(S[i]/2^i, i=1..6));
end:


#########################################
# Section 4.1: Helper functions

# Binomial with symbolic n
# Try: Binom(n,3);
Binom := proc(n,k) option remember;
local i;
mul(n-i,i=0..k-1)/k!;
end:

# Output: Multinomial
# Try: MulNo(5,[2,1,2]);
# MulNo(n,[2,n-2]);
MulNo := proc(n,S)
local i,T,s;

T := [seq(S[i], i=1..nops(S)-1)];

if add(s, s in S) <> n then
	ERROR();
fi:

for s in T do
	if s < 0 then
		return(0);
	fi:
od:

mul(n-i,i=0..n-S[nops(S)]-1)/mul(s!,s in T);
end:


# Falling factor
# Try: Fall(a,3);
Fall := proc(a,k) option remember;
local i;
mul(a-i,i=0..k-1);
end:


######################################
# Section 5: Work with independence



# Input: numeric k and n (k<=n)
# and symbolic p,q
# Output: generating functions
# Try: MyH(1,3,p,q);
MyH := proc(k,n,p,q) option remember;
local m;

1/2^(2^n)*add(  binomial(2^n,m)*(p*q+1-p)^binomial(m,2^k)  
,m=0..2^n);
end:


# Try: Myp(1,3);
Myp := proc(k,n) option remember;
binomial(n,k)*(2^k)!/2^(k*2^k)/2^((n-k)*(2^k-1));

end:

# Try: Mean1(1,3);
Mean1 := proc(k,n) option remember;
binomial(n,k)*2^(n-k)/2^(2^k);
end:







######################################
# Part 2 
# Dominos on a Boolean Board
######################################

# Section 1: Empirical Counting

# Input: Board of width m and length n
# Output: List L = [a_i],
# a_i is the number of boards that contained
# exactly i dominos.
# Try: ListBB(2,1);
# [seq(ListBB(2,i),i=1..7)];
ListBB := proc(m,n) option remember;
local i,B,C,S;

S := [0$(2*m*n-m-n+1)];

for i from 0 to 2^(m*n)-1 do
	B := BBoard(i,m,n); 
	C := CountD(B,m,n);
	S[C+1] := S[C+1]+1;
od:

S;
end:


# Input: number s, size of the board m.n
# Output: the s_th boolean board of size m.n
# Try: [seq(BBoard(s,2,1),s=0..3)];
BBoard := proc(s,m,n) option remember;
local i,j,ns,T;

if s >= 2^(m*n) then
	ERROR("s is out of range");
fi:
	
ns := s;
T := [];

for i from m*n-1 by -1 to 0 do
	T := [op(T), floor(ns/(2^i))];
	ns := ns - T[nops(T)]*2^i;
od:

Matrix( [seq([seq(T[i+m*j],i=1..m)],j=0..n-1)]);
end:



# Input: Board T of size m-by-n
# Output: the number [1,1] or [0,0]
# (next to each other vertically, horizontally)
# in the board
# Try: T := BBoard(34,3,2);
# CountD(T,3,2);
CountD :=proc(T,m,n) option remember;
local i,j,count;

count := 0;

for i from 1 to n-1 do
	for j from 1 to m do
		if T[i,j]=T[i+1,j] then
			count := count+1;
		fi:
	od:
od:

for i from 1 to n do
	for j from 1 to m-1 do
		if T[i,j]=T[i,j+1] then
			count := count+1;
		fi:
	od:
od:

count:	
end:



##############################################
# Section2: Counting the rth-moment

# Output: r-th moment of number of domino
# from the boolean board of size m-by-n
# Try: [seq(MoDoBo(1,2,i),i=1..7)];
# [seq(2^(2*i-1)*(3*i-2)/2^(2*i),i=1..7)];
MoDoBo := proc(r,m,n) option remember;
local i,tol,LBB;

LBB := ListBB(m,n);

tol := add(LBB[i],i=1..nops(LBB));
add(i^r*LBB[i+1],i=0..nops(LBB)-1)/tol;
end:


# Output: r-th moment about the mean 
# of number of dominos from the 
# boolean board of size m-by-n
# Try: [seq(MoMeDoBo(2,2,i),i=1..7)];
# [seq( MoDoBo(1,2,i)/2,i=1..7)];
# [seq( MoMeDoBo(2,1,n),n=1..12)];
MoMeDoBo := proc(r,m,n) option remember;
local i,tol,me,LBB;

LBB := ListBB(m,n);

tol := add(LBB[i],i=1..nops(LBB));
me := MoDoBo(1,m,n);
add((i-me)^r*LBB[i+1],i=0..nops(LBB)-1)/tol;
end:


# Output: Factorial moment about the mean 
# of number of dominos from the 
# boolean board of size m-by-n
# Try: [seq(FactMeDoBo(2,2,i),i=1..7)];
FactMeDoBo := proc(r,m,n) option remember;
local i,tol,me,LBB;

LBB := ListBB(m,n);

tol := add(LBB[i],i=1..nops(LBB));
me := MoDoBo(1,m,n);
add( mul(i-me-j,j=0..r-1)*LBB[i+1],i=0..nops(LBB)-1)/tol;
end:


# Output: Generating function from
# the board of size m-by-n
# Try: seq(print(GenBB(1,i,q)),i=1..10);
GenBB := proc(m,n,q) option remember;
local i,tol,A;

A :=ListBB(m,n);
tol := add(A[i],i=1..nops(A));

factor(add( A[i]*q^(i-1),i=1..nops(A)))/tol;
end:


# Fomula of E[X^r] in term of mu
# Try: [seq(ForMoBB(m,r),r=1..5)];
ForMoBB := proc(m,r) option remember;
local i,a,A;

A:=add(Fall(a,i)/2^i*stirling2(r,i),i=1..r);
factor(subs(a=2*m,A));

end:

##############################################
# Section 3: Formula for moments 1-by-n case
# Functions:
# CenGenBB1(n,q), RecCenGenBB1(n,q),
# MyTaylor(r,n), LeadTermBB1(r),
# SecLeadTermBB1(r),

# Try: seq(CenGenBB1(n,q),n=1..5);
CenGenBB1 := proc(n,q) option remember;
factor(GenBB(1,n,q)/q^((n-1)/2));
end:


# Try: seq(RecCenGenBB1(n,q),n=1..5);
RecCenGenBB1 := proc(n,q) option remember;
if n=1 then
	return(1);
fi:	

factor((q^(1/2)+q^(-1/2))/2*RecCenGenBB1(n-1,q));
end:


# Try: MyTaylor(6,3);
MyTaylor := proc(r,n) option remember;
local s,t,A;

A:=taylor(RecCenGenBB1(n,exp(t/s)),t=0,r+1);
coeff(A,t^r)*s^r;

end:


# Try: LeadTermBB1(10);
LeadTermBB1 := proc(r) option remember;
local i,t,n,A,B;

A:=[seq(MyTaylor(r,i),i=1..r+5)];
B := expand(GuessPol(A,1,n));

coeff(B,n,degree(B,n));
end:


# Try: SecLeadTermBB1(10,n);
SecLeadTermBB1 := proc(r) option remember;
local i,t,n,A,B;

A:=[seq(MyTaylor(r,i),i=1..r+5)];
B := expand(GuessPol(A,1,n));

coeff(B,n,degree(B,n)-1);
end:


##############################################
# Section 4: Case 2-by-n board

# Generating function when ends
# with the same last vertical numbers 
# Try: [seq(GenA(n,q),n=1..5)];
# [seq(expand(GenBB(2,n,q))*2^(2*n),n=5)];
# [seq(GenA(n,q)+GenB(n,q),n=5)];
GenA := proc(n,q) option remember;
local A;
if n=1 then 
	return(2*q);
fi:

A:=GenA(n-1,q)*q*(1+q^2)+GenB(n-1,q)*2*q^2;	
expand(A);
end:

# Generating function when ends with 
# the different last vertical numbers 
GenB := proc(n,q) option remember;
local B;
if n=1 then 
	return(2);
fi:

B:=GenA(n-1,q)*2*q+GenB(n-1,q)*(1+q^2);	
expand(B);
end: