######################################################################
## HorizTilings.txt: Save this file as  HorizTilings.txt to use it,  #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
#then type: read `HorizTilings.txt`: <Enter>                         #
## Then follow the instructions given there                          #
##                                                                   #
## Written by George Spahn and Doron Zeilberger, Rutgers University ,#
##  DoronZeil@gmail.com      .                                       # 
######################################################################

with(combinat):

Digits:=100:

print(`First Written: May 2021: tested for Maple 18  `):
print():
print(`This is HorizTilings.txt, A Maple package`):
print(`accompanying George Spahn and Doron Zeilberger's  article: `):
print(` Automatic Counting of Generalized Latin Rectangles and Trapezoids `):
print(`-------------------------------`):

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

print(`-------------------------------`):
print(`For a list of the supporting functions type: ezra1();`):
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print():
 print(`-------------------------------`):


ezra1:=proc()
if args=NULL then

print(` The supporting procedures are:  CheckGFzX, CP, EQSVAR, FindTiles,IsGoodState, KidsReg, KidsState,  NuTilings, Rect,   RegToState, Sch, StateToReg, Tilings,  `):
print(` `):

else
 ezra(args):
fi:
end:

ezra:=proc()

if args=NULL then

print(` HorizTILINGS.txt: A Maple package for weighted counting of arbitary tilinds in general domains`):
print(`The Horizontal Tiling procedures are: GFzX, WtTilings `):




elif nargs=1 and args[1]=CheckGFzX then
print(`CheckGFzX(k,TIL,N): Checks that the first N coefficients (in X) of GFzX(k,TIL,z,X) are consistent with the output of WtTilings. Try:`):
print(`CheckGFzX(2,{ {[0,0]}, {[0,1]},{[0,0],[0,1]}} ,10);`):
print(``):
print(`CheckGFzX(2,{{[0, 0]}, {[0, 1]}, {[0, 0], [1, 1]}, {[0, 1], [1, 0]}},10);`):

elif nargs=1 and args[1]=CP then
print(`CP(S): given a set of lattice points, S, finds the point [i,j] such that i is as small as possible and for that i, j is as small as possible . try:`):
print(`CP({[0,0],[1,0],[2,0],[2,2],[-2,-5]});`):


elif nargs=1 and args[1]=EQSVAR then
print(`EQSVAR(k,TIL,z,X,Y): inputs a  positive integer k, a set of tiles TIL (obeying the conditions, see the paper) and symbols z, X,Y, outputs`):
print(`the set of equations, followed by the set of variables that enables one to find the generating function for the weight-enumerator of tilings`):
print(`of Rect(n,k) by the tiles of TIL. Try:`):
print(`EQSVAR(2,{ {[0,0]}, {[0,1]},{[0,0],[0,1]}} ,z,X,Y);`):
print(`EQSVAR(2,{{[0, 0]}, {[0, 1]}, {[0, 0], [1, 1]}, {[0, 1], [1, 0]}},z,X,Y);`):

elif nargs=1 and args[1]=FindTiles then
print(`FindTiles(T,y): inputs a set of tiles T and a height y, outputs the subset of T that contains a point with coordinates y. It also outputs the x-coordinate of those with that y.`):
print(`So the output is a set of pairs. Try:`):
print(`FindTiles(Tiles3(),1);`):

elif nargs=1 and args[1]=GFzX then
print(`GFzX(k,TIL,z,X): The rational function in X, where the MacLaurin coefficients of x^n, that is polynomial in z[tile] for tile in TIL that is the weight-enumerator of tilings with the tiles from TIL of`):
print(`of Rect(n,k). In other words it is WtTilings(Rect(n,k),TIL,z). Try:`):
print(` GFzX(2,{ {[0,0]}, {[0,1]},{[0,0],[0,1]}} ,z,X);`):
print(`GFzX(2,{{[0, 0]}, {[0, 1]}, {[0, 0], [1, 1]}, {[0, 1], [1, 0]}},z,X);`):


elif nargs=1 and args[1]=IsGoodState then
print(`IsGoodState(k,S): Is the state S=[Ome,n1] a good state for Rect(n,k)? Try:`):
print(`IsGoodState(2,[{[0,0],[1,1]},2]);`):
print(`IsGoodState(2,[{[0,0],[1,1]},3]);`):
print(`IsGoodState(2,[{[1,1],[1,2]},3]);`):

elif nargs=1 and args[1]=KidsReg then
print(`KidsReg(Reg,TIL): Given a concrete region, and a set of Tiles, outputs all its children, the pairs [tile,NewRegion] for all the possibilities of`):
print(`placing a tile legally contained in the canonical point, and NewRegion is the new region where the points belonging to the tile have been removed.`):
print(`Try:`):
print(`KidsReg(Rect(6,2),{{[0,0]},{[0,1]}, {[0,0],[0,1]}, {[0,0],[1,1]},{[0,1],[1,0]}});`):
print(`KidsReg(Rect(6,3),{{[0,0]},{[0,1]},{[0,2]}, {[0,0],[0,1]}, {[0,0],[1,1]},{[0,1],[1,0]}, {[0,0],[2,2]},{[0,1],[1,2]} });`):



elif nargs=1 and args[1]=KidsState then

print(`KidsState(k,S,TIL): Given a state S, for the width-k problem `):
print(`and a set of Tiles, TIL, outputs all its children-states, the  [tile,NewState,shift] for all the possibilities of`):
print(`placing a tile legally contained in the canonical point, and NewRegion is the new region where the points belonging to the tile have been joined and asjusted`):
print(`KidsState(2,[{},0], {{[0,0]},{[0,1]}, {[0,0],[0,1]}, {[0,0],[1,1]},{[0,1],[1,0]}});`):
print(`KidsState(3,[{},0],{{[0,0]},{[0,1]},{[0,2]}, {[0,0],[0,1]}, {[0,0],[1,1]},{[0,1],[1,0]}, {[0,0],[2,2]},{[0,1],[1,2]} });`):



elif nargs=1 and args[1]=Rect then
print(`Rect(a,b): The rectangle [0,a-1]x[0,b-1] of a*b points. Try:`):
print(`Rect(2,3);`):


elif nargs=1 and args[1]=Sch then
print(`Sch(k,TIL): inputs a positive integer k and a set of horizontal Tiles, TIL (that lie in 0<=y<=k-1) and outputs a scheme for weight-enumerating`):
print(`Rect(n,k) with the tiles of TIL by the weight product of z[tile] over all tiles that participate. Try:`):
print(`Sch(2,{{[0,0]},{[0,1]},{[0,0],[0,1]}});`):

elif nargs=1 and args[1]=StateToReg then
print(`StateToReg(k,S,N1): translates the abstract state S=[ome,n1] into a concrete region a subset of Rect(N1,k). Try:`):
print(`StateToReg(3,[{},0],10);`):
print(`StateToReg(3,[{[0,0],[2,2]},3],10);`):



elif nargs=1 and args[1]=RegToState then
print(`RegToState(R,k,n1): Given a region R that is a subset of Rect(n1,k) for some n1, finds the corresponding state. Try:`):
print(`RegToState(Rect(10,3),3,10);`):

elif nargs=1 and args[1]=Tilings then
print(`Tilings(S,T): Given a set of lattice points S, and a set of tiles (where one can do horizontal translation only, so tiles that are vertical translations are distinct), outputs the set of all tilings as a set of `):
print(`[RegionCoveredByTile,Tile] `):
print(`Tilings(Rect(2,2),{ {[0,0],[0,1]},{[0,0]},{[0,1]} });`):


elif nargs=1 and args[1]=WtTilings then
print(`WtTilings(S,T,z): Given a set of lattice points S, and a set of tiles, outputs the weight-enumerator`):
print(`according to the index variables z[tile], by clever Dynmaical programming.`):
print(`WtTilings(Rect(3,3),{{[0,0]}, {[0,1]}, {[0,2]}, {[0,0],[1,1],[2,2]}},z);`):


elif nargs=1 and args[1]=WtTilingsBF then
print(`WtTilingsBF(S,T,z): Given a set of lattice points S, and a set of tiles, outputs the weight-enumerator`):
print(`according to the index variables z[tile], by BRUTE FORCE. For checking purposes only. Should be the`):
print(`same as WtTiling(S,T,z). Try:`):
print(`WtTilingsBF(Rect(2,2),{ {[0,0],[0,1]},{[0,0]},{[0,1]}},z);`):


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

fi:



end:




#Rect(a,b): The rectangle [0,a-1]x[0,b-1] of a*b points. Try:
#Rect(2,3);
Rect:=proc(a,b) local i,j:
{seq(seq([i,j],i=0..a-1),j=0..b-1)}:
end:


#Tilings(S,T): Given a set of lattice points S, and a set of tiles, outputs the set of all tilings as a set of 
#[RegionCoveredByTile,Tile] 
#Tilings(Rect(2,2),{ {[0,0],[0,1]},{[0,0]},{[0,1]} });
Tilings:=proc(S,T) local gu,pt,T1,t,pt1,S1,mu1,mu,tS:
option remember:

if {seq(IsLegalTile(t), t in T)}<>{true} then
print(`Not all tiles in `, T, `are legal `):
 RETURN(FAIL):
fi:

if S={} then
 RETURN({{}}):
fi:

pt:=CP(S):

T1:=FindTiles(T,pt[2]):

if T1={} then
 RETURN({}):
fi:

gu:={}:

for t in T1 do
tS:={seq([pt[1]-t[2],0]+pt1,pt1 in t[1])}:

if tS minus S={} then
 S1:=S minus tS:
 mu:=Tilings(S1,T):
 gu:=gu union {seq({op(mu1),[t[1],tS]},mu1 in mu)}:
fi:

od:

gu:

end:

#CP(S): given a set of lattice points, S, finds the point [i,j] such that i is as small as possible and for that i, j is as small as possible . try:
#CP({[0,0],[1,0],[2,0],[2,2],[-2,-5]});
CP:=proc(S) local j0,i0,s,Sb:
i0:=min(seq(s[1],s in S)):
Sb:={}:
for s in S do

if s[1]=i0 then
 Sb:=Sb union {s}:
fi:

od:
j0:=min(seq(s[2], s in Sb)):

[i0,j0]:

end:


#FindTiles(T,y): inputs a set of tiles T and a height y, outputs the set of pairs [t,x] where
#t contains a point of the form [x,y]. Try:
#FindTiles(Tiles3(),1);
FindTiles:=proc(T,y) local t,t1,gu,i:

gu:={}:

for t in T do
 if member(y,{seq(t1[2], t1 in t)}) then
  for i from 1 to nops(t) while t[i][2]<>y do od:
  gu:=gu union {[t,t[i][1]]}:
 fi:
od:

gu:

end:


#IsLegalTile(S): inputs a set of lattice points in the plane and decides whether it is a legal tile.
#(i) The intersection with every vertical line y=i has at most one element
#(ii): The smallest x-coordinate is 0. Try:
#IsLegalTile({[0,1],[1,0]});
IsLegalTile:=proc(S) local s:
option remember:
if nops([seq(s[2],s in S)]) <>nops({seq(s[2],s in S)}) then
 RETURN(false):
fi:


true:

end:


#NuTilings(S,T): Given a set of lattice points S, and a set of tiles, outputs the NUMBER of tilings in terms of the indexed variables z[tile]
#Try:
#NuTilings(Rect(2,2),{ {[0,0],[1,0]},{[0,0],[0,1]},{[0,0]}});

NuTilings:=proc(S,T) local gu,pt,T1,t,pt1,S1,tS:
option remember:

if {seq(IsLegalTile(t), t in T)}<>{true} then
print(`Not all tiles in `, T, `are legal `):
 RETURN(FAIL):
fi:

if S={} then
 RETURN(1):
fi:

pt:=CP(S):

T1:=FindTiles(T,pt[2]):

if T1={} then
 RETURN(0):
fi:

gu:=0:

for t in T1 do
tS:={seq([pt[1]-t[2],0]+pt1,pt1 in t[1])}:

if tS minus S={} then
 S1:=S minus tS:

 gu:=gu+NuTilings(S1,T):
fi:

od:

gu:

end:



#WtTilings(S,T,z): Given a set of lattice points S, and a set of tiles, outputs the weight-enumerator set of all tilings in terms of the indexed variables z[tile]
#Try:
#WtTilings(Rect(2,2),{ {[0,0],[1,0]},{[0,0],[0,1]},{[0,0]}},z);
WtTilings:=proc(S,T,z) local gu,pt,T1,t,pt1,S1,tS:
option remember:

if {seq(IsLegalTile(t), t in T)}<>{true} then
print(`Not all tiles in `, T, `are legal `):
 RETURN(FAIL):
fi:

if S={} then
 RETURN(1):
fi:


pt:=CP(S):

T1:=FindTiles(T,pt[2]):

if T1={} then
 RETURN(0):
fi:

gu:=0:

for t in T1 do
tS:={seq([pt[1]-t[2],0]+pt1,pt1 in t[1])}:

if tS minus S={} then
 S1:=S minus tS:

 gu:=expand(gu+z[t[1]]*WtTilings(S1,T,z)):
fi:

od:

gu:

end:

#KidsReg(Reg,TIL): Given a concrete region, and a set of Tiles, outputs all its children, the pairs [tile,NewRegion] for all the possibilities of
#placing a tile legally contained in the canonical point, and NewRegion is the new region where the points belonging to the tile have been removed.
#Try:
#KidsReg(Rect(6,2),{{[0,0]},{[0,1]}, {[0,0],[0,1]}, {[0,0],[1,1]},{[0,1],[1,0]}});
#KidsReg(Rect(6,3),{{[0,0]},{[0,1]},{[0,2]}, {[0,0],[0,1]}, {[0,0],[1,1]},{[0,1],[1,0]}, {[0,0],[2,2]},{[0,1],[1,2]} });
KidsReg:=proc(Reg,TIL) local gu, Reg1, TIL1, pt, t, tS,pt1:

if {seq(IsLegalTile(t), t in TIL)}<>{true} then
print(`Not all tiles in `, TIL, `are legal `):
 RETURN(FAIL):
fi:

if Reg={} then
 RETURN({}):
fi:


pt:=CP(Reg):

TIL1:=FindTiles(TIL,pt[2]):

if TIL1={} then
 RETURN({}):
fi:

gu:={}:

for t in TIL1 do
tS:={seq([pt[1]-t[2],0]+pt1,pt1 in t[1])}:

if tS minus Reg={} then
 Reg1:=Reg minus tS:

gu:=gu union {[t[1],Reg1]}:
fi:

od:

gu:


end:


#IsGoodState(k,S): Is the state S=[Ome,n1] a good state for Rect(n,k)? Try:
#IsGoodState(2,{[0,0],[1,1]},2]);
#IsGoodState(2,{[0,0],[1,1]},3]);
#IsGoodState(3,{[1,1],[1,2]},3]);
IsGoodState:=proc(k,S) local Ome,n1,pt,i:


if not type(S,list) and nops(S)=2 then
 RETURN(false):
fi:

Ome:=S[1]:
n1:=S[2]:

if Ome={} and n1=0 then
 RETURN(true):
fi:

if not (type(Ome,set) and type(n1,integer) and n1>=0) then
 RETURN(false):
fi:

if {seq(pt[2],pt in Ome)} minus {seq(i,i=0..k-1)}<>{} then
 RETURN(false):
fi:

if  min(seq(pt[1],pt in Ome))<>0 then
 RETURN(false):
fi:

if {seq([n1-1,i],i=0..k-1)} minus Ome={} then
 RETURN(false):
fi:

true:

end:


#StateToReg(k,S,N1): translates the abstract state S=[ome,n1] into a concrete region a subset of Rect(N1,k). Try:
#StateToReg(3,[{},0],10);
#StateToReg(3,[{[0,0],[2,2]},3],10);
StateToReg:=proc(k,S,N1) local i,j,Ome,n1:

if not IsGoodState(k,S) then
 print(S, ` is not a good state for k=`,k):
 RETURN(FAIL):
fi:

if not type(S,list) then
 RETURN(FAIL):
fi:


Ome:=S[1]:
n1:=S[2]:

Ome union {seq(seq([i,j],j=0..k-1),i=n1..N1-1)}:


end:


#RegToState(R,k,N1): Given a region R that is a subset of Rect(N1,k) for some N1, finds the corresponding state and shift. Try:
#RegToState(Rect(10,3),3,10);
RegToState:=proc(R,k,N1) local j,n1,Ome,katan,J,i,pt:
if R minus Rect(N1,k)<>{} then
 RETURN(FAIL):
fi:

for J from N1-1 by -1 while {seq(seq([i,j],j=0..k-1),i=J..N1-1)} minus R={} do od:

n1:=J+1:

if n1=0 then
 RETURN([{},0]):
fi:

Ome:=R minus  {seq(seq([i,j],j=0..k-1),i=n1..N1-1)}:


if Ome={} then
 RETURN([{},0],n1):
fi:

katan:=min(seq(pt[1],pt in Ome)):

Ome:={seq(pt-[katan,0],pt in Ome)}:

n1:=n1-katan:

[Ome,n1],katan:
end:






#KidsState(k,S,TIL): Given a state S, for the width-k problem 
#and a set of Tiles, TIL, outputs all its children-states, the  [tile,NewState,shift] for all the possibilities of
#placing a tile legally contained in the canonical point, and NewRegion is the new region where the points belonging to the tile have been joined and asjusted
#KidsState(2,[{},0], {{[0,0]},{[0,1]}, {[0,0],[0,1]}, {[0,0],[1,1]},{[0,1],[1,0]}});
#KidsState(3,[{},0],{{[0,0]},{[0,1]},{[0,2]}, {[0,0],[0,1]}, {[0,0],[1,1]},{[0,1],[1,0]}, {[0,0],[2,2]},{[0,1],[1,2]} });
KidsState:=proc(k,S,TIL) local ti,N1,R,mu,mu1,gu,R1,S1,sh,pt:

N1:=S[2]+max(seq(seq(pt[1],pt in ti),ti in TIL))+10:

R:=StateToReg(k,S,N1):

mu:=KidsReg(R,TIL):

gu:={}:

for mu1 in mu do
 ti:=mu1[1]:
 R1:=mu1[2]:
 S1:=RegToState(R1,k,N1):
 sh:=S1[2]:
 S1:=S1[1]:
 gu:=gu union {[ti,S1,sh]}:
od:

gu:

end:


#Sch(k,TIL): inputs a positive integer k and a set of horizontal Tiles, TIL (that lie in 0<=y<=k-1) and a symbol z, and outputs a scheme for weight-enumerating
#Rect(n,k) with the tiles of TIL by the weight product of z[tile] over all tiles that participate. Try:
#Sch(2,{{[0,0]},{[0,1]},{[0,0],[0,1]}});
Sch:=proc(k,TIL) local gu,t,pt,i,StillToDo,AlreadyDone,S,lu:


if {seq(IsLegalTile(t), t in TIL)}<>{true} then
print(`Not all tiles in `, TIL, `are legal `):
 RETURN(FAIL):
fi:

for t in TIL do
 if {seq(pt[2],pt in t)} minus {seq(i,i=0..k-1)}<>{} then
 print(`The tiles must have points with y-coordonate between 0 and`, k-1):
 RETURN(FAIL):
 fi:
od:

gu:=[]:

StillToDo:={[{},0]}:
AlreadyDone:={}:

 while StillToDo<>{} do
  S:=StillToDo[1]:
  AlreadyDone:= AlreadyDone union {S}:
  lu:=KidsState(k,S,TIL):
  gu:=[op(gu),[S,lu]]:
  StillToDo:=StillToDo union {seq(lu[i][2],i=1..nops(lu))} minus AlreadyDone:
od:

gu:

end:



#EQSVAR(k,TIL,z,X,Y): inputs a  positive integer k, a set of tiles TIL (obeying the conditions, see the paper) and symbols z, X,Y, outputs
#the set of equations, followed by the set of variables that enables one to find the generating function for the weight-enumerator of tilings
#of Rect(n,k) by the tiles of TIL. Try:
#EQSVAR(2,{ {[0,0]}, {[0,1]},{[0,0],[0,1]}} ,z,X,Y);
EQSVAR:=proc(k,TIL,z,X,Y) local Sc,eq,var,Sta,i,eq1,j,lu:
Sc:=Sch(k,TIL):

var:={}:
eq:={}:

for i from 1 to nops(Sc) do
 Sta:=Sc[i][1]:

 var:=var union {Y[Sta]}:
 eq1:=Y[Sta]:

if Sta=[{},0] then
 eq1:=eq1-1:
fi:


for j from 1 to nops(Sc[i][2]) do
 lu:=Sc[i][2][j]:
 eq1:=eq1-z[lu[1]]*X^lu[3]*Y[lu[2]]:

od:

eq:=eq union {eq1}:
od:

eq,var:


end:

#GFzX(k,TIL,z,X): The rational function in X, where the MacLaurin coefficients of x^n, that is polynomial in z[tile] for tile in TIL that is the weight-enumerator of tilings with the tiles from TIL of
#of Rect(n,k). In other words it is WtTilings(Rect(n,k),TIL,z). Try:
#GFzX(2,{ {[0,0]}, {[0,1]},{[0,0],[0,1]}} ,z,X);
GFzX:=proc(k,TIL,z,X) local Y:
normal(subs(solve(EQSVAR(k,TIL,z,X,Y)),Y[[{},0]])):
end:


#CheckGFzX(k,TIL,N): Checks that the first N coefficients (in X) of GFzX(k,TIL,z,X) are consistent with the output of WtTilings. Try:
#CheckGFzX(2,{ {[0,0]}, {[0,1]},{[0,0],[0,1]}} ,10);
CheckGFzX:=proc(k,TIL,N)  local gu,z,X,i:

gu:=GFzX(k,TIL,z,X):

gu:=taylor(gu,X=0,N+1):

evalb({seq(expand(WtTilings(Rect(i,k),TIL,z) -expand(coeff(gu,X,i))),i=1..N)}={0}):

end: