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

print(`Welcome to Domino, a Maple package`):
print(`written by Thotsaporn (Aek) Thanatipanonda`):
print(`Version: June 5, 2018 `):
print(``):
print(`The program is a implement of the paper`):
print(`"Statistics of Domino Tilings`):
print(`on a Rectangular Board"`):
print(``):


####################################
# Section 1:  
# Number of Vertical dominos 
# on tiling of an 2-by-n board
# (numerically mostly)
#
# Functions:
# GenV2(n,v), SumV(r,n), MoV(r,n),  
# SumMeanV(r,n), MoMeanV(r,n), 
#
# CenGenV2(n,v), RecCenGenV2(n,v),
# NewMoMeanV(r,n)
####################################

# Input: the length of the board n
# and variable v
# Output: Generating function 
# of V on an 2-by-n board
# 
# Try: 
# seq(GenV2(n,v),n=1..10);
# seq(GenHV(2,n,v,1),n=1..10);
GenV2 := proc(n,v) option remember;
local A;

if n = 0 then
	return(1);
elif n=1 then
	return(v);
fi:	

A := v*GenV2(n-1,v)+GenV2(n-2,v);
factor(expand(A));
end:


# Input: the r-th Sum  and 
# the length of the board n
# Output: the r-th Sum  
# of V on 2-by-n board
# Try: [seq(SumV(2,n),n=1..10)];
# [seq(SumPowHV(2,0,2,n),n=1..10)];
SumV := proc(r,n) option remember;
local i,v,A;

A := GenV2(n,v);

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

subs(v=1,A);
end:


# Input: the r-th moment 
# the length of the board n
# Output: the r-th moment 
# of V on 2-by-n board
# Try: seq(MoV(1,n),n=1..10);
MoV := proc(r,n) option remember;
SumV(r,n)/SumV(0,n);
end:


# Input: the r-th Sum 
# about the mean and the 
# length of the board n
# Output: the r-th Sum about the 
# mean of V on 2-by-n board
# computed from straight Sum.
# Try: [seq(SumMeanV(2,n),n=1..10)];
SumMeanV := proc(r,n) option remember;
local i,mu,A;

mu := SumV(1,n)/SumV(0,n);

add(  binomial(r,i)*SumV(i,n)*(-mu)^(r-i) ,i=0..r);
end:



# Input: the r-th moment about
# the mean and the length 
# of the board n
# Output: the r-th moment about
# the mean of V on 2-by-n board
# computed directly from the
# straight moments.
# Try: [seq(MoMeanV(2,n),n=1..10)];
# [seq(SumMeanV(2,n)/SumV(0,n),n=1..10)];
MoMeanV := proc(r,n) option remember;
local i,mu,A;

mu := MoV(1,n);
add(  binomial(r,i)*MoV(i,n)*(-mu)^(r-i) ,i=0..r);
end:

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

# Input: the length of the board n
# and variable v
# Output: Centralized Generating 
# function of V on an 2-by-n board
# computed directly from GenV2
#
# Try: 
# seq(CenGenV2(n,v),n=0..5);
CenGenV2 := proc(n,v) option remember;
GenV2(n,v)/SumV(0,n)/v^MoV(1,n);
end:


# Input: the length of the board n
# and variable v
# Output: Centralized Generating 
# function of V on an 2-by-n board
# computed from the recurrence
# Try: 
# seq(RecCenGenV2(n,v),n=0..5);
# seq(CenGenV2(n,v),n=0..5);
RecCenGenV2 := proc(n,v) option remember;
local A;
if n = 0 or n=1 then
	return(1);
fi:	

A:= 1/SumV(0,n)*(
v*SumV(0,n-1)*RecCenGenV2(n-1,v)/v^(MoV(1,n)-MoV(1,n-1))
+SumV(0,n-2)*RecCenGenV2(n-2,v)/v^(MoV(1,n)-MoV(1,n-2)));

factor(A);
end:


# Comment: Double check MoMeanV 
# with another one from 
# generating function

# Input: the r-th moment about
# the mean and the length 
# of the board n
# Output: the r-th moment about
# the mean of V on 2-by-n board
# computed from generating function

# Try: [seq(NewMoMeanV(2,n),n=1..10)];
# [seq(MoMeanV(2,n),n=1..10)];
NewMoMeanV := proc(r,n) option remember;
local i,v,A;

A := CenGenV2(n,v);

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

subs(v=1,A);
end:



#################################
# Section 2: 2-by-n board
# Grand Generating Functions, 
# Conjectures of moments
# 
# Functions:
# GrandGen2(r,N,t), FastSumV(r,N),
# ConjSumV(r,n,F), ConjMoV(r,n,F),
# BigConjMoV(k,n,r),
# ConjMoMeanV(r,n),
#
# ForMoV(r,n)
#################################

# Input: non-negative integer r and N
# symbolic t
# Output: the first N terms of 
# the r-th moment computed from 
# the grand generating function
# of 2-by-n board
# Try: 
# GrandGen2(2,10,t);
# [seq(SumV(3,n),n=1..10)];
GrandGen2 := proc(r,N,t) option remember;
local i,v,f,A;

f := 1/(1-v*t-t^2);

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

A:=taylor(f,t,N+1);
subs(v=1,A);
end:


# Input: the r-th Sum  and 
# the length of the board N
# Output: the first N term 
# of the r-th Sum of V 
# on 2-by-n board
# Try: 
# FastSumV(3,80);
# [seq(SumV(3,n),n=0..80)];
FastSumV := proc(r,N) option remember;
local t,T;

T := GrandGen2(r,N,t);

[seq(coeff(T,t,i),i=0..N)];
end:


# Input: non-negative integer r,
# symbolic n and F
# Output: The conjectured formula
# of SumV in term of n and F_n
# Try: 
# ConjSumV(2,n,F);
ConjSumV := proc(r,n,F) option remember;
local nn,j,deg,S,T,eq,var;

T := FastSumV(r,2*r+5);
S := FastSumV(0,2*r+5);
deg := r;

eq := { seq( T[nn] =  
add(add( a[i,j]*(nn-1)^i,i=0..deg)
*S[nn-j],j=0..1 )
,nn=3..2*r+5)};
var := {seq(seq(a[i,j],i=0..deg),j=0..1)};

var:=solve(eq,var);

add(factor(subs(var,add( a[i,j]*n^i,i=0..deg)))
*F[n-j],j=0..1 );

end:


# Input: non-negative integer r,
# symbolic n 
# Output: The conjectured formula
# of the r-th moment of V in term of n 
# Try: 
# ConjMoV(2,n);
ConjMoV := proc(r,n) option remember;
local i,F,A,S;

A := expand(ConjSumV(r,n,F)/F[n]);

S := {F[n-1]=1,F[n]=(1+sqrt(5))/2};
A := simplify(expand(subs(S,A)));

add( factor(coeff(A,n,i))*n^i,i=0..degree(A));
end:


# Input: non-negative integer k,
# symbolic n and r
# Output: the k-th biggest term
# of Moment V-formula as variables
# in n and r
# Try: BigConjMoV(2,n,r);
BigConjMoV  := proc(k,n,r) option remember;
local i,F,S,A,B,C;

A := [];
B := [];
for i from 0 to 2*k+2 do
	S := ConjSumV(i,n,F);
	A := [op(A),
	coeff(coeff(S,n,degree(S,n)-k+1),F[n],1)];
	B := [op(B), 
	coeff(coeff(S,n,degree(S,n)-k+1),F[n-1],1)];
od:

C := [seq(factor(A[i]+B[i]*2/(1+sqrt(5)))
*5^((i-1)/2),i=1..2*k+3)];
GuessPol(C,0,r)*n^(r-k+1)/5^(r/2);

end:


# Input: non-negative integer r,
# symbolic n
# Output: The conjectured formula
# of Moment about the mean of V 
# in term of n
# Try: ConjMoMeanV(2,n);
ConjMoMeanV := proc(r,n) option remember;
local i,F,S,A;

A := add( (-1)^i*binomial(r,i)
*ConjSumV(r-i,n,F)*ConjSumV(1,n,F)^i/F[n]^(i+1)
, i=0..r);

S := {F[n-1]=1,F[n]=(1+sqrt(5))/2};

A := simplify(expand(subs(S,A)));

add( factor(coeff(A,n,i))*n^i,i=0..r);
end:


######################################################
# Input: symbolic r and n
# Output: Pre-compute for 
# the conjecture of E[V^r] 
# in term of moment r and 
# board of length n
# Try: ForMoV(r,n);
ForMoV := proc(r,n) option remember;
n^r/5^(r/2)*(1+r*(2*r+sqrt(5)-3)/sqrt(5)/n
+1/15*r*(r-1)*(6*r^2-32*r+6*r*5^(1/2)+37-9*5^(1/2))/n^2
+2/75*5^(1/2)*r*(r-1)*(r-2)*(2*r^3+3*r^2*5^(1/2)-23*r^2
-16*r*5^(1/2)+77*r+16*5^(1/2)-73)/n^3
+1/225*r*(r-1)*(r-2)*(r-3)*(6*r^4-120*r^3+12*r^3*5^(1/2)
-138*r^2*5^(1/2)+800*r^2-2100*r+432*r*5^(1/2)
-393*5^(1/2)+1849)/n^4);
end: 






######################################
# Section 3: 
# Proof of conjecture of
# moment about the mean
# 
# Functions:
# TaylorV(r,n), MainRelation(r,n),
# MyId1(), MyId2(), MyId3()
#
######################################

# Comment: Checking the Maclaurin
# series of CenGenV2 with the 
# moment about the mean

# Input: the r-th moment about
# the mean and the length of 
# the board n
# Output: coeff of t^r of 
# the generating function
# CenGenV2(e^t)

# Try: seq(TaylorV(r,3),r=0..8);
# seq(MoMeanV(r,3)/r!,r=0..8);
TaylorV :=proc(r,n) option remember;
local s,t,A;

A:=taylor(CenGenV2(n,exp(t)),t=0,r+1);
coeff(simplify(A),t,r);
end:


# Test recurrence with 
# E[(V-mu)^r], numeric r

# Input: non-negative integer r,
# symbolic n
# Output: The difference values
# on both sides of the relation

# Try: 
# seq(MainRelation(i,n),i=1..11);
MainRelation := proc(r,n) option remember;
local i,p,a,b,L,R;

p := (1+sqrt(5))/2;
a := 1-1/sqrt(5);
b := -2/sqrt(5);

L := ConjMoMeanV(r,n)/r!
-1/p*ConjMoMeanV(r,n-1)/r!-1/p^2*ConjMoMeanV(r,n-2)/r!;

R:= add(a^i/i!/p*ConjMoMeanV(r-i,n-1)/(r-i)!,i=1..r)
+add(b^i/i!/p^2*ConjMoMeanV(r-i,n-2)/(r-i)!,i=1..r);

factor(L-R);
end:




# Comment: Test of coeff of n^r

# Input: NULL
# Output: Verify the first 
# identity of an even case.
# Try: MyId1();
MyId1 := proc() option remember;
local a,b,p;

a:= 1-1/sqrt(5);
b:= -2/sqrt(5);
p := (1+sqrt(5))/2;

simplify(a/p+b/p^2);
end:


# Input: Null
# Output: Verify the second 
# identity of an even case.
# Try: MyId2();
MyId2 := proc() option remember;
local a,b,p,L,R;

a:= 1-1/sqrt(5);
b:= -2/sqrt(5);
p := (1+sqrt(5))/2;

L:= (1/p+2/p^2)*2/5/sqrt(5);
R := a^2/(2*p)+b^2/(2*p^2);

simplify(R/L);
end:

# Input: NULL
# Output: Verify the identity 
# of an odd case.
# Try: MyId3();
MyId3 := proc() option remember;
local a,b,x,y,p,r,L,R;

x := 2/15;
y := 2/5/sqrt(5);
a := 1-1/sqrt(5);
b := -2/sqrt(5);
p := (1+sqrt(5))/2;


L:=(1/p+2/p^2)*x*y*r;

R :=-(a/p+2*b/p^2)*y
+(a^2/(2*p)+b^2/(2*p^2))*x*(r-1)
+(a^3/6/p+b^3/6/p^2);

factor(simplify(R/L));
end:




#######################################
# Section 3.1: 
# Proof of conjecture of
# the straight moment! 
#
# Functions:
# GenPrV(n,v), SRelation(r,n),
# MyId4(r)
########################################

# Input: the length of the board n
# and variable v
# Output: Probability generating 
# function of V on an 2-by-n board

# Try: seq(GenPrV(n,v),n=0..5);
# seq(factor(GenV2(n,v)/SumV(0,n)),n=0..5);
GenPrV := proc(n,v) option remember;
local A;

if n=0 then
	return(1);
elif n=1 then
	return(v);
fi:	

A:= v*SumV(0,n-1)/SumV(0,n)*GenPrV(n-1,v)
+SumV(0,n-2)/SumV(0,n)*GenPrV(n-2,v);

factor(A);
end:


# Test recurrence with 
# E[V^r], numeric r

# Input: non-negative integer r,
# symbolic n
# Output: The difference values
# on both sides of the relation
# Try: 
# seq(SRelation(i,n),i=1..9);
SRelation := proc(r,n) option remember;
local i,p,L,R;

p := (1+sqrt(5))/2;

L := ConjMoV(r,n)/r!
-1/p*ConjMoV(r,n-1)/r!-1/p^2*ConjMoV(r,n-2)/r!;

R:= add(1/i!/p*ConjMoV(r-i,n-1)/(r-i)!,i=1..r);

factor(L-R);
end:


# Input: symbolic r 
# Output: Verify the identity 
# of coeff of n^(r-1) in the 
# relation of straight moment
# Try: MyId4(r);
MyId4 := proc(r) option remember;
local p,L,R;

p := (1+sqrt(5))/2;

L := (2*r+sqrt(5)-3)/5*(1/p+2/p^2)-(1/p+4/p^2)/2/sqrt(5);
R := 1/p*(-1+(2*r+sqrt(5)-5)/sqrt(5)+sqrt(5)/2);

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




########################################
# Section 4: General board size,
# Counting Moment of horizontal
# or vertical Dominos
# 
# Functions:
# ListHV(m,n), CountB(T),
# 
# SumPowHV(a,b,m,n), MoHV(a,b,m,n), 
# GenHV(m,n,v,t),
#
########################################

# Input: Board of width m and length n
# Output: List L = [a_i],
# a_i is the number of possilbe boards 
# that contained exactly i vertical dominos.
#
# Try: ListHV(2,1);
ListHV := proc(m,n) option remember;
local T;

T := [[0$n]$m];
CountB(T,m,n); 
end:


# Input: Board T of size m-by-n
# Output: the distribution of 
# total number of H of each tiling 
# from the rest of the board
# Try: CountB([[1,0],[1,0],[0,0]],3,2);
CountB :=proc(T,m,n) option remember;
local i,j,k,a,L,NT,R,S;

# Find the next spot
i := 1:
j := 1:

while T[i][j] = 1 do
	if i = m then
		if j = n then
			return([[1]]);
		else	
			i := 1;
			j := j+1;
		fi:
	else
		i := i+1;
	fi:
od:


# put H or V and return
R := [[]];
S := [[]];

if i < m and T[i+1,j]=0 then
	NT := T;
	NT[i][j] := 1:
	NT[i+1][j] := 1:
	L := CountB(NT,m,n);
	
	if L <> [[]] then	
		R := [seq([0,op(L[k])],k=1..nops(L)),
			[0$(1+nops(L[1]))]];
	fi:
fi:
	
if j < n and T[i,j+1]=0 then
	NT := T;
	NT[i][j] := 1:
	NT[i][j+1] := 1:
	L := CountB(NT,m,n);
	
	if L <> [[]] then	
		S := [[0$(1+nops(L[1]))],
			seq([op(L[k]),0],k=1..nops(L))];
	fi:	
	
fi:

if R <> [[]] and S <> [[]] then
	return(R+S);
elif R=[[]] then	
	return(S);
else
	return(R);
fi:

end:


# Input: numeric a,b,m and n
# Output: sum of the V^a*H^b 
# of tiling by Vertical-Horizontal dominos 
# from the board of size m-by-n
# Try: A:=[seq(SumPowHV(0,0,2,n),n=1..7)];
# 
SumPowHV := proc(a,b,m,n) option remember;
local i,j,LHV;

LHV := ListHV(m,n);

add(add(i^a*j^b*LHV[j+1,i+1]
,j=0..nops(LHV)-1),i=0..nops(LHV[1])-1);
end:


# Input: numeric a,b,m and n
# Output: moment of the V^a*H^b (E[V^a*H^b]) 
# of tiling by Vertical-Horizontal dominos 
# from the board of size m-by-n

# Try: [seq(MoHV(1,1,2,n),n=1..7)];
# [seq(evalf(MoHV(1,0,2,n)),n=1..17)];
MoHV := proc(a,b,m,n) option remember;
local i,j,tol,LHV;

LHV := ListHV(m,n);

tol := add(add(LHV[j+1,i+1]
,j=0..nops(LHV)-1),i=0..nops(LHV[1])-1);

SumPowHV(a,b,m,n)/tol;
end:

# Input: numeric m,n and 
# symbolic v and t
# Output: Bivariate Generating 
# function (v,t) from empirical 
# result of the board of size m-by-n
# Try: 
# seq(GenHV(4,n,v,1),n=1..10);
GenHV := proc(m,n,v,t) option remember;
local i,j,tol,A,M;

if n=0 or m=0 then
	return(1);
fi:
	
A := ListHV(m,n);

#tol := add(add(A[j+1,i+1]
#,j=0..nops(A)-1),i=0..nops(A[1])-1);
tol := 1;

M := add(add(A[j+1,i+1]*v^i*t^j
,j=0..nops(A)-1),i=0..nops(A[1])-1)/tol;

factor(M);
end:


##################################
# section 4.1: 
# Recurrences of generating function
# on an m-by-n board
#
# Functions:
# GenVm(m,n,q), 
# GuessGenRec(m,p,t,N,v)
# 
###################################


# Input: numeric m,n and symbolic v
# Output: Generating function of v
# (vertical tile) from
# the board of size m-by-n computed
# from empirical data
# Try: 
# seq(GenVm(4,n,v),n=1..10);
GenVm := proc(m,n,v) option remember;
GenHV(m,n,v,1);
end:

# Input: numeric m,p,t and symbolic n,v
# Output: t-th order recurrences with
# polynomial degree p of on an m-by-n board.
# Try: GuessGenRec(2,1,2,N,v);
# GuessGenRec(3,2,2,N,v);
# GuessGenRec(4,2,4,N,v);
# GuessGenRec(5,4,4,N,v);
# GuessGenRec(6,6,8,N,v);
# GuessGenRec(7,12,8,N,v);
GuessGenRec := proc(m,p,t,N,v) option remember;
local i,j,k,c,n,a,P,T,eq,var;

if m mod 2 = 0 then
	c:=1; 
else
	c:=2:
fi:

n := 2;
while degree(GenVm(m,n,v),v) < (p+1)*t+6 do
	n := n+c;
od:	

n := n+c;
P := GenVm(m,n,v);

T := add(add(a[i,j]*v^i,i=0..p)
*GenVm(m,n-c*j,v),j=1..min(t,n/c));
T := expand(T);

eq := {seq( coeff(P,v,k)=coeff(T,v,k)
,k=0..degree(P,v))};
var := {seq(seq(a[i,j],i=0..p),j=1..t)};

var := solve(eq,var);

if var = op({}) then
	return("No Solution");
fi:
	
H(N)=subs(var,add(add(a[i,j]*v^i,i=0..p)
*H(N-c*j),j=1..t));
end: