############################################################################
## LatinTrapezoids.txt: Save this file as LatinTrapezoids.txt to use it,   #
# stay in the                                                              #
## same directory, get into Maple (by typing: maple <Enter> )              #
#then type: read `LatinTrapezoids.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: May 2021: tested for Maple 18  `):
print():
print(`This is LatinTrapezoids.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/LatinTrapezoids.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(`-------------------------------`):

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

 print(`-------------------------------`):



print(`For a list of the Symbolic Dynamical Porgramming functions type: ezraS();`):
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print():

 print(`-------------------------------`):



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

 print(`-------------------------------`):


ezraR:=proc()
if args=NULL then

print(` The Rectangles  procedures are: `):
print(`Bnz, IsGoodR1, LR, R3, R3G, R3bf , R3GSeq,  R3GSeqF, R3Seq, R3SeqF, SchR3, TilesRk, WtRec3, WtRec3G `):

else
 ezra(args):
fi:
end:


ezraS:=proc()
if args=NULL then

print(` The Symbolic Trapezoid  procedures are: `):
print(` EqSta, EvalSch, RegionToState , Sch, SolveSch, StateToRegion `):

else
 ezra(args):
fi:
end:

ezraT:=proc()
if args=NULL then

print(` The Trapezoid  procedures are: `):
print(`AntOld, Ant, IsGood1, LT, SchT3, T3, T3Seq, T3SeqBF , WtTrap3 `):

else
 ezra(args):
fi:
end:

ezra1:=proc()
if args=NULL then

print(` The supporting procedures are: `):
print(` AveAndMoms, AveAndMomsS, CP, ContsT, FindRandT, FindTiles, GR, IsLegalTile, PL, PrintLT, Rect, SeqFromRec, Tiles3, Trap , UV, Wt, WtTilingsBF`):

else
 ezra(args):
fi:
end:

ezra:=proc()

if args=NULL then

print(` LatinTrapezoids.txt: A Maple package for weighted counting of arbitary tilinds in general domains`):
print(`The MAIN procedures are: NuTilings,Tilings, WtTilings `):




elif nargs=1 and args[1]=almostLSgf then
print(`almostLSgf(n,k,x)`):
print(` Takes input positive integers n,k and variable name x `):
print(`  Outputs the generating function where the coefficient of `):
print(`  x^t gives the number of latin squares such that each row is a  `):
print(`  permutation of n, the bottom row is [1..n], and there are  `):
print(` exactly t instances of repeated numbers in columns (t violations). Try:`):
print(`almostLSgf(3,3,x);`):


elif nargs=1 and args[1]=Alpha then
print(`Alpha(f,x,N): Given a probability generating function f(x) (a polynomial, rational function and possibly beyond)`):
print(`returns a list, of length N, whose  (i) First entry is the average`):
print(`(ii): Second entry is the variance`):
print(`for i=3...N, the i-th entry is the so-called alpha-coefficients`):
print(`that is the i-th moment about the mean divided by the`):
print(`variance to the power i/2 (For example, i=3 is the skewness`):
print(`and i=4 is the Kurtosis)`):
print(`If f(1) is not 1, than it first normalizes it by dividing`):
print(`by f(1) (if it is not 0) .`):
print(`For example, try:`):
print(`Alpha(((1+x)/2)^100,x,4);`):

elif nargs=1 and args[1]=AntOld then
print(`AntOld(n,t): The weight-enumerator of the trapezoid Trap(n+2,3) with respect to the indexed variable z[tile] where the tiles are in Tiles3() (q.v.). Try:`):
print(`AntOld(3,t);`):

elif nargs=1 and args[1]=Ant then
print(`Ant(n0,x23,x2,x3,t1): The weight-enumerator of the trapezoid Trap(n+2,3) with respect to the active variables x23,x2,x3,t1 used in SchT3(z,x23,x2,x3,t1) (q.v.). Try:`):
print(`Ant(3,x23,x2,x3,t1);`):



elif nargs=1 and args[1]=AveAndMoms then
print(`AveAndMoms(f,x,N): Given a probability generating function`):
print(`f(x) (a polynomial, rational function and possibly beyond)`):
print(`returns a list whose first entry is the average  (alias expectation, alias mean)`):
print(`followed by the variance, the third moment (about the mean) and so on, until the N-th moment (about the mean).`):
print(`If f(1) is not 1, than it first normalizes it by dividing by f(1) (if it is not 0) .`):
print(`For example, try:`):
print(`AveAndMoms(((1+x)/2)^100,x,4);`):

elif nargs=1 and args[1]=AveAndMomsS then
print(`AveAndMomsS(L,x,K,n): inputs a list of polynomials L giving the weight-enumerator according to some random variable defined on a sequene of related sets`):
print(`conjectures explicit expressions, as rational functions in n for the expectation, varaiance, and the 3rd through the K-th moments about the mean . Try:`):
print(`AveAndMomsS(R3GSeqF(30,x), x,4,n); `):


elif nargs=1 and args[1]=Bnz then
print(`Bnz(n,z): The weight-enumerator of the Rectangle Rect(n,3) with respect to the indexed variable z[tile] where the tiles are in TilesRk(3) (q.v.). Try:`):
print(`Bnz(3,z);`):



elif nargs=1 and args[1]=Cnz then
print(`Cnz(n,z): The weight-enumerator of the Rectangle Rect(n,3) with respect to the indexed variable z[tile] where the tiles are in TilesSL(3) (q.v.). Try:`):
print(`Cnz(3,z);`):

elif nargs=1 and args[1]=ContsT then
print(`ContsT(old1): The SET OF legal continuations of the trapezoind old1`):
print(`Try:`):
print(`Conts([[1,2,3,4]]);`):

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]=EqSta then
print(`EqSta(SS,n,c,T,z,X,A): inputs an abstract state,SS, phrased in terms of the symbol n, (the width k is nops(SS)), and spacing c, with a list of tiles T, and symbol z (for indexed  variables z[tile]) and a variable X`):
print(`and a symbol A (such that A[SS] represents the weight-enumerator of the tilings of that state), ouputs the equation for A[SS] in terms of other A[state], followed by`):
print(`the set of "children states", i.e. all state such that A[state] shows up.`):
print(`outputs. Try:`):
print(`EqSta([[[], [0, n]], [[], [0, n]], [[], [0, n]]],n,1,TilesRk(3),z,X,A);`):
print(`EqSta([[[], [0, 2*n]], [[], [1, -1 + 2*n]], [[], [2, -2 + 2*n]]],n,2,Tiles3(),z,X,A);`):

elif nargs=1 and args[1]=EvalSch then
print(`EvalSch(Sc,z,i,n0): inputs a a scheme Sc, phrased in terms of the variable z, an integer i between 1 and nops(SC[1]) (alias nops(Sc[2])) and a positive integer n0`):
print(`outputs the value of the i-th component at n=n0. Try:`):
print(`Sc:=Sch([[[], [0, n]], [[], [0, n]], [[], [0, n]]],n,1,TilesRk(3),t,z)[1]; EvalSch(Sc,z,1,4);`):
print(`Sc:=Sch([[[], [0, 2*n]], [[], [1, -1 + 2*n]], [[], [2, -2 + 2*n]]],n,2,Tiles3(),t,z)[1];   EvalSch(Sc,z,1,2); `):

elif nargs=1 and args[1]=FindRandT then
print(`FindRandT(n,k,K): finds a random Latin trapezoid with k rows and base n, trying K times. Returns FAIl of nothing found. Try:`):
print(`FindRandT(5,5,100);`):

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]=GR then
print(`GR(L,n): inputs set  {[pt,value]} tries to find a rational function . try:`):
print(`GR({seq([i,(i+1)/(3*i+4)],i=1..10)},n);`):


elif nargs=1 and args[1]=IsGood1 then
print(`IsGood1(T): Given a potential trapezoid (in the right-angled format)`):
print(` that you know that its horiz. ,  vertical`):
print(`and its NW->SE diagonals for the truncated version (where you chopped`):
print(`the last row), returns true iff the new one is also OK`):
print(`Try:`):
print(`IsGood1([[1,2,3],[2,1]]);`):


elif nargs=1 and args[1]=IsGoodR1 then
print(`IsGoodR1(T): Given a potential rectangle`):
print(` that you know that except for the last row is Latin`):
print(` returns true iff the new one is also OK`):
print(`Try:`):
print(`IsGoodR1([[1,2,3],[2,3,1]]);`):

elif nargs=1 and args[1]=IsLegalTile then
print(`IsLegalTile(S): inputs a set of lattice points in the plane and decides whether it is a legal tile.`):
print(`(i) The intersection with every vertical line y=i has at most one element`):
print(`(ii): The smallest x-coordinate is 0. Try:`):
print(`IsLegalTile({[0,1],[1,0]});`):

elif nargs=1 and args[1]=LR then
print(`LR(n,k): The SET OF REDUCED k by n REDUCED LATIN  RECTANGLES. Try:`):
print(` LR(4,3);`):

elif nargs=1 and args[1]=LT then
print(`LT(n,k): The SET OF REDUCED n by k LATIN TRAPEZOIDS`):
print(`all different. Try`):
print(` LT(4,3);`):

elif nargs=1 and args[1]=NuTilings then
print(`NuTilings(S,T): Given a set of lattice points S, and a set of tiles, outputs the NUMBER of  tilings`):
print(`of S by the tiles of T. Try:`):
print(`NuTilings(Rect(2,2),{ {[0,0],[0,1]},{[0,0]},{[0,1]}});`):

elif nargs=1 and args[1]=PL then
print(`PL(S): plots the set of lattice points S. Try:`):
print(`PL(Rect(5,3));`):

elif nargs=1 and args[1]=PrintLT then
print(`PrintLT(T): Prints the latin trapezoid T. Try:`):
print(`PrintLT(LT(5,5)[1]);`):


elif nargs=1 and args[1]=R3 then
print(`R3(n): The number of 3 by n latin rectangles done "cleverly". Try:`):
print(`R3(4);`):

elif nargs=1 and args[1]=R3G then
print(`The weight-enumerator of 3 by n matrices whose bottom row is 1...n and the middle and top rows are permutions of [1,...,n] according to the weight x^(NumberOfViolations). Try:`):
print(`R3G(4,x);`):



elif nargs=1 and args[1]=R3GSeq then
print(`R3GSeq(N,x): The first N terms of the sequence weight-enumerator of 3 by n matrices where the bottom is [1, ..., n] and the middle and top rows are permutations of [1,...,n]`):
print(`according to the weight number of occurences of two distinct entries in the same column are equal. When x=0 it is the same as R3Seq(N), and when x=1 it is the sequence {n!^2}.`):
print(`Warning: for a much faster version use R3GSeqF(M,x)`):
print(` Try:`):
print(`R3GSeq(10,x);`):


elif nargs=1 and args[1]=R3GSeqF then
print(`R3GSeqF(N,x): The first N terms of the sequence "number of 3 by n latin rectangles" very fast, using the fifth-order recurrence gotten from R3GSeq(80,x). Try:`):
print(`R3GSeqF(30,x);`):

elif nargs=1 and args[1]=R3Seq then
print(`R3Seq(N): The first N terms of the sequence "number of 3 by n latin rectangles" done "cleverly", but directly.`):
print(`Warning: for a much faster version use R3SeqF(N)`):
print(`Try: `):
print(`R3Seq(30);`):


elif nargs=1 and args[1]=R3SeqF then
print(`R3SeqF(N): The first N terms of the sequence "number of 3 by n latin rectangles" very fast, using the fifth-order recurrence gotten from R3Seq(50). Try:`):
print(`R3SeqF(30);`):


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]=RegionToState then
print(`RegionToState(S,k,c,n): Given a set of lattice points in the region 0<=y<=k-1 where the default-spacing in each horizontal line is c, and a symbol n`):
print(`codes it as a list`):
print(`of length k, where row i+1, consits of two lists [L1,L2] where L1 is the list of x-coordinates for the "sporadic" points at the start`):
print(`whose y-coordinate is i, and L2 is the pair`):
print(`[start1,end1] where you have a continuum [[start1,i],[start1+c,i], ..., [end1,i]], and end1 is SYMBOLIC, in terms of n. It also returns the numeric n0 such that n=n0 gives that specifc region`):
print(` For example, try (it should give [[[],[0,18]],[[],[1,17]],[[],[2,16]]]):`):
print(`RegionToState(Trap(10,3),3,2,n);`):

elif nargs=1 and args[1]=R3bf then
print(`R3bf(n): The sequence of the number of 3 by n latin rectanles done "stupidly", by direct counting. ONLY FOR CHECMKING.. Try:`):
print(`R3bf(4);`):


elif nargs=1 and args[1]=Sch then
print(`Sch(SS,n,c,T,z,X): The scheme for computing the tilings g.f. for the state SS. Try:`):
print(`Sch([[[], [0, n]], [[], [0, n]], [[], [0, n]]],n,1,TilesRk(3),t,z);`):
print(`Sch([[[], [0, 2*n]], [[], [1, -1 + 2*n]], [[], [2, -2 + 2*n]]],n,2,Tiles3(),t,z);`):


elif nargs=1 and args[1]=SchR3 then
print(`SchR3(z,x23,x2,x3,t1): The  scheme needed to compute  R3(n) (q.v.), only using the active variables x23 (joining rows 2 and 3) t1 (all tiles touching the first row), x2 (singletons in 2nd row),`):
print(`x3 (singletons in 3rd row),. try:`):
print(`SchR3(z,x23,x2,x3,t1); `):

elif nargs=1 and args[1]=SchT3 then
print(`SchT3(z,x23,x2,x3,t1): The 16 by 16 scheme needed to compute  T3(n) (q.v.), only using the active variables x23 (joining rows 2 and 3) t1 (all tiles touching the first row), x2 (singletons in 2nd row),`):
print(`x3 (singletons in 3rd row). Try:`):
print(`SchT3(z,x23,x2,x3,t1); `):

elif nargs=1 and args[1]=SchSL3 then
print(`SchSL3(z,x23,x2,x3,t1): The  scheme needed to compute  SL3Seq(N) (q.v.), only using the active variables x23 (joining rows 2 and 3) t1 (all tiles touching the first row), x2 (singletons in 2nd row),`):
print(`x3 (singletons in 3rd row),. Try:`):
print(`SchSL3(z,x23,x2,x3,t1);`):
 
elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends it to the first K values. Try:`):
print(`SeqFromRec(N-n-1,n,N,[1],10);`):
print(`SeqFromRec(N^2-N-1,n,N,[1,1],10);`):


elif nargs=1 and args[1]=SL3 then
print(`SL3(n): The number of super Latin 3 by n retangles. Try:`):
print(`SL3(4);`):

elif nargs=1 and args[1]=SL3Seq then
print(`$L3Seq(N): The first N terms of the sequence enumerating super Latin 3 by n retangles. Try:`):
print(`SL3Seq(10);`):

elif nargs=1 and args[1]=SLR then
print(`SLR(n,k): The SET OF REDUCED n by k SUPER LATIN RETANGLES`):
print(`The bottom line is 1,...,n, `):
print(`all other rows are permutations of {1,...,n} and all columns have distinct entries. and also the NW and NE diagonals are different. Try`):
print(`SLR(4,3);`):

elif nargs=1 and args[1]=SolveSch then
print(`SolveSch(Sc,z): inputs a a scheme Sc, phrased in terms of the variable z, outputs the generating function, in z, for the first component. Try:`):
print(`Sc:=Sch([[[], [0, n]], [[], [0, n]], [[], [0, n]]],n,1,TilesRk(3),t,z)[1]; SolveSch(Sc,z);`):
print(` Sc:=SchT3(z,x23,x2,x3,t1): SolveSch(Sc,z); `):

elif nargs=1 and args[1]=StateToRegion then
print(`StateToRegion(SS,k,c,n,n0): Given an abstract state of width k and spacing c, and a concrete n0, finds the corresponding concrete region. Try:`):
print(`S:=Trap(10,3):evalb(S=StateToRegion(RegionToState(S,3,2,n)[1],3,2,n,9));`):
print(`S:=Rect(10,3):evalb(S=StateToRegion(RegionToState(S,3,1,n)[1],3,1,n,9));`):

elif nargs=1 and args[1]=T3 then
print(`T3(n): The number of latin trapezoids with n+2 cells at the bottom, n+1 in the middle, and n on the top, done "cleverly". Try:`):
print(`T3(2);`):


elif nargs=1 and args[1]=T3Seq then
print(`T3Seq(n): The sequence of the number of latin trapezoids with n+2 cells at the bottom, n+1 in the middle, and n on the top, done "cleverly", using inclusion-exclusion. Try:`):
print(`T3Seq(4);`):


elif nargs=1 and args[1]=T3SeqBF then
print(`T3SeqBF(n): The sequence of the number of latin trapezoids with n+2 cells at the bottom, n+1 in the middle, and n on the top, done "stupidly", by direct counting. ONLY FOR CHECMKING.. Try:`):
print(`T3SeqBF(4);`):

elif nargs=1 and args[1]=Tiles3 then
print(`Tiles3(): The set of tiles for the problem of counting Latin trapezoids with 3 rows. try:`):
print(`Tiles3();`):


elif nargs=1 and args[1]=TilesSL3 then
print(`TIlesSL3(): The tiles for a super-Latin 3 by n rectangle. Try:`):
print(`TilesSL3();`):


elif nargs=1 and args[1]=TilesRk then
print(`TilesRk(k): The set of tiles for the problem of counting Latin RECTANGLES with k rows. try:`):
print(`TilesRk(3);`):

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]=Trap then
print(`Trap(a,b): The trapezoid with a lattice points at the bottom, a-1, at the second row,...a-b+1 dots at the top (b-th) row such that the points interlace i.e.`):
print(`the points on y=i are [i,i],[i+2,i],..., [i+2*(a-i-1),i]`):
print(`Trap(6,4);`):

elif nargs=1 and args[1]=UV then
print(`UV(n,k): The Universal SET OF REDUCED matrices with the first  1,...,n, and the other rows permutations of {1,...n}`):
print(`UV(4,3);`):

elif nargs=1 and args[1]=Wt then
print(`Wt(T,z): The weight of the tiling T according to the indexed variables z[tile]. try:`):
print(`Wt({[[0,0],{[0,0],[1,0]}], [[0,1],{[0,0],[1,0]}]} ,z);`):


elif nargs=1 and args[1]=WtRec3 then
print(`WtRec3(M,n,z): The weight of a monomial in z[tiles] where tiles are in TilesRk(3)`):

elif nargs=1 and args[1]=WtRecSL3 then
print(`WtRecSL3(M,n,z): The weight of a monomial in z[tiles] where tiles are in TileSL3().`):
print(`Try: gu:=Cnz(2,t); WtRecSL3(op(1,gu),2,t);`):

elif nargs=1 and args[1]=WtRec3G then
print(`WtRec3G(M,n,z,x): The GENERALIZED weight, according to x^#Violations, of a monomial in z[tiles] where tiles are in TilesRk(3)`):

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),TilesRk(3),z);`):
print(`WtTilings(Trap(4,3),Tiles3(),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);`):

elif nargs=1 and args[1]=WtTrap3 then
print(`WtTrap3(M,n,z): The weight of a monomial in z[tiles] where tiles are in Tiles3()`):

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

fi:



end:




#START GuessRat
#GR1(L,n,d1,d2): inputs set  {[pt,value]} tries to find a rational function with degree d1 on top and d2 at bottom
GR1:=proc(L,n,d1,d2) local eq,var,a,b,T,B,i,f:

if nops(L)<d1+d2+5 then
 print(`The set of data points should be at least of length`, d1+d2+5):
RETURN(FAIL):
fi:

T:=add(a[i]*n^i,i=0..d1-1)+n^d1:
B:=add(b[i]*n^i,i=0..d2):
var:={seq(a[i],i=0..d1-1),seq(b[i],i=0..d2)}:
f:=T/B:
eq:={seq(numer(subs(n=L[i][1],f)-L[i][2]),i=1..nops(L))}:
var:=solve(eq,var):
if var=NULL then
 RETURN(FAIL):
fi:
factor(subs(var,f)):
end:


#GR(L,n): inputs set  {[pt,value]} tries to find a rational function . try:
#GR({seq([i,(i+1)/(3*i+4)],i=1..10)},n);
GR:=proc(L,n) local d,gu,d1:

 for d from 1 to nops(L)-6 do
  for d1 from 0 to d do
  gu:=GR1(L,n,d1,d-d1):
  if gu<>FAIL then
   RETURN(gu):
 fi:
 od:
od:
FAIL:
end:



#End GuessRat

#Start From HISTBRUT

#AveAndMomsSt(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list whose first entry is the average 
#(alias expectation, alias mean)
#followed by the moments
#so on, until the N-th moment (about the mean).
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#AveAndMomsSt(((1+x)/2)^100,x,4);

AveAndMomsSt:=proc(f,x,N) local mu,F,memu1,gu,i:
mu:=simplify(subs(x=1,f)):

if mu=0 then
print(f, `is neither a prob. generating function nor can it be made so`):
RETURN(FAIL):
fi:

F:=f/mu:


memu1:=simplify(subs(x=1,diff(F,x))):

gu:=[memu1]:


F:=x*diff(F,x):
for i from 2 to N do
 F:=x*diff(F,x):
 gu:=[op(gu),simplify(subs(x=1,F))]:
od:

gu:

end:

#AveAndMoms(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list whose first entry is the average 
#(alias expectation, alias mean)
#followed by the variance, the third moment (about the mean) and
#so on, until the N-th moment (about the mean).
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#AveAndMoms(((1+x)/2)^100,x,4);

AveAndMoms:=proc(f,x,N) local mu,F,memu1,gu,i:
mu:=simplify(subs(x=1,f)):

if mu=0 then
print(f, `is neither a prob. generating function nor can it be made so`):
RETURN(FAIL):
fi:

F:=f/mu:


memu1:=simplify(subs(x=1,diff(F,x))):

gu:=[memu1]:

F:=F/x^memu1:

F:=x*diff(F,x):
for i from 2 to N do
 F:=x*diff(F,x):
 gu:=[op(gu),simplify(subs(x=1,F))]:
od:

gu:

end:


#Alpha(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list, of length N, whose 
#(i) First entry is the average
#(ii): Second entry is the variance
#for i=3...N, the i-th entry is the so-called alpha-coefficients
#that is the i-th moment about the mean divided by the
#variance to the power i/2 (For example, i=3 is the skewness
#and i=4 is the Kurtosis)
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#Alpha(((1+x)/2)^100,x,4);
Alpha:=proc(f,x,N) local gu,i:
 gu:=AveAndMoms(f,x,N):
 if gu=FAIL then
  RETURN(gu):
 fi:

if gu[2]=0 then
print(`The variance is 0`):
  RETURN(FAIL):
fi:

[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..N)]:

end:


#End From HISTBRUT

###from FindRec.txt

Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:
 
#SeqFromRec(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values. Try:
#SeqFromRec(N-n-1,n,N,[1],10);
#SeqFromRec(N^2-N-1,n,N,[1,1],10);

SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), normal(-add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),j1=1..L))]:
od:

gu:

end:
###End rom FindRec.txt

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


#Trap(a,b): The trapezoid with a lattice points at the bottom, a-1, at the second row,...a-b+1 dots at the top (b-th) row such that the points interlace i.e.
#the points on y=i are [i,i],[i+2,i],..., [i+2*(a-i-1),i]
#Trap(6,4);
Trap:=proc(a,b) local j,k:
{seq(seq([j+2*k,j],k=0..a-1-j),j=0..b-1)}:
end:

#PL(S): plots the set of lattice points S. Try:
#PL(Rect(5,3));
PL:=proc(S)
plot(S,axes=none,style=point):
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:

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


#Wt(T,z): The weight of the tiling T according to the indexed variables x[tile]. try:
#Wt({[[0,0],{[0,0],[1,0]}], [[0,1],{[0,0],[1,0]}],x);
Wt:=proc(T,z) local i:
mul(z[T[i][1]],i=1..nops(T)):
end:



#WtTilingsBF(S,T,z): Given a set of lattice points S, and a set of tiles, outputs the weight-enumerator
#according to the index variables x[tile], by BRUTE FORCE. For checking purposes only. Should be the
#same as WtTiling(S,T,z). Try
#WtTilingsBF(Rect(2,2),{ {[0,0],[1,0]},{[0,0],[0,1]},{[0,0]}});
WtTilingsBF:=proc(S,T,z) local gu,i:
gu:=Tilings(S,T) :
add(Wt(gu[i],z),i=1..nops(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:




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





#Tiles3(): The set of tiles for the problem of counting Latin trapezoids with 3 rows. try:
#Tiles3():
Tiles3:=proc():
{ 
{[0,0]},{[0,1]},{[0,2]},
{[0,0],[1,1],[2,2]},
{[0,0],[1,1]},
{[0,0],[2,2]},
{[0,1],[1,2]},
{[0,2],[1,1],[2,0]},
{[0,2],[1,1]},
{[0,2],[2,0]},
{[0,1],[1,0]},
{[0,2],[2,0],[3,1]},
{[0,1],[1,0],[3,2]},
{[0,1],[1,2],[3,0]},
{[0,0],[2,2],[3,1]},
{[0,1],[1,0],[1,2]},
{[0,0],[0,2],[1,1]}
}:

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:

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

true:

end:

#WtTrap3(M,n,z): The weight of a monomial in z[tiles] where tiles are in Tiles3()
WtTrap3:=proc(M,n,z) local T,d,w,a2,a3,t,ava,T1:

T:=Tiles3():

T1:={}:

for t in T do
 if (member([0,0],t)  or member([1,0],t) or member([2,0],t) or member([3,0],t) ) and nops(t)>1 then
  T1:=T1 union {t}:
 fi:
od:


w:=normal(M/mul(z[t]^degree(M,z[t]),t in T)):


if not type(w,integer) then
 RETURN(FAIL):
fi:


#d is the number of tiles that join line second and third lines, i.e. the tiles {[0,1],[1,2]}, {[0,2],[1,1]}
d:=degree(M, z[{[0,1],[1,2]}] )+degree(M,   z[{[0,2],[1,1]}]  ):



#a2 is the number of cells left empty in the second row
a2:=degree(M,z[{[0,1]}]):


#a3 is the number of  cells left empty in the third (top) row
a3:=degree(M,z[{[0,2]}]):

ava:=n+2-add(degree(M,z[t]),t in T1):

w:=w*binomial(ava,d)*d!:

w:=w*(a2+1)!*(a3+2)!/2:


w:=w*(-1)^( degree(M,z[{[0,0], [1,1] }]) +degree(M,z[{[0,0], [2,2]}]) +degree(M,z[{[0,1] ,[1,2]}]) 
+ degree(M,z[{[0,2], [1,1]}])+degree(M,z[{[0,2],[2,0]}])+degree(M,z[{[0,1] ,[1,0]}]))*
2^(degree(M,z[{[0,0], [1,1] ,[2,2]}]) +  degree(M,z[{[0,2], [1,1] ,[2,0]}]) ):

w:

end:


#AntOldOld(n,t): The weight-enumerator of the trapezoid Trap(n+2,3) with respect to the indexed variable z[tile] where the tiles are in Tiles3() (q.v.). Try:
#AntOldOld(3,t);
AntOldOld:=proc(n,t)
WtTilings(Trap(n+2,3),Tiles3(),t):
end:

#T3OldOld(n): The number of latin trapezoids with n+2 cells at the bottom, n+1 in the middle, and n on the top, done "cleverly". Try:
#T3OldOld(2);
T3OldOld:=proc(n) local t,gu,i:
gu:=AntOldOld(n,t):
add(WtTrap3(op(i,gu),n,t),i=1..nops(gu)):
end:


#AntOld(n0,t): The weight-enumerator of the trapezoid Trap(n+2,3) with respect to the indexed variable z[tile] where the tiles are in Tiles3() (q.v.). Try:
#AntOld(3,t);
AntOld:=proc(n0,t) local n,z,Sc;
Sc:=Sch([[[], [0, 2*n]], [[], [1, -1 + 2*n]], [[], [2, -2 + 2*n]]],n,2,Tiles3(),t,z)[1]:
EvalSch(Sc,z,1,n0+1):
end:


#T3Old(n): The number of latin trapezoids with n+2 cells at the bottom, n+1 in the middle, and n on the top, done "cleverly". Try:
#T3Old(2);
T3Old:=proc(n) local t,gu,i:
gu:=AntOld(n,t):
add(WtTrap3(op(i,gu),n,t),i=1..nops(gu)):
end:


T3SeqOldOld:=proc(N) local n,gu:
gu:=[]:
for n from 1 to N do
 gu:=[op(gu),T3OldOld(n)]:
print(gu):
od:

gu:
end:


T3SeqOld:=proc(N) local n,gu:
gu:=[]:
for n from 1 to N do
 gu:=[op(gu),T3Old(n)]:
print(gu):
od:

gu:
end:



#DIRECT CONSTRUCTION OF Latin Trapezoids
with(combinat):

#IsGood1(T): Given a potential trapezoid (in the right-angled format)
# that you know that its horiz. ,  vertical
#and its NW->SE diagonals for the truncated version (where you chopped
#the last row), returns true iff the new one is also OK
#Try:
#IsGood1([[1,2,3],[2,1]]);
IsGood1:=proc(T) local k,i,j:
k:=nops(T):
for i from 1 to nops(T[k]) do
#checking that all NW->SW diagonals are OK

 if member(T[k][i],{seq(T[k-j][i+j],j=1..k-1)}) then
   RETURN(false):
 fi:

#Checking that all vertical lines are distinct
if member(T[k][i],{seq(T[j][i],j=1..k-1)}) then
   RETURN(false):
 fi:
od:

true:

end:

#LT(n,k): The SET OF REDUCED n by k LATIN TRAPEZOIDS
#The bottom line is 1,...,n, and each floor above it
#of length n-1,n-2, .., n-k+1 has all DIFFERENT entries
#also VERTICALLY all different also A[i][i],A[i+1][i+1], ... are
#all different
LT:=proc(n,k) local S, Sold,PNF,old1,f,new1,i:
option remember:

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

Sold:=LT(n,k-1):


PNF:=permute(n,n-k+1):
S:={}:
for old1 in Sold do
  for f in PNF do
    new1:=[op(old1),f]:


   if IsGood1(new1) then
      S:=S union {new1}:
   fi:
  od:
od:
S:
end:

#PT(T): prints a trapezoid
PT:=proc(T) local i:

for i from 1 to nops(T) do
 print(op(T[i])):
od:

end:


#T3SeqBF(n): The sequence of the first n number of latin trapezoids with n+2 cells at the bottom, n+1 in the middle, and n on the top, done "stupidly", by direct counting. 
#ONLY FOR CHECKING!
#Try:
#T3SeqBF(1);
T3SeqBF:=proc(n)  local n1: [seq(nops(LT(n1+2,3)),n1=1..n)]: end:

#END DIRECT CONSTRUCTION OF Latin Trapezoids

###start Latin Rectangles



#IsGoodR1(T): Given a potential rectangle (in the right-angled format)
# that you know that except for the bottom row is is Latin
#returns true iff the new one is also OK
#Try:
#IsGoodR1([[1,2,3],[2,1,3]]);
IsGoodR1:=proc(T) local n,k,i,j:
k:=nops(T):
n:=nops(T[1]):

#Checking that all colums are distinct
for i from 1 to n do
if member(T[k][i],{seq(T[j][i],j=1..k-1)}) then
   RETURN(false):
 fi:
od:

true:

end:

#LR(n,k): The SET OF REDUCED n by k LATIN RETANGLES
#The bottom line is 1,...,n, 
#all other rows are permutations of {1,...,n} and all columns have distinct entries. Try
#LR(4,3);
LR:=proc(n,k) local S, Sold,PNF,old1,f,new1,i:
option remember:

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

Sold:=LR(n,k-1):


PNF:=permute(n):
S:={}:
for old1 in Sold do
  for f in PNF do
    new1:=[op(old1),f]:


   if IsGoodR1(new1) then
      S:=S union {new1}:
   fi:
  od:
od:
S:
end:


#R3bf(n); The first n terms of the sequence enumerating 3 by Latin rectanles, doe by BRUCE formce
R3bf:=proc(n)  local n1: [seq(nops(LR(n1,3)),n1=1..n)]: end:



#NuViolsR(M): inputs a matrix M and outputs the number of pairs such that M[i][j1]=M[i][j2]
NuViolsR:=proc(M) local i,j1,j2,co:
co:=0:

for i from 1 to nops(M[1]) do
for j1 from 1 to nops(M) do
 for j2 from j1+1 to nops(M) do
  if M[j1][i]=M[j2][i] then
   co:=co+1:
  fi:
od:
od:
od:
  
co:

end:

#R3Gbf(n,x); The first n terms of the sequence R3G(n,x) by brute force. Only for checking. Try:
#R3Gbf(4,x);
R3Gbf:=proc(n,x) local mu,S,i,j,i1:
mu:=permute(n):
S:={seq(seq([[seq(i1,i1=1..n)],mu[i],mu[j]],j=1..nops(mu)),i=1..nops(mu))}:
sort(add(x^NuViolsR(S[i]),i=1..nops(S))):
end:







#TilesRk(k): The set of tiles for the problem of counting Latin rectangles with k rows. try:
#TilesRk(3):
TilesRk:=proc(k) local S,i:
S:={seq([0,i],i=0..k-1)}:
powerset(S) minus {{}}:
end:

#WtRec3(M,n,z): The weight of a monomial in z[tiles] where tiles are in TilesRk(3)
WtRec3:=proc(M,n,z) local T,d,w,a2,a3,t:


T:=TilesRk(3):

w:=normal(M/mul(z[t]^degree(M,z[t]),t in T)):


if not type(w,integer) then
 RETURN(FAIL):
fi:

#d is the number of tiles that join line second and third rows, i.e. the tiles {[0,1],[1,2]}, {[0,2],[1,1]}
d:=degree(M, z[{[0,1],[0,2]}] ):
w:=w*binomial(n-degree(M,z[{[0,0],[0,2]}]) -degree(M,z[{[0,0],[0,1]}]) -degree(M,z[{[0,0],[0,1],[0,2]}]) ,d)*d!:



#a2 is the number of cells left empty in the second row
a2:=degree(M,z[{[0,1]}]):


#a3 is the number of  cells left empty in the third (top) row
a3:=degree(M,z[{[0,2]}]):

w:=w*a2!*a3!:


w:=w*(-1)^( degree(M,z[{[0,0],[0,1]}])+degree(M,z[{[0,0],[0,2]}])+degree(M,z[{[0,1],[0,2]}]) )*2^degree(M,z[{[0,0],[0,1],[0,2]}]):

w:

end:



#WtRec3G(M,n,z,x): The Generalized weight, according to x^(# of violations) of a monomial in z[tiles] where tiles are in TilesRk(3)
WtRec3G:=proc(M,n,z,x) local T,d,w,a2,a3,t:


T:=TilesRk(3):

w:=normal(M/mul(z[t]^degree(M,z[t]),t in T)):


if not type(w,integer) then
 RETURN(FAIL):
fi:

#d is the number of tiles that join line second and third rows, i.e. the tiles {[0,1],[1,2]}, {[0,2],[1,1]}
d:=degree(M, z[{[0,1],[0,2]}] ):
w:=w*binomial(n-degree(M,z[{[0,0],[0,2]}]) -degree(M,z[{[0,0],[0,1]}]) -degree(M,z[{[0,0],[0,1],[0,2]}]) ,d)*d!:



#a2 is the number of cells left empty in the second row
a2:=degree(M,z[{[0,1]}]):


#a3 is the number of  cells left empty in the third (top) row
a3:=degree(M,z[{[0,2]}]):

w:=w*a2!*a3!:


w:=expand(w*(x-1)^( degree(M,z[{[0,0],[0,1]}])+degree(M,z[{[0,0],[0,2]}])+degree(M,z[{[0,1],[0,2]}]) )*((x-1)^2*(2+x))^degree(M,z[{[0,0],[0,1],[0,2]}])):

w:

end:


#Bnz(n,z): The weight-enumerator of the Rectangle Rect(n,3) with respect to the indexed variable z[tile] where the tiles are in TilesRk(3) (q.v.). Try:
#Bnz(3,z);
Bnz:=proc(n,z)
WtTilings(Rect(n,3),TilesRk(3),z):
end:



#R3(n): The number of 3 by n latin rectangles 
#R3(3);
R3:=proc(n) local z,gu,i:
gu:=Bnz(n,z):
add(WtRec3(op(i,gu),n,z),i=1..nops(gu)):
end:


#R3G(n,x): The weight-enumerator of 3 by n matrices whose bottom row is 1...n and the middle and top rows are permutions of [1,...,n] according to the weight x^(NumberOfViolations). Try:
#R3G(3,x);
R3G:=proc(n,x) local z,gu,i:
gu:=Bnz(n,z):
add(WtRec3G(op(i,gu),n,z,x),i=1..nops(gu)):
end:

#R3SeqOldOld(N): The first N terms, starting at n=1, of  "The number of 3 by n latin rectangles". Try:
#R3SeqOldOld(10);
R3SeqOldOld:=proc(N) local n:
[seq(R3(n),n=1..N)]:
end:



#R3SeqOld(N): The first N terms, starting at n=1, of  "The number of 3 by n latin rectangles", done directly (but cleverly). Try:
#R3SeqOld(10);
R3SeqOld:=proc(N) local Sc,n,z,t,gu,gu1,L,n0,i:
Sc:=Sch([[[], [0, n]], [[], [0, n]], [[], [0, n]]],n,1,TilesRk(3),t,z)[1]; 
gu:=SolveSch(Sc,z):
gu:=taylor(gu,z=0,N+2):

L:=[0]:

for n0 from 1 to N-1 do
 gu1:=expand(coeff(gu,z,n0)):
 L:=[op(L),add(WtRec3(op(i,gu1),n0+1,t ) , i=1..nops(gu1))]:
od:
L:

end:

#R3SeqF(N): The first N terms, starting at n=1, of  "The number of 3 by n latin rectangles", very fast, using the recurrence. try:
#R3SeqF(40);
R3SeqF:=proc(K) local n,N,ope,INI,i:
ope:=OpeR3(n,N,0):
INI:=[0, 0, 2, 24, 552]:
SeqFromRec(ope,n,N,INI,K):
end:


R3GSeq:=proc(N,x) local Sc,n,z,t,gu,gu1,L,n0,i:
Sc:=Sch([[[], [0, n]], [[], [0, n]], [[], [0, n]]],n,1,TilesRk(3),t,z)[1]; 
gu:=SolveSch(Sc,z):
gu:=taylor(gu,z=0,N+2):

L:=[x^3]:

for n0 from 1 to N-1 do
 gu1:=expand(coeff(gu,z,n0)):
 L:=[op(L),add(WtRec3G(op(i,gu1),n0+1,t,x ) , i=1..nops(gu1))]:
od:
L:

end:



#R3GSeqF(N,x): A fast version of R3GSeq(N,x) (q.v.) using the pre-computed recurrence, gotten from R3SeqG(80,x). Try:
#R3GSeqF(40,x);
R3GSeqF:=proc(K,x) local n,N,ope,INI:
ope:=OpeR3(n,N,x):
INI:=[x^3, x^6+3*x^2, x^9+9*x^5+6*x^3+18*x^2+2, x^12+18*x^8+24*x^6+72*x^5+45*x^4+176*x^3+144*x^2+72*x+24, x^15+30*x^11+60*x^9+180*x^8+225*x^7+860*x^6+972*x^5+2340*x^4+3720*x^3+3300*x^2+2160*x+552]:
SeqFromRec(ope,n,N,INI,K):
end:

OpeR3:=proc(n,N,x):
2*(x+2)*(x-1)^6*(n+4)*(3+n)*(n+2)*(n+1)*(x^2+2*n+x+8)/(x^2+2*n+x+6)-(x-1)^4*(n+4)*(3+n)*(n+2)*(-x^5-2*n*x^3+x^4+6*n*x^2-x^3+4*n^2+12*n*x+25*x^2+20*n+44*x+12)/(x^2+2*n+x+6)*N-(x-1)^3*(n+4)*(3+n)*(x^5+4*n*x^3+x^4+4*n^2*x+n*x^2+13*x^3-2*n^2+23*n*x-2*x^2-16*n+23*x-26)/(x^2+2*n
+x+6)*N^2+(x-1)*(n+4)*(2*n*x^4+x^5+5*n^2*x^2+3*n*x^3+10*x^4+2*n^3-n^2*x+39*n*x^2+8*x^3+22*n^2-14*n*x+72*x^2+82*n-33*x+102)/(x^2+2*n+x+6)*N^3-(x^5+n^2*x^2+3*n*x^3+x^4+2*n^3+3*n^2*x+10*n*x^2+10*x^3+24*n^2+19*n*x+24*x^2+98*n+28*x+136)/(x^2+2*n+x+6)*N^4+N^5:
end:

###end Latin Rectangles


#RegionToState(S,k,c,n): Given a set of lattice points in the region 0<=y<=k-1 where the default-spacing in each horizontal line is c, 
#and a symbol n
#codes it as a list
#of length k, where row i+1, consits of two lists [L1,L2] where L1 is the list of x-coordinates for the "sporadic" points at the start
#whose y-coordinate is i, and L2 is the pair
#[start1,end1], wnd1 is symbolic, where you have a continuum [[start1,i],[start1+c,i], ..., [end1,i]]. For example, try (it should give [[[],[0,18]],[[],[1,17]],[[],[2,16]]]):
#RegionState(Trap(10,3),3,2);

RegionToState:=proc(S,k,c,n) local T,gu,y1,pt,en1,st1,i,n0,k0,lu:


for y1 from 0 to k-1 do
 T[y1]:={}:
od:


for pt in S do
 T[pt[2]]:= T[pt[2]] union {pt[1]}:
od:

for y1 from 0 to k-1 do
T[y1]:=sort(convert(T[y1],list)):
od:

gu:=[]:

for y1 from 0 to k-1 do
if T[y1]=[] then
 gu:=[op(gu),[[],[]] ]:
fi:


en1:=T[y1][-1]:

if nops(T[y1])=1 then
 gu:=[op(gu),[[],[en1,en1]]]:
fi:

for i from nops(T[y1]) by -1 to 2 while T[y1][i]-T[y1][i-1]=c do od:

st1:=T[y1][i]:

gu:=[op(gu),[[op(1..i-1,T[y1])],[st1,en1]] ]:

 od:
gu:

k0:=gu[1][2][2] mod c:
n0:=(gu[1][2][2]-k0)/c:

lu:=[]:

for i from 1 to nops(gu) do
 en1:=k0+n*c+gu[i][2][2]-gu[1][2][2]:

lu:=[op(lu),[gu[i][1],[gu[i][2][1],en1]]]:
od:

lu,n0:


end:


#StateToRegion(SS,k,c,n,n0): Given an abstract state of width k and spacing c, and a concrete n0, finds the corresponding concrete region. Try:
#S:=Trap(10,3):evalb(S=StateToRegion(RegionToState(S,3,2,n),3,2,n,10));
#S:=Rect(10,3):evalb(S=StateToRegion(RegionToState(S,3,1,n),3,1,n,10));

StateToRegion:=proc(SS,k,c,n,n0) local S,i,mu,j,kama,i1:

S:={}:

for i from 0 to k-1 do
 mu:=SS[i+1]:
 S:=S union {seq([mu[1][i1],i],i1=1..nops(mu[1]) )}:
 kama:=(subs(n=n0,mu[2][2])-mu[2][1])/c+1:
 S:=S union {seq([mu[2][1]+(j-1)*c,i],j=1..kama)}:
od:

S:

end:



#EqSta(SS,n,c,T,z,X,A): inputs an abstract state,SS, phrased in terms of the symbol n, (the width k is nops(SS)), and spacing c, with a list of tiles T, and symbol z (for indexed  variables z[tile]) and a variable X
#and a symbol A (such that A[SS] represents the weight-enumerator of the tilings of that state), ouputs the equation for A[SS] in terms of other A[state], followed by
#the set of "children states", i.e. all state such that A[state] shows up.
#outputs. Try:
#EqSta([[[], [0, n]], [[], [0, n]], [[], [0, n]]],n,1,TilesRk(3),t,z);
#EqSta([[[], [0, 2*n]], [[], [1, -1 + 2*n]], [[], [2, -2 + 2*n]]],n,2,Tiles3(),t,z,A);
EqSta:=proc(SS,n,c,T,z,X,A) local k,eq,gu,t,L1,n0,n0a,pt1,S,S1,pt,T1,tS,j,TA,n0Min,smol,INI1:
#k is the width
k:=nops(SS):



#L1 is the maximal length of the tiles in T
L1:= max(seq(max(seq(pt[1],pt in t))-min(seq(pt[1],pt in t)),t in T)):

n0:=L1+2:
S:=StateToRegion(SS,k,c,n,n0):


pt:=CP(S):





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


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

eq:=0:
gu:={}:


for t in T1 do


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


if tS minus S={} then


for j from n0 by -1  to 0 while tS minus StateToRegion(SS,k,c,n,j)={} do od:

j:=max(0,j):


 S1:=S minus tS:


smol:=min(seq(pt[1],pt in S1)):
S1:={seq(pt-[smol,0],pt in S1)}:

S1:=RegionToState(S1,k,c,n):



 n0a:=S1[2]:
 S1:=S1[1]:


TA[S1]:=j+1:


 eq:=eq+A[S1]*z[t[1]]*X^(n0-n0a):
 gu:=gu union {S1}:
fi:

od:

n0Min:=max(seq(TA[S1], S1 in gu)):

INI1:=[seq(WtTilings(StateToRegion(SS,k,c,n,j), T,z),j=1..n0Min)]:

[eq,INI1],gu:


end:
                                                                                                               
                                                                                                                    


#Sch(SS,n,c,T,z,X): The scheme for computing the tilings g.f. for the state SS, followed by the legend.. Try:
#Sch([[[], [0, 2*n]], [[], [1, -1 + 2*n]], [[], [2, -2 + 2*n]]],n,2,Tiles3(),t,z);
Sch:=proc(SS,n,c,T,z,X) local A,gu,StillToDo, AlreadyDone,mu,SS1,STATES,INIS,i,EQS,i1,TA,B,j:

StillToDo:={SS}:

AlreadyDone:={}:

gu:=[]:

while StillToDo<>{} do
 SS1:=StillToDo[1]:
 mu:=EqSta(SS1,n,c,T,z,X,A):
 StillToDo:=StillToDo minus {SS1}:
 AlreadyDone:=AlreadyDone union {SS1}:
 StillToDo:= StillToDo union (mu[2] minus AlreadyDone):

  gu:=[op(gu),[SS1,mu[1]]]:
od:

gu:

STATES:=[seq(gu[i][1],i=1..nops(gu))]:
INIS:=[seq(gu[i][2][2],i=1..nops(gu))]:

STATES, INIS:

EQS:=[]:

for i from 1 to nops(gu) do
EQS:=[op(EQS),subs({seq(A[STATES[i1]]=B[i1],i1=1..nops(STATES))},gu[i][2][1])]:
od:

EQS, INIS:

for i from 1 to nops(STATES) do
 for j from 1 to nops(STATES) do
  TA[i,j]:=coeff(EQS[i],B[j]):
 od:
od:

[[seq([seq(TA[i,j],j=1..nops(STATES))],i=1..nops(STATES))],INIS],STATES:

end:





#EvalSch(Sc,z,i,n0): inputs a a scheme Sc, phrased in terms of the variable z, an integer i between 1 and nops(SC[1]) (alias nops(Sc[2])) and a positive integer n0
#outputs the value of the i-th component at n=n0. Try:
#Sc:=Sch([[[], [0, n]], [[], [0, n]], [[], [0, n]]],n,1,TilesRk(3),t,z)[1]; EvalSch(Sc,z,1,4);
EvalSch:=proc(Sc,z,i,n0) local su,j,mu,k:
option remember:

if n0<=0 then
 RETURN(0):
fi:

if n0<=nops(Sc[2][i]) then
 RETURN(Sc[2][i][n0]):
fi:

su:=0:

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

 if coeff(mu,z,0)<>0 then
   su:=expand(su+coeff(mu,z,0)*EvalSch(Sc,z,j,n0)):
 fi:

 for k from 1 to degree(mu,z) do
  su:=expand(su+coeff(mu,z,k)*EvalSch(Sc,z,j,n0-k)):
od:
od:
su:
end:


#SchT3(z,x23,x2,x3,t1): The 16 by 16 scheme needed to compute  T3(n) (q.v.), only using the active variables x23 (joining rows 2 and 3) t1 (all tiles touching the first row), x2 (singletons in 2nd row),
#x3 (singletons in 3rd row),
SchT3:=proc(z,x23,x2,x3,t1) local n,TIL,TIL1, Sc, a,t:
option remember:
TIL:=Tiles3():


TIL1:={}:

for a in TIL do
 if (member([0,0],a)  or member([1,0],a) or member([2,0],a) or member([3,0],a) ) and nops(a)>1 then
  TIL1:=TIL1 union {a}:
 fi:
od:



Sc:=Sch([[[], [0, 2*n]], [[], [1, -1 + 2*n]], [[], [2, -2 + 2*n]]],n,2,Tiles3(),t,z)[1]:


for  a in TIL do
 if nops(a)=2 then
  Sc:=subs(t[a]=-t[a],Sc):
 fi:

if (a={[0,0], [1,1] ,[2,2]} or a={[0,2], [1,1] ,[2,0]} ) then
   Sc:=subs(t[a]=2*t[a],Sc):
 fi:

od:


for  a in TIL1 do
 Sc:=subs(t[a]=t1,Sc):
od:

Sc:=subs({t[{[0,1],[1,2]}]=x23,t[{[0,2],[1,1]}]=x23},Sc):

Sc:=subs(t[{[0,1]}]=x2,Sc):
Sc:=subs(t[{[0,2]}]=x3,Sc):
Sc:=subs(t[{[0,0]}]=1,Sc):

Sc:

end:




#Ant(n0,x23,x2,x3,t1): The weight-enumerator of the trapezoid Trap(n+2,3) with respect to the active variables x23,x2,x3,t1 used in SchT3(z,x23,x2,x3,t1) (q.v.). Try
#Ant(3,x23,x2,x3,t1);
Ant:=proc(n0,x23,x2,x3,t1) local z,Sc;
Sc:=SchT3(z,x23,x2,x3,t1):
EvalSch(Sc,z,1,n0+1):
end:




#WtTrap3New(M,n,x23,x2,x3,t1 ): The weight of a monomial in z[tiles] where tiles are in Tiles3()
WtTrap3New:=proc(M,n,x23,x2,x3,t1 ) local T,d,w,a2,a3,t,ava,T1:



w:=normal(M/(x23^degree(M,x23)*x2^degree(M,x2)*x3^degree(M,x3)*t1^degree(M,t1))):


if not type(w,integer) then
 RETURN(FAIL):
fi:


#d is the number of tiles that join line second and third lines, i.e. the tiles {[0,1],[1,2]}, {[0,2],[1,1]}
d:=degree(M, x23 ):



#a2 is the number of cells left empty in the second row
a2:=degree(M,x2):


#a3 is the number of  cells left empty in the third (top) row
a3:=degree(M,x3):

ava:=n+2-degree(M,t1):

w:=w*binomial(ava,d)*d!:

w:=w*(a2+1)!*(a3+2)!/2:
w:

end:

#T3(n): The number of latin trapezoids with n+2 cells at the bottom, n+1 in the middle, and n on the top, done "cleverly". Try:
#T3(2);
T3:=proc(n) local gu,i, x23,x2,x3,t1:
gu:=Ant(n,x23,x2,x3,t1):

add(WtTrap3New(op(i,gu),n,x23,x2,x3,t1 ) , i=1..nops(gu)):
end:


T3SeqOld:=proc(N) local n,gu:
gu:=[]:
for n from 1 to N do
 gu:=[op(gu),T3(n)]:
od:

gu:
end:




#SolSch(Sc,z): inputs a a scheme Sc, phrased in terms of the variable z, outputs the generating function for the first (most important) component
#Try:
#Sc:=Sch([[[], [0, n]], [[], [0, n]], [[], [0, n]]],n,1,TilesRk(3),t,z)[1]; SolveSch(Sc,z);
SolveSch:=proc(Sc,z) local X,i1,eq,var,yemin,n0,j,START,i,lu,eq1:
option remember:

var:={seq(X[i1],i1=1..nops(Sc[1]))}:

n0:=max(seq(nops(Sc[2][i1]),i1=1..nops(Sc[2][1])))+3:

for i1 from 1 to nops(Sc[1]) do
START[i1]:=add(EvalSch(Sc,z,i1,j)*z^j,j=0..n0):
od:

for i from 1 to nops(Sc[1]) do
  lu:=expand(START[i]-add(Sc[1][i][j]*START[j],j=1..nops(Sc[1][i]))):
  lu:=add(coeff(lu,z,i1)*z^i1,i1=0..n0):
  yemin[i]:=lu:
od:


eq:={}:

for  i from 1 to nops(Sc[1]) do
 eq1:=X[i]-yemin[i]-add(Sc[1][i][j]*X[j],j=1..nops(Sc[1][i])):
 eq:=eq union {eq1}:
od:

var:=solve(eq,var):

normal(subs(var,X[1])):


end:






T3Seq:=proc(N) local Sc,gu,z,gu1,n0,i,L,t1,x23,x2,x3:


 Sc:=SchT3(z,x23,x2,x3,t1): 
 gu:=SolveSch(Sc,z):
gu:=taylor(gu,z=0,N+2):

L:=[]:

for n0 from 1 to N do
 gu1:=expand(coeff(gu,z,n0+1)):
 L:=[op(L),add(WtTrap3New(op(i,gu1),n0,x23,x2,x3,t1 ) , i=1..nops(gu1))]:
od:
L:

end:






#AveAndMomsS(L,x,K,n): inputs a list of polynomials L giving the weight-enumerator according to some random variable defined on a sequene of related seta
#conjectures explicit expressions, as rational functions in n for the expectation, varaiance, and the 3rd through the K-th moments about the mean . Try:
#AveAndMomsS(R3GSeqF(30,x) x,4,n):

AveAndMomsS:=proc(L,x,K,n) local mu,i,j,gu,L1,gu1:

if nops(L)<20 then
 print(`List must have at least 20 members`):
 RETURN(FAIL):
fi:

mu:=[seq(AveAndMoms(L[i],x,K),i=1..nops(L))]:

gu:=[]:

for j from 1 to K do
 L1:=[seq([i,mu[i][j]],i=5..nops(mu))]:
 gu1:=GR(L1,n):
 if gu1=FAIL then
 print(`We only managed to get up to`, nops(gu), `moments `):
   RETURN(gu):
 else
 gu:=[op(gu),gu1]:
 fi:
od:

gu:

end:


#AveAndMomsStS(L,x,K,n): inputs a list of polynomials L giving the weight-enumerator according to some random variable defined on a sequene of related seta
#conjectures explicit expressions, as rational functions in n for the expectation, and 3rd through the K-th moments . Try:
#AveAndMomsStS(R3GSeqF(30,x) x,4,n):

AveAndMomsStS:=proc(L,x,K,n) local mu,i,j,gu,L1,gu1:

if nops(L)<20 then
 print(`List must have at least 20 members`):
 RETURN(FAIL):
fi:

mu:=[seq(AveAndMomsSt(L[i],x,K),i=1..nops(L))]:

gu:=[]:

for j from 1 to K do
 L1:=[seq([i,mu[i][j]],i=5..nops(mu))]:
 gu1:=GR(L1,n):
 if gu1=FAIL then
 print(`We only managed to get up to`, nops(gu), `moments `):
   RETURN(gu):
 else
 gu:=[op(gu),gu1]:
 fi:
od:

gu:

end:










#SchR3(z,x23,x2,x3,t1): The  scheme needed to compute  R3(n) (q.v.), only using the active variables x23 (joining rows 2 and 3) t1 (all tiles touching the first row), x2 (singletons in 2nd row),
#x3 (singletons in 3rd row),
SchR3:=proc(z,x23,x2,x3,t1) local n,TIL,TIL1, Sc, a,t:
option remember:
TIL:=TilesRk(3):

TIL1:={}:

for a in TIL do
 if member([0,0],a)  and nops(a)>1 then
  TIL1:=TIL1 union {a}:
 fi:
od:


Sc:=Sch([[[], [0, n]], [[], [0, n]], [[], [0, n]]],n,1,TilesRk(3),t,z)[1]; 

for  a in TIL do
 if nops(a)=2 then
  Sc:=subs(t[a]=-t[a],Sc):
 fi:

if a={[0,0], [0,1] ,[0,2]}  then
   Sc:=subs(t[a]=2*t[a],Sc):
 fi:

od:


for  a in TIL1 do
 Sc:=subs(t[a]=t1,Sc):
od:

Sc:=subs(t[{[0,1],[0,2]}]=x23,Sc):

Sc:=subs(t[{[0,1]}]=x2,Sc):
Sc:=subs(t[{[0,2]}]=x3,Sc):
Sc:=subs(t[{[0,0]}]=1,Sc):

Sc:

end:









#WtRec3New(M,n,x23,x2,x3,t1 ): The weight of a monomial in z[tiles] where tiles are in TilesR3()
WtRec3New:=proc(M,n,x23,x2,x3,t1 ) local T,d,w,a2,a3,t,ava,T1:



w:=normal(M/(x23^degree(M,x23)*x2^degree(M,x2)*x3^degree(M,x3)*t1^degree(M,t1))):


if not type(w,integer) then
 RETURN(FAIL):
fi:


#d is the number of tiles that join line second and third lines, i.e. the tiles {[0,1],[1,2]}, {[0,2],[1,1]}
d:=degree(M, x23 ):



#a2 is the number of cells left empty in the second row
a2:=degree(M,x2):


#a3 is the number of  cells left empty in the third (top) row
a3:=degree(M,x3):

ava:=n-degree(M,t1):

w:=w*binomial(ava,d)*d!:

w:=w*a2!*a3!:
w:

end:

R3Seq:=proc(N) local Sc,gu,z,gu1,n0,i,L,t1,x23,x2,x3:


 Sc:=SchR3(z,x23,x2,x3,t1): 

 gu:=z*SolveSch(Sc,z):
gu:=taylor(gu,z=0,N+1):

L:=[0]:

for n0 from 2 to N do
 gu1:=expand(coeff(gu,z,n0)):


 L:=[op(L),add(WtRec3New(op(i,gu1),n0,x23,x2,x3,t1 ) , i=1..nops(gu1))]:
od:
L:

end:





###start Super Latin Rectangles



#IsGoodSLR1(T): Given a potential rectangle (in the right-angled format)
# that you know that except for the bottom row is is Latin
#returns true iff the new one is also OK
#Try:
#IsGoodSR1([[1,2,3],[2,1,3]]);
IsGoodSLR1:=proc(T) local n,k,i,j,a:
k:=nops(T):
n:=nops(T[1]):

#Checking that all colums are distinct
for i from 1 to n do
if member(T[k][i],{seq(T[j][i],j=1..k-1)}) then
   RETURN(false):
 fi:
od:


#Checking that all y=x direction diagonals  are distinct

for i from k to n do
 if member(T[k][i],{seq(T[k-a][i-a],a=1..k-1)}) then
   RETURN(false):
 fi:
od:



for i from 1 to n-k+1 do
 if member(T[k][i],{seq(T[k-a][i+a],a=1..k-1)}) then
   RETURN(false):
 fi:
od:


true:

end:


#SLR(n,k): The SET OF REDUCED n by k SUPER LATIN RETANGLES
#The bottom line is 1,...,n, 
#all other rows are permutations of {1,...,n} and all columns have distinct entries. 
#and also the NE and NW diagonals are different. Try
#SLR(4,3);
SLR:=proc(n,k)  local ret,rec,ps,r,p,x,y,ok,i: 
option remember:
if k=1 then 
 RETURN({[[seq(i,i=1..n)]]}):
fi:
ret:={}:
rec:= SLR(n,k-1):
ps:=permute(n):

for r in rec do
 for p in ps do
  ok:=true:
  for x from 1 to n do
   for y from 1 to k-1 do
    if p[x] = r[y][x] then
     ok:=false:
     break:
    fi:
   od:
   if ok = false then break: fi:

   for y from max(1,k-x+1) to k-1 do
    if p[x] = r[y][x-(k-y)] then
     ok:=false:
     break:
    fi:
   od:
   if ok = false then break: fi:

   for y from max(1,k-(n-x)) to k-1 do
    if p[x] = r[y][x+(k-y)] then
     ok:=false:
     break:
    fi:
   od:
   if ok = false then break: fi:
   
  od:
  if ok=true then
   ret:=ret union {[op(r),p]}:
  fi:
 od:
od:

ret:
end:


#SLR3bf(n); The first n terms of the sequence enumerating 3 by n Super Latin rectanles, doe by BRUCE formce
SLR3bf:=proc(n)  local n1: [seq(nops(SLR(n1,3)),n1=1..n)]: end:



#NuViolsR(M): inputs a matrix M and outputs the number of pairs such that M[i][j1]=M[i][j2]
NuViolsR:=proc(M) local i,j1,j2,co:
co:=0:

for i from 1 to nops(M[1]) do
for j1 from 1 to nops(M) do
 for j2 from j1+1 to nops(M) do
  if M[j1][i]=M[j2][i] then
   co:=co+1:
  fi:
od:
od:
od:
  
co:

end:

#R3Gbf(n,x); The first n terms of the sequence R3G(n,x) by brute force. Only for checking. Try:
#R3Gbf(4,x);
R3Gbf:=proc(n,x) local mu,S,i,j,i1:
mu:=permute(n):
S:={seq(seq([[seq(i1,i1=1..n)],mu[i],mu[j]],j=1..nops(mu)),i=1..nops(mu))}:
sort(add(x^NuViolsR(S[i]),i=1..nops(S))):
end:



#TilesSL3(): The tiles for a super-Latin 3 by n rectangle. Try:
#TilesSL3();
TilesSL3:=proc() local gu:

#singleton tiles
gu:={{[0,0]},{[0,1]},{[0,2]}}:

#tiles from row 1 and row 2 only
gu:=gu union {  {[0,1],[1,0]}, {[0,0],[0,1]},{[0,0],[1,1]} }:

#tiles from row 2 and row 3 only
gu:=gu union {  {[0,2],[1,1]}, {[0,1],[0,2]},{[0,1],[1,2]} }:

#tiles from row 1 and row 3 only
gu:=gu union {  {[0,2],[2,0]}, {[0,0],[0,2]},{[0,0],[2,2]} }:

#tiles from all three rows (that contain three edges, hence weight 2)
gu:=gu union {{[0,2],[1,1],[2,0]},{[0,0],[0,1],[0,2]},{[0,0],[1,1],[2,2]}}:

#tiles from all three rows (hybrid vertical and diagonal)
gu:=gu union {
{[0,1],[0,2],[1,0]},
{[0,2],[1,0],[1,1]},
{[0,0],[0,1],[1,2]},
{[0,0],[2,1],[2,2]},
{[0,2],[2,0],[2,1]},
{[0,0],[0,1],[2,2]},
{[0,1],[0,2],[2,0]}
}:

#tiles from all three rows (hybrid both diagonals)
gu:=gu union {
{[0,2],[2,0],[3,1]}, {[0,1],[1,0],[3,2]},
{[0,1],[1,2],[3,0]},{[0,0],[2,2],[3,1]},
{[0,0],[1,1],[0,2]},{[0,1],[1,0],[1,2]},
{[0,0],[1,1],[1,2]}
}


end:





#SchSL3(z,x23,x2,x3,t1): The  scheme needed to compute  SL3Seq(N) (q.v.), only using the active variables x23 (joining rows 2 and 3) t1 (all tiles touching the first row), x2 (singletons in 2nd row),
#x3 (singletons in 3rd row),. TryL
#SchSL3(z,x23,x2,x3,t1):
SchSL3:=proc(z,x23,x2,x3,t1) local n,TIL,TIL1, TIL23, TIL3edges, Sc, a,t:
option remember:
TIL:=TilesSL3():


#TIL1 is the set of non-trivial tiles that touch the x-axis
TIL1:={}:

for a in TIL do
 if (member([0,0],a)  or member([1,0],a) or member([2,0],a) or member([3,0],a) ) and nops(a)>1 then
  TIL1:=TIL1 union {a}:
 fi:
od:

#TIL23 is the set of non-trivial tiles that touch the 2nd and 3rd rows but not the first
TIL23:={{[0,1],[0,2]},{[0,2],[1,1]},{[0,1],[1,2]}}:



#TIL3edges is the set of tiles with 3 violations (so they give weight 2)
TIL3edges:={
{[0,2],[1,1],[2,0]},
{[0,0],[0,1],[0,2]},
{[0,0],[1,1],[2,2]},
{[0,0],[1,1],[0,2]},
{[0,1],[1,0],[1,2]}
}:


Sc:=Sch([[[], [0, n]], [[], [0, n]], [[], [0, n]]],n,1,TIL,t,z)[1]:


for  a in TIL do
 if nops(a)=2 then
  Sc:=subs(t[a]=-t[a],Sc):
 fi:

if member(a,TIL3edges) then
   Sc:=subs(t[a]=2*t[a],Sc):
fi:

od:

for  a in TIL1 do
 Sc:=subs(t[a]=t1,Sc):
od:

Sc:=subs({seq(t[a]=x23,a in TIL23)},Sc):


Sc:=subs(t[{[0,1]}]=x2,Sc):
Sc:=subs(t[{[0,2]}]=x3,Sc):
Sc:=subs(t[{[0,0]}]=1,Sc):

Sc:

end:



#SL3(n): The number of super Latin 3 by n retangles. Try:
#SL3(4);
SL3new:=proc(n) local z,z23,x2,x3,t1,Sc,gu1,i:
Sc:=SchSL3(z,x23,x2,x3,t1):
gu1:=EvalSch(Sc,z,1,n):
add(WtRec3New(op(i,gu1),n,x23,x2,x3,t1 ),i=1..nops(gu1)):
end:



#SL3Seq(N): The first N terms of the sequence enumerating super Latin 3 by n retangles. Try:
#SL3Seq(10);
SL3Seq:=proc(N) local Sc,z,x23,x2,x3,t1,n0,L,gu,gu1,i:

Sc:=SchSL3(z,x23,x2,x3,t1):


L:=[0]:

for n0 from 2 to N do
 gu1:=EvalSch(Sc,z,1,n0):
 L:=[op(L),add(WtRec3New(op(i,gu1),n0,x23,x2,x3,t1 ) , i=1..nops(gu1))]:
od:
L:
end:

#SL3SeqF(N): The first N terms of the sequence enumerating super Latin 3 by n retangles. Try:
#SL3SeqF(10);
SL3SeqF:=proc(N) local Sc,z,x23,x2,x3,t1,n0,L,gu,gu1,i:

Sc:=SchSL3(z,x23,x2,x3,t1):

gu:=SolveSch(Sc,z):


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



L:=[0]:

for n0 from 2 to N do
 gu1:=expand(coeff(gu,z,n0)):
 L:=[op(L),add(WtRec3New(op(i,gu1),n0,x23,x2,x3,t1 ) , i=1..nops(gu1))]:
od:
L:

end:






#WtRecSL3(M,n,t): The weight of a monomial in z[tiles] where tiles are in TileSL3().
WtRecSL3:=proc(M,n,t) local TIL,a,M1,d,w,a2,a3,TIL1,TIL23, TIL3edges,t1,x2,x3,x23,ava:


TIL:=TilesSL3():



#TIL1 is the set of non-trivial tiles that touch the x-axis
TIL1:={}:

for a in TIL do
 if (member([0,0],a)  or member([1,0],a) or member([2,0],a) or member([3,0],a) ) and nops(a)>1 then
  TIL1:=TIL1 union {a}:
 fi:
od:


#TIL23 is the set of non-trivial tiles that touch the 2nd and 3rd rows but not the first
TIL23:={{[0,1],[0,2]},{[0,2],[1,1]},{[0,1],[1,2]}}:


#TIL3edges is the set of tiles with 3 violations (so they give weight 2)
TIL3edges:={
{[0,2],[1,1],[2,0]},
{[0,0],[0,1],[0,2]},
{[0,0],[1,1],[2,2]},
{[0,0],[1,1],[0,2]},
{[0,1],[1,0],[1,2]}
}:

M1:=M:

for  a in TIL do
 if nops(a)=2 then
  M1:=subs(t[a]=-t[a],M1):
 fi:

if member(a,TIL3edges) then
   M1:=subs(t[a]=2*t[a],M1):
fi:

od:

for  a in TIL1 do
 M1:=subs(t[a]=t1,M1):
od:

M1:=subs({seq(t[a]=x23,a in TIL23)},M1):


M1:=subs(t[{[0,1]}]=x2,M1):
M1:=subs(t[{[0,2]}]=x3,M1):
M1:=subs(t[{[0,0]}]=1,M1):

w:=normal(M1/(x2^degree(M1,x2)*x3^degree(M1,x3)*x23^degree(M1,x23)*t1^degree(M1,t1))):

if not type(w,integer) then
 print(w):
 RETURN(FAIL):
fi:

w:



ava:=n-degree(M1,t1):


d:=degree(M1,x23):
a2:=degree(M1,x2):
a3:=degree(M1,x3):

w:=w*binomial(ava,d)*d!*a2!*a3!:
w:

end:


#Cnz(n,z): The weight-enumerator of the Rectangle Rect(n,3) with respect to the indexed variable z[tile] where the tiles are in TilesSL3() (q.v.). Try:
#Cnz(3,z);
Cnz:=proc(n,z):
WtTilings(Rect(n,3),TilesSL3(),z):
end:


#SL3Future(n): The number of super Latin 3 by n retangles. Try:
#SL3(4);
SL3Future:=proc(n) local z,x23,x2,x3,t1,Sc,gu1,i:
Sc:=SchSL3(z,x23,x2,x3,t1):
gu1:=EvalSch(Sc,z,1,n):
add(WtRec3New(op(i,gu1),n,x23,x2,x3,t1 ),i=1..nops(gu1)):
end:



#SL3(n): The number of super Latin 3 by n retangles. Try:
#SL3(4);
SL3:=proc(n) local t,gu,i:
gu:=Cnz(n,t):
add(WtRecSL3(op(i,gu),n,t) ,i=1..nops(gu)):
end:



#Cnz(n,z): The weight-enumerator of the Rectangle Rect(n,3) with respect to the indexed variable z[tile] where the tiles are in TilesSL3() (q.v.). by brute force
#Cnz(3,z);
Cnz:=proc(n,z):
WtTilings(Rect(n,3),TilesSL3(),z):

end:


#CnzBF(n,z): The weight-enumerator of the Rectangle Rect(n,3) with respect to the indexed variable z[tile] where the tiles are in TilesSL3() (q.v.). by brute force
#CnzBF(3,z);
CnzBF:=proc(n,z) local gu:
gu:=WtTilingsBF(Rect(n,3),TilesSL3(),z):
gu:
end:





#UV(n,k): The Universal SET OF REDUCED matrices with the first  1,...,n, and the other rows permutations of {1,...n}
#UV(4,3);
UV:=proc(n,k) local S, Sold,PNF,old1,f,new1,i:
option remember:

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

Sold:=UV(n,k-1):


PNF:=permute(n):
S:={}:
for old1 in Sold do
  for f in PNF do
    new1:=[op(old1),f]:
      S:=S union {new1}:
 od:
od:
S:
end:


#IsVPSL3(pt1,pt2): it the pair {pt1,pt2} such that A[pt1]=A[pt2] intruduces a violation for 3 by n Super Latin Rectangles. Try IsVPSL3,1],[1,0]);
#IsVPSL3 ([0,1],[2,2]);
#Is
IsVPSL3:=proc(pt1,pt2) local F:
F:={[0,1],[0,-1],[0,2],[0,-2],[1,1],[-1,-1],[-1,1],[1,-1],[2,2],[-2,-2],[2,-2],[-2,2]}:
member(pt1-pt2,F):
end:



#NuGPSL3(S): Given a set S finds the number of good pairs
NuGPSL3:=proc(S) local co,i,j:
co:=0:
for i from 1 to nops(S) do
 for j from i+1 to nops(S) do
  if IsVPSL3(S[i],S[j]) then
    co:=co+1:
  fi:
 od:
od:

co:

end:



#CHT(S): Given a set of lattice points finds the translate where the smallest x coordinate is 0. 
CHT:=proc(S) local s,c:
c:=min(seq(s[1],s in S)):
{seq(s-[c,0],s in S)}:
end:

#TilesSL3N(): The tiles for a super-Latin 3 by n rectangle. 
#Returns a list where
#the first entry is the set of singleton tiles
#the second entry is the set of tiles with two points (and one violation)
#the third entry is the set of tiles with three points and two violation
#the fourth entry is the set of tiles with three points and three violation
#Try:
#TilesSL3N();
TilesSL3N:=proc() local S,i,j,k,gu2,gu32,gu33,P:
S:={seq(seq([i,j],i=0..3),j=0..2)}:

gu2:={}:

for i from 1 to nops(S) do
 for j from i+1 to nops(S) do

    if S[i][2]<>S[j][2] then
   P:={S[i],S[j]}:
    if NuGPSL3(P)=1 then
        
      gu2:=gu2 union {CHT(P)}:
   fi:
  fi:

 od:
od:

gu32:={}:
gu33:={}:


for i from 1 to nops(S) do
 for j from i+1 to nops(S) do
  for k from j+1 to nops(S) do

     if (S[i][2]<>S[j][2] and S[i][2]<>S[k][2] and S[j][2]<>S[k][2] ) then
   P:={S[i],S[j],S[k]}:

    if NuGPSL3(P)=2 then
      gu32:=gu32 union {CHT(P)}:
   fi:

    if NuGPSL3(P)=3 then
      gu33:=gu33 union {CHT(P)}:
   fi:
fi:


 od:
od:
od:


[{{[0,0]},{[0,1]},{[0,2]}},gu2,gu32,gu33]
  

end:


#PrintLT(T): Prints the latin trapezoid T. Try:
#PrintLT(LT(5,5)[1]);
PrintLT:=proc(T) local i:

for i from 1 to nops(T) do
 print(op(T[nops(T)-i+1])):
od:

end:




#RandPreT(n,k): A random pre-trpezoid with k rows and base n. Try
#TrandPreT(5,5);
RandPreT:=proc(n,k) local T,i,i1,c,pi:
T:=[[seq(i,i=1..n)]]:

for i from 2 to k do
pi:=randperm(n-i+1):
c:=convert(randcomb(n,n-i+1),list):
c:=[seq(c[pi[i1]],i1=1..nops(pi))]:
T:=[op(T),c]:
od:

T:

end:


#IsGood(T): Given a potential trapezoid (in the right-angled format) decides whether it is indeed Latin
#IsGood([[1,2,3],[2,1]]);
IsGood:=proc(T) local k,i,j,T1:
k:=nops(T):

for i from 2 to k do
T1:=[op(1..i,T)]:
if not IsGood1(T1) then
 RETURN(false):
fi:
od:

true:

end:


#FindRandT(n,k,K): finds a random Latin trapezoid with k rows and base n, trying K times. Returns FAIl of nothing found. Try:
#FindRandT(5,5,100);
FindRandT:=proc(n,k,K) local T,i:

for i from 1 to K do
T:=RandPreT(n,k):
if IsGood(T) then
 RETURN(T):
fi:
od:

FAIL:

end: