######################################################################
##SD.txt: Save this file as  SD.txt                                  #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read SD.txt<Enter>                                 #
##Then follow the instructions given there                           #
##                                                                   #
##Written by George Spahn and                                        #
#Doron Zeilberger, Rutgers University ,                              #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: Fall 2023/ Winter 2024
 
print(`Created:  Fall  2023/ Winter 2024`):
print(` This is  SD.txt `):
print(`It is the main Maple package that accompanies the article `):
print(`Enumerating Seating Arrangments that Obey Social Distancing `):
print(`by George Spahn and Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):



print(`---------------------------------------`):
 print(`For a list of the Pre-computed procedures type ezraPC();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Story procedures type ezraS();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):


print(`---------------------------------------`):
 print(`For a list of the R22 procedures type ezraR22();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):


print(`---------------------------------------`):
 print(`For a list of the Checking procedures type ezraC();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

#Digits:=50:
with(combinat):
with(LinearAlgebra):







ezraR22:=proc()

if args=NULL then
 print(` The two procedures for avoiding 11 in a line are: GFu, fn `):

 print(``):

else
ezra(args):
fi:

end:

ezraPC:=proc()

if args=NULL then
 print(` The Pre-computed procedures are: GFrPC, GFrDimerPC, GFrKingsPC `):

 print(``):

else
ezra(args):
fi:

end:

ezraRSA:=proc()

if args=NULL then
 print(` The RSA procedures are: CanBHg, RSA,  RSAcollect, RSAcollectG, RSAcollectGall, RSA1g, RSAg, RSAgGF, RSAsimu, RSAsimuG, TransPol, TransPol1, TransPolS `):

 print(``):

else
ezra(args):
fi:

end:

ezraC:=proc()

if args=NULL then
 print(` The Checking procedures are: CheckGFr, CheckGFrN, IsBad1, IsGood, IsMaxGood, MGmats, WtMGmats  `):

 print(``):

else
ezra(args):
fi:

end:


ezraS:=proc()

if args=NULL then
 print(` The story procedures are: KingStory5, Paper1RowA, Paper1RowB, PaperC, PaperCdimer, PaperCkings,  PaperCsp, `):

 print(``):

else
ezra(args):
fi:

end:




ezraOld:=proc()
 
if args=NULL then
 print(` The old procedures are:  GFrT, TMold,  WCSold `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()
 
if args=NULL then
 print(` The supporting procedures are:  Alpha, AsyAv, CanBH, GFrN, RSA1, SetToMat, TM3`):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:  DensStat, DensStatG, GF1rowA, GF1rowB, GFr, GFrPC, LimAvDen,  SeqsR, `):
 print(` `):



elif nargs=1 and args[1]=Alpha then
print(`Alpha(f,x,N): Given a probability generating function`):
print(`f(x) (a polynomial, rational function and possibly beyond)`):
print(`returns a list, of length N, whose `):
print(`(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]=AsyAv then
print(`AsyAv(f,x,z): Given a rational generating  function f of x and z signifying weight enumerators of some statistic, find the`):
print(`constant alpha such that the (average of the statistics)/n converges to it for items of size n. Try:`):
print(`AsyAv(1/(1-x*z-x^2),x,z);`):


elif nargs=1 and args[1]=CanBH then
print(`CanBH(pt,B) can you build a house on pt if there houses at B. Try:`):
print(`CanBH([1,1],{[3,1]});`):

elif nargs=1 and args[1]=CanBHg then
print(`CanBHg(pt,SP,B,r,s) can you build a house on pt if there are houses at B (a subset of [1,s]x[q,r] and still avoid the patterns in the set of patterns SP. Try`):
print(`CanBHg([1,1],{{[0,0],[0,1],[0,2],[1,1]}}, {[3,1]},5,5 );`):

elif nargs=1 and args[1]=CheckGFr then
print(`CheckGFr(r,S,C): checks GFr(r,S,z,x) up to the coeff. x^C .`):
print(`It should give you [0,0,..,0].`):
print(` Try:`):
print(`CheckGFr(3,{{[0,0],[0,1],[0,2],[1,1]}},4);`):
print(`So far it does not work for polyonimos (or set of them) with two rows. for example for avoiding a cross, it does not work. `):
print(``):
print(`CheckGFr(3,{{[1,0],[1,1],[1,2],[0,1],[2,1]}},4);`): 


elif nargs=1 and args[1]=CheckGFrN then
print(`CheckGFrN(r,S,C): checks GFr(r,S,z,x) up to the coeff. x^C .`):
print(`It should give you [0,0,..,0].`):
print(` Try:`):
print(`CheckGFrN(3,{{[0,0],[0,1],[0,2],[1,1]}},4,3);`):
print(`CheckGFrN(3,{{[1,0],[1,1],[1,2],[0,1],[2,1]}},4,3);`): 

elif nargs=1 and args[1]=DensStat then
print(`DensStat(r,s,K): the density statistics for maximal configurations of an r by s building lot, up to the K-th moment. It also returns the actual distribution. `):
print(`Try:`):
print(`DensStat(5,1000,6);`):


elif nargs=1 and args[1]=DensStatG then
print(`DensStatG(r,s,SP,K): the density statistics for maximal configurations of an r by s lot avoiding the patterns of SP, up to the K-th moment.`):
print(`Try:`):
print(`DensStatG(3,50, {{[0,0],[0,1],[0,2],[1,1]}},6);`):

elif nargs=1 and args[1]=GFrT then
print(`GFrT(r,z,x): The generating function in x for the number of maximal T-avoiding 0-1 matrices with r colums, according to z^Number of ones. The`):
print(`coefficient of x^c is the weight-enumerator of c by  r such matrices.  Using the original old approach much slower.`):
print(`Here only for old-time-sake, Same as`):
print(`GFr(r,{[0,0],[0,1],[0,2],[1,1]},z,x); Try:`):
print(`GFrT(2,z,x);`):


elif nargs=1 and args[1]=fn then
print(`fn(n,z): The prob. generating function for the number of ones in maximal 11-avoiding 0-1 vectors using RSA. Try:`):
print(`fn(100,z);`):

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

print(`GFr(r,S,z,x): `):

print(`Inputs a  set of patterns, S, and a positive integer r, and variables z and x,`):
print(`Outputs The generating function in x for the number of maximal S-avoiding 0-1 matrices with r columns, according to z^Number of ones.`):
print(`The coefficient of x^c is the weight-enumerator of maximal c by r matrices avoiding the configurations in S.`):
print(``):
print(`For example for the original T-avoiding maximal configurations, with 3 rows, type:`):
print(`GFr(3,{{[0,0],[0,1],[0,2],[1,1]}},z,x);`): 
print(`or by using a reverse T `):
print(`GFr(3,{{[1,0],[1,1],[1,2],[0,1]}},z,x);`): 
print(`These should be the same `):
print(``):
print(`For avoiding both of them type`):
print(``):
print(`GFr(3,{{[0,0],[0,1],[0,2],[1,1]},{[1,0],[1,1],[1,2],[0,1]}},z,x);`): 
print(``):
print(`For avoiding dimers `):
print(``):
print(`GFr(3,{{[0,0],[0,1]},{[0,0],[1,0]}},z,x);`): 
print(`For avoiding a cross, type: `):
print(``):
print(`GFr(3,{{[1,0],[1,1],[1,2],[0,1],[2,1]}},z,x);`): 



elif nargs=1 and args[1]=GF1rowA then
print(`GF1rowA(k,z,x): The bi-variate generating function whose coefficient of x^n is the weight-enumerator, according to z^NumberOfOccupiedSeats in maximal seating arrangements`):
print(`with n seats in one row, obeying the social-distancing requirement that you can't have a block of k continguous occupied seats, in other words the`):
print(`weight enumerator of 0-1 vectors of length n avoiding the factor 1...1 (1 repeated k times), and such that if you change any 0 to a 1 you would`):
print(`get a violation. Try:`):
print(`GF1rowA(3,z,x);`):


elif nargs=1 and args[1]=GF1rowB then
print(`GF1rowB(k,z,x): The bi-variate generating function whose coefficient of x^n is the weight-enumerator, according to z^NumberOfOccupiedSeats in maximal seating arrangements`):
print(`with n seats in one row, obeying the social-distancing requirement that any two people must have at least k empty seats between them. In other words, the`):
print(`weight enumerator of 0-1 vectors of length n avoiding the patterns 11, 1*1, 1**1, ... 1*^(k-1)1`):
print(`get a violation. Try:`):
print(`GF1rowB(3,z,x);`):


elif nargs=1 and args[1]=GF1rowBconj then
print(`GF1rowB(k,z,x): The conjectured expression for GB1rowB(k,z,x). Try:`):
print(`evalb([seq(GF1rowB(k1,z,x),k1=1..5)]=[seq(GF1rowBconj(k1,z,x),k1=1..5)]);`):


elif nargs=1 and args[1]=GFrN then
print(`GFrN(r,S,z,x,a): `):
print(`Like GFr(r,S,z,x) but with another parameter a to get the off-set right`):
print(`Inputs a  set of patterns, S, and a positive integer r, and variables z and x,`):
print(`Outputs The generating function in x for the number of maximal S-avoiding 0-1 matrices with r columns, according to z^Number of ones.`):
print(`The coefficient of x^c is the weight-enumerator of maximal c by r matrices avoiding the configurations in S.`):
print(``):
print(`For example for the original T-avoiding maximal configurations, with 3 rows, type:`):
print(`GFrN(3,{{[0,0],[0,1],[0,2],[1,1]}},z,x,3);`): 
print(`or by using a reverse T `):
print(`GFrN(3,{{[1,0],[1,1],[1,2],[0,1]}},z,x,3);`): 
print(`These should be the same `):
print(``):
print(`For avoiding both of them type`):
print(``):
print(`GFrN(3,{{[0,0],[0,1],[0,2],[1,1]},{[1,0],[1,1],[1,2],[0,1]}},z,x,3);`): 
print(``):
print(`For avoiding dimers `):
print(``):
print(`GFrN(3,{{[0,0],[0,1]},{[0,0],[1,0]}},z,x,3);`): 
print(`For avoiding a cross, type: `):
print(``):
print(`GFrN(3,{{[1,0],[1,1],[1,2],[0,1],[2,1]}},z,x,3);`): 

elif nargs=1 and args[1]=GFrDimerPC then
print(`GFrDimerPC(C,z,x): the precumputed generating function for the weight-enumearor of r by C (c from 1 to 5) 0-1 matrices that aviod dimers ( i.e. avoid two consectutive horizontal 1's and `):
print(`avoid two consecutive vertical 1s). Try:`):
print(`GFrDimerPC(5,z,x)`):

elif nargs=1 and args[1]=GFrKingsPC then
print(`GFrKingsPC(C,z,x): The pre-comuted generating function according to the weight x^r*z^NumberOfKings in all maximal ways of placing non-attacking kings on an r by C chess board.`):
print(`If C s between 1 and 5 it is very fast (pre-comuted), but if C>=6 is is very slow. Try:`):
print(`GFrKingsPC(5,z,x);`):

elif nargs=1 and args[1]=GFrPC then
print(`GFrPC(r,z,x): The Pre-computed, for r<=5, generating function in x for the number of maximal T-avoiding 0-1 matrices with r columns, according to z^Number of ones. The`):
print(`coefficient of x^c is the weight-enumerator of c by  r such matrices. r should be between 2 and 6. If r>6 then it returns the very slow GFr(r,{[0,0],[0,1],[0,2],[1,1]},z,x). Try:`):
print(`GFrPC(6,z,x);`):

elif nargs=1 and args[1]=GFu then
print(`GFu(n,z): The prob. generating function for max. 11 avoiding vectors of length n for the r.v. number of ones. Try: `):
print(`GFu(10,z);`):

elif nargs=1 and args[1]=IsBad1 then
print(`IsBad1(M,P): Given a pattern  P and a 0-1 matrix decides whether it contains a pattern of ones that is a translate of P. Try:`):
print(`IsBad1([[0,1,0],[1,1,1],[0,1,0]],{[0,1],[1,0],[1,1],[1,2],[2,1]});`):

elif nargs=1 and args[1]=IsGood then
print(`IsGood(M,S) : Is the matrix M good with repsect to the set of patterns S?`):
print(`Try:`):
print(`IsGood([[0,1,0],[1,1,1],[0,1,0]],{{[0,1],[1,0],[1,1],[1,2],[2,1]}});`):

elif nargs=1 and args[1]=IsMaxGood then
print(`IsMaxGood(M,S) : Is the matrix M maximally good with repsect to the set of patterns  S?`):
print(`Try:`):
print(`IsMaxGood([[0,1,0],[1,1,1],[0,0,0]],{{[0,1],[1,0],[1,1],[1,2],[2,1]}});`):


elif nargs=1 and args[1]=KingStory5 then
print(`KingStory5(n): The story for non-attacking kinds for a 5 times n board. Try:`):
print(`KingStory5(100);`):

elif nargs=1 and args[1]=LimAvDen then
print(`LimAvDen(C): The limit of the average density of maximal T-avoiding r by C 0-1 matrices as r goes to infinity. Try:`):
print(`LimAvDen(5);`):


elif nargs=1 and args[1]=LimAvDenG then
print(`LimAvDenG(C): The limit of the average density of maximal SP-avoiding r by C 0-1 matrices as r goes to infinity. Try:`):
print(`LimAvDenG(5,{{[0,0],[1,0],[2,0],[1,1]}});`):

elif nargs=1 and args[1]=MGmats then
print(`MGmats(r,c,S): all the r by c 0-1 matrices that  are maximally good with respect to the set of patterns S. Try:`):
print(`MGmats(3,3,{{[0,1],[1,0],[1,1],[1,2],[2,1]}});`):
print(`MGmats(3,3,{{[0,0],[0,1],[0,2],[1,1]}});`):


elif nargs=1 and args[1]=Paper1RowA then
print(`Paper1RowA(B,z,x): Inputs a positive integer B and outputs a paper with the bi-variate generating function of maximal sittings on one row`):
print(`with at most b consecutive seats,  for b from 2 to B.Try:`):
print(`Paper1RowA(5,z,x);`):

elif nargs=1 and args[1]=Paper1RowB then
print(`Paper1RowB(B,z,x): Inputs a positive integer B and outputs a paper with the bi-variate generating function of maximal sittings on one row`):
print(`with minimimal distance>=b for b from 1 to B. It also confirms the general conjecture (probably fairly easy to prove) Try:`):
print(`Paper1RowB(5,z,x);`):

elif nargs=1 and args[1]=PaperC then
print(`PaperC(C,z,x,K1,K2,K3,K4): Inputs a positive integer C, and variables  z and x, and positive intgegers K1,K2,K3,K4`):
print(`Outputs A paper about maximal T-avoiding 0-1 r by C matrices, for general r. x is the varibale corresponding to the number of rows, and z to the number of 1s. Try:`):
print(`PaperC(3,z,x,30,100,6,1000):`):

elif nargs=1 and args[1]=PaperCdimer then
print(`PaperCdimer(C,z,x,K1,K2,K3,K4): A paper about maximal Dimer-avoiding 0-1 r by C matrices, for general r. Try:`):
print(`PaperCdimer(4,z,x,30,100,4,10):`):


elif nargs=1 and args[1]=PaperCkings then
print(`PaperCkings(C,z,x,K1,K2,K3,K4): A paper about placing non-attacking kings  on a C by r chessboard for general r. Try:`):
print(`PaperCkings(4,z,x,30,100,4,100):`):


elif nargs=1 and args[1]=PaperCsp then
print(`PaperCsp(C,SP,z,x,K1,K2,K3,K4): Inputs a positive integer C, a set of patterns, SP, and variables  z and x, and positive intgegers K1,K2,K3,K4`):
print(`Outputs A paper about maximal SP-avoiding 0-1 r by C matrices, for general r. x is the varibale corresponding to the number of rows, and z to the number of 1s. Try:`):
print(`PaperCsp(3,{{[0,0],[0,1],[0,2],[1,1]}},z,x,30,100,6,1000):`):

elif nargs=1 and args[1]=RSA then
print(`RSA(r,s): all maximal r by s {0,1} matrices that are T-avoiding. Try:`):
print(`RSA(3,3);`):


elif nargs=1 and args[1]=RSA1 then
print(`RSA1(r,s,pi): random sequential adsorption of an r by s lots avoiding T. Try:`):
print(`RSA1(3,3,[1,2,3,4,5,6,7,8,9]);`):

elif nargs=1 and args[1]=RSA1g then
print(`RSA1g(r,s,SP,pi): random sequential adsorption of an r by s lot avoiding the patterms in SP, using the permutation pi. Try:`):
print(`RSA1g(3,3,{{[0,0],[0,1],[0,2],[1,1]}}, [1,2,3,4,5,6,7,8,9]);`):

elif nargs=1 and args[1]=RSAg then
print(`RSAg(SP,r,s): all maximal r by s {0,1} matrices that avoide the polyonlinoes in SP. Followed by the frequency table. Try`):
print(`RSAg({{[0,0],[0,1],[0,2],[1,1]}},3,3);`):

elif nargs=1 and args[1]=RSAgGF then
print(`RSAgGF(SP,r,s,z): the prob. generating function of r by s in z, of {0,1} matrices that avoide the polyonlinoes in SP. Try`):
print(`RSAgGF({{[0,0],[0,1],[0,2],[1,1]}},3,2,z);`):
print(`RSAgGF({{[0,0],[1,0]}},5,1,z);`):

elif nargs=1 and args[1]=RSAcollect then
print(`RSAcollect(r,s,K): collects many maximal r by s 0-1 matrices avoiding T using K random permutations.`):
print(`Also returns how many are missing. Try:`):
print(`RSAcollect(3,3,1000); `):




elif nargs=1 and args[1]=RSAcollectG then
print(`RSAcollectG(r,s,SP,K): collects many maximal r by s 0-1 matrices avoiding the patterns in SP,  using K random permutations.`):
print(`RSAcollectG(3,3,{{[0,0],[0,1],[0,2],[1,1]}},1000);`):

elif nargs=1 and args[1]=RSAcollectGall then
print(`RSAcollectGall(r,s,SP,K): collects many maximal r by s 0-1 matrices avoiding the patterns in SP,  using K random permutations.`):
print(`Also returns how many are missing. Try:`):
print(`RSAcollectGall(3,3,{{[0,0],[0,1],[0,2],[1,1]}},1000):`):

elif nargs=1 and args[1]=RSAsimu then
print(`RSAsimu(r,s,K1,K2): the analog of  DensStat(r,s,K1) but using RSA simulation with K2 random permutations. Try:`):
print(`RSAsimu(5,20,6,1000);`):

elif nargs=1 and args[1]=RSAsimuG then
print(`RSAsimuG(r,s,SP,K1,K2): the analog of  DensStatG(r,s,SP,K1) but using RSA simulation with K2 random permutations. Try:`):
print(`RSAsimuG(3,10,{{[0,0],[0,1],[0,2],[1,1]}},6,1000);`):

elif nargs=1 and args[1]=SeqsR then
print(`SeqsR(r,z,K): The first K terms of the weight-enumerator according to the number of ones of`):
print(`number of MAXIMAL T-avoiding r by n 0-1 matrices, foloowed by their total number,`):
print(`followed by the first K terms of the number of EFFICIENT ones. Try:`):
print(`followed by those almost efficient ones`):
print(`Note that for (r=3), i.e. 3 by c matrices, the total number is A347725, and the number of efficient ones is a(2*i)=3*i-2,  a(2*i-1)=1`):
print(`for (r=4), i.e. 4 by c matrices, the total number is A189735, and the number of efficient ones is A15518`):
print(`for (r=5), i.e. 5 by c matrices, none of them were in the OEIS (Dec. 14, 2023)`):
print(`Try:SeqsR(3,z,30);`):

elif nargs=1 and args[1]=SetToMat then
print(`SetToMat(B,r,s): converts the subset of [1,r]x[1,s] to a 0-1 matrix. Try:`):
print(`SetToMat({[1,1],[2,2]},2,2);`):

elif nargs=1 and args[1]=TMold then
print(`TMold(r,z): The transfer matrix for maximal T-avoiding 0-1 matrices with r colums, according to the weight z^NumberOfOnes. Try:`):
print(`TMold(r,z);`):


elif nargs=1 and args[1]=TM then
print(`TM(r,S, z): The transfer matrix for maximal S-avoiding 0-1 matrices with r colums, avoiding the configurations in S, a set of patterns, according to the weight z^NumberOfOnes. Try:`):
print(`TM(3,{ {[0,0],[0,1],[0,2],[1,1]}},z);`):


elif nargs=1 and args[1]=TM3 then
print(`TM3(r,polys,z)` ):
print(` is the rows, polys is a set of sets of matrix coordinates, z is a weight variable. Try: `):
print(` TM3(3,{{[0,0],[1,1],[1,0],[2,0]}},z); `): 

elif nargs=1 and args[1]=TransPol then
print(`TransPol(P,pt,r,s): inputs a pattern P  a lattice point pt=[i,j] (1<=i<=r,1<=j<=s)`):
print(`outputs the set of translated pattern if pt is identified with any of the members of P`):
print(`and it still in [1,r]x[1,s]. Try:`):
print(`Try:`):
print(`TransPol({[0,0],[0,1],[0,2],[1,1]},[3,5],6,7);`):

elif nargs=1 and args[1]=TransPol1 then
print(`TransPol1(P,pt1,pt,r,s): inputs a pattern P a member of  P, pt1, and  a lattice point pt=[i,j] (1<=i<=r,1<=j<=s)`):
print(`outputs the translated pattern if pt is identified with pt1 and it is contained in [1,r]x[1,s] (if it does not it returns FAIL)`):
print(`Try:`):
print(`TransPol1({[0,0],[0,1],[0,2],[1,1]},[1,1],[3,5],6,7);`):

elif nargs=1 and args[1]=TransPolS then
print(`TransPolS(SP,pt,r,s): inputs a set of pattern SP  a lattice point pt=[i,j] (1<=i<=r,1<=j<=s)`):
print(`outputs the set of translated pattern if pt is identified with any of the members of any member P of SP`):
print(`and it still in [1,r]x[1,s]. Try:`):
print(`TransPolS({{[0,0],[0,1],[0,2],[1,1]}},[3,5],6,7);`):

elif nargs=1 and args[1]=WCSold then
print(`WCSold(rows,col,z): The weight-enumerator according to z^number of 1's of rows time columns 0,1 matrices that avoid a T, i.e a submatrix of the form [[1,1,1],[*,1,*]] with`):
print(` a maximal number of ones (i.e. maximal configurations according to the paper "Packing density of COmbinatorial Settlement Planning Models", by`):
print(`Mate Puljiz, Stjepan Sebek, and Josop Zubrinic, Amer. Math. Monthly v. 130, No. 10 (Dec. 2023). Try:`):
print(`WCSold(3,3,z);`):


elif nargs=1 and args[1]=WtMGmats then
print(`WtMGmats(r,c,S,z): the weight-enumerator of all r by c 0-1 matrices that maximally good with respect to the set of patterns S. Try:`):
print(`WtMGmats(3,3,{{[0,0],[0,1],[0,2],[1,1]}},z);`):

else
print(`There is no ezra for`,args):
fi:
 
end:




nextcolumns := proc(c,r)
local L,T,c3,flag,i,j:
L := []:
T:=combinat:-cartprod([seq([seq(i,i=0..1)],j=1..r)]):
while not T[finished] do
c3:=T[nextvalue]() :
flag := false:


for i from 1 to r-2 do
        if {c[2][i] , c[2][i+1], c[2][i+2],c3[i+1]} = {1} then
                flag := true:
                break:
        fi:
od:
if flag then
        next:
fi:
for i from 1 to r do
        if c[2][i] = 1 then next: fi:
        if i > 1 and i < r and ({c[1][i] , c[1][i+1], c[1][i-1]} = {1}) then next: fi:

        if i > 1 and i < r and ({c[2][i-1] , c[2][i+1], c3[i]} = {1}) then next: fi:
        if i > 2 and ({c[2][i-1] , c[2][i-2], c3[i-1]} = {1}) then next: fi:
        if i < r-1 and ({c[2][i+2] , c[2][i+1], c3[i+1]} = {1}) then next: fi:
        flag := true:
        break:
od:
if  flag then next: fi:
L := [op(L),[c[2],c3]]
od:
L:
end:

gen_states := proc(r):
local L,T,c3,i,j:
L := []:
T:=combinat:-cartprod([seq([seq(i,i=0..1)],j=1..2*r)]):
while not T[finished] do
c3:=T[nextvalue]() :
L:=[op(L),[[seq(c3[i],i=1..r)],[seq(c3[i],i=r+1..2*r)]]]:
od:
L:
end:

#TMold(r,z): The transfer matrix for maximal T-avoiding 0-1 matrices with r colums, according to the weight z^NumberOfOnes. Try:
#TMold(r,z);
TMold:=proc(r,z)  local S,M,i,t,w:
option remember:
S := gen_states(r):
M:=Matrix(nops(S)+2,nops(S)+2):
for i from 1 to nops(S) do

        for t in nextcolumns(S[i],r) do
                #print(t,ListTools:-Search(t,S)):
                M[i,ListTools:-Search(t,S)] := weight(t,r):
        od:
        if final_state(S[i],r) then
                M[i,nops(S)+2] := 1:
        fi:
        if initial_state(S[i],r) then
                M[nops(S)+1,i] := weight(S[i],r):
        fi:
od:
M:
end:


initial_state:=proc(c,r)
if c[1] = [0$r] then return(true): fi:
return(false):
end:

final_state := proc(c,r)
local i:
if c[2][1] = 0 or c[2][r] = 0 then return(false): fi:
for i from 2 to r-1 do
if c[2][i] = 0 and (0 in  {c[1][i] , c[1][i+1], c[1][i-1]}) then
return(false)
fi:
od:
return(true):
end:



weight := proc(c,r)

z^ListTools:-Occurrences(1,c[2]):
end:

#WCSold(rows,col,z): The weight-enumerator according to z^number of 1's of rows time columns 0,1 matrices that avoid a T, i.e a submatrix of the form [[1,1,1],[*,1,*]] with
# a maximal number of ones (i.e. maximal configurations according to the paper "Packing density of COmbinatorial Settlement Planning Models", by
#Mate Puljiz, Stjepan Sebek, and Josop Zubrinic, Amer. Math. Monthly v. 130, No. 10 (Dec. 2023). Try:
#WCSold(3,3,z);
WCSold:=proc(rows, cols,z)
local p:
p := (TMold(cols,z)^(rows+1))[2^(2*cols)+1][2^(2*cols)+2]:
expand(p):
end:


#GFrT(r,z,x): The generating function in x for the number of maximal T-avoiding 0-1 matrices with r rows, according to z^Number of ones. The
#coefficient of x^c is the weight-enumerator of r by  c such matrices. Same as WCSold(r,c,z); Try:
#GFrT(2,z,x);
GFrT:=proc(r,z,x) local A,d:
A:=TMold(r,z):
d:=Dimension(A)[1]:
normal(evalm(((IdentityMatrix(d)-x*A))^(-1))[d-1,d]/x):
end:


#GFrPC(r,z,x): The pre-computed generating function in x for the number of maximal T-avoiding 0-1 matrices with r rows, according to z^Number of ones. The
#coefficient of x^c is the weight-enumerator of r by  c such matrices. Same as WCSold(r,c,z); r should be <=6. Try:
#GFrPC(3,z,x);
GFrPC:=proc(r,z,x) local L:


L:=
[-1/(x*z-1), -1/(x*z^2-1), -(x^2*z^5+x*z^3-2*x*z^2+1)/(x^3*z^7+x^2*z^5+2*x*z^2-\
1), -(2*x^2*z^5+x*z^4-2*x*z^3-x*z^2+1)/(3*x^2*z^6-2*x^2*z^5+2*x*z^3+x*z^2-1), -
(2*x^4*z^15-2*x^4*z^14-x^3*z^12+2*x^3*z^11-x^3*z^10+x^2*z^9-5*x^2*z^8+4*x^2*z^7
-x*z^5+2*x*z^4+x*z^3-1)/(8*x^5*z^17-2*x^4*z^14+4*x^4*z^13+3*x^3*z^11-7*x^3*z^10
+x^2*z^8-8*x^2*z^7-2*x*z^4-x*z^3+1), (2*x^6*z^25+6*x^6*z^24-2*x^5*z^23+6*x^5*z^
22+4*x^5*z^21-5*x^5*z^20-4*x^4*z^19+21*x^4*z^18-20*x^4*z^17-7*x^4*z^16-2*x^3*z^
15+19*x^3*z^14-28*x^3*z^13+12*x^3*z^12+3*x^2*z^10-2*x^2*z^9-2*x^2*z^8+x*z^6-5*x
*z^4+1)/(10*x^6*z^25+14*x^6*z^24+12*x^5*z^21-17*x^5*z^20-12*x^4*z^17-35*x^4*z^
16-20*x^3*z^13+8*x^3*z^12-6*x^2*z^9-6*x^2*z^8-5*x*z^4+1)]:
if r>6 then
GFr(r,z,x):
else
L[r]:
fi:
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:



#DensStat(r,s,K): the density statistics for maximal configurations of an r by s lot, up to the K-th moment.
#Try:
#DensStat(5,1000,6);
DensStat:=proc(r,s,K) local f,z,x,lu,Dist,i:
f:=GFrPC(r,z,x):
f:=coeff(taylor(f,x=0,s+1),x,s):
f:=f/subs(z=1,f):

Dist:=[]:

for i from ldegree(f,z) to degree(f,z) do
 Dist:=[op(Dist),[i/(r*s),evalf(coeff(f,z,i))]]:
od:

lu:=evalf(Alpha(f,z,K),20):
lu:=[lu[1]/(r*s),lu[2]/(r*s),op(3..K,lu)]:
lu,Dist:

end:



#CanBH(pt,B) can you build a house on pt if there are houses at B. Try:
#CanBH([1,1],{[3,1]});
CanBH:=proc(pt,B)

 if  {pt+[0,1],pt+[0,2],pt+[1,1]} minus B={} then
  RETURN(false):
 fi:

 if  {pt+[0,-1],pt+[0,1],pt+[1,0]} minus B={} then
  RETURN(false):
 fi:

 if  {pt+[0,-2],pt+[0,-1],pt+[1,-1]} minus B={} then
  RETURN(false):
 fi:


 if  {pt+[-1,-1],pt+[-1,0],pt+[-1,1]} minus B={} then
  RETURN(false):
 fi:

true:
end:

#RSA1(r,s,pi): random sequential adsorption of an r by s lots avoiding T. Try:
#RSA1(3,3,[1,2,3,4,5,6,7,8,9]);
RSA1:=proc(r,s,pi) local L,i,j,B:
if nops(pi)<>r*s then
 RETURN(FAIL):
fi:

L:=[ seq(seq([i,j],j=1..s),i=1..r)]:


B:={}:

for i from 1 to nops(L) do
 if CanBH(L[pi[i]],B) then
  B:=B union {L[pi[i]]}:
 fi:
od:
B:
end:

#SetToMat(B,r,s): converts the subset of [1,r]x[1,s] to a 0-1 matrix. Try:
#SetToMat({[1,1],[2,2]},2,2);
SetToMat:=proc(B,r,s) local i,j,b,T:


for i from 1 to r do
 for j from 1 to s do
  T[i,j]:=0:
 od:
od:

for b in B do
 T[op(b)]:=1:
od:


[seq([seq(T[i,j],j=1..s)],i=1..r)]:

end:



#RSA(r,s): all maximal r by s {0,1} matrices that are T-avoiding. Try:
#RSA(3,3);
RSA:=proc(r,s) local lu,gu,T,pi,M,i:
lu:=permute(r*s):

gu:={}:

for i from 1 to nops(lu) do
 M:=SetToMat(RSA1(r,s,lu[i]),r,s):
 if not member(M,gu) then
   gu:=gu union {M}:
   T[M]:=1:
  else
  T[M]:=T[M]+1:
 fi:
od:

[gu,op(T)]:

end:



#RSAcollect(r,s,K): collects many maximal r by s 0-1 matrices avoiding T using K random permutations.
#Also returns how many are missing. Try:
#RSAcollect(3,3,1000):
RSAcollect:=proc(r,s,K) local gu,f,co,gu1,i:
gu:={seq(RSA1(r,s,randperm(r*s)),i=1..K)}:

f:=subs(z=1,GFrPC(s,z,x)):
co:=coeff(taylor(f,x=0,r+1),x,r):

gu:={seq(SetToMat(gu1,r,s),gu1 in gu)}:
[gu,co-nops(gu)]:

end:

#RSAsimu(r,s,K1,K2): the analog of  DensStat(r,s,K1) but using RSA simulation with K2 random permutations. Try:
#RSAsimu(5,20,6,1000);
RSAsimu:=proc(r,s,K1,K2) local gu,z,f,i,Dist,lu:
f:=add(z^nops(RSA1(r,s,randperm(r*s))),i=1..K2):
f:=f/subs(z=1,f):

Dist:=[]:
for i from ldegree(f,z) to degree(f,z) do
 Dist:=[op(Dist),[i/(r*s),evalf(coeff(f,z,i))]]:
od:

lu:=evalf(Alpha(f,z,K1),20):
lu:=[lu[1]/(r*s),lu[2]/(r*s),op(3..K1,lu)]:
lu,Dist:
end:


#SeqsR(r,z,K): The first K terms of the weight-enumerator according to the number of ones of
#number of MAXIMAL T-avoiding r by n 0-1 matrices, foloowed by their total number,
#followed by the first K terms of the number of EFFICIENT ones. 
#followed by those almost efficient ones
#Try:
#SeqsR(3,z,30);
SeqsR:=proc(r,z,K) local f,f1,gu,gu1,gu2,gu3,i:
f:=GFrPC(r,z,x):
f1:=taylor(f,x=0,K+1):
gu:=expand([seq(coeff(f1,x,i),i=1..K)]):
gu1:=subs(z=1,gu):
gu2:=[seq(lcoeff(gu[i],z),i=1..nops(gu))]:
gu3:=[seq(coeff(gu[i],z,degree(gu[i],z)-1),i=1..nops(gu))]:
[gu,gu1,gu2,gu3]:
end:


#PaperC(C,z,x,K1,K2,K3,K4): A paper about maximal T-avoiding 0-1 r by C matrices, for general r. Try:
#PaperC(4,z,x,30,100,6,10000):

PaperC:=proc(C,z,x,K1,K2,K3,K4) local f,A,r,f1,lu,mu,t0,avden:
t0:=time():
f:=GFrPC(C,z,x):

print(`Investigating Maximal T-avoiding configurations of 0-1 matrices with`, C , `columns and a general number of rows`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Let `, A[C](r)(z), ` be the weight-enumerator of MAXIMAL r by `, C, ` 0-1 matrices avoiding a T pattern `):
print(``):
print(matrix([[1,1,1],[`*`,1, `*`]])):
print(``):
print(`according to the weight z^NumberOfOnes (or equivalently, z^SumOfEntries)`):
print(``):
print(`Then `):
print(``):
print(Sum(A[C](r)(z)*x^r,r=0..infinity)=f):
print(``):
print(`and in Maple notation `):
print(``):
lprint(f):
print(``):
f1:=subs(z=1,f):
print(`It  follows that the generating function for the total number is`):
print(``):
print(Sum(A[C](r)(1)*x^r,r=0..infinity)=f1):
print(``):
lu:=SeqsR(C,z,K1):
print(`For the sake of the OEIS, the first`, K1, `terms are `):
print(``):
print(lu[2]):
print(``):
if C=3 then
 print(`This seems to be OEIS sequence A347725 `):
elif C=4 then
 print(`This seems to be OEIS sequence A189735 `):
elif C=5 then
 print(`This is not yet in the OEIS `):
fi:
print(``):

print(`The first`, K1, `terms of the efficient matrices (those with a maximal number of ones) are`):
print(``):
print(lu[3]):
print(``):
if C=4 then
 print(`This seems to be OEIS sequence A15518 `):
fi:


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

#
avden:=LimAvDen(C):

print(`First let us note that the EXACT limiting average density as the number or rows goes to infinity is`, avden):
print(``):
print(`and in Maple notation`):
print(``):
lprint(avden):
print(``):
print(`and in floating-point `):
print(``):
lprint(evalf(avden)):

print(``):

print(`Now let's specialize and look at the statistical distribution of`, K2, ` by `, C , `such matrices, and look at the staistics of the densitiy`):
print(``):

print(``):
lu:=DensStat(C,K2,K3):

print(``):
print(`The average, standard-deviation, and third through the`, K3, `scaled momenets are `):
print(``):
print(lu[1]):
print(``):
print(`Here is a plot of the distribution`):
print(``):
print(plot(lu[2])):
print(``):
print(`Now let's compare it to RSA simulation with`, K4, `random permutations of`, C*K2):
print(``):
mu:=RSAsimu(K2,C,K3,K4):
print(``):
print(`The average, s.d., and first few scaled moments up to the`, K3, `are, for this particular simulation `):
print(``):
print(mu[1]):
print(``):
print(`and a plot of the densitiy distribution is`):
print(``):
print(plot(mu[2])):
print(``):
print(`The ratio of the exact average (under the uniform distribution) and the average of the RSA simulations is`):
print(``):
print(mu[1][1]/lu[1][1]):
print(``):
print(`-------------------------------------`):
print(``):
print(`This ends this atricle that took`, time()-t0, `seconds to generate most of the time was spent by the RSA simulation`):
print(``):


end:

#AsyAv(f,x,z): Given a rational generating  function f of x and z signifying weight enumerators of some statistic, find the
#constant alpha such that the (average of the statistics)/n converges to it for items of size n.

#Try:
#AsyAv(1/(1-x*z-x^2),x,z);
AsyAv:=proc(f,x,z) local a,f0,f1,A1,A2,lu,i,aluf,si:

f0:=normal(subs(z=1,f)):

lu:=[solve(denom(f0),x)]:


aluf:=lu[1]:
si:=evalf(abs(lu[1])):

for i from 2 to nops(lu) do
 if evalf(abs(lu[i]))<si then
   aluf:=lu[i]:
   si:=evalf(abs(lu[i])):
 fi:
od:

a:=aluf:

f1:=normal(subs(z=1,diff(f,z))):



A1:=subs(x=a,numer(f0))/subs(x=a,diff(denom(f0),x)):

A2:=subs(x=a,numer(f1))/subs(x=a,diff(denom(f0),x))^2:

-A2/(A1*a):

end:


#LimAvDen(R): The limit of the average density of maximal T-avoiding R by c 0-1 matrices as c goes to infinity. Try:
#LimAvDen(5);
LimAvDen:=proc(R) local f,x,z:
f:=GFrPC(R,z,x):
AsyAv(f,x,z)/R:
end:



#LimAvDenG(C,SP): The limit of the average density of maximal SP-avoiding r by C 0-1 matrices as c goes to infinity. Try:
#LimAvDenG(5,{{[0,0],[0,1],[0,2],[1,1]});
LimAvDenG:=proc(R,SP) local f,x,z:
f:=GFr(R,SP,z,x):
AsyAv(f,x,z)/R:
end:


#New stuff

#checks whether everything except ref is 1s. always return -1 when out of bounds.
checkones := proc(c,xy,ref,poly)
local flag,p,newx,newy:
flag := false:
for p in poly do
	if p = ref then next: fi:
	newx := xy[1]+p[1]-ref[1]:
	newy := xy[2]+p[2]-ref[2]:
	if not(newx > 0 and newx <= nops(c) and newy > 0 and newy <= nops(c[1])) then
		return(-1):
	fi:

	if c[newx][newy] = 0 then  
		flag := true:
	fi:
od:
if flag then return(0): fi:
return(1):
end:





nextcolumns2 := proc(c,r,wid,polys)
local L,T,c3,flag,i,j,p,flag2,poly:
L := []:
T:=combinat:-cartprod([seq([seq(i,i=0..1)],j=1..r)]):
while not T[finished] do 
c3:=T[nextvalue]() :
if not(is_good2([op(2..nops(c),c),c3],polys)) then next: fi:
flag := false:		


for i from 1 to r do
	if c3[i] = 0 then next: fi:
	for poly in polys do
	for p in poly do
		if checkones([op(c),c3],[2*wid-1,i],p,poly) = 1 then
			flag := true:
			break:	
		fi:
	
		
	od:
	if flag then break: fi:
	od:
	if flag then break: fi:
od:	
if flag then next: fi:
flag := false:
for i from 1 to r do
	if c[wid][i] = 1 then next: fi:
	flag2 := false:
	for poly in polys do
	for p in poly do 
		if checkones([op(c),c3],[wid,i],p,poly) = 1 then 
			flag2 := true:
			break:
		fi:
	od:
	if flag2 then break: fi:
	od:
	if flag2=false then  
		flag := true:
		break:
	fi:
	
	
od:
if flag then next: fi:
L := [op(L),[op(2..nops(c),c),c3]]:
od: 
L:
end:

gen_states2 := proc(r,wid,polys):
local L,T,c3,i,j,m1,m2,p,cols:
L := []:

T:=combinat:-cartprod([seq([seq(i,i=0..1)],j=1..(2*wid-2)*r)]):
while not T[finished] do 
c3:=T[nextvalue]() :
cols := [seq([seq(c3[i],i=((j-1)*r+1)..(j*r))], j=1..(2*wid-2))]:
if is_good2(cols,polys) then  
	L:=[op(L),cols]:
fi:
od:
L:
end:

wi := proc(poly)
local m1,m2,p:
m1 := -100:
m2 := 100:

for p in poly do
if p[1] < m2 then m2 := p[1]: fi:
if p[1] > m1 then m1 := p[1]: fi:
od:
m1 - m2+1:
end:


TM:=proc(r,polys,z):
local S,M,i,t,w,m1,m2,p,wid,poly:

wid := -100:
for poly in polys do
if wi(poly) > wid then
	wid := wi(poly):
fi:
od:
S := gen_states2(r,wid,polys):
M:=Matrix(nops(S)+2,nops(S)+2):
for i from 1 to nops(S) do

	for t in nextcolumns2(S[i],r,wid,polys) do
		#print(t,ListTools:-Search(t,S)):
		M[i,ListTools:-Search(t,S)] := weight2(t,z):
	od:
	if final_state2(S[i],r,wid,polys) then 
		M[i,nops(S)+2] := 1:
	fi:
	if initial_state2(S[i],r,wid) then 
		M[nops(S)+1,i] := init_weight(S[i],z):
	fi:
od: 
M:
end:


initial_state2 := proc(c,r,wid)
local i:
for i from 1 to wid-1 do
if not(c[i] = [0$r]) then return(false): fi:
od:
return(true):
end:

final_state2 := proc(c,r,wid,polys)
local i,j,flag,p,poly:
for i from wid to nops(c) do
for j from 1 to r do

	if c[i][j] = 0 then
		flag := false:
		for poly in polys do
		for p in poly do
			if checkones(c,[i,j],p,poly) = 1 then
				flag := true:
				break:
			fi:
		od:
		if flag then break: fi:
		od:
		if flag = false then return(false): fi:
	fi: 

od:
od:

return(true):
end:


init_weight := proc(c,z) 
local s,ci:
s := 0:
for ci in c do
s := s + ListTools:-Occurrences(1,ci):
od:
z^s:
end:


weight2 := proc(c,z):
#1:
z^ListTools:-Occurrences(1,c[nops(c)]):
end:

count_sols2 := proc(rows, cols,poly)
local p:
p := (TM(rows,poly)^(cols+1))[-2][-1]:
print(p):
#print(subs(z=1,p)):
end:

is_good2 := proc(c, polys)
local co,i,j,flag,p,poly:
	co := 0:
	for i from 1 to nops(c) do
		if c[i] <> [0$(nops(c[1]))] then
			break:
		fi:
		co := co + 1:
	od:
	for i from (co + 1) to nops(c) do:
	for j from 1 to nops(c[1]) do
		flag := false:
		for poly in polys do
		for p in poly do
			if c[i][j] = 1 and checkones(c,[i,j],p,poly) = 1 then
				return(false):
			fi:
			if checkonesflex(c,[i,j],p,poly) = 1 then
				flag := true:
			fi:
		od:
		od:
		if c[i][j]=0 and not(flag) then
			return(false):
		fi:
	od:
	od:
	return(true):
end:

checkonesflex := proc(c,xy,ref,poly)
local flag,p,newx,newy:
flag := false:
for p in poly do
	if p = ref then next: fi:
	newx := xy[1]+p[1]-ref[1]:
	newy := xy[2]+p[2]-ref[2]:
	if newy <= 0 or newy > nops(c[1]) then
		return(-1):
	fi:
	if newx <= 0 or newx > nops(c) then 
		next:
	fi:

	if c[newx][newy] = 0 then  
		flag := true:
	fi:
od:
if flag then return(0): fi:
return(1):
end:



#GFrOld(r,S,z,x): The generating function in x for the number of maximal S-avoiding 0-1 matrices with r columns, according to z^Number of ones., where S is any set of polyonimoes. The
#coefficient of x^c is the weight-enumerator of c by  r such matrices. 

GFrOld:=proc(r,S,z,x) local A,d,f:
option remember:


A:=TM(r,S,z):
d:=Dimension(A)[1]:

f:=factor(evalm(((IdentityMatrix(d)-x*A))^(-1))[d-1,d]):




end:


#GFr(r,S,z,x): The generating function in x for the number of maximal S-avoiding 0-1 matrices with r columns, according to z^Number of ones., where S is any set of polyonimoes. The
#coefficient of x^c is the weight-enumerator of c by  r such matrices. 

GFr:=proc(r,S,z,x) local A,d,f,a,i,s:
option remember:

a:=max(seq(seq(s[1], s in S[i]),i=1..nops(S)))+1:
GFrN(r,S,z,x,a):


end:



#IsBad1(M,P): Given a pattern P and a 0-1 matrix decides whether it contains a pattern of ones that is a translate of P. Try:
#IsBad1([[0,1,0],[1,1,1],[0,1,0]],{[0,1],[1,0],[1,1],[1,2],[2,1]});
IsBad1:=proc(M,P) local i,j,P1,r,c,S,p1,p2:

r:=nops(M):
c:=nops(M[1]):

S:={seq(seq([i,j],j=1..c),i=1..r)}:
for i from 1 to r do
 for j from 1 to c do
  P1:={seq([i,j]+p1,p1 in P)}:
    if P1 minus S={} and {seq(M[p2[1]][p2[2]],p2 in P1)}={1} then
       RETURN(true):
    fi:
 od:
od:
false:
end:

 
#IsGood(M,S) : Is the matrix M good with repsect to the set of patterns S?
#Try:
#IsGood([[0,1,0],[1,1,1],[0,1,0]],{{[0,1],[1,0],[1,1],[1,2],[2,1]}});
IsGood:=proc(M,S)  local P:

for P in S do
 if IsBad1(M,P) then
  RETURN(false):
 fi:
od:
true:

end:



 
#IsMaxGood(M,S) : Is the matrix M maximally good with repsect to the set of patterns S?
#Try:
#IsMaxGood([[0,1,0],[1,1,1],[0,0,0]],{{[0,1],[1,0],[1,1],[1,2],[2,1]}});
IsMaxGood:=proc(M,S)  local P,i,j,M1:

if not IsGood(M,S) then
  RETURN(false):
fi:

for i from 1 to nops(M) do
 for j from 1 to nops(M[1]) do
  
  if M[i][j]=0 then

    M1:=[op(1..i-1,M),[op(1..j-1,M[i]),1,op(j+1..nops(M[i]),M[i])],op(i+1..nops(M),M)]:
    if  IsGood(M1,S) then
       RETURN(false):
    fi:
  fi:
 od:
od:
true:
end:

#Vecs(n): all 0-1 vectors of length n.
Vecs:=proc(n) local gu,gu1:
option remember:
if n=0 then
 RETURN({[]}):
fi:

gu:=Vecs(n-1):

{seq([op(gu1),0],gu1 in gu), seq([op(gu1),1],gu1 in gu)}:

end:

#Mats(r,c): all r by c 0-1 matrices. Try:
#Mats(3,3);
Mats:=proc(r,c) local gu,gu1,lu,lu1,ku:
option remember:
lu:=Vecs(c):
if r=1 then
  RETURN({seq([lu1],lu1 in lu)}):
fi:
gu:=Mats(r-1,c):

{seq(seq([op(gu1),lu1],gu1 in gu),lu1 in lu)}:
end:


#MGmats(r,c,S): all the r by c 0-1 matrices that  are maximally good with respect to the set of patterns  S. Try:
#MGmats(3,3,{{[0,1],[1,0],[1,1],[1,2],[2,1]}});
MGmats:=proc(r,c,S) local lu,lu1,gu:
option rememeber:
lu:=Mats(r,c):
gu:={}:
for lu1 in lu do
 if IsMaxGood(lu1,S) then
  gu:=gu union {lu1}:
 fi:
od:
gu:
end:


#WtMGmats(r,c,S,z): the weight-enumerator of all r by c 0-1 matrices that maximally good with respect to the set of patterns S. Try:
#WtMGmats(3,3,{{[0,1],[1,0],[1,1],[1,2],[2,1]}},z);
WtMGmats:=proc(r,c,S,z) local gu,gu1,i,j:
gu:=MGmats(r,c,S):

add(z^(add(add(gu1[i][j],j=1..c),i=1..r)),gu1 in gu):

end:



#CheckGFr(r,S,C): checks GFr(r,S,z,x) up to the coeff. x^C . Try:
#CheckGFr(3,{{[0,0],[0,1],[0,2],[1,1]}},4);
CheckGFr:=proc(r,S,C) local z,x,f,c,gu1,gu2:
f:=GFr(r,S,z,x):
f:=taylor(f,x=0,C+1):

gu1:=[seq(expand(coeff(f,x,c)),c=1..C)]:
gu2:= [seq(WtMGmats(c,r,S,z),c=1..C)]:

evalb(gu1=gu2):

end:


###start generalized RSA

#TransPol1(P,pt1,pt,r,s): inputs a pattern P a member of  P, pt1, and  a lattice point pt=[i,j] (1<=i<=r,1<=j<=s)
#outputs the translated pattern if pt is identified with pt1 and it is contained in [1,r]x[1,s] (if it does not it returns FAIL)
#Try:
#TransPol1({[0,0],[0,1],[0,2],[1,1]},[1,1],[3,5],6,7);
TransPol1:=proc(P,pt1,pt,r,s) local x,i,j,S:
if not member(pt1,P) then
 print(pt1, `should be a member of`, P):
 RETURN(FAIL):
fi:

if not (1<=pt[1] and pt[1]<=r and 1<=pt[2] and pt[2]<=s) then
 print(pt , `is not in `, [1,r]*[1,s]):
 RETURN(FAIL):
fi:

S:={seq(pt+x-pt1,x in P)}:

if S minus {seq(seq([i,j],j=1..s),i=1..r)}<>{} then
 FAIL:
else
 S:
fi:


end:

 



#TransPol(P,pt,r,s): inputs a pattern P  a lattice point pt=[i,j] (1<=i<=r,1<=j<=s)
#outputs the set of translated pattern if pt is identified with any of the members of P
#and it still in [1,r]x[1,s]. Try:
#Try:
#TransPol({[0,0],[0,1],[0,2],[1,1]},[3,5],6,7);
TransPol:=proc(P,pt,r,s) local pt1,S,lu:
S:={}:
for pt1 in P do
lu:=TransPol1(P,pt1,pt,r,s):
 if lu<>FAIL then
  S:=S union {lu}:
 fi:
od:
S:
end:



#TransPolS(SP,pt,r,s): inputs a set of patterns  SP  a lattice point pt=[i,j] (1<=i<=r,1<=j<=s)
#outputs the set of translated pattern if pt is identified with any of the members of any member P of SP
#and it still in [1,r]x[1,s]. Try:
#TransPolS({{[0,0],[0,1],[0,2],[1,1]}},[3,5],6,7);
TransPolS:=proc(SP,pt,r,s) local P:
{seq(op(TransPol(P,pt,r,s)), P in SP)}:
end:




#CanBHg(pt,SP,B,r,s) can you build a house on pt if there are houses at B (a subset of [1,s]x[q,r] and still avoid the patterns in the set of patterns SP. Try
#CanBHg([1,1],{{[0,0],[0,1],[0,2],[1,1]}}, {[3,1]},5,5 );
CanBHg:=proc(pt,SP,B,r,s) local B1,SP1,P1:

SP1:=TransPolS(SP,pt,r,s):

B1:=B union {pt}:

for P1 in SP1 do
 if P1 minus B1={} then
  RETURN(false):
 fi:
od:


true:
end:



#RSA1g(r,s,SP,pi): random sequential adsorption of an r by s lot avoiding the paterns in SP, using the permutation pi. Try:
#RSA1g(3,3,{{[0,0],[0,1],[0,2],[1,1]}}, [1,2,3,4,5,6,7,8,9]);
RSA1g:=proc(r,s,SP,pi) local L,i,j,B:
if nops(pi)<>r*s then
 RETURN(FAIL):
fi:

L:=[ seq(seq([i,j],j=1..s),i=1..r)]:


B:={}:

for i from 1 to nops(L) do
 if CanBHg(L[pi[i]],SP,B,r,s) then
  B:=B union {L[pi[i]]}:
 fi:
od:
B:
end:



#RSAg(SP,r,s): all maximal r by s {0,1} matrices that avoide the polyonlinoes in SP. Followed by the frequency table. Try
#RSAg({{[0,0],[0,1],[0,2],[1,1]}},3,3);
RSAg:=proc(SP,r,s) local lu,gu,T,pi,M,i:
lu:=permute(r*s):

gu:={}:

for i from 1 to nops(lu) do
 M:=SetToMat(RSA1g(r,s,SP,lu[i]),r,s):
 if not member(M,gu) then
   gu:=gu union {M}:
   T[M]:=1:
  else
  T[M]:=T[M]+1:
 fi:
od:

[gu,op(T)]:

end:


#RSAgGF(SP,r,s,z): the prob. generating function of r by s in z, of {0,1} matrices that avoide the polyonlinoes in SP. Try
#RSAgGF({{[0,0],[0,1],[0,2],[1,1]}},3,3,z);
RSAgGF:=proc(SP,r,s,z) local lu,pi:
lu:=permute(r*s):

add(z^nops(RSA1g(r,s,SP,pi,r,s)), pi in lu)/nops(lu):

end:


#RSAsimuG(r,s,SP,K1,K2): the analog of  DensStatG(r,s,SP,K1) but using RSA simulation with K2 random permutations. Try:
#RSAsimuG(3,10,{{[0,0],[0,1],[0,2],[1,1]}},6,1000);
RSAsimuG:=proc(r,s,SP,K1,K2) local gu,z,f,i,Dist,lu:
f:=add(z^nops(RSA1g(r,s,SP,randperm(r*s))),i=1..K2):
f:=f/subs(z=1,f):

Dist:=[]:
for i from ldegree(f,z) to degree(f,z) do
 Dist:=[op(Dist),[i/(r*s),evalf(coeff(f,z,i))]]:
od:

lu:=evalf(Alpha(f,z,K1),20):
lu:=[lu[1]/(r*s),lu[2]/(r*s),op(3..K1,lu)]:
lu,Dist:
end:



#DensStatG(r,s,SP,K): the density statistics for maximal configurations of an r by s lot avoiding the patterns of SP, up to the K-th moment.
#Try:
#DensStatG(3,50, {{[0,0],[0,1],[0,2],[1,1]}},6);
DensStatG:=proc(r,s,SP,K) local f,z,x,lu,Dist,i:
f:=GFr(r,SP,z,x):
f:=coeff(taylor(f,x=0,s+1),x,s):
f:=f/subs(z=1,f):

Dist:=[]:

for i from ldegree(f,z) to degree(f,z) do
 Dist:=[op(Dist),[i/(r*s),evalf(coeff(f,z,i))]]:
od:

lu:=evalf(Alpha(f,z,K),20):
lu:=[lu[1]/(r*s),lu[2]/(r*s),op(3..K,lu)]:
lu,Dist:

end:





#RSAcollectGall(r,s,SP,K): collects many maximal r by s 0-1 matrices avoiding the patterns in SP,  using K random permutations.
#Also returns how many are missing. Try:
#RSAcollectG(3,3,{{[0,0],[0,1],[0,2],[1,1]}},1000):
RSAcollectGall:=proc(r,s,SP,K) local gu,f,co,gu1,i:

gu:={seq(RSA1g(r,s,SP,randperm(r*s)),i=1..K)}:

f:=subs(z=1,GFr(s,SP,z,x)):
co:=coeff(taylor(f,x=0,r+1),x,r):

gu:={seq(SetToMat(gu1,r,s),gu1 in gu)}:

[gu,co-nops(gu)]:

end:


#RSAcollectG(r,s,SP,K): collects many maximal r by s 0-1 matrices avoiding the patterns in SP,  using K random permutations.
#RSAcollect(3,3,{{[0,0],[0,1],[0,2],[1,1]}},1000):
RSAcollectG:=proc(r,s,SP,K) local gu,gu1,i:

gu:={seq(RSA1g(r,s,SP,randperm(r*s)),i=1..K)}:
{seq(SetToMat(gu1,r,s),gu1 in gu)}:
end:




#PaperCsp(C,SP,z,x,K1,K2,K3,K4): A paper about maximal SP-avoiding 0-1 r by C matrices, for general x. Try:
#PaperCsp(C,z,x,30,100,6,10000):

PaperCsp:=proc(C,SP,z,x,K1,K2,K3,K4) local f,A,r,f1,lu,mu,t0,avden,i:
t0:=time():
f:=GFr(C,SP,z,x):

print(`Investigating Maximal`, SP, `avoiding configurations of 0-1 matrices with`, C , `columns and a general number of rows`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Let `, A[C](r)(z), ` be the weight-enumerator of MAXIMAL r by `, C, ` 0-1 matrices avoiding a  pattern in the set`):
print(``):
print(SP):
print(``):
print(`according to the weight z^NumberOfOnes (or equivalently, z^SumOfEntries)`):
print(``):
print(`Then `):
print(``):
print(Sum(A[C](r)(z)*x^r,r=0..infinity)=f):
print(``):
print(`and in Maple notation `):
print(``):
lprint(f):
print(``):
f1:=subs(z=1,f):
print(`It  follows that the generating function for the total number is`):
print(``):
print(Sum(A[C](r)(1)*x^r,r=0..infinity)=f1):
print(``):
lu:=[seq(coeff(taylor(f1,x=0,K1+1),x,i),i=1..K1)]:
print(`For the sake of the OEIS, the first`, K1, `terms are `):
print(``):
print(lu):
print(``):




avden:=evalf(AsyAv(f,x,z)/C):

print(`First let us note that the EXACT (In floating point) limiting average density as the number or rows goes to infinity is`, avden):


print(``):

print(`Now let's specialize and look at the statistical distribution of`, K2, ` by `, C , `such matrices, and look at the staistics of the densitiy`):
print(``):

print(``):

lu:=DensStatG(C,K2,SP,K3):

print(``):
print(`The average, standard-deviation, and third through the`, K3, `scaled momenets are `):
print(``):
print(lu[1]):
print(``):
print(`Here is a plot of the distribution`):
print(``):
print(plot(lu[2])):
print(``):
print(`Now let's compare it to RSA simulation with`, K4, `random permutations of`, C*K2):
print(``):
mu:=RSAsimuG(K2,C,SP,K3,K4):
print(``):
print(`The average, s.d., and first few scaled moments up to the`, K3, `are, for this particular simulation `):
print(``):
print(mu[1]):
print(``):
print(`and a plot of the densitiy distribution is`):
print(``):
print(plot(mu[2])):
print(``):
print(`The ratio of the exact average (under the uniform distribution) and the average of the RSA simulations is`):
print(``):
print(mu[1][1]/lu[1][1]):
print(``):
print(`-------------------------------------`):
print(``):
print(`This ends this atricle that took`, time()-t0, `seconds to generate most of the time was spent by the RSA simulation`):
print(``):


end:



#GFrDimerPC(C,z,x): the precumputed generating function for the weight-enumearor of r by C (c from 1 to 5) 0-1 matrices that aviod dimers ( i.e. avoid two consectutive horizontal 1's and 
#avoid two consecutive vertical 1s). Try:
#GFrDimerPC
GFrDimerPC:=proc(C,z,x) local L:
L:=
[-x*z*(x+1)^2/(x^3*z+x^2*z-1), -2*x*z/(x^2*z+x*z-1), -(2*x^5*z^4+2*x^3*z^3+2*x^
2*z^3-6*x^2*z^2-x*z^2-x*z-z-1)*x*z/(x^5*z^4+2*x^4*z^4-x^4*z^3+x^3*z^4-4*x^3*z^3
-x^2*z^3-x*z+1), (x^6*z^6-x^5*z^6+x^5*z^5-2*x^5*z^4-3*x^4*z^4+2*x^3*z^4-7*x^3*z
^3+2*x^3*z^2-4*x^2*z^3+7*x^2*z^2-x*z^2+4*x*z+3)*x*z^2/(x^6*z^6+x^5*z^6+x^5*z^5+
2*x^4*z^5-x^4*z^4+2*x^3*z^5-4*x^3*z^4-x^3*z^3-2*x^2*z^3-x^2*z^2-x*z^2+1), -(3*x
^13*z^17-x^13*z^16-4*x^12*z^17+10*x^12*z^16+2*x^11*z^17-9*x^11*z^16-9*x^11*z^15
+17*x^11*z^14+18*x^10*z^15+x^11*z^13-18*x^10*z^14-4*x^9*z^15-22*x^10*z^13-8*x^9
*z^14+4*x^10*z^12+47*x^9*z^13+2*x^8*z^14+4*x^9*z^12-x^9*z^11-43*x^8*z^12-12*x^8
*z^11+9*x^7*z^12+3*x^8*z^10+2*x^7*z^11+40*x^7*z^10-10*x^6*z^11-30*x^7*z^9+36*x^
6*z^10-x^7*z^8-74*x^6*z^9+2*x^5*z^10+33*x^6*z^8-19*x^5*z^9-5*x^6*z^7+45*x^5*z^8
-25*x^5*z^7+13*x^4*z^8-x^5*z^6-33*x^4*z^7+11*x^4*z^6-2*x^3*z^7-9*x^4*z^5+7*x^3*
z^6-26*x^3*z^5+11*x^3*z^4-5*x^2*z^5+7*x^2*z^4+9*x^2*z^3+x^2*z^2+x*z^3+3*x*z^2+2
*x*z+z+3)*x*z^2/(x^13*z^17-2*x^12*z^17+2*x^12*z^16+x^11*z^17+2*x^11*z^16-5*x^11
*z^15-4*x^10*z^16+2*x^11*z^14+10*x^10*z^15+x^9*z^16-6*x^10*z^14-3*x^9*z^15+4*x^
10*z^13+5*x^9*z^14-14*x^9*z^13-2*x^8*z^14-4*x^9*z^12+10*x^8*z^13-x^9*z^11+3*x^8
*z^12-3*x^7*z^13-7*x^8*z^11+10*x^7*z^12-2*x^8*z^10-26*x^7*z^11+x^6*z^12+19*x^7*
z^10-2*x^6*z^11-2*x^7*z^9+12*x^6*z^10-4*x^6*z^9+3*x^5*z^10-6*x^6*z^8-3*x^5*z^9-\
22*x^5*z^8-x^4*z^9+x^5*z^7+x^5*z^6-8*x^4*z^7+10*x^4*z^6-x^3*z^7+3*x^4*z^5+11*x^
3*z^5+x^2*z^5+x^2*z^4+2*x^2*z^3+x*z^2-1)]:
if C>=1 and C<=5 then
 RETURN(L[C]):
else
GFr(C,{{[0,0],[0,1]},{[0,0],[1,0]}},z,x):
fi:
end:




#GFrKingsPC(C,z,x): The pre-comuted generating function according to the weight x^r*z^NumberOfKings in all maximal ways of placing non-attacking kings on an r by C chess board.
#If C s between 1 and 5 it is very fast (pre-comuted), but if C>=6 is is very slow. Try:
#GFrKingsPC(5,z,x);
GFrKingsPC:=proc(C,z,x) local L:
L:=
[-x*z*(x+1)^2/(x^3*z+x^2*z-1), -2*x*z*(x+1)^2/(2*x^3*z+2*x^2*z-1), -(x^5*z^3+x^
5*z^2-x^3*z^3+x^3*z+2*x^2*z^2+x^2*z+x*z^2-x^2-3*x*z-2*x-z-1)*x*z/(x^6*z^4+x^6*z
^3-x^5*z^4-x^5*z^3+x^4*z^3+x^4*z^2+x^3*z^3-x^3*z^2-x^3*z-x^2*z^2-x^2*z-x*z+1),
-(6*x^6*z^3+9*x^5*z^2-6*x^4*z^3+3*x^4*z^2-3*x^3*z^2+3*x^3*z+3*x^2*z^2+2*x^2*z-3
*x^2+3*x*z-12*x-3)*x*z^2/(6*x^7*z^5-6*x^6*z^5+9*x^6*z^4-6*x^5*z^4+3*x^4*z^4+x^4
*z^3+3*x^3*z^3-6*x^3*z^2-4*x^2*z^2-x*z+1), -(x^18*z^17+4*x^18*z^16+2*x^17*z^17+
3*x^18*z^15+7*x^17*z^16+x^16*z^17+x^17*z^15+3*x^16*z^16-7*x^17*z^14-3*x^16*z^15
-2*x^15*z^16-3*x^17*z^13-11*x^16*z^14-2*x^15*z^15-2*x^14*z^16-6*x^16*z^13+10*x^
15*z^14-2*x^14*z^15-10*x^15*z^13+12*x^14*z^14-12*x^15*z^12-6*x^14*z^13+x^13*z^
14+x^12*z^15-15*x^14*z^12-4*x^13*z^13+3*x^12*z^14+12*x^14*z^11-30*x^13*z^12-2*x
^12*z^13+9*x^14*z^10-29*x^13*z^11-19*x^12*z^12+2*x^11*z^13-6*x^13*z^10-65*x^12*
z^11+4*x^11*z^12-4*x^10*z^13+4*x^12*z^10+12*x^11*z^11+2*x^10*z^12+12*x^12*z^9+
31*x^11*z^10+25*x^10*z^11+27*x^11*z^9-21*x^10*z^10-3*x^9*z^11-27*x^11*z^8+18*x^
10*z^9-22*x^9*z^10+6*x^8*z^11-9*x^11*z^7+31*x^10*z^8-30*x^9*z^9-6*x^8*z^10+9*x^
10*z^7+71*x^9*z^8-22*x^8*z^9-12*x^9*z^7+36*x^8*z^8+x^7*z^9+20*x^8*z^7+22*x^7*z^
8-4*x^6*z^9+14*x^8*z^6+52*x^7*z^7+2*x^6*z^8+31*x^8*z^5+8*x^7*z^6-7*x^6*z^7+3*x^
8*z^4+6*x^7*z^5-45*x^6*z^6-x^5*z^7-3*x^7*z^4-53*x^6*z^5-6*x^5*z^6+x^4*z^7+2*x^6
*z^4+17*x^5*z^5+3*x^4*z^6-3*x^6*z^3+10*x^5*z^4+18*x^4*z^5-15*x^5*z^3-x^4*z^4+x^
3*z^5-9*x^5*z^2-34*x^4*z^3-4*x^3*z^4-5*x^4*z^2-16*x^3*z^3-2*x^2*z^4+4*x^3*z^2-4
*x^2*z^3-3*x^3*z-4*x^2*z^2-x^2*z-x*z^2+3*x^2+5*x*z+12*x+z+3)*x*z^2/(x^19*z^19+4
*x^19*z^18+x^18*z^19+3*x^19*z^17+3*x^18*z^18-2*x^18*z^17-7*x^18*z^16-x^17*z^17-\
2*x^16*z^18-3*x^18*z^15-4*x^17*z^16-x^16*z^17-3*x^17*z^15+13*x^16*z^16-x^15*z^
17-8*x^16*z^15-2*x^15*z^16+x^14*z^17-12*x^16*z^14+2*x^15*z^15+3*x^14*z^16-x^15*
z^14-6*x^14*z^15+13*x^15*z^13-28*x^14*z^14+6*x^13*z^15+9*x^15*z^12-41*x^14*z^13
+4*x^13*z^14-4*x^12*z^15-15*x^14*z^12-23*x^13*z^13+x^12*z^14+23*x^13*z^12+34*x^
12*z^13+12*x^13*z^11+10*x^12*z^12-9*x^11*z^13+12*x^12*z^11-26*x^11*z^12+6*x^10*
z^13-30*x^12*z^10+18*x^11*z^11-2*x^10*z^12-9*x^12*z^9+66*x^11*z^10-38*x^10*z^11
+18*x^11*z^9+5*x^9*z^11-34*x^10*z^9+35*x^9*z^10-4*x^8*z^11+46*x^9*z^9-4*x^8*z^
10+22*x^9*z^8+10*x^8*z^9+34*x^9*z^7-32*x^8*z^8-2*x^7*z^9+3*x^9*z^6-34*x^8*z^7-\
18*x^7*z^8+x^6*z^9-6*x^8*z^6-9*x^7*z^7+7*x^6*z^8+8*x^7*z^6+12*x^6*z^7-3*x^7*z^5
+x^6*z^6+x^5*z^7-22*x^6*z^5+x^5*z^6-10*x^6*z^4-3*x^5*z^5-3*x^4*z^6+7*x^5*z^4-7*
x^4*z^5-6*x^4*z^4-2*x^4*z^3-2*x^3*z^4+4*x^3*z^3+6*x^3*z^2+3*x^2*z^3+4*x^2*z^2+x
*z-1)]:

if (C>=1 and C<=5) then
 L[C]:
else
GFr(C,{{[0,0],[0,1]},{[0,0],[1,0]},{[0,0],[1,1]},{[0,1],[1,0]}},z,x):
fi:
end:






#PaperCdimer(C,z,x,K1,K2,K3,K4): A paper about maximal Dimer-avoiding 0-1 r by C matrices, for general r. Try:
#PaperCdimer(4,z,x,30,100,6,10000):

PaperCdimer:=proc(C,z,x,K1,K2,K3,K4) local f,A,r,f1,lu,mu,t0,avden,i:
t0:=time():

f:=GFrDimerPC(C,z,x):


print(`Investigating Maximal Dimer-avoiding configurations of 0-1 matrices with`, C , `columns and a general number of rows`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Let `, A[C](r)(z), ` be the weight-enumerator of MAXIMAL r by `, C, ` 0-1 matrices avoiding Dimers, i.e. two consecutive horizontal 1s and also avoiding two consecutive vertical ones `):
print(``):
print(`according to the weight z^NumberOfOnes (or equivalently, z^SumOfEntries)`):
print(``):
print(`Then `):
print(``):
print(Sum(A[C](r)(z)*x^r,r=0..infinity)=f):
print(``):
print(`and in Maple notation `):
print(``):
lprint(f):
print(``):
f1:=subs(z=1,f):
print(`It  follows that the generating function for the total number is`):
print(``):
print(Sum(A[C](r)(1)*x^r,r=0..infinity)=f1):
print(``):
lu:=[seq(coeff(taylor(f1,x=0,K1+2),x,i),i=1..K1)]:
print(`For the sake of the OEIS, the first`, K1, `terms are `):
print(``):
print(lu):
print(``):

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


avden:=evalf(AsyAv(f,x,z)/C):
print(`First let us note that the EXACT limiting average density as the number or rows goes to infinity is`, avden):
print(``):
print(``):

print(`Now let's specialize and look at the statistical distribution of`, K2, ` by `, C , `such matrices, and look at the staistics of the densitiy`):
print(``):

print(``):

lu:=DensStatG(C,K2,{{[0,0],[0,1]} , {[0,0],[1,0]}},K3):



print(``):
print(`The average, standard-deviation, and third through the`, K3, `scaled momenets are `):
print(``):
print(lu[1]):
print(``):
print(`Here is a plot of the distribution`):
print(``):
print(plot(lu[2])):
print(``):
print(`Now let's compare it to RSA simulation with`, K4, `random permutations of`, C*K2):
print(``):
mu:=RSAsimuG(K2,C,{{[0,0],[0,1]},{[0,0],[1,0]}}, K3,K4):
print(``):
print(`The average, s.d., and first few scaled moments up to the`, K3, `are, for this particular simulation `):
print(``):
print(mu[1]):
print(``):
print(`and a plot of the densitiy distribution is`):
print(``):
print(plot(mu[2])):
print(``):
print(`The ratio of the exact average (under the uniform distribution) and the average of the RSA simulations is`):
print(``):
print(mu[1][1]/lu[1][1]):
print(``):
print(`-------------------------------------`):
print(``):
print(`This ends this atricle that took`, time()-t0, `seconds to generate most of the time was spent by the RSA simulation`):
print(``):


end:





#PaperCkings(C,z,x,K1,K2,K3,K4): A paper about maximal  Non-attacking kinds 0-1 r by C matrices, for general r. Try:
#PaperCkings(4,z,x,30,100,6,10000):

PaperCkings:=proc(C,z,x,K1,K2,K3,K4) local f,A,r,f1,lu,mu,t0,avden,i:
t0:=time():

f:=GFrKingsPC(C,z,x):

#KAKI
print(`Investigating Maximal  configurations of Placing non-attacking kings on a rectangular chessboard with`, C , `columns and a general number of rows`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Let `, A[C](r)(z), ` be the weight-enumerator of MAXIMAL r by `, C, ` 0-1 matrices avoiding adjacent Ones either horizontally, vertically, or diagonally `):
print(``):
print(`according to the weight z^NumberOfOnes (or equivalently, z^SumOfEntries`):
print(``):
print(`Then `):
print(``):
print(Sum(A[C](r)(z)*x^r,r=0..infinity)=f):
print(``):
print(`and in Maple notation `):
print(``):
lprint(f):
print(``):
f1:=subs(z=1,f):
print(`It  follows that the generating function for the total number is`):
print(``):
print(Sum(A[C](r)(1)*x^r,r=0..infinity)=f1):
print(``):
lu:=[seq(coeff(taylor(f1,x=0,K1+2),x,i),i=1..K1)]:
print(`For the sake of the OEIS, the first`, K1, `terms are `):
print(``):
print(lu):
print(``):

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


avden:=evalf(AsyAv(f,x,z)/C):
print(`First let us note that the EXACT limiting average density as the number or rows goes to infinity is`, avden):
print(``):
print(``):

print(`Now let's specialize and look at the statistical distribution of`, K2, ` by `, C , `such matrices, and look at the staistics of the densitiy`):
print(``):

print(``):

lu:=DensStatG(C,K2,{{[0,0],[0,1]} , {[0,0],[1,0]}, {[0,0],[1,1]},{[0,1],[1,0]}},K3):



print(``):
print(`The average, standard-deviation, and third through the`, K3, `scaled momenets are `):
print(``):
print(lu[1]):
print(``):
print(`Here is a plot of the distribution`):
print(``):
print(plot(lu[2])):
print(``):
print(`Now let's compare it to RSA simulation with`, K4, `random permutations of`, C*K2):
print(``):
mu:=RSAsimuG(K2,C,{{[0,0],[0,1]},{[0,0],[1,0]}, {[0,0],[1,1]},{[0,1],[1,0]}}, K3,K4):
print(``):
print(`The average, s.d., and first few scaled moments up to the`, K3, `are, for this particular simulation `):
print(``):
print(mu[1]):
print(``):
print(`and a plot of the densitiy distribution is`):
print(``):
print(plot(mu[2])):
print(``):
print(`The ratio of the exact average (under the uniform distribution) and the average of the RSA simulations is`):
print(``):
print(mu[1][1]/lu[1][1]):
print(``):
print(`-------------------------------------`):
print(``):
print(`This ends this atricle that took`, time()-t0, `seconds to generate most of the time was spent by the RSA simulation`):
print(``):


end:




#GFrN(r,S,z,x,a): The generating function in x for the number of maximal S-avoiding 0-1 matrices with r columns, according to z^Number of ones., where S is any set of patterns. The
#coefficient of x^c is the weight-enumerator of c by  r such matrices. 

GFrN:=proc(r,S,z,x,a) local A,d,f,a1,C1,C2,lu1,lu2,ka,i,vu,vu1,c:
option remember:


f:=GFrOld(r,S,z,x):

lu1:=expand([seq(coeff(taylor(f,x=0,a+1),x,a1),a1=1..a)]):
lu2:=[seq(WtMGmats(c,r,S,z),c=1..a)]:
vu:=add(WtMGmats(c,r,S,z)*x^c,c=1..a):

 ka:=lu2[nops(lu2)]:

for i from 1 to nops(lu1) while lu1[i]<>ka do od:

if i<=nops(lu1) then
f:=normal(x^(nops(lu1)-i)*f):
vu1:=add(expand(coeff(taylor(f,x=0,nops(lu1)+1),x,c))*x^c,c=1..a):
RETURN(normal(1+f+vu-vu1)):

fi:

ka:=lu1[nops(lu1)]:

for i from 1 to nops(lu2) while lu2[i]<>ka do od:
if i=nops(lu2)+1 then
 RETURN(FAIL):
fi:

f:=normal(f/x^(nops(lu2)-i)):

vu1:=add(expand(coeff(taylor(f,x=0,nops(lu1)+1),x,c))*x^c,c=1..a):
RETURN(normal(1+f+vu-vu1)):


end:


#CheckGFrN(r,S,a,C): checks GFNr(r,S,z,x,a) up to the coeff. x^C . Try:
#CheckGFrN(3,{{[0,0],[0,1],[0,2],[1,1]}},4);
CheckGFrN:=proc(r,S,a,C) local z,x,f,c,gu1,gu2:
f:=GFrN(r,S,z,x,a):

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

f:=taylor(f,x=0,C+1):

gu1:=[seq(expand(coeff(f,x,c)),c=1..C)]:
gu2:= [seq(WtMGmats(c,r,S,z),c=1..C)]:

evalb(gu1=gu2):

#evalb([seq(expand(coeff(f,x,c))-WtMGmats(c,r,S,z),c=2..C)]=[0$(C-1)]):
end:


#GF1rowA(k,z,x): The bi-variate generating function whose coefficient of x^n is the weight-enumerator, according to z^NumberOfOccupiedSeats in maximal seating arrangements
#with n seats in one row, obeying the social-distancing requirement that you can't have a block of k continguous occupied seats, in other words the
#weight enumerator of 0-1 vectors of length n avoiding the factor 1...1 (1 repeated k times), and such that if you change any 0 to a 1 you would
#get a violation. Try:
#GF1rowA(3,z,x);
GF1rowA:=proc(k,z,x) local i,S:
 S:={{seq([i,0],i=0..k-1)}}:
 GFr(1,S,z,x) :
end:


#GF1rowB(k,z,x): The bi-variate generating function whose coefficient of x^n is the weight-enumerator, according to z^NumberOfOccupiedSeats in maximal seating arrangements
#with n seats in one row, obeying the social-distancing requirement that any two people must have at least k empty seats between them. In other words
#weight enumerator of 0-1 vectors of length n avoiding the patterns 11, 1*1, 1**1, ... 1*^(k-1)1
#get a violation. Try:
#GF1rowB(3,z,x);
GF1rowB:=proc(k,z,x) local i,S:
 S:={seq({[0,0],[i,0]},i=1..k)}:
 GFr(1,S,z,x) :
end:

#GF1rowBconj(k,z,x): The conjectured expression for GF1rowB(k,z,x). Try:
#evalb([seq(GF1rowB(k1,z,x),k1=1..5)]=[seq(GF1rowBconj(k1,z,x),k1=1..5)]);
GF1rowBconj:=proc(k,z,x):
simplify(normal((x^k-1)/(1-x)/x^k+(1-x^(k+1))/x^k/((1-x)-x^(k+1)*(1-x^(k+1))*z))):
end:




#Paper1RowA(B,z,x): Inputs a positive integer B and outputs a paper with the bi-variate generating function of maximal sittings on one row
#with at most b consecutive occupied seats for b from 2 to B
#Paper1RowA(5,z,x);
Paper1RowA:=proc(B,z,x) local b,gu,mu,gu1,mu1,t0:

t0:=time():

print(`Generating functions and Limiting Densitiy for maximal sitting arrangements in 1 row where the Social Distancing  prohibits b adjacent seats for b from 2 to`, B):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let A[b](m,n) be the number of ways that m people can be seated in a row of n chairs such that is is forbidden to have more than b consecutive occupied seatgs and such that`):
print(`no emtpy seat can be occupied without violating this restriction. Let f[b](z,x) be the bi-variate generating function`):
print(``):
print(f[b](z,x)= Sum(Sum(A[b](m,n)*z^m,m=0..n)*x^n,n=0..infinity) ):
print(``):
gu:=[]:
mu:=[]:
print(`We have`):

for b from 2 to B do
 gu1:=GF1rowA(b,z,x):
print(``):
print(f[b](z,x)=gu1):
print(``):
print(`and in Maple notation`):
print(``):
lprint(f[b](z,x)=gu1):
print(``):
gu:=[op(gu),gu1]:


mu1:=evalf(AsyAv(gu1,x,z)):
print(``):
print(`The limiting average density of occupied chairs, as n goes to infinity  is`):
print(``):
lprint(mu1):
print(``):
mu:=[op(mu),mu1]:

od:
print(``):
print(`-------------------------------`):
print(``):
print(`To sum up the list of generating functions from b=1 to b=`, B, ` is `):
print(``):
lprint(gu):
print(``):
print(`and the list of limiting average densities from b=1 to b=`, B, ` is `):
print(``):
lprint(mu):
print(``):
print(`------------------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `second to generate `):
print(``):

end:




#Paper1RowB(B,z,x): Inputs a positive integer B and outputs a paper with the bi-variate generating function of maximal sittings on one row
#with minimimal distance>=b for b from 1 to B. Try:
#Paper1RowB(5,z,x);
Paper1RowB:=proc(B,z,x) local b,gu,mu,gu1,mu1,t0,CON,b1:

t0:=time():

print(`Generating functions and Limiting Densitiy for maximal sitting arrangements in 1 row where the Social Distancing requires at least b empty seats between two poeple for b from 1 to`, B):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let A[b](m,n) be the number of ways that m people can be seated in a row of n chairs such that  you need at least b empty seats between people, for b from 1 to`, B):
print(`no emtpy seat can be occupied without violating this. Let f[b](z,x) be the bi-variate generating function`):
print(``):
print(f[b](z,x)= Sum(Sum(A[b](m,n)*z^m,m=0..n)*x^n,n=0..infinity) ):
print(``):
gu:=[]:
mu:=[]:
print(`It seems that a closed-form expression for SYMBOLIC (i.e. ALL) b is`):
print(``):
CON:=GF1rowBconj(b,z,x):
print(CON):
print(``):
print(` and in Maple notation `):
print(``):
lprint(CON):

print(``):

print(`We have:`):

for b1 from 1 to B do
 gu1:=GF1rowB(b1,z,x):
print(``):
print(f[b1](z,x)=gu1):
print(``):
print(`and in Maple notation`):
print(``):
lprint(f[b1](z,x)=gu1):
print(``):
gu:=[op(gu),gu1]:

if normal(subs(b=b1,CON)-gu1)<>0 then
 print(`The conjecture is wrong`):
else
print(`The conjecture is confirmed for b=`, b1):
fi:


mu1:=evalf(AsyAv(gu1,x,z)):
print(``):
print(`The limiting average density of occupied chairs, as n goes to infinity  is`):
print(``):
lprint(mu1):
print(``):
mu:=[op(mu),mu1]:

od:
print(``):
print(`-------------------------------`):
print(``):
print(`To sum up the list of generating functions from b=1 to b=`, B, ` is `):
print(``):
lprint(gu):
print(``):
print(`and the list of limiting average densities from b=1 to b=`, B, ` is `):
print(``):
lprint(mu):
print(``):
print(`and we also confimed the conjecture for all b from 1 to`, B):
print(``):
print(`------------------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `second to generate `):
print(``):

end:


####Start the R22 case


#GFu(n,z): The prob. generating function for max. 11 avoiding vectors of lentgth n for the r.v. number of ones. TryL 
#GFu(10,z);
GFu:=proc(n,z) local x,f:
f:=expand(coeff(taylor((x^2*z+x*z+1)/(x^3*z+x^2*z-1),x=0,n+1),x,n)):
f/subs(z=1,f):
end:

#fn(n,z): The prob. generating function for the number of ones in maximal 11-avoiding 0-1 vectors using RSA. Try:
#fn(100,z);
fn:=proc(n,z) local i:
option remember:

if n=0 then
 1:
elif n=1 then
 z:
elif n=2 then
 z:
elif n=3 then
 2/3*z^2+1/3*z:
else
expand(z/n*(2*fn(n-2,z)+add(fn(i-2,z)*fn(n-i-1,z),i=2..n-1))):
fi:



end:
####

#KingStory5(n): The story for non-attacking kinds for a 5 times n board
KingStory5:=proc(n) local f,x,z,gu,gu1,mu1:
f:=
-(x^18*z^17+4*x^18*z^16+2*x^17*z^17+3*x^18*z^15+7*x^17*z^16+x^16*z^17+x^17*z^15
+3*x^16*z^16-7*x^17*z^14-3*x^16*z^15-2*x^15*z^16-3*x^17*z^13-11*x^16*z^14-2*x^
15*z^15-2*x^14*z^16-6*x^16*z^13+10*x^15*z^14-2*x^14*z^15-10*x^15*z^13+12*x^14*z
^14-12*x^15*z^12-6*x^14*z^13+x^13*z^14+x^12*z^15-15*x^14*z^12-4*x^13*z^13+3*x^
12*z^14+12*x^14*z^11-30*x^13*z^12-2*x^12*z^13+9*x^14*z^10-29*x^13*z^11-19*x^12*
z^12+2*x^11*z^13-6*x^13*z^10-65*x^12*z^11+4*x^11*z^12-4*x^10*z^13+4*x^12*z^10+
12*x^11*z^11+2*x^10*z^12+12*x^12*z^9+31*x^11*z^10+25*x^10*z^11+27*x^11*z^9-21*x
^10*z^10-3*x^9*z^11-27*x^11*z^8+18*x^10*z^9-22*x^9*z^10+6*x^8*z^11-9*x^11*z^7+
31*x^10*z^8-30*x^9*z^9-6*x^8*z^10+9*x^10*z^7+71*x^9*z^8-22*x^8*z^9-12*x^9*z^7+
36*x^8*z^8+x^7*z^9+20*x^8*z^7+22*x^7*z^8-4*x^6*z^9+14*x^8*z^6+52*x^7*z^7+2*x^6*
z^8+31*x^8*z^5+8*x^7*z^6-7*x^6*z^7+3*x^8*z^4+6*x^7*z^5-45*x^6*z^6-x^5*z^7-3*x^7
*z^4-53*x^6*z^5-6*x^5*z^6+x^4*z^7+2*x^6*z^4+17*x^5*z^5+3*x^4*z^6-3*x^6*z^3+10*x
^5*z^4+18*x^4*z^5-15*x^5*z^3-x^4*z^4+x^3*z^5-9*x^5*z^2-34*x^4*z^3-4*x^3*z^4-5*x
^4*z^2-16*x^3*z^3-2*x^2*z^4+4*x^3*z^2-4*x^2*z^3-3*x^3*z-4*x^2*z^2-x^2*z-x*z^2+3
*x^2+5*x*z+12*x+z+3)*x*z^2/(x^19*z^19+4*x^19*z^18+x^18*z^19+3*x^19*z^17+3*x^18*
z^18-2*x^18*z^17-7*x^18*z^16-x^17*z^17-2*x^16*z^18-3*x^18*z^15-4*x^17*z^16-x^16
*z^17-3*x^17*z^15+13*x^16*z^16-x^15*z^17-8*x^16*z^15-2*x^15*z^16+x^14*z^17-12*x
^16*z^14+2*x^15*z^15+3*x^14*z^16-x^15*z^14-6*x^14*z^15+13*x^15*z^13-28*x^14*z^
14+6*x^13*z^15+9*x^15*z^12-41*x^14*z^13+4*x^13*z^14-4*x^12*z^15-15*x^14*z^12-23
*x^13*z^13+x^12*z^14+23*x^13*z^12+34*x^12*z^13+12*x^13*z^11+10*x^12*z^12-9*x^11
*z^13+12*x^12*z^11-26*x^11*z^12+6*x^10*z^13-30*x^12*z^10+18*x^11*z^11-2*x^10*z^
12-9*x^12*z^9+66*x^11*z^10-38*x^10*z^11+18*x^11*z^9+5*x^9*z^11-34*x^10*z^9+35*x
^9*z^10-4*x^8*z^11+46*x^9*z^9-4*x^8*z^10+22*x^9*z^8+10*x^8*z^9+34*x^9*z^7-32*x^
8*z^8-2*x^7*z^9+3*x^9*z^6-34*x^8*z^7-18*x^7*z^8+x^6*z^9-6*x^8*z^6-9*x^7*z^7+7*x
^6*z^8+8*x^7*z^6+12*x^6*z^7-3*x^7*z^5+x^6*z^6+x^5*z^7-22*x^6*z^5+x^5*z^6-10*x^6
*z^4-3*x^5*z^5-3*x^4*z^6+7*x^5*z^4-7*x^4*z^5-6*x^4*z^4-2*x^4*z^3-2*x^3*z^4+4*x^
3*z^3+6*x^3*z^2+3*x^2*z^3+4*x^2*z^2+x*z-1):

print(`The number of Maximal Non-Attacking Kings Configurations on a 5 times`, n, `chessboard `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`The number of ways of placing m non-attacking kings on a 5 times n chessboard is the coefficient of z^m x^n in the Taylor expansion of`):
print(``):
print(f):
print(``):
print(`and in Maple notation`):
print(``):
lprint(f):
print(``):
gu:=expand(coeff(taylor(f,x=0,n+1),x,n)):
gu1:=subs(z=1,gu):
print(`In particular the total number on a 5 times`, n, `chessboard is`):
print(``):
lprint(gu1):
print(``):
print(`The smallest number of kings possible is`, ldegree(gu,z), `and the number of ways of doing it is`):
print(``):
mu1:=coeff(gu,z,ldegree(gu,z)):
print(``):
print(mu1):
print(``):
print(`The largest number of kings possible is`, degree(gu,z), `and the number of ways of doing it is`):
print(``):
mu1:=coeff(gu,z,degree(gu,z)):
print(``):
lprint(mu1):
print(``):
end:

gen_states3 := proc(r,wid,polys):
local L,T,c3,i,j,k,m1,m2,p,cols:
L := []:
for k from 0 to 2*wid-2 do 
T:=combinat:-cartprod([seq([seq(i,i=0..1)],j=1..k*r)]):
while not T[finished] do 
c3:=T[nextvalue]() :
cols := [([-1$r])$(2*wid-2-k),seq([seq(c3[i],i=((j-1)*r+1)..(j*r))], j=1..k)]:
if is_good3(cols,polys) then  
	L:=[op(L),cols]:
fi:
od:
od:
L:
end:


is_good3 := proc(c, polys)
local co,i,j,flag,p,poly:
	co := 0:
	for i from 1 to nops(c) do
		if c[i] <> [-1$(nops(c[1]))] then
			break:
		fi:
		co := co + 1:
	od:
	for i from (co + 1) to nops(c) do:
	for j from 1 to nops(c[1]) do
		flag := false:
		for poly in polys do
		for p in poly do
			if c[i][j] = 1 and checkones3(c,[i,j],p,poly) = 1 then
				return(false):
			fi:
			if checkonesflex3(c,[i,j],p,poly) = 1 then
				flag := true:
			fi:
		od:
		od:
		if c[i][j]=0 and not(flag) then
			return(false):
		fi:
	od:
	od:
	return(true):
end:



 	
#checks whether everything except ref is 1s. always return -1 when out of bounds.
checkones3 := proc(c,xy,ref,poly)
local flag,p,newx,newy:
flag := false:
for p in poly do
	if p = ref then next: fi:
	newx := xy[1]+p[1]-ref[1]:
	newy := xy[2]+p[2]-ref[2]:
	if not(newx > 0 and newx <= nops(c) and newy > 0 and newy <= nops(c[1])) then
		return(-1):
	fi:

	if c[newx][newy] <> 1 then  
		flag := true:
	fi:
od:
if flag then return(0): fi:
return(1):
end:



checkonesflex3 := proc(c,xy,ref,poly)
local flag,p,newx,newy:
flag := false:
for p in poly do
	if p = ref then next: fi:
	newx := xy[1]+p[1]-ref[1]:
	newy := xy[2]+p[2]-ref[2]:
	if newy <= 0 or newy > nops(c[1]) then
		return(-1):
	fi:
	if newx <= 0 or newx > nops(c) then 
		next:
	fi:

	if c[newx][newy] <> 1 then  
		flag := true:
	fi:
od:
if flag then return(0): fi:
return(1):
end:













nextcolumns3 := proc(c,r,wid,polys)
local L,T,c3,flag,i,j,p,flag2,poly:
L := []:
T:=combinat:-cartprod([seq([seq(i,i=0..1)],j=1..r)]):
while not T[finished] do 
c3:=T[nextvalue]() :
if not(is_good3([op(2..nops(c),c),c3],polys)) then next: fi:
flag := false:		


for i from 1 to r do
	if c3[i] = 0 then next: fi:
	for poly in polys do
	for p in poly do
		if checkones3([op(c),c3],[2*wid-1,i],p,poly) = 1 then
			flag := true:
			break:	
		fi:
	
		
	od:
	if flag then break: fi:
	od:
	if flag then break: fi:
od:	
if flag then next: fi:
flag := false:
if c[wid][1] <> -1 then
for i from 1 to r do
	if c[wid][i] = 1 then next: fi:
	flag2 := false:
	for poly in polys do
	for p in poly do 
		if checkones3([op(c),c3],[wid,i],p,poly) = 1 then 
			flag2 := true:
			break:
		fi:
	od:
	if flag2 then break: fi:
	od:
	if flag2=false then  
		flag := true:
		break:
	fi:
	
	
od:
if flag then next: fi:
fi:
L := [op(L),[op(2..nops(c),c),c3]]:
od: 
L:
end:


wi := proc(poly)
local m1,m2,p:
m1 := -100:
m2 := 100:

for p in poly do
if p[1] < m2 then m2 := p[1]: fi:
if p[1] > m1 then m1 := p[1]: fi:
od:
m1 - m2+1:
end:


#printf("the relevant matrix entry has index M[1][-1]"):
TM3 := proc(r,polys,z):
local S,M,i,t,w,m1,m2,p,wid,poly:

wid := -100:
for poly in polys do
if wi(poly) > wid then
	wid := wi(poly):
fi:
od:
S := gen_states3(r,wid,polys):
M:=Matrix(nops(S)+1,nops(S)+1):
for i from 1 to nops(S) do

	for t in nextcolumns3(S[i],r,wid,polys) do
		#print(t,ListTools:-Search(t,S)):
		M[i,ListTools:-Search(t,S)] := weight2(t,z):
	od:
	if final_state3(S[i],r,wid,polys) then 
		M[i,nops(S)+1] := 1:
	fi:
	
od: 
M:
end:



final_state3 := proc(c,r,wid,polys)
local i,j,flag,p,poly:
for i from wid to nops(c) do
for j from 1 to r do

	if c[i][j] = 0 then
		flag := false:
		for poly in polys do
		for p in poly do
			if checkones3(c,[i,j],p,poly) = 1 then
				flag := true:
				break:
			fi:
		od:
		if flag then break: fi:
		od:
		if flag = false then return(false): fi:
	fi: 

od:
od:

return(true):
end:



weight2 := proc(c,z)
#1:
z^ListTools:-Occurrences(1,c[nops(c)]):
end:




GFrGeorge:=proc(r,S,z,x) local A,A2,d:

A:=TM3(r,S,z):
d:=(LinearAlgebra:-Dimension(A))[1]:
A2 := ((LinearAlgebra:-IdentityMatrix(d))-x*A)^(-1):
#print(A,A2,A2[1][d]):
normal(simplify((A2[1][d])/x)):




end: