######################################################################
## RecTILINGS: Save this file as  RecTILINGS to use it,              #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read RecTILINGS: <Enter>                          #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  zeilberg at math dot rutgers dot edu.                            # 
######################################################################

with(combinat):

Digits:=100:

print(`First Written: Dec. 2005: tested for Maple 10  `):
print(`Version of Dec. , 2005:  `):
print():
print(`This is RecTILINGS, A Maple package`):
print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's  article: `):
print(`"Automatic CounTilings" `):

print():
print(`The most current version is available on WWW at:`):
print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/RecTILINGS .`):
print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezra1();`):
print():
 

ezra1:=proc()
if args=NULL then

print(` The supporting procedures are: Asy1`):
print(` DS, EqG, EqGt, Growth1, KmnF, KZmn, Red, Sefer, SeferW, Zmn , Zinn`):
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
else
 ezra(args):
fi:
end:

ezra:=proc()
if args=NULL then

print(` RecTILINGS: A Maple package for counting tilings `):
print(`of rectangles by any set of rectangular tiles`):
print(`of arbitrary sizes`):
print(`The MAIN procedures are`):
print(` AaveT, Astat,AstatV, AstatVs, GF1, GF1t, Sidra, SidraW ,  Sipur, SipurW `):

elif nargs=1 and args[1]=AaveT then
print(`AaveT(k,ListOfTiles): Given a postive integer k and a symbol n`):
print(`The list of asymptotic averages (divided by n) of the number of tiles`):
print(`of tiling a k by n rectangular board by the  tiles in`):
print(`the list of rectangular tiles Tiles ([a,b] means a b by a rectangle)`):
print(`followed by the average total number of total number of tiles (divided by n)`):
print(`For example, try: AaveT(2,[[1,2],[2,1]]) `):


elif nargs=1 and args[1]=Aitken then
print(`Aitken(L): The conjectured limit of the sequence L`):
print(`from its three last elements, using Aitken's fromula`):
print(`For example, try: Aitken([3.1,3.01,3.001]);`):

elif nargs=1 and args[1]=Astat then
print(`Astat(k,ListOfTiles,n): Given a postive integer k and a symbol n`):
print(` and given a List of rectangular Tiles`):
print(`ListOfTiles`):
print(`returns: (1) The lists of for the`):
print(`asymptotics (in n) for pairs [average,variance] for the`):
print(`r.v. number of tiles of type ListOfTiles[i], in a random tiling by the rectangular tiles`):
print(`of Tiles (in order), followed by the pair [av,variance] for the`):
print(`Total number of tiles, followed by the asympto. correlation`):
print(`matrix for example try:`):
print(`Astat(4,[[1,2],[2,1]],n):`):

elif nargs=1 and args[1]=AstatV then
print(`AstatV(k,ListOfTiles,n): Verbose version of Astat (q.v.) `):
print(`For example, try:`):
print(`AstatV(4,[[1,2],[2,1]],n):`):

elif nargs=1 and args[1]=AstatVs then
print(`AstatVs(K,ListOfTiles,n): AstatV(k,ListOfTiles,n) for k from 1 to K `):
print(`For example, try:`):
print(`AstatVs(4,[[1,2],[2,1]],n):`):

elif nargs=1 and args[1]=Asy1 then
print(`Asy1(f,t,n): the Asymptotics of the coeff. of the rational function f`):
print(`try for example: Asy1(1/(1-t-t^2),t,n);`):

elif nargs=1 and args[1]=DS then
print(`D(v): the distingushed Square followed by the number of`):
print(`available squares below`):
print(`For example, try: DS([0,0,0,0]);`):

elif nargs=1 and args[1]=EqG then
print(`EqG(v,A,t,H,Tiles): the equation for A[v] in terms of H[Tiles]`):
print(`where Tiles is the set of rectangular tiles`):
print(`For example, try: EqG([0,0],A,t,H,{[1,2],[2,1]});`):

elif nargs=1 and args[1]=EqtG then
print(`EqtG(v,A,t,Tiles): the equation for A[v]`):
print(`where Tiles is the set of rectangular tiles`):
print(`For example, try: EqtG([0,0],A,t,{[1,2],[2,1]});`):

elif nargs=1 and args[1]=GF1 then
print(`GF1(k,t,H,Tiles): Given a positive integer k, a variable t`):
print(`A letter H, and a set of pos. integer-pairs Tiles`):
print(`outputs the rational function in t`):
print(`whose general coeff. (of t^n) is the weight-enumeator of the set`):
print(`of tilings of an n by k rectangle by `):
print(`tiles drawn from the set of rectangular tiles, Tiles, given as`):
print(`a set of pairs [a,b] (representing an a by b rectangle)`):
print(` a is the vertical dimension and b is the horizontal dimension `):
print(`The weight of a tiling is the product of H[a,b]'s `):
print(`For example, try: GF1(2,t,H,{[1,2],[2,1]});`):

elif nargs=1 and args[1]=GF1t then
print(`GF1t(k,t,Tiles): Given a positive integer k, a variable t`):
print(`and a set of pos. integer-pairs Tiles`):
print(`outputs the rational function in t`):
print(`whose general coeff. (of t^n) is the number `):
print(`of tilings of an n by k rectangle by `):
print(`tiles drawn from the set of rectangular tiles, Tiles, given as`):
print(`a set of pairs [a,b] (representing an a by b rectangle)`):
print(` a is the vertical dimension and b is the horizontal dimension `):
print(`For example, try: GF1t(2,t,{[1,2],[2,1]});`):

elif nargs=1 and args[1]=Growth1 then
print(`Growth1(f,t): the growth of the coeff. of the rational function f`):
print(`try for example: Growth1(1/(1-t-t^2),t);`):


elif nargs=1 and args[1]=KmnF then
print(`KmnF(m,n,z): The Kasteleyn polynomial for tiling`):
print(`an m by n board by domino tiles computed`):
print(`according to Kasteleyn's fromula as it appears in`):
print(`his paper "The statistics of Dimers on a Lattice I...."`):
print(`Physica 27 (1961) Eq. (13) p. 1215`):
print(`Using Floating point`):
print(`For example, try: KmnF(3,4,z);`):

elif nargs=1 and args[1]=KZmn then
print(`KZmn(m,n,z): The Kasteleyn polynomial for an m by n board `):
print(`For example, try: KZmn(3,2,z);`):

elif nargs=1 and args[1]=Red then
print(`Red(v): the reduction of a vector`):
print(`For example, try: Red([-1,-1,-2,-2] `);


elif nargs=1 and args[1]=Sefer then
print(`Sefer(t,K,N,L1,L2) : All the Sipurs for Tiles`):
print(`ALL the non-empty subsets of tiles of dimensions<=L1`):
print(`of cardinality  L2 `):
print(`For example, try: Sefer(t,2,10,2,2);`):

elif nargs=1 and args[1]=SeferW then
print(`SeferW(t,H,K,N,L1,L2) : All the SipurWs for Tiles`):
print(`ALL the non-empty subsets of tiles of dimensions<=L1`):
print(`of cardinality  L2 `):
print(`For example, try: SeferW(t,H,2,10,2,2);`):

elif nargs=1 and args[1]=Sidra then
print(`Sidra(k,Tiles,N): the first N+1 terms in the`):
print(`enumerating sequence for the number of tilings `):
print(`of an n by k rectangle using tiles from the set`):
print(`or Rectangular tiles Tiles {[a,b]} where`):
print(`[a,b] stands for an a by b rectangular tile`):
print(` For example, try: Sidra(2,{[1,2],[2,1]},20); `):

elif nargs=1 and args[1]=SidraW then
print(`SidraW(k,H,Tiles,N): the first N+1 terms in the`):
print(`weight-enumerating sequence for the number of tilings `):
print(`of an n by k rectangle using tiles from the set`):
print(`or Rectangular tiles Tiles {[a,b]} where`):
print(`[a,b] stands for an a by b rectangular tile`):
print(` and the corresponding weight is H[a,b]`):
print(` For example, try: SidraW(2,H,{[1,2],[2,1]},20); `):

elif nargs=1 and args[1]=Sipur then
print(`Sipur(Tiles,t,K,N): Tells the story for the`);
print(`enumeration of rectangular boards of width up to K`):
print(`using the rectangular tiles Tiles. t is the variable for`):
print(`the length of the board `):
print(`N is a positive integer for the lengths of the`):
print(`desired beginnings for the enumetaring sequence`):
print(`Sipur({[1,2],[2,1]},t,4,10);`):

elif nargs=1 and args[1]=SipurW then
print(`SipurW(Tiles,t,H,K,N): Tells the story for the`);
print(`enumeration of rectangular boards of width up to K`):
print(`using the rectangular tiles Tiles. t is the variable for`):
print(`the length of the board and H is a symbol for the`):
print(`indexed variables H[a,b] accounting for the tile-types`):
print(`N is a positive integer for the lengths of the`):
print(`desired beginnings for the enumetaring sequence`):
print(`SipurW({[1,2],[2,1]},t,H,4,10);`):

elif nargs=1 and args[1]=Zinn then
print(`Zinn(resh,n): Given a list resh, conjectures mu and theta`):
print(`such that resh[i]=mu^n*n^theta`):
print(`For example, try: Zinn(Sidra(6,{[1,2],[2,1]},40),39);`):

elif nargs=1 and args[1]=Zmn then
print(`Zmn(m,n,H,Tiles): the weight-enumerator of all tilings of`):
print(`an m by n rectangle with rectangular tiles from the`):
print(`set Tiles, for example, try:`):
print(`Zmn(2,2,H,{[1,2],[2,1]});`):

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

fi:


end:

###Begin General Rectangular Tiles

#Red(v): the reduction of a vector
Red:=proc(v) local i,m: m:=max(op(v)):
[seq(v[i]-m,i=1..nops(v))]:
end:

#DS(v): the distinguished square, and the number of squares
#below it
DS:=proc(v) local m,i,k:
m:=min(op(v)):
for i from 1 while v[i]<>m do od:
for k from 1 to nops(v)-i while v[i+k]=m do od:
i,k:
end:

#GF1t(k,t,Tiles): the generating function for tilings of an n by k
#rectangle by  rectangular tiles from the set Tiles, given as
#pairs [a,b] (an a by b rectangle)
#For example, try: GF1t(2,t,{2},{2});
GF1t:=proc(k,t,Tiles)
local eq,OldVar,A,NewVar,gu,var,i:
option remember:
OldVar:={}:
NewVar:={[0$k]}:
eq:={}:

while NewVar<>{} do
 gu:=EqtG(NewVar[1],A,t,Tiles):
 OldVar:=OldVar union {NewVar[1]}:
 NewVar:=NewVar minus {NewVar[1]}:
 eq:=eq union {gu[1]}:
 NewVar:=NewVar union (gu[2] minus OldVar):
od:

var:={seq(A[op(OldVar[i])],i=1..nops(OldVar))}:

var:=solve(eq,var):
normal(subs(t=t^(1/k),normal(subs(var,A[0$k])))):
end:


#EqtG(v,A,t,Tiles): the equation for A[v] in terms of t
#where Tiles is the set of rectangular tiles
#For example, try: EqtG([0,0],A,t,{[1,2],[2,1]});
EqtG:=proc(v,A,t,Tiles) local gu,i,k,i1,eq,v1,i2,var,a,b,tile:
gu:=DS(v):i:=gu[1]: k:=gu[2]:
eq:=A[op(v)]:
var:={}:

for i1 from 1 to nops(Tiles) do
tile:=Tiles[i1]:
a:=tile[2]: b:=tile[1]:

if b<=k then
  v1:=Red([op(1..i-1,v),seq(v[i2]+a,i2=i..i+b-1),
 op(i+b..nops(v),v)]):
 var:=var union {v1}:
 eq:=eq-t^(a*b)*A[op(v1)]:
 var:=var union {v1}:
fi:

od:


if v=[0$nops(v)] then eq:=eq-1: fi:
eq,var:
end:


#EqG(v,A,t,H,Tiles): the equation for A[v] in terms of H[Tiles]
#where Tiles is the set of rectangular tiles
#For example, try: EqG([0,0],A,t,H,{[1,2],[2,1]});
EqG:=proc(v,A,t,H,Tiles) local gu,i,k,i1,eq,v1,i2,var,a,b,tile:
gu:=DS(v):i:=gu[1]: k:=gu[2]:
eq:=A[op(v)]:
var:={}:
for i1 from 1 to nops(Tiles) do
tile:=Tiles[i1]:
a:=tile[2]: b:=tile[1]:

if b<=k then
  v1:=Red([op(1..i-1,v),seq(v[i2]+a,i2=i..i+b-1),
 op(i+b..nops(v),v)]):
 var:=var union {v1}:
 eq:=eq-t^(a*b)*H[op(tile)]*A[op(v1)]:
 var:=var union {v1}:
fi:

od:


if v=[0$nops(v)] then eq:=eq-1: fi:
eq,var:
end:







#GF1(k,t,H,Tiles): the generating function for tilings of an n by k
#rectangle by  rectangular tiles from the set Tiles, given as
#pairs [a,b] (an a by b rectangle)
#The weight of a tiling is the product of H[op(tile)]'s and t^n
#For example, try: GF1(2,t,H,{[1,2],[2,1]});
GF1:=proc(k,t,H,Tiles)
local eq,OldVar,A,NewVar,gu,var,i:
option remember:
OldVar:={}:
NewVar:={[0$k]}:
eq:={}:

while NewVar<>{} do
 gu:=EqG(NewVar[1],A,t,H,Tiles):
 OldVar:=OldVar union {NewVar[1]}:
 NewVar:=NewVar minus {NewVar[1]}:
 eq:=eq union {gu[1]}:
 NewVar:=NewVar union (gu[2] minus OldVar):
od:

var:={seq(A[op(OldVar[i])],i=1..nops(OldVar))}:

var:=solve(eq,var):
normal(subs(t=t^(1/k),normal(subs(var,A[0$k])))):
end:

#Sidra(k,Tiles,N): the first N+1 terms in the
#enumerating sequence for the number of tilings using
#of an n by k rectangle using tiles from the set
#or Rectangular tiles Tiles {[a,b]} where
#[a,b] stands for an a by b rectangular tile
#For example, try: Sidra(2,{[1,2],[2,1]},20);
Sidra:=proc(k,Tiles,N) local gu,t,i:
gu:=GF1t(k,t,Tiles):
[seq(coeff(taylor(gu,t=0,N+1),t,i),i=0..N)]:
end:


#Zmn(m,n,H,Tiles): the weight-enumerator of all tilings of
#an m by n rectangle with rectangular tiles from the
#set Tiles, for example, try:
#Zmn(2,2,H,{[1,2],[2,1]});
Zmn:=proc(m,n,H,Tiles) local lu,t:
lu:=GF1(m,t,H,Tiles):
lu:=taylor(lu,t=0,n+1):
lu:=sort(expand(coeff(lu,t,n))):
end:


#KZmn(m,n,z): The Kastelyan polynomial
#For example, try: KZmn(2,4,z);
KZmn:=proc(m,n,z) local H,gu:
gu:=subs({H[1,2]=1,H[2,1]=z},Zmn(m,n,H,{[1,2],[2,1]})):
gu:=subs(z=sqrt(z),gu);
end:

#Growth1(f,t): the growth of the coeff. of the rational function f
#try for example: Growth1(1/(1-t-t^2),t);
Growth1:=proc(f,t) local i,aluf,ku,shi,P:
P:=denom(normal(f)):
if degree(P,t)=0 then
 RETURN(FAIL):
fi:
ku:={fsolve(P,t)}:
aluf:=ku[1]:
shi:=abs(aluf):

for i from 2 to nops(ku) do
 if abs(ku[i])<shi then
    aluf:=ku[i]:
    shi:=abs(ku[i]):
 fi:
od:
1/shi:
end:

#FindMult1(P): Finds the highest multiplicity
#of the polynomial P in t
FindMult1:=proc(P,t)
local k:
for k from 1 to degree(P,t) while degree(gcd(P,diff(P,t$k)),t)>0 do od:
k:
end:


#Asy1(f,t,n): the Asymptotics of the coeff. of the rational function f
#try for example: Asy1(1/(1-t-t^2),t,n);
Asy1:=proc(f,t,n) local i,aluf,ku,shi,P,Q,f1:
f1:=normal(f):
Q:=numer(f1):
P:=denom(f1):
ku:={fsolve(P,t)}:
aluf:=ku[1]:
shi:=abs(aluf):

for i from 2 to nops(ku) do
 if abs(ku[i])<shi then
    aluf:=ku[i]:
    shi:=abs(ku[i]):
 fi:
od:
-(1/shi)^n*subs(t=shi,Q)/subs(t=shi,diff(P,t))/shi:
end:



#SidraW(k,H,Tiles,N): the first N+1 terms in the
#weight-enumerating sequence for the number of tilings using
#of an n by k rectangle using tiles from the set
#or Rectangular tiles Tiles {[a,b]} where
#[a,b] stands for an a by b rectangular tile
#and the weight corresponding to [a,b] is H[a,b]
#For example, try: SidraW(2,H,{[1,2],[2,1]},20);
SidraW:=proc(k,H,Tiles,N) local gu,t,i:
gu:=GF1(k,t,H,Tiles):
[seq(expand(coeff(taylor(gu,t=0,N+1),t,i)),i=0..N)]:
end:



#KmnF(m,n,z): The Kasteleyn polynomial for tiling
#an m by n board by domino tiles computed
#according to Kasteleyn's fromula as it appears in
#his paper "The statistics of Dimers on a Lattice I...."
#Physica 27 (1961) Eq. (13) p. 1215
#Using Floating point
#For example, try: KmnF(3,4,z);
KmnF:=proc(m,n,z) local gu,K,L,i:
Digits:=100:
gu:=1:


if n mod 2=1 then
 RETURN(FAIL):
fi:

gu:=2^(n*trunc(m/2)):


 for K from 1 to m/2 do
 for L from 1 to n/2 do
   gu:=gu*(z*cos(K*Pi/(m+1))^2+cos(L*Pi/(n+1))^2):
 od:
 od:

gu:=add(round(evalf(coeff(gu,z,i)))*z^i,i=0..degree(gu,z)):

RETURN(gu):

end:

#Pashet1Old(P,w,K): Given a polynomial in the primitive
#Kth root of unity, simplifies it
Pashet1Old:=proc(P,w,K) local gu1,i,gu2,gu3:
gu1:=0:
for i from ldegree(P,w) to degree(P,w) do
gu1:=gu1+coeff(P,w,i)*w^(i mod 2*K):
od:

gu2:=0:
for i from 0 to K-1 do
 gu2:=gu2+coeff(gu1,w,i)*w^i:
od:

for i from K to 2*K-1 do
 gu2:=gu2-coeff(gu1,w,i)*w^(i-K):
od:

if K mod 2 =0 then
 RETURN(gu2):
fi:

if K mod 2=1 then
 gu3:=0:
 for i from 0 to K-2 do
   gu3:=gu3+coeff(gu2,w,i)*w^i:
 od:

gu3:=gu3+coeff(gu2,w,K-1)*add(w^i*(-1)^(i+1),i=0..K-2):
gu3:=expand(gu3):
RETURN(gu3):
fi:
end:

#Pashet1(P,w,K): Given a polynomial in the primitive
#Kth root of unity, simplifies it
Pashet1:=proc(P,w,K) :
rem(numer(normal(P)),cyclotomic(2*K,w),w):
end:


cos1:=proc(z): (z+1/z)/2:end:
#Kmn(m,n,z): The Kasteleyn polynomial for tiling
#an m by n board by domino tiles computed
#according to Kasteleyn's fromula as it appears in
#his paper "The statistics of Dimers on a Lattice I...."
#Physica 27 (1961) Eq. (13) p. 1215
#Using Floating point
#For example, try: Kmn(3,4,z);
Kmn:=proc(m,n,z,wm,wn) 
#local gu,K,L,wm,wn:
local gu,K,L:
Digits:=100:
gu:=1:


if n mod 2=1 then
 RETURN(FAIL):
fi:

gu:=2^(n*trunc(m/2)):


 for K from 1 to m/2 do
 for L from 1 to n/2 do
   gu:=gu*(z*cos1(wm^K)^2+cos1(wn^L)^2):
 od:
 od:
gu:=expand(gu):

gu:=Pashet1(gu,wm,m+1):
gu:=Pashet1(gu,wn,n+1):

#gu:=subs({wm=exp(Pi*I/(m+1)),wn=exp(Pi*I/(n+1))},gu):
#gu:=expand(gu):
#gu:=evalc(gu):
#gu:=evalf(gu):
#gu:=add(round(evalf(coeff(gu,z,i)))*z^i,i=0..degree(gu,z)):
RETURN(gu):

end:







#Zinn(resh,n): Given a list resh, conjectures mu and theta
#such that resh[i]=mu^n*n^theta
Zinn:=proc(resh,n)
local s1,s2:
s1:=sn(resh,n):
s2:=sn(resh,n-1):
evalf(2*(s1+s2)/(s1-s2)^2),
evalf(sqrt(op(n+1,resh)/op(n-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))):
end:
 
 
sn:=proc(resh,n):
-1/log(op(n+1,resh)*op(n-1,resh)/op(n,resh)^2):
#evalf(%):
end:
 
#Aitken(L): The conjectured limit of the sequence L
#from its three last elements, using Aitken's fromula
Aitken:=proc(L)
local n,a,b,c:
n:=nops(L):

if n<3 then
 RETURN(FAIL):
fi:

a:=L[n-2]: b:=L[n-1]: c:=L[n]:

c-(c-b)^2/(c-2*b+a):
end:


#SipurW(Tiles,t,H,K,N): Tells the story for the
#enumeration of rectangular boards of width up to K
#using the rectangular tiles Tiles. t is the variable for
#the length of the board and H is a symbol for the
#indexed variables H[a,b] accounting for the tile-types
#N is a positive integer for the lengths of the
#desired beginnings for the enumetaring sequence 
#For example, try:
#SipurW({[1,2],[2,1]},t,H,4,10);

SipurW:=proc(Tiles,t,H,K,N) local k,lu,gu,vu:
print(`This is the story of the number of tilings of  rectangular`):
print(`boards of width <=`,K, `using rectangular tiles of dimensions`):
print(Tiles):
print(`where [a,b] is the rectngle of width a and length b`):
print(`In the generating functions for weight-enumerators given below`):
print(` H[a,b] corresponds to an a by b  rect. tile (a vertical dim. b horiz. dim.) is`):
vu:=[]:

for k from 1 to K do
 print(`The number of tilings using the tiles`):
 print(Tiles):
 print(`of a `, k, ` by n rectangular board is the coeff. of t^n in the rational function`):
 gu:=GF1t(k,t,Tiles):
 lprint(gu):
 print(`The first `, N+1, `terms are `):
 print(Sidra(k,Tiles,N)):
 lu:=Growth1(gu,t):
 print(`The asymptotic rate of growth is`,lu):
 print(`and adjusted it is`, lu^(1/k)):
 vu:=[op(vu),lu^(1/k)]:
 print(``):
 print(`The weight-enumerator of the set of the above-mentioed set of tilings`):
 print(`is the coeff. of t^n in the rational function`):


 gu:=GF1(k,t,H,Tiles):
 lprint(gu):
 print(``):
od:
 print(`To summarize, the sequence of adjusted growth-rates `):
 print(`For tilings with the set of tiles`, Tiles , ` is `):
 print(vu):
end:

#SeferW(t,H,K,N,L1,L2) : All the SipurWs for Tiles
#ALL the non-empty subsets of tiles of dimensions<=L1
#of cardinality L2 
#For example, try: SeferW(t,H,4,10,2,2);
SeferW :=proc(t,H,K,N,L1,L2)  local gu,a,b,i,mu:

mu:={seq(seq([a,b],a=1..L1),b=1..L1)};

print(`This is a book about enumerating all tilings with a set`):
print(` of rectangular tiles {[a,b]} with a and b between 1 and `, L1):
print(`with `, L2, `tiles  `):


print(`The set of all tiles is`, mu):

gu:=choose(mu,L2):
 for i from 1 to nops(gu) do
  SipurW(gu[i],t,H,K,N):
 od:

end:


Sipur:=proc(Tiles,t,K,N) local k,lu,gu,vu:
print(`This is the story of the number of tilings of  rectangular`):
print(`boards of width <=`,K, `using rectangular tiles of dimensions`):
print(Tiles):
print(`where [a,b] is the rectngle of width a and length b`):
vu:=[]:

for k from 1 to K do
 print(`The number of tilings using the tiles`):
 print(Tiles):
 print(`of a `, k, ` by n rectangular board is the coeff. of t^n in the rational function`):
 gu:=GF1t(k,t,Tiles):
 lprint(gu):
 print(`The first `, N+1, `terms are `):
 print(Sidra(k,Tiles,N)):
 lu:=Growth1(gu,t):
 print(`The asymptotic rate of growth is`,lu):
 print(`and adjusted it is`, lu^(1/k)):
 vu:=[op(vu),lu^(1/k)]:
 print(``):
 print(``):
od:
 print(`To summarize, the sequence of adjusted growth-rates `):
 print(`For tilings with the set of tiles`, Tiles , ` is `):
 print(vu):
vu:
end:

#Sefer(t,H,K,N,L1,L2) : All the SipurWs for Tiles
#ALL the non-empty subsets of tiles of dimensions<=L1
#of cardinality L2 
#For example, try: SeferW(t,H,4,10,2,2);
Sefer :=proc(t,H,K,N,L1,L2)  local gu,a,b,i,mu:

mu:={seq(seq([a,b],a=1..L1),b=1..L1)};

print(`This is a book about enumerating all tilings with a set`):
print(` of rectangular tiles {[a,b]} with a and b between 1 and `, L1):
print(`with `, L2, `tiles  `):


print(`The set of all tiles is`, mu):

gu:=choose(mu,L2):
 for i from 1 to nops(gu) do
  Sipur(gu[i],t,H,K,N):
 od:

end:


###############Begin Stuff from StatGF


Digits:=100:
#Asy(f,t,n): Given a rational functionf of the variab;e t, finds the
#asymptotics (in floating point) if the coeff. of t^n
#For example, try: Asy(1/(1-t-^2),t,n);

Asy:=proc(f,t,n) local gu,gu1,ku,i,alufim,shi,muam,lu1,gu2:

if f=0 then
 RETURN(0):
fi:

gu:=convert(f,parfrac,t,complex):



if not type(gu,`+`) then
 RETURN(Asy1(gu,t,n)):
fi:

ku:={seq({solve(denom(op(i,gu)),t)}[1],i=1..nops(gu))}:



alufim:={ku[1]}:
shi:=abs(ku[1]):

for i from 2 to nops(ku) do
muam:=abs(ku[i]):

if muam<shi then
  alufim:={ku[i]}:
  shi:=muam:

elif muam=shi then
 alufim:=alufim union {ku[i]}:
fi:
od:

gu1:=0:

for i from 1 to nops(gu) do
lu1:=op(i,gu):

if member({solve(denom(lu1),t)}[1],alufim) then

gu1:=gu1+ lu1:
fi:

od:


if not type(gu1,`+`)  then
 RETURN(Asy1(gu1,t,n)):
fi:

gu2:=0:


for i from 1 to nops(gu1) do

 gu2:=gu2+Asy1(op(i,gu1),t,n):
od:

gu2:

end:


#Asy1(f,t,n): The asymptotics of a single term
Asy1:=proc(f1,t,n) local top1,bot1,a,m:


top1:=numer(f1):
bot1:=denom(f1):


a:={solve(bot1,t)}[1]:



if type(bot1,`^`) then
 m:=op(2,bot1):
top1:=top1/coeff(op(1,bot1),t,1):
else
 m:=1:
top1:=top1/coeff(bot1,t,1):
fi:


top1*(-1)^m/a^m*simplify((m+n-1)!/(m-1)!/n!)*(1/a)^n:
end:

#CheckAsy(f,t): checks the correctness of Asy(f,t,n)
#For example, try: CheckAsy(1/(1-t),t);
CheckAsy:=proc(f,t)
local gu,lu,n,i:
gu:=taylor(f,t=0,100):
lu:=Asy(f,t,n):

{seq(coeff(gu,t,i)/subs(n=i,lu),i=80..99)}:

end:


#AsyM(F,t,VarList,ExpList,n): Given F, a rational function of t whose
#coeffs. of t^n are polynomials in the variables of VarList
#corresponding to r.v.
#copmutes the asymptotic value of the Expectation of 
#the product of VarList[i]^ExpList[i]
#For example, to find the asymptotic second moment
#of the number of Heads in tossing n coins type:
#AsyM(1/(1-x*t-y*t^2)),t,[x,y],[2,0],n);
AsyM:=proc(F,t,VarList,ExpList,n)
local gu,i,j,lu,ku:
gu:=F:

for i from 1 to nops(VarList) do
 for j from 1 to ExpList[i] do
   gu:=VarList[i]*diff(gu,VarList[i]):
   gu:=normal(gu):
 od:
od:


gu:=subs({seq(VarList[i]=1,i=1..nops(VarList))},gu):


lu:=subs({seq(VarList[i]=1,i=1..nops(VarList))},F):


ku:=normal(Asy(gu,t,n)/Asy(lu,t,n)):

if degree(ku,n)=FAIL then
 RETURN(FAIL):
fi:

ku:
end:


#Trim2(f,n,L): trims a floating-point expression to L-digit appx.
Trim2:=proc(f,n,L)
local gu,i:

if f=0 then
RETURN(0):
fi:
if f=FAIL then
 RETURN(FAIL):
fi:

gu:=0:

for i from 0 to degree(f,n) do
if abs(coeff(f,n,i))>10^(-L-1) then
 gu:=gu+evalf(coeff(f,n,i),L)*n^i:
fi:
od:
gu:
end:

Trim1:=proc(f,n,L) Trim2(coeff(f,I,1),n,L)*I+ Trim2(coeff(f,I,0),n,L): end:
#Amoment(F,t,VarList,i,k,n): The asymptotic k^th moment about the
#mean  of the r.v. corresponding to the i^th entry in VarList
#where F is a rational function in t whose coeffs. of t^n are
#weight-enumerators according to statistics represented by VarList
#For example (assuming that you have CountTilings), to get
#the asymtotic (in n) variance of the number of horizontal tiles
#in a domino-tiling of a 4 by n board, type:
#Amoment(GF1(4,t,H,{[1,2],[2,1]}),t,[H[1,2],H[2,1]],1,2,n):
Amoment:=proc(F,t,VarList,i,k,n) local gu,j,lu,mu,i1:
gu:=[seq(AsyM(F,t,VarList,[0$(i-1),j,0$(nops(VarList)-i)],n),j=1..k)]:

mu:=gu[1]:

lu:=(-1)^k*mu^k:

for i1 from 1 to k do
 lu:=lu+(-1)^(k-i1)*binomial(k,i1)*mu^(k-i1)*gu[i1]:
od:

Trim1(expand(lu),n,40):


end:


#Akurtois(F,t,VarList,i,k,n): The asymptotic Kurtosis 
#of the r.v. corresponding to the i^th entry in VarList
#where F is a rational function in t whose coeffs. of t^n are
#weight-enumerators according to statistics represented by VarList
#For example (assuming that you have CountTilings), to get
#the asymtotic (in n) Kurtosis of the number of horizontal tiles
#in a domino-tiling of a 4 by n board, type:
#AKurtosis(GF1(4,t,H,{[1,2],[2,1]}),t,[H[1,2],H[2,1]],1,n):
Akurtosis:=proc(F,t,VarList,i,n) local m2,m4:
m2:=Amoment(F,t,VarList,i,2,n):
m4:=Amoment(F,t,VarList,i,4,n):

Trim1(normal(coeff(m4,n,degree(m4,n))*n^degree(m4,n)/
(coeff(m2,n,degree(m2,n))*n^degree(m2,n))^2),n,40):

end:





#Aaverage(F,t,VarList,i,n): The asymptotic expectation
#of the r.v. corresponding to the i^th entry in VarList
#where F is a rational function in t whose coeffs. of t^n are
#weight-enumerators according to statistics represented by VarList
#For example (assuming that you have CountTilings), to get
#the asymptotic (in n) average of the number of horizontal tiles
#in a domino-tiling of a 4 by n board, type:
#Aaverage(GF1(4,t,H,{[1,2],[2,1]}),t,[H[1,2],H[2,1]],1,n):
Aaverage:=proc(F,t,VarList,i,n) :
Trim1(AsyM(F,t,VarList,[0$(i-1),1,0$(nops(VarList)-i)],n),n,40):
end:


#Cov(F,VarList,i,j): Given a polynomial F in the variable-list VarList
#finds the covariance between the r.v. associated with VarList[i] and
#VarList[j]. For example, try:
#Cov(x+y^2,[x,y],1,1):
Cov:=proc(F,VarList,i,j) local x,y,gu,av1,av2,var1,var2,lu,i1:
x:=VarList[i]: y:=VarList[j]:
lu:={seq(VarList[i1]=1,i1=1..nops(VarList))}:
gu:=subs(lu,y*diff(x*diff(F,x),y))/subs(lu,F):
av1:=subs(lu,x*diff(F,x))/subs(lu,F):
av2:=subs(lu,y*diff(F,y))/subs(lu,F):
var1:=expand(subs(lu,x*diff(x*diff(F,x),x))/subs(lu,F)-av1^2):
var2:=expand(subs(lu,y*diff(y*diff(F,y),y))/subs(lu,F)-av2^2):
normal((gu-av1*av2)/sqrt(var1*var2)):
end:

#Acovariance(F,t,VarList,i,j,n): The asymptotic covariance
#of the r.v. corresponding to the i^th entry in VarList
#and the j^th entry (1<=i<j<=nops(varList))
#where F is a rational function in t whose coeffs. of t^n are
#weight-enumerators according to statistics represented by VarList
#For example (assuming that you have CountTilings), to get
#the asymptotic (in n) covariance of the number of horizontal tiles
#and the number of vertical tiles
#in a domino-tiling of a 4 by n board, type:
#Acovariance(GF1(4,t,H,{[1,2],[2,1]}),t,[H[1,2],H[2,1]],1,2,n):
Acovariance:=proc(F,t,VarList,i,j,n) local gu,gu1,gu2:
gu:=AsyM(F,t,VarList,[0$(i-1),1,0$(j-i-1),1,0$(nops(VarList)-j)],n):

gu1:=Aaverage(F,t,VarList,i,n) :
gu2:=Aaverage(F,t,VarList,j,n) :
Trim1(expand(gu-gu1*gu2),n,40):
end:


#Acorrelation(F,t,VarList,i,j): The asymptotic correlation
#of the r.v. corresponding to the i^th entry in VarList
#and the j^th entry (1<=i<j<=nops(varList))
#where F is a rational function in t whose coeffs. of t^n are
#weight-enumerators according to statistics represented by VarList
#For example (assuming that you have CountTilings), to get
#the asymptotic  correlation of the number of horizontal tiles
#and the number of vertical tiles
#in a domino-tiling of a 4 by n board, type:
#Acorrelation(GF1(4,t,H,{[1,2],[2,1]}),t,[H[1,2],H[2,1]],1,2):
Acorrelation:=proc(F,t,VarList,i,j) local gu,gu1,gu2,n:
gu:=Acovariance(F,t,VarList,i,j,n):

gu1:=Amoment(F,t,VarList,i,2,n):
gu2:=Amoment(F,t,VarList,j,2,n):

gu:=coeff(gu,n,degree(gu,n)):

gu1:=sqrt(coeff(gu1,n,degree(gu1,n))):
gu2:=sqrt(coeff(gu2,n,degree(gu2,n))):

if gu1=0 or gu2=0 then
 RETURN(FAIL):
fi:

Trim1(evalf(gu/gu1/gu2,20),n,20):
end:


#aCorMatrix(F,t,VarList): The asymptotic correlation Matrix
#of the random variables corresponding to the r.v. corresponding
#to those in  VarList
#where F is a rational function in t whose coeffs. of t^n are
#weight-enumerators according to statistics represented by VarList
#For example (assuming that you have CountTilings), to get
#the asymptotic correlation of the number of horizontal tiles
#and the number of vertical tiles
#in a domino-tiling of a 4 by n board, type:
#aCorMatrirx(GF1(4,t,H,{[1,2],[2,1]}),t,[H[1,2],H[2,1]],n):
aCorMatrix:=proc(F,t,VarList) local T,i,j:
for i from 1 to nops(VarList) do
 for j from i+1 to nops(VarList) do
  T[i,j]:=Acorrelation(F,t,VarList,i,j):
 od:
od:
[seq(  
[seq(T[j,i],j=1..i-1),1,seq(T[i,j],j=i+1..nops(VarList))]
 , i=1..nops(VarList))]:

end:




###############End Stuff from StatGF

#Astat(k,ListOfTiles,n): Given a List of rectangular Tiles
#ListOfTiles
#returns: (1) The lists of for the
#asymptotics (in n) for pairs [average,variance] for a random
#tiling by the rectangular tiles
#of Tiles (in order), 
#followed by the av. and variance for the
#Total number of tiles, followed by the asympto. correlation
#matrix for example try:
#Astat(4,[[1,2],[2,1]],n):
Astat:=proc(k,Tiles,n) local t,H,F,VarList,i,s,Fs:


VarList:=[seq(H[op(Tiles[i])],i=1..nops(Tiles))]:

F:=GF1(k,t,H,Tiles):


Fs:=subs({seq(H[op(Tiles[i])]=s,i=1..nops(Tiles))},F):

[[seq([Aaverage(F,t,VarList,i,n) ,Amoment(F,t,VarList,i,2,n)],i=1..nops(VarList))],
[Aaverage(Fs,t,[s],1,n) ,Amoment(Fs,t,[s],1,2,n) ],aCorMatrix(F,t,VarList)]:

end:

#The list of asymptotic averages of the number of tiles
#of tiling a k by n rectangular board by the  tiles in
#the list of rectangular tiles Tiles ([a,b] means a b by a rectangle)
#For example, try: AaveT(2,[[1,2],[2,1]]) 
#followed by the total number of tiles
AaveT:=proc(k,Tiles) local t,H,F,VarList,i,s,Fs,n:

VarList:=[seq(H[op(Tiles[i])],i=1..nops(Tiles))]:

F:=GF1(k,t,H,Tiles):


Fs:=subs({seq(H[op(Tiles[i])]=s,i=1..nops(Tiles))},F):

[seq(coeff(Aaverage(F,t,VarList,i,n),n,1) ,i=1..nops(VarList)), coeff(Aaverage(Fs,t,[s],1,n),n,1) ]:

end:


#AstatV(k,ListOfTiles,n): verbose version of Astat
#For example try:
#AstatV(4,[[1,2],[2,1]],n):
AstatV:=proc(k,Tiles,n) local gu,i:
gu:=Astat(k,Tiles,n):
print(`Consider a random tiling of the rectangular board of dimension`, k, ` by n `):
print(k, `  is the vertical dimension and n the horizontal one`):
print(`Using the tiles in the list of tiles`):
print(Tiles):
print(`here [a,b] means an a by b  rectangular tile `):
print(` a is the vertical dimension and b the horizontal one`):
print(`We are interested in the asymptotics as n goes to infinity`):
print(`of the statistics for the number of tiles`):

for i from 1 to nops(Tiles) do
 print(`The asymptotic (in n)  average number of tiles of type`, Tiles[i] , `in such a random tiling is`):
 print(evalf(coeff(gu[1][i][1],n,1 ),10)*n, `[variance=`, evalf(coeff(gu[1][i][2],n,1),10)*n, `]` ):
# print(`and the variance is`, evalf(coeff(gu[1][i][2],n,1),10)*n):
od:

print(`The average total number of tiles is`):
print(evalf(coeff(gu[2][1],n,1),10)*n, `[variance=`, evalf(coeff(gu[2][2],n,1),10)*n, `]` ):
#print(`and the variance is`, evalf(coeff(gu[2][2],n,1),10)*n ):
print(`Finally, the asymptotic correlation matrix, in the order of appearace of`, Tiles, ` is `):
print(evalf(gu[3],10)):
end:



AstatVs:=proc(K,Tiles,n) local i:
print(`This is info about statistics for radnom tilings of rectangular board`):
print(`of width from 2 to`, K):

for i from 2 to K do
 AstatV(i,Tiles,n):
od:
end: