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

with(combinat):

Digits:=100:

print(`First Written: July 2021: tested for Maple 18  `):
print():
print(`This is GenLatinRecs.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(`-------------------------------`):
print(`For  help with the main procedures from HorizTilings.txt, and a list of the MAIN functions there,`):
print(` type "ezraH();". For specific help type "ezraH(procedure_name);" `):

print(`-------------------------------`):
print(`For  a list of the supporing procedures from HorizTilings.txt`):
print(` type "ezraH1();". For specific help type "ezraH(procedure_name);" `):


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:  EQS2, EQS3, GatherTiles, GLR, GLRslow, GLRseqBF,  IsGoodR, IsTile, Matkn `):
print(` `):

else
 ezra(args):
fi:
end:

ezra:=proc()

if args=NULL then

print(` GenLatinRecs.txt: A Maple package for counting generalized Latin rectangles`):
print(`The MAIN procedures are:  Kernel2, Kernel3, GenDerSeq, GenDerSeqI, GLR3seq, GLR3seqI, GLRseq ,  Paper2, Paper3`):



elif nargs=1 and args[1]=EQS2 then
print(`EQS2(R,x,X,Z): Given a set of integers R, and variables x and X outputs the rational function f(x,X) such that the number of permutations pi of length n`):
print(`such that pi[i]-i is NOT in R is the coefficient of X^n in Int(exp(-x)*R(x,X),x=0..infinity), or equivalenty applyting the umbra x^n->n! to it. `):
print(`To get the classical derangement numbers, type:`):
print(`EQS2({0},x,X,Z);`):

elif nargs=1 and args[1]=EQS3 then
print(`EQS3(R12,R13,R23,x1,x2,x3,x23,X,Z): Given three sets of integers R12,R13,R23,  and variables x1,x2,x3,x23 and X outputs the system of equations that enables`):
print(`the computation of  Kernel3(R12,R13,R23,x1,x2,x3,x23,X) and GLR3seq(R12,R13,R23,N);.  Try:`):
print(`EQS3({0},{0},{0},x1,x2,x3,x23,X,Z);`):

elif nargs=1 and args[1]=GatherTiles then
print(`GatherTiles(k,R): inputs a positive integer k and a set of restictions {[[a,b],S]} finds the corresponding set of tiles. try:`):
print(`GatherTiles(3,{[[1,2],{-1,0,1}], [[2,3],{-1,0,1}], [[1,3],{-2,0,2}]  });`):


elif nargs=1 and args[1]=GenDerSeq then
print(`GenDerSeq(R,N): Given a set of integers R (0 and negatives are OK) and a positive integer N, outputs the first N terms of the sequence enumerating`):
print(`the number of permutations of pi, {1,...,n} such that pi[i]-i is NOT in R. `):
print(`For the first 30 terms of the number of derangements, type:`):
print(`GenDerSeq({0},30); `):
print(`For the first 30 terms of the number of (straight) Menage numbers, type:`):
print(`GenDerSeq({0,1},30); `):


elif nargs=1 and args[1]=GenDerSeqI then
print(`GenDerSeqI(R,N): same as GenDerSeq(R,N) (q.v.), but using iteration rather than solving`):
print(`For the first 30 terms of the number of (straight) Menage numbers, type:`):
print(`GenDerSeqI({0,1},30); `):

elif nargs=1 and args[1]=GLR then
print(`GLR(n,k,R): Given a set of pairs,R {[[a,b],S]} where 1<=a<b<=k and S is a set of integers`):
print(`outputs the set of reduced k by n matrices where every row is a permutation of {1,...,n}, the first row is [1,...,n] and obeying the condition`):
print(`M[a][i]<>M[b][i+s] for each s in whenver both i and i+s are between 1 and n . For example, to get the set of reduced Latin`):
print(`3 by 5 rectangles, type`):
print(`GLR(5,3,{ [[1,2],{0}], [[1,3],{0}], [[2,3],{0}] }); `):
print(`To get the set of derangements of length 5, type:`):
print(`GLR(5,2,{[[1,2],{0}]});`):
print(`To get the set of Menage permutations of length 5, type:`):
print(`GLR(5,2,{[[1,2],{0,1}]});`):
print(`To get the set of 3 by 5 super Latin Rectangles, type:`):
print(`GLR(5,3,{[[1,2],{-1,0,1}], [[2,3],{-1,0,1}], [[1,3],{-2,0,2}]  });`):


elif nargs=1 and args[1]=GLRseq then
print(`GLRseq(k,R,N): Given a set of pairs,R {[[a,b],S]} where 1<=a<b<=k and S is a set of integers`):
print(`outputs the sequence enumerating reduced k by n matrices where every row is a permutation of {1,...,n}, the first row is [1,...,n] and obeying the condition`):
print(`M[a][i]<>M[b][i+s] for each s in S, in whenver both i and i+s are between 1 and n .`):

print(``):
print(` For example, to get the first 10 terms of the sequence enumerating reduced SUPER Latin`):
print(`3 by n rectangles, type`):
print(`GLRseq(3,{ [[1,2],{0,1,-1}], [[1,3],{0,2,-2}], [[2,3],{0,1,-1}] },7); `):
print(``):
print(` For example, to get the first 10 terms of the sequence enumerating reduced Latin`):
print(`3 by n rectangles, type`):
print(`GLRseq(3,{ [[1,2],{0}], [[1,3],{0}], [[2,3],{0}] },10); `):

print(` To get the first 10 terms of the sequence enumerating  super Latin`):
print(`3 by n rectangles, type`):
print(`GLRseq(3,{ [[1,2],{0}], [[1,3],{0}], [[2,3],{0}] },10); `):

print(`GLRseq(3,{[[1,2],{-1,0,1}], [[2,3],{-1,0,1}], [[1,3],{-2,0,2}]  },5);`):

print(``):
print(`To get the first 7 terms  of derangements of length n, type:`):
print(`GLRseq(2,{[[1,2],{0}]},7);`):
print(`To get the first seven terms of the sequence enumeratingMenage permutations of length n, type:`):
print(`GLRseq(2,{[[1,2],{0,1}]},7);`):



elif nargs=1 and args[1]=GLR3seq then
print(`GLR3seq(R12,R13,R23,N): Given  three sets of integers R12, R13,R23  and a positive integer N`):
print(`outputs the first N terms of the sequence enumerating 3 by n matrices`):
print(`where each row is a permutation of [1,...,n] and the first row is [1,...,n] (to get the unreduced number simply multiply by n!)`):
print(` M[1,i]<>M[2,j] whenever j-i is in R12`):
print(` M[1,i]<>M[3,j] whenever j-i is in R13`):
print(` M[2,i]<>M[3,j] whenever j-i is in R23`):
print(`It is the same output as`):
print(`GLRseq(3,{ [[1,2],R12], [[1,3],R13], [[2,3],R23]},N);`):
print(`but hopefully faster.`):
print(`For example, to get the first 30 terms of the sequence enuerating reduced Latin Rectangles with 3 rows, try:`):
print(`GLR3seq({0},{0},{0},30);`):


elif nargs=1 and args[1]=GLR3seqI then
print(`GLR3seqI(R12,R13,R23,N): same as GLR3seq(R12,R13,R23,N) (q.v.) but using iteration rather than solving (for cases where solving takes too much time and/or space`):
print(`For example, to get the first 30 terms of the sequence enuerating reduced Latin Rectangles with 3 rows, try:`):
print(`GLR3seqI({0},{0},{0},30);`):


elif nargs=1 and args[1]=GLRseqBF then
print(`GLRseqBF(k,R,N): Same as GLRseq(k,R,N) (q.v.), but by brute force. FOR CHECKING PURPOSES ONLY.`):

print(`3 by n rectangles, type`):
print(`GLRseqBF(3,{ [[1,2],{0}], [[1,3],{0}], [[2,3],{0}] },7); `):
print(`To get the first 7 terms  of derangements of length n, type:`):
print(`GLRseqBF(2,{[[1,2],{0}]},7);`):
print(`To get the first seven terms of the sequence enumeratingMenage permutations of length n, type:`):
print(`GLRseqBF(2,{[[1,2],{0,1}]},7);`):

elif nargs=1 and args[1]=GLRslow then
print(`GLRslow(n,k,R): Like GLR(n,k,R) but the slow way. For checking purposes only.`):
print(`outputs the set of reduced k by n matrices where every row is a permutation of {1,...,n}, the first row is [1,...,n] and obeying the condition`):
print(`M[a][i]<>M[b][i+s] for each s in S, whenver both i and i+s are between 1 and n . For example, to get the set of reduced Latin`):
print(`3 by 5 rectangles, type`):
print(`GLRslow(5,3,{ [[1,2],{0}], [[1,3],{0}], [[2,3],{0}] }); `):
print(`To get the set of derangements of length 5, type:`):
print(`GLRslow(5,2,{[[1,2],{0}]});`):
print(`To get the set of Menage permutations of length 5, type:`):
print(`GLRslow(5,2,{[[1,2],{0,1}]});`):


elif nargs=1 and args[1]=IsGoodR then
print(`IsGoodR(M,R) : Given a matrix of integers M say of dimensions k by N and a set of restrictions R, decides whether M is good`):
print(`Try: `):
print(`IsGoodR([[1,2,3],[2,3,1],[3,1,2]],{ [[1,2],{0}],[[1,3],{0}],[[2,3],{0}] });`):


elif nargs=1 and args[1]=IsTile then
print(`IsTile(T,k,R): Given a set of lattice points T with y-coordinates all different and in {0,...,k-1} decides whether T is a legal tile in the enumeration of`):
print(`generakized Latin rectangles with condition set  R.It also returns the number of edges.  Try:`):
print(`IsTile({[0,0],[0,2]},3,{ [[1,2],{0}], [[1,3],{0}], [[2,3],{0}] }); `):
print(`IsTile({[0,0],[1,1],[2,2]},3,{[[1,2],{-1,0,1}], [[2,3],{-1,0,1}], [[1,3],{-2,0,2}]  }); `):


elif nargs=1 and args[1]=Kernel2 then
print(`Kernel2(R,x,X): Given a set of integers R, and variables x and X outputs the rational function f(x,X) such that the number of permutations pi of length n`):
print(`such that pi[i]-i is NOT in R is the coefficient of X^n in Int(exp(-x)*R(x,X),x=0..infinity), or equivalenty applyting the umbra x^n->n! to it. `):
print(`To get the classical derangement numbers, type:`):
print(`Kernel2({0},x,X);`):



elif nargs=1 and args[1]=Kernel3 then
print(`Kernel3(R12,R13,R23,x1,x2,x3,x23,X): Given three sets of integers R12, R13, R23, and variables x2, x3, x23, and X outputs the rational function f(x2,x3,x23,X) such that `):
print(`enables the computation of the number of 3 by n matrices`):
print(`where each row is a permutation of [1,...,n] and the first row is [1,...,n] (to get the unreduced number simply multiply by n!)`):
print(` M[1,i]<>M[2,j] whenever j-i is in R12`):
print(` M[1,i]<>M[3,j] whenever j-i is in R13`):
print(` M[2,i]<>M[3,j] whenever j-i is in R23`):
print(`via the umbral operation`):
print(` x1^a1*x2^a2*x3^a2*x23^a23 ->binomial(n-a1,a23)*a23!*a2!*a3! `):
print(`To get the Kernel for enumerating 3 by reduced Latin rectangles. Try:`):
print(`Kernel3({0},{0},{0},x1,x2,x3,x23,X);`):

elif nargs=1 and args[1]=Matkn then
print(`Matkn(k,n): All the k by n matrices all whose rows are permutations of {1,...,n} and first row is [1,...,n] . Try:`):
print(`Matkn(3,4);`):


elif nargs=1 and args[1]=Paper2 then
print(`Paper2(K,N): inputs positive integers K and N and outputs an article with the first N terms of the sequenes enumerating generalized derangelents`):
print(`i.e. permutations pi of {1,...n} such that pi[i]-i is not in S for all subsets S of {-K,...K}. try:`):
print(`Paper2(1,50);`):


elif nargs=1 and args[1]=Paper3 then
print(`Paper3(K,N): inputs positive integers K and N and outputs an article with the first N terms of the sequenes enumerating 3-rowed matrices such that`):
print(`the bottom row is [1,...,n] and M[1,i]<>M[1,j], M[1,i]<>M[3,j], M[2,i]<>M[3,j] whenever j-i belongs to the set S`):
print(`for all subsets S of {-K,...K}. try:`):
print(`Paper3(1,20);`):

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

fi:



end:




###Start From HotizTlings.txt

ezraH1:=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:

ezraH:=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:



#EQSVARold(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:
#EQSVARold(2,{ {[0,0]}, {[0,1]},{[0,0],[0,1]}} ,z,X,Y);
EQSVARold:=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:



#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:=0:

if Sta=[{},0] then
 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 {Y[Sta]=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:


###End From HotizTlings.txt


#IsGoodR(M,R) : Given a matrix of integers M say of dimensions k by N and a set of restrictions R, decides whether M is good
#Try:
#IsGoodR([[1,2,3],[2,3,1],[3,1,2]],{ [[1,2],{0}],[[1,3],{0}],[[2,3],{0}] });
IsGoodR:=proc(M,R) local n,k,r,a,b,S,s,i1,i:

k:=nops(M):
n:=nops(M[1]):


if {seq(convert(M[i],set),i=1..k)}<>{{seq(i1,i1=1..n)}} then
print(`not all rows are permutations`):
 RETURN(false):
fi:

if k=1 then
 RETURN(true):
fi:


for r in R do
a:=r[1][1]:
b:=r[1][2]:
if not a>=1 and a<b and b<=k then
 print(r, ` is a wrong resriction`):
 RETURN(FAIL):
fi:

 S:=r[2]:
for i from 1 to n do
 for s in S do
   if i+s>=1 and i+s<=n and M[a][i]=M[b][i+s] then
    RETURN(false):
  fi:
 od:
od:
od:
true:
end:

 
#Matkn(k,n): All the k by n matrices all whose rows are permutations of {1,...,n} and first row is [1,...,n]
Matkn:=proc(k,n) local mu,mu1,lu1,lu,i1:
option remember:

lu:=permute(n):

if k=1 then
 RETURN({seq([[seq(i1,i1=1..n)]],lu1 in lu)}):
fi:

mu:= Matkn(k-1,n):

{seq(seq([op(mu1),lu1],mu1 in mu),lu1 in lu)}:

end:

 



#GLRslow(n,k,R): Given a set of pairs,R {[[a,b],S]} where 1<=a<b<=k and S is a set of integers
#outputs the set of reduced k by n matrices where every row is a permutation of {1,...,n}, the first row is [1,...,n] and obeying the condition
#M[a][i]<>M[b][i+s] for each s in S, whenver both i and i+s are between 1 and n . For example, to get the set of reduced Latin
#3 by 5 rectangles, type
#GLR(5,3,{ [[1,2],{0}], [[1,3],{0}], [[2,3],{0}] }); Done the SLOW way
#To get the set of derangements of length 5, type:
#GLRslow(5,2,{[[1,2],{0}]);
#To get the set of Menage permutations of length 5, type:
#GLRsslow(5,2,{[[1,2],{0,1}]);
GLRslow:=proc(n,k,R) local mu,mu1,gu:
option remember:


mu:=Matkn(k,n):

gu:={}:

for mu1 in mu do
if IsGoodR(mu1,R) then
  gu:=gu union {mu1}:
fi:
od:
gu:

end:



#GLR(n,k,R): Given a set of pairs,R {[[a,b],S]} where 1<=a<b<=k and S is a set of integers
#outputs the set of reduced k by n matrices where every row is a permutation of {1,...,n}, the first row is [1,...,n] and obeying the condition
#M[a][i]<>M[b][i+s] for each s in whenver both i and i+s are between 1 and n . For example, to get the set of reduced Latin
#3 by 5 rectangles, type
#GLR(5,3,{ [[1,2],{0}], [[1,3],{0}], [[2,3],{0}] }); 
#To get the set of derangements of length 5, type:
#GLR(5,2,{[[1,2],{0}]);
#To get the set of Menage permutations of length 5, type:
#GLR(5,2,{[[1,2],{0,1}]);
GLR:=proc(n,k,R) local R1,NR,i,pa,mu,mu1,gu,M,lu,lu1:
option remember:

if k=1 then
 RETURN({[[seq(i,i=1..n)]]}):
fi:

R1:={}:

NR:={}:

for pa in R do
 if pa[1][2]<k then
  R1:=R1 union {pa}:
 else
  NR:=NR union {pa}:
 fi:
od:


mu:=GLR(n,k-1,R1):
lu:=permute(n):

gu:={}:

for mu1 in mu do
 for lu1 in lu do
  M:=[op(mu1),lu1]:
  if IsGoodR(M,NR) then
    gu:=gu union {M}:
  fi:
od:
od:
gu:


end:



#Neighs(V,ED1,v): Given a set of vertices V and a set of edges ED1 and a vertex v in V oututs all its neighbors, Try:
#Neighs({1,2,3,4,5},{{1,2},{1,3},{4,5}},1);
Neighs:=proc(V,ED1,v)  local u,N:

if not member(v,V) then
 RETURN(FAIL):
fi:

N:={}:

for u in V minus {v} do
 if member({u,v},ED1) then
  N:=N union {u}:
 fi:
od:
N:
end:

  
#CC(V,ED1,v): The connected component of v in the graph (V,ED1). Try:
#CC({1,2,3,4,5},{{1,2},{1,3},{4,5}},1);
CC:=proc(V,ED1,v) local C,C1,C2,v1:

C:={v}:

C1:=C union {seq(op(Neighs(V,ED1,v1)),v1 in C)}:


while C<>C1 do
C2:=C1 union{seq(op(Neighs(V,ED1,v1)),v1 in C1)}:
C:=C1:
C1:=C2:
od:

C1:

end:

#IsCon(V,ED1): Is the graph (V,ED1) connected?
IsCon:=proc(V,ED1) 

if  V={} then
 RETURN(true):
fi:

evalb(CC(V,ED1,V[1])=V):
end:

#IsTile(T,k,R): Given a set of lattice points T with y-coordinates all different and in {0,...,k-1} decides whether T is a legal tile in the enumeration of
#generakized Latin rectangles with condition set  R. Try:
#IsTile({[0,0],[0,2]},3,{ [[1,2],{0}], [[1,3],{0}], [[2,3],{0}] });  It also returns the number of edges
IsTile:=proc(T,k,R) local i,j,t,ED1,ed1,a,b,TA,katan,gadol:
if {seq(t[2],t in T)} minus {seq(i,i=0..k-1)}<>{} then
 print(`The y-coordinates are not all between 0 and`, k-1):
 RETURN(false):
fi:

for i from 1 to nops(T) do
for j from i+1 to nops(T) do
 if T[i][2]=T[j][2] then
  print(T[i],T[j], `have the same y-coordinate `):
  RETURN(false):
fi:
od:
od:

for i from 1 to k do
 for j from i+1 to k do
  TA[[i,j]]:={}:
 od:
od:

for i from 1 to nops(R) do
TA[R[i][1]]:=R[i][2]:
od:


ED1:={}:

#finds the set of edges in the graph whose vertices are the members of the tile T

for i from 1 to nops(T) do
 for j from i+1 to nops(T) do
  ed1:={T[i],T[j]}:
  a:=min(T[i][2],T[j][2]):
  b:=max(T[i][2],T[j][2]):
  if T[i][2]=a then
    katan:=T[i]:
    gadol:=T[j]:
  else
    katan:=T[j]:
    gadol:=T[i]:
  fi:

 if member(gadol[1]-katan[1],TA[[a+1,b+1]]) then
  ED1:=ED1 union {ed1}:
 fi:
od:
od:

[IsCon(T,ED1),nops(ED1)]:

end:



#CanTile(T) : The translation ot the tile T such that the smallest x-coordinate is 0. Try:
#CanTile({[-2,0],[1,1]});
CanTile:=proc(T) local pt,katan:
katan:=min(seq(pt[1],pt in T)):
{seq(pt-[katan,0],pt in T)}:
end:



#GatherTiles(k,R): inputs a positive integer k and a set of restictions {[[a,b],S]} finds the corresponding set of tiles. try:
#FindTiles(3,{[[1,2],{-1,0,1}], [[2,3],{-1,0,1}], [[1,3],{-2,0,2}]  });
GatherTiles:=proc(k,R) local TA,gu12,gu13,gu23, gu1213,gu,i,j,s,s1,s2,gu1223,gu1323,gu1,mu,ti,yofee:

if not type(k, integer) and k>=1 then
print(` bad input`):
 RETURN(FAIL):
fi:

if k>=4 then
 print(`Not yet implemented`):
fi:


for i from 1 to k do
 for j from i+1 to k do
  TA[[i,j]]:={}:
 od:
od:


for i from 1 to nops(R) do
TA[R[i][1]]:=R[i][2]:
od:

if {seq(R[i][1],i=1..nops(R))} minus {seq(seq([i,j],j=i+1..k),i=1..k)}<>{} then
 print(R, `is bad input`):
 RETURN(FAIL):
fi:

if k=2 then
yofee:={seq({[0,0],[s,1]},s in R[1][2]) ,{[0,0]},{[0,1]}  }:
RETURN({seq(CanTile(ti), ti in yofee)} ):
fi:


if k=3 then

#tiles with one edge
gu12:={seq({[0,0],[s,1]},s in TA[[1,2]])}:
gu13:={seq({[0,0],[s,2]},s in TA[[1,3]])}:
gu23:={seq({[0,1],[s,2]},s in TA[[2,3]] )}:

#tiles with two or three edges

gu1213:={seq(seq({[0,0],[s1,1],[s2,2]},s1 in TA[[1,2]]),s2 in TA[[1,3]])}:
gu1223:={seq(seq({[0,1],[-s1,0],[s2,2]},s1 in TA[[1,2]]),s2 in TA[[2,3]])}:
gu1323:={seq(seq({[0,2],[-s1,1],[-s2,0]},s1 in TA[[2,3]]),s2 in TA[[1,3]])}:


gu:=gu12 union gu13 union gu23 union gu1213 union gu1223 union gu1323:

gu:={seq(CanTile(gu1) ,gu1 in gu)}:


mu:={ {[0,0]},{[0,1]},{[0,2]} }:

for gu1 in gu do

if IsLegalTile(gu1)  and IsTile(gu1,k,R)[1] then
  mu:=mu union {gu1}:
fi:
od:
RETURN(mu):

fi:



FAIL:



end:





#GLRseqBF(k,R,N): Same as GLRseq(k,R,N) but by brute force. For checking purposes only.
#3 by n rectangles, type
#GLRseqBF(3,{ [[1,2],{0}], [[1,3],{0}], [[2,3],{0}] },5); 
#To get the first 7 terms  of derangements of length n, type:
#GLRseqBF(2,{[[1,2],{0}],7);
#To get the first seven terms of the sequence enumeratingMenage permutations of length n, type:
#GLRseqBF(2,{[[1,2],{0,1}]});

GLRseqBF:=proc(k,R,N) local n1:
[seq(nops(GLR(n1,k,R)),n1=1..N)]:
end:




#GLRseq(k,R,N): Given a set of pairs,R {[[a,b],S]} where 1<=a<b<=k and S is a set of integers
#outputs the sequence enumerating reduced k by n matrices where every row is a permutation of {1,...,n}, the first row is [1,...,n] and obeying the condition
#M[a][i]<>M[b][i+s] for each s in S,  whenver both i and i+s are between 1 and n . For example, to get the first 5 terms of the sequence enumerating reduced Latin
#3 by n rectangles, type
#GLRseq(3,{ [[1,2],{0}], [[1,3],{0}], [[2,3],{0}] },5); 
#To get the first 7 terms  of derangements of length n, type:
#GLRseq(2,{[[1,2],{0}],7);
#To get the first seven terms of the sequence enumeratingMenage permutations of length n, type:
#GLRseq(2,{[[1,2],{0,1}]});

GLRseq:=proc(k,R,N) local TIL,z,i,x23,ti1,n0,gu,x1,x2,x3,mu1,mu,gu1,gu11,i1,pt,c:


if not type(k, integer) and k>=2 then
print(k, `must be an integer >=2`):
RETURN(FAIL):
fi:

if k>3 then
 print(`k>=4 are not yet implemented`):
 RETURN(FAIL):
fi:

TIL:=GatherTiles(k,R):

gu:=[seq(WtTilings(Rect(i,k),TIL,z),i=1..N)]:


if k=2 then

for ti1 in TIL do

 if ti1={[0,0]} then
  gu:=subs(z[ti1]=1,gu):

elif  ti1={[0,1]} then
  gu:=subs(z[ti1]=x2,gu):

else 
  gu:=subs(z[ti1]=-1,gu):
fi:

od:

fi:


if k=3 then


for ti1 in TIL do



if ti1={[0,0]} then
  gu:=subs(z[ti1]=1,gu):

elif ti1={[0,1]} then
  gu:=subs(z[ti1]=x2,gu):

elif ti1={[0,2]} then

  gu:=subs(z[ti1]=x3,gu):

elif nops(ti1)=2 and member(0,{seq(pt[2],pt in ti1)}) then

  gu:=subs(z[ti1]=-x1,gu):

elif nops(ti1)=2 and not member(0,{seq(pt[2],pt in ti1)}) then

  gu:=subs(z[ti1]=-x23,gu):

elif nops(ti1)=3 and IsTile(ti1,k,R)[2]=2 then

  gu:=subs(z[ti1]=x1,gu):

elif nops(ti1)=3 and IsTile(ti1,k,R)[2]=3 then

  gu:=subs(z[ti1]=2*x1,gu):


 fi:


od:


fi:




if k=2 then
mu:=[]:
for n0 from 1 to N do
gu1:=gu[n0]:
 mu1:=add(coeff(gu1,x2,i)*i!,i=0..degree(gu1,x2)):
mu:=[op(mu),mu1]:
od:


RETURN(mu):
fi:



if k=3 then

mu:=[]:



for n0 from 1 to N do

gu1:=expand(gu[n0]):


if type(gu1,`+`) then

mu1:=0:

for i1 from 1 to nops(gu1) do

gu11:=op(i1,gu1):

c:=normal(gu11/(x1^degree(gu11,x1)*x2^degree(gu11,x2)*x3^degree(gu11,x3)*x23^degree(gu11,x23))):

if not type(c, integer) then
 print(`Something bad happened`):
 RETURN(FAIL):
fi:

mu1:=mu1+c*binomial(n0-degree(gu11,x1),degree(gu11,x23))*degree(gu11,x23)!*degree(gu11,x2)!*degree(gu11,x3)!:
od:

mu:=[op(mu),mu1]:

else


c:=normal(gu1/(x1^degree(gu1,x1)*x2^degree(gu1,x2)*x3^degree(gu1,x3)*x23^degree(gu1,x23))):

if not type(c, integer) then
 print(`Something bad happened`):
 RETURN(FAIL):
fi:

mu1:=c*binomial(n0-degree(gu1,x1),degree(gu1,x23))*degree(gu1,x23)!*degree(gu1,x2)!*degree(gu1,x3)!:


mu:=[op(mu),mu1]:

fi:


od:

RETURN(mu):
fi:


end:





#Kernel2(R,x,X): Given a set of integers R, and variables x and X outputs the rational function f(x,X) such that the number of permutations pi of length n
#such that pi[i]-i is NOT in R is the coefficient of X^n in Int(exp(-x)*R(x,X),x=0..infinity), or equivalenty applyting the umbra x^n->n! to it. 
#To get the classical derangement numbers, type:
#Kernel2({0},x,X);
Kernel2:=proc(R,x,X) local Z,gu:
option remember:

gu:=EQS2(R,x,X,Z):


subs(solve(gu[1],convert(gu[2],set)),Z[1]):

end:



#GenDerSeq(R,N): Given a set of integers R (0 and negatives are OK) and a positive integer N, outputs the first N terms of the sequence enumerating
#the number of permutations of pi, {1,...,n} such that pi[i]-i is NOT in R. 
#For the first 30 terms of the number of derangements, type:
#GenDerSeq({0},30); 
#For the first 30 terms of the number of (straight) Menage numbers, type:
#GenDerSeq({0,1},30); 

GenDerSeq:=proc(R,N) local gu,x,X,lu,L,i,j:

 gu:=Kernel2(R,x,X):

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

L:=[]:

for i from 1 to N do
lu:=expand(coeff(gu,X,i)):
L:=[op(L),add(j!*coeff(lu,x,j),j=0..degree(lu,x))]:
od:
L:

end:








#Kernel3(R12,R13,R23,x1,x2,x3,x23,X): Given three sets of integers R12, R13, R23, and variables x2, x3, x23, and X outputs the rational function f(x2,x3,x23,X) such that 
#enables the computation of the number of 3 by n matrices
#where each row is a permutation of [1,...,n] and the first row is [1,...,n] (to get the unreduced number simply multiply by n!)
# M[1,i]<>M[2,j] whenever j-i is in R12
# M[1,i]<>M[3,j] whenever j-i is in R13
# M[2,i]<>M[3,j] whenever j-i is in R23
#via the umbral operation
#x1^a1*x2^a2*x3^a2*x23^a23 ->binomial(n-a1,a23)*a23!*a2!*a3!
#To get the Kernel for enumerating 3 by reduced Latin rectangles. Try:
#Kernel3({0},{0},{0},x1,x2,x3,x23,X);
Kernel3:=proc(R12,R13,R23,x1,x2,x3,x23,X) local Z,gu:
option remember:

gu:=EQS3(R12,R13,R23,x1,x2,x3,x23,X,Z):

normal(subs(solve(gu[1],convert(gu[2],set)),Z[1])):

end:





#GLR3seq(R12,R13,R23,N): Given  three sets of integers R12, R13,R23  and a positive integer N
#outputs the first N terms of the sequence enumerating 3 by n matrices
#where each row is a permutation of [1,...,n] and the first row is [1,...,n] (to get the unreduced number simply multiply by n!)
# M[1,i]<>M[2,j] whenever j-i is in R12
# M[1,i]<>M[3,j] whenever j-i is in R13
# M[2,i]<>M[3,j] whenever j-i is in R23
#It is the same output as
#GLRseq(3,{ [[1,2],R12], [[1,3],R13], [[2,3],R23]});
#but hopefeully faster.
#For example, to get the first 30 terms of the sequence enuerating reduced Latin Rectangles with 3 rows, try:
#GLR3seq({0},{0},{0},30);
GLR3seq:=proc(R12,R13,R23,N) local gu,x1,x2,x3,x23,X,gu1,gu11,mu,mu1,n0,c,i1:

 gu:=Kernel3(R12,R13,R23,x1,x2,x3,x23,X):

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

mu:=[]:


for n0 from 1 to N do

gu1:=expand(coeff(gu,X,n0)):

mu1:=0:

for i1 from 1 to nops(gu1) do

gu11:=op(i1,gu1):

c:=normal(gu11/(x1^degree(gu11,x1)*x2^degree(gu11,x2)*x3^degree(gu11,x3)*x23^degree(gu11,x23))):

if not type(c, integer) then
 print(`Something bad happened`):
 RETURN(FAIL):
fi:

mu1:=mu1+c*binomial(n0-degree(gu11,x1),degree(gu11,x23))*degree(gu11,x23)!*degree(gu11,x2)!*degree(gu11,x3)!:
od:

mu:=[op(mu),mu1]:
od:

mu:

end:

#Paper2(K,N): inputs positive integers K and N and outputs an article with the first N terms of the sequenes enumerating generalized derangelents
#i.e. permutations pi of {1,...n} such that pi[i]-i is not in S for all subsets S of {-K,...K}. try:
#Paper2(1,50);

Paper2:=proc(K,N) local gu,i,gu1,co,t0:

t0:=time():
gu:=powerset({seq(i,i=-K..K)}) minus {{}}:

co:=0:

print(`The first N terms of the sequences enumerating permutations that that pi[i]-i is never in the set S for all subets of`, {seq(i,i=-K..K)}):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

for gu1 in gu do
co:=co+1:
print(``):
print(`Prop. number`, co):
print(``):
print(`The first`, N, `terms, starting at n=1 of the number of permutations of {1,...,n} such that pi[i]-i is never in`, gu1, `are `):
print(``):
lprint(GenDerSeq(gu1,N)):
print(``):
print(`---------------------------------------------`):
print(``):

od:


print(`This concludes this paper that took`, time()-t0, `seconds to generate . `):

end:



#Paper3(K,N): inputs positive integers K and N and outputs an article with the first N terms of the sequenes enumerating 3-rowed matrices such that
#the bottom row is [1,...,n] and M[1,i]<>M[1,j], M[1,i]<>M[3,j], M[2,i]<>M[3,j] whenever j-i belongs to the set S
#for all subsets S of {-K,...K}. try:
#Paper3(1,30);

Paper3:=proc(K,N) local gu,i,gu1,co,t0:

t0:=time():
gu:=powerset({seq(i,i=-K..K)}) minus {{}}:

co:=0:

print(`The first N terms of the sequences enumerating  Rduced Generalized 3 by n Latin Rectangles with respect to all subests of`, {seq(i,i=-K..K)}):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

for gu1 in gu do
co:=co+1:
print(``):
print(`Prop. number`, co):
print(``):
print(`The first`, N, `terms, starting at n=1 of the number of 3-rowed matrices where each row is a permutation of {1....,n}, the first  row is [1,...n]`):
print(`and M[1,i]<>M[1,j], M[1,i]<>M[3,j], M[2,i]<>M[3,j] whenever j-i belongs to`, gu1, `are `):
print(``):
lprint(GLR3seq(gu1,gu1,gu1,N)):
print(``):
print(`---------------------------------------------`):
print(``):

od:


print(`This concludes this paper that took`, time()-t0, `seconds to generate . `):

end:





#EQS2(R,x,X,Z): Given a set of integers R, and variables x and X outputs the rational function f(x,X) such that the number of permutations pi of length n
#such that pi[i]-i is NOT in R is the coefficient of X^n in Int(exp(-x)*R(x,X),x=0..infinity), or equivalenty applyting the umbra x^n->n! to it. 
#To get the classical derangement numbers, type:
#EQS2({0},x,X,Z);
EQS2:=proc(R,x,X,Z) local TIL,eq,var,z,gu,ti1,Y,i:
option remember:
TIL:=GatherTiles(2,{[[1,2],R]}):

gu:=EQSVAR(2,TIL,z,X,Y):

eq:=gu[1]:

var:=gu[2]:

for ti1 in TIL do

 if ti1={[0,0]} then
  eq:=subs(z[ti1]=1,eq):

elif  ti1={[0,1]} then
  eq:=subs(z[ti1]=x,eq):

else 
  eq:=subs(z[ti1]=-1,eq):
fi:

od:

var:=var minus {Y[[{},0]]}:

var:=[Y[[{},0]],op(var)]:


eq:=subs({seq(var[i]=Z[i],i=1..nops(var))},eq):

eq,[seq(Z[i],i=1..nops(var))]:

end:




Chop1:=proc(f,x,N) local i:add(coeff(f,x,i)*x^i,i=0..N):end:
Chop:=proc(L,x,N) local i: [seq(Chop1(L[i],x,N),i=1..nops(L))]:end:


#GenDerSeqI(R,N): Same  as GenDerSeq(R,N) (q.v.) but using iterations
#For the first 30 terms of the number of (straight) Menage numbers, type:
#GenDerSeqI({0,1},30); 

GenDerSeqI:=proc(R,N) local gu,x,X,Z,eq,j,i,M,v0,v1,v2,L,lu:

gu:=EQS2(R,x,X,Z):

eq:=gu[1]:

gu:=[seq(subs(eq,Z[i]),i=1..nops(eq))]:


M:=[seq([seq(coeff(gu[i],Z[j],1),j=1..nops(gu))],i=1..nops(gu))]:


v0:=Chop([1,0$(nops(M)-1)],X,N):

v1:=Chop(expand([1+add(M[1][j]*v0[j],j=1..nops(M[1])),seq(add(M[i][j]*v0[j],j=1..nops(M[i])),i=2..nops(M))]),X,N):

while v0<>v1 do
v2:=Chop(expand([1+add(M[1][j]*v1[j],j=1..nops(M[1])),seq(add(M[i][j]*v1[j],j=1..nops(M[i])),i=2..nops(M))]),X,N):
v0:=v1:
v1:=v2:
od:

gu:=Chop(v2,X,N)[1]:


L:=[]:

for i from 1 to N do

lu:=coeff(gu,X,i):

L:=[op(L),add(j!*coeff(lu,x,j),j=0..degree(lu,x))]:
od:
L:

end:




#EQS3(R12,R13,R23,x1,x2,x3,x23,X,Z): Given three sets of integers R12,R13,R23,  and variables x1,x2,x3,x23 and X outputs the system of equations that enables
#the computation of  Kernel3(R12,R13,R23,x1,x2,x3,x23,X) and GLR3seq(R12,R13,R23,N);.  Try:
#EQS3({0},{0},{0},x1,x2,x3,x23,X,Z);
EQS3:=proc(R12,R13,R23,x1,x2,x3,x23,X,Z) local R,TIL,eq,var,z,gu,ti1,Y,i,pt:
option remember:

R:={[[1,2],R12], [[1,3],R13], [[2,3],R23] }:
TIL:=GatherTiles(3,R):

gu:=EQSVAR(3,TIL,z,X,Y):



eq:=gu[1]:

var:=gu[2]:

for ti1 in TIL do



if ti1={[0,0]} then
  eq:=subs(z[ti1]=1,eq):

elif ti1={[0,1]} then
  eq:=subs(z[ti1]=x2,eq):

elif ti1={[0,2]} then

  eq:=subs(z[ti1]=x3,eq):

elif nops(ti1)=2 and member(0,{seq(pt[2],pt in ti1)}) then

  eq:=subs(z[ti1]=-x1,eq):

elif nops(ti1)=2 and not member(0,{seq(pt[2],pt in ti1)}) then

  eq:=subs(z[ti1]=-x23,eq):

elif nops(ti1)=3 and IsTile(ti1,3,R)[2]=2 then

  eq:=subs(z[ti1]=x1,eq):

elif nops(ti1)=3 and IsTile(ti1,3,R)[2]=3 then

  eq:=subs(z[ti1]=2*x1,eq):


 fi:


od:

eq:=subs({seq(var[i]=Z[i],i=1..nops(var))},eq):


eq,[seq(Z[i],i=1..nops(var))]:
end:




#GLR3seqI(R12,R13,R23,N):  Like GLR3seq(R12,R13,R23,N)  but using iteration. Try:
#GLR3seqI({0},{0},{0},30);
GLR3seqI:=proc(R12,R13,R23,N) local gu,x1,x2,x3,x23,X,gu1,gu11,mu,mu1,n0,c,i1,Z, i, j, x, M, eq, v0, v1, v2:



gu:=EQS3(R12,R13,R23,x1,x2,x3,x23,X,Z):

eq:=gu[1]:

gu:=[seq(subs(eq,Z[i]),i=1..nops(eq))]:


M:=[seq([seq(coeff(gu[i],Z[j],1),j=1..nops(gu))],i=1..nops(gu))]:


v0:=Chop([1,0$(nops(M)-1)],X,N):

v1:=Chop(expand([1+add(M[1][j]*v0[j],j=1..nops(M[1])),seq(add(M[i][j]*v0[j],j=1..nops(M[i])),i=2..nops(M))]),X,N):

while v0<>v1 do
v2:=Chop(expand([1+add(M[1][j]*v1[j],j=1..nops(M[1])),seq(add(M[i][j]*v1[j],j=1..nops(M[i])),i=2..nops(M))]),X,N):
v0:=v1:
v1:=v2:
od:

gu:=Chop(v2,X,N)[1]:



mu:=[]:


for n0 from 1 to N do

gu1:=coeff(gu,X,n0):

mu1:=0:

for i1 from 1 to nops(gu1) do

gu11:=op(i1,gu1):

c:=normal(gu11/(x1^degree(gu11,x1)*x2^degree(gu11,x2)*x3^degree(gu11,x3)*x23^degree(gu11,x23))):

if not type(c, integer) then
 print(`Something bad happened`):
 RETURN(FAIL):
fi:

mu1:=mu1+c*binomial(n0-degree(gu11,x1),degree(gu11,x23))*degree(gu11,x23)!*degree(gu11,x2)!*degree(gu11,x3)!:
od:

mu:=[op(mu),mu1]:
od:

mu:

end: