######################################################################
##SD1row.txt: Save this file as  SD1row.txt                          #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read SD1row.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
 
print(`Created:  Fall  2023`):
print(` This is  SD1row.txt `):
print(`It is one of the packages that accompany the article `):
print(`Enumerating Seating Arrangments that Obey Social Distancing Restrictions of Any Kind`):
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 Story procedures type ezraS();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):


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 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(`---------------------------------------`):


with(combinat):


#with(LinearAlgebra):



ezraC:=proc()

if args=NULL then
 print(` The checking procedures are: `):

print(` CheckGFav, CheckGFmav `):

else
ezra(args):
fi:

end:

ezraPC:=proc()

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

 print(`  GFmavApc `):

else
ezra(args):
fi:

end:

ezraS:=proc()

if args=NULL then
 print(` The story procedures are: PaperS, PaperSpc `):

 print(``):

else
ezra(args):
fi:

end:



ezra1:=proc()
 
if args=NULL then
 print(` The supporting procedures are: ABC, AF, Alpha, AsyAv, ConjGR, ConjGR1, DensStat, v, GW , GFgr, GGav, GGmav, GWs, IsBad1, IsBad, IsMax, Lang, Lang1, LegalC,  MGW, MGWrsa, RSA1, RSAex, WDS  `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:  GFav, GFmav, GFmavP, LimAvDen, UNI, RSAsimu`):
 print(` `):



elif nargs=1 and args[1]=ABC then
print(`ABC(W,k): Given a set of words W finds all factors of length k. Try:`):
print(`W:=MGW(10,{[1,1,1]}); ABC(W,3);`):


elif nargs=1 and args[1]=AF then
print(`AF(w,k): all the factors of length k in the word w. Try:`):
print(`AF([1,0,0,1],2);`):

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,t): Given a rational generating  function f of x and t signifying weight enumerators of some statistic, find the`):
print(`(t for length, x for number of 1s) signifying weight enumerators of some statistic, find the number of 1s`):
print(`constant alpha such that the (average of the statistics)/n converges to it for items of size n.`):
print(`Try:`):
print(`AsyAv(1/(1-t*x-t^2),x,t);`):

elif nargs=1 and args[1]=CheckGFav then
print(`CheckGFav(S,K,N1,N2,N3): checks GFav(S,K,N1,N2,x,t) for the coefficients up t^N3. Try`):
print(`CheckGFav({[1,1]},3,10,14,16);`):


elif nargs=1 and args[1]=CheckGFmav then
print(`CheckGFmav(S,K,N1,N2,N3): checks GFmav(S,K,N1,N2,x,t) for the coefficients up t^N3. Try`):
print(`CheckGFmav({[1,1]},3,10,14,16);`):

elif nargs=1 and args[1]=ConjGR then
print(`ConjGR(W,K): A conjectured grammar for the set of words in {0,1} all of the same length. w with memory<=K `):
print(`it rertuns the memory-size, followed by the grammar ([k,G])`):
print(`Try:`):
print(`Try: W:=MGW(10,{[1,1]}); ConjGR(W,2); `):

elif nargs=1 and args[1]=ConjGR1 then
print(`ConjGR1(W,k): The conjectured abbreviated grammar for the set of words in {0,1} with segments of length k`):
print(`It returns`):
print(`the alphabet, the starters,  the enders, and the follower table, if there is one, or returns FAIL. Try`):
print(`Try: W:=GW(10,{[1,1]}); ConjGR1(W,2); `):



elif nargs=1 and args[1]=DensStat then
print(`DensStat(S,k,K): the density statistics for maximal  words in {0,1} avoiding S, up to the K-th moment. It also returns the actual distribution. `):
print(`Try:`):
print(`DensStat({[1,1]},1000,6);`):


elif nargs=1 and args[1]=GFav then
print(`GFav(S,K,N1,N2,x,t): the generating function for words avoiding S. K, N1,N2 are guessing parameters. Try:`):
print(`GFav({[1,1]},3,10,14,x,t);`):



elif nargs=1 and args[1]=GFgr then
print(`GFgr(G,x,t): Given a grammar finds the weight enumerator of the set of words generated by it with weight t^length*x^(#1). Try:`):
print(`G:=GGmav({[1,1]},2,10,14); GFgr(G,x,t);`):

elif nargs=1 and args[1]=GFmav then
print(`GFmav(S,K,N1,N2,x,t): the generating function for MAXIMAL words avoiding S. K, N1,N2 are guessing parameters. Try:`):
print(`GFmav({[1,1]},3,10,14,x,t);`):



elif nargs=1 and args[1]=GFmavApc then
print(`GFmavApc(i,x,t): The precomputed values of GFmavP({[1$i]},x,t);  In other words the bi-variate generating function in x and t whose coefficient of x^a*t^b `):
print(`is the number of maximal words in {0,1} of length a avoiding i consective 1's with b 1s. i must be from 1 to 7.`):
print(`Try:`):
print(`GFmavApc(6,x,t);`):

elif nargs=1 and args[1]=GGav then
print(`GGav(S,K,N1,N2): Given a set S of words in {0,1} and a positive integers K, and N, guesses a type-3 grammer for`):
print(`of memory <=K for the`):
print(`the set of binary words avoiding S as facors (consecutive substrings). It tests it up to words of length N. Try:`):
print(`GGav({[1,1]},2,8,12);`):


elif nargs=1 and args[1]=GGmav then
print(`GGmav(S,K,N1,N2): Given a set S of MAXIMAL words in {0,1} and a positive integers K, and N, guesses a type-3 grammer for`):
print(`of memory <=K for the`):
print(`the set of binary words avoiding S as facors (consecutive substrings). It tests it up to words of length N. Try:`):
print(`GGmav({[1,1]},2,8,12);`):

elif nargs=1 and args[1]=GFmavP then
print(`GFmavP(S,K,x,t): same as GFmav(S,6,12,15,z,t): (q.v.). Try: `):
print(`GFmavP({[1,1]},x,t);`):

elif nargs=1 and args[1]=GW then
print(`GW(k,S): The set of words in {0,1} of length k avoiding the mistakes in S. Try: `):
print(`GW(4,{[1,1]});`):


elif nargs=1 and args[1]=GWs then
print(`GWs(k,S): The set of words in {0,1} of length k avoiding the mistakes in S. The STUPID WAY. Only for checking. Try: `):
print(`GWs(4,{[1,1]});`):

elif nargs=1 and args[1]=IsBad then
print(`IsBad(w,S): Given a word w and a set of mistakes S finds out whether it has at least one of them. Try:`):
print(`IsBad([1,0,1,1,0],{[0,1,1]}); `);

elif nargs=1 and args[1]=IsBad1 then
print(`IsBad1(w,m): does the word w  contain m as factor?. Try:`):
print(`IsBad1([1,0,1,1,0],[0,1,1]); `);

elif nargs=1 and args[1]=IsMax then
print(`IsMax(w,S): inputs a word  w in the alphabet {0,1} and a set of "mistakes" S, such that w avoids the words in S as factors (consecutive subwords)`):
print(`and decides whether it is maximal, i.e. whether changing any of the 0's to 1 will create a mistake.`):
print(`Try:`):
print(`IsMax([1,1,0,1,1],{[1,1,1]});`):


elif nargs=1 and args[1]=Lang then
print(`Lang(G,n): Given a type-3 grammar   G=[k,[St,En,A,CT]] where k is the length of the letters (windows), St are the starters, En the enders,A, the alphabet, and  and n a positive integer n`):
print(`outputs all the  words of length n obeying this grammar. Try:`):
print(`W:=GW(12,{[1,1]}); G:=ConjGR(W,1); Lang(G,12);`):

elif nargs=1 and args[1]=Lang1 then
print(`Lang1(k,G,n): Given a type-3 grammar  with alphabet of length k, G=[St,En,A,CT] where St are the starters, En the enders,A, the alphabet, and  and n a positive integer n`):
print(`outputs all the prefixes of words of length n obeying this grammar. Try:`):
print(`W:=GW(10,{[1,1]}); G:=ConjGR1(W,2); Lang1(2,G,10);`):

elif nargs=1 and args[1]=LegalC then
print(`LegalC(k,G,w): Given a positive integer k (the length of the windows) and type-3 grammar, and a binary words, outputs the set of all its legal continuations`):
print(`Try:W:=GW(10,{[1,1]}); G:=ConjGR1(W,2); LegalC(2,G,[1,0,0,1]);`):


elif nargs=1 and args[1]=LimAvDen then
print(`LimAvDen(S): The limit of the average density of maximal words in {0,1} avoiding the patterns in S. Try:`):
print(`LimAvDen({[1,1]}); `):

elif nargs=1 and args[1]=MGW then
print(`MGW(k,S): The set of words in {0,1} of length k maximally avoiding the mistakes in S. `):
print(`MGW(4,{[1,1]});`):

elif nargs=1 and args[1]=MGWrsa then
print(`MGWrsa(k,S): MGW(k,S) via RSA, VERY slow. Only for checking. Try:`):
print(`MGWrsa(6,{[1,1]});`):



elif nargs=1 and args[1]=PaperS then
print(`PaperS(S,x,t,K1,K2,K3,K4): A paper about maximal S-avoiding maximal words in {0,1} comparing the uniform distribution of the density with`):
print(`RSA simulation K4 times. K1 is the number of terms for the OEIS. K2 is the length of the row taken both for the uniform distribution and`):
print(`the RSA simulation, K3 is the number of scaled moments to take. Here is t is the variable responsible for the length and x for the number of ones. Try:`):
print(`PaperS({[1,1]},x,t,30,100,6,1000):`):


elif nargs=1 and args[1]=PaperSpc then
print(`PaperSpc(b,x,t,K1,K2,K3,K4): Inputs a positive integer b, from 2 to 7, Outputs A paper about maximal [1$b] -avoiding maximal words in {0,1} comparing the uniform distribution of the density with`):
print(`RSA simulation K4 times. K1 is the number of terms for the OEIS. K2 is the length of the row taken both for the uniform distribution and`):
print(`the RSA simulation, K3 is the number of scaled moments to take. `):
print(`Here is t is the variable responsible for the length and x for the number of ones. Try:`):
print(`PaperSpc(2,x,t,30,100,6,1000):`):

elif nargs=1 and args[1]=RSA1 then
print(`RSA1(S,pi): Inputs a set of forbidden factors S, and a permutation pi, uses Random Sequential Adsorption to generate  the corresponding maximal S avoding word of length  N. try:`):
print(`RSA1({[1,1]},[1,2,3,4,5]);`):

elif nargs=1 and args[1]=RSAex then
print(`RSAex(k,S,z): The exact prob. generating function for words of length k avoding S in the variable z. Try:`):
print(`RSAex(7,{[1,1]},z);`):

elif nargs=1 and args[1]=RSAsimu then
print(`RSAsimu(k,S,K1,K2): The approximate RSA statistics for maximal words in {0,1} according to the number of ones, simulationg K1 times. It reurns the average, s.d. followed by the scaled moments moments up to K2. Try:`):
print(`RSAsimu(10,{[1,1]},1000,4);`):

elif nargs=1 and args[1]=UNI then
print(`UNI(k,S,z): The distribution of the number of ones according to the uniform distribution of S-avoiding maximal words. Try:`):
print(`UNI(100,{[1,1]},z);`):

elif nargs=1 and args[1]=WDS then
print(`WDS(k): all the words in {0,1} of length k. Try:`):
print(`WDS(5);`):


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





#Start from Histabrut
#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:

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

end:


#End from Histabrut
#WDS(k): all the words in {0,1} of length k. Try:
#WDS(5);
WDS:=proc(k) local gu,gu1:
option remember:
if k=0 then
 RETURN({[]}):
fi:
gu:=WDS(k-1):

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

#IsBad1(w,m): does the word w  contain m as factor?
IsBad1:=proc(w,m) local i:

if nops(w)<nops(m) then
 RETURN(false):
fi:

for i from 1 to nops(w)-nops(m)+1 do
 if [op(i..i+nops(m)-1,w)]=m then
   RETURN(true):
 fi:
od:
false:
end:


#IsBad1e(w,m): does the word w  end m as factor?
IsBad1e:=proc(w,m):

if nops(w)<nops(m) then
 RETURN(false):
fi:


 if [op(nops(w)-nops(m)+1..nops(w),w)]=m then
   RETURN(true):
 fi:

false:
end:

#IsBad(w,S): Given a word w and a set of mistakes S finds out whether it has at least one of them
IsBad:=proc(w,S) local s:
for s in S do
 if IsBad1(w,s) then
   RETURN(true):
 fi:
od:
false:
end:


#IsBadE(w,S): Given a word w and a set of mistakes S finds out whether it ends at least one of them
IsBadE:=proc(w,S) local s:
for s in S do
 if IsBad1e(w,s) then
   RETURN(true):
 fi:
od:
false:
end:

#GWs(k,S): The set of words in {0,1} of length k avoiding the mistakes in S. The stupid way. Try:
#GWs(4,{[1,1]});
GWs:=proc(k,S) local gu,lu,lu1:
option remember:
lu:=WDS(k):
gu:={}:
for lu1 in lu do
 if not IsBad(lu1,S) then
   gu:=gu union {lu1}:
 fi:
od:
gu:

end:


#GW(k,S): The set of words in {0,1} of length k avoiding the mistakes in S. 
#GW(4,{[1,1]});
GW:=proc(k,S) local gu,lu,lu1,w:
option remember:

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

lu:=GW(k-1,S):

gu:={}:

for lu1 in lu do
w:=[op(lu1),0]:

if not IsBadE(w,S) then
 gu:=gu  union {w}:
fi:

w:=[op(lu1),1]:
if not IsBadE(w,S) then
 gu:=gu  union {w}:
fi:

od:

gu:
end:

#IsMax(w,S): inputs a word  w in the alphabet {0,1} and a set of "mistakes" S, such that w avoids the words in S as factors (consecutive subwords)
#and decides whether it is maximal, i.e. whether changing any of the 0's to 1 will create a mistake.
#Try:
#IsMax([1,1,0,1,1],{[1,1,1]});
IsMax:=proc(w,S) local n,i,w1:

n:=nops(w):

if IsBad(w,S) then
 print(w, `is already bad`):
 RETURN(FAIL):
fi:

for i from 1 to n do
if w[i]=0 then
 w1:=[op(1..i-1,w),1,op(i+1..n,w)]:
 if not IsBad(w1,S) then
  RETURN(false):
 fi:
fi:
od:
true:
 
end:



#MGW(k,S): The set of words in {0,1} of length k maximally avoiding the mistakes in S. 
#MGW(4,{[1,1]});
MGW:=proc(k,S) local lu,gu,lu1:
 lu:=GW(k,S):
gu:={}:

 for lu1 in lu do
 if IsMax(lu1,S) then
   gu:=gu union {lu1}:
 fi:
od:
gu:
end:



#AF(w,k): all the factors of length k in the word w. Try:
#AF([1,0,0,1],2);
AF:=proc(w,k) local i:
{seq([op(i..i+k-1,w)],i=1..nops(w)-k+1)}:
end:

#ABC(W,k): Given a set of words W finds all factors of length k. Try:
#W:=MGW(10,{[1,1,1]}); ABC(W,3);
ABC:=proc(W,k) local w:
{seq(op(AF(w,k)), w in W)}:
end:




#ConjGR1(W,k): The conjectured abbreviated grammar for the set of words in {0,1} with segments of length k
#It returns
#the alphabet, the starters,  the enders, and the follower table, if there is one, or returns FAIL. Try
#Try: W:=GW(10,{[1,1]}); ConjGR(W,2); 
ConjGR1:=proc(W,k) local A,St,En,T,w,a,i:
A:=ABC(W,k):

St:={seq([op(1..k,w)], w in W)}:
En:={seq([op(nops(w)-k+1..nops(w),w)], w in W)}:


for a in A do
T[a]:={}:
od:

for w in W do
for i from 1 to nops(w)-k do
 T[[op(i..i+k-1,w)]]:=T[[op(i..i+k-1,w)]] union {w[i+k]}
od:
od:


[St,En,A, op(T)]:

end:

#LegalC(k,G,w): Given a positive integer k (the length of the windows) and type-3 grammar, and a binary words, outputs the set of all its legal continuations
#Try:W:=GW(10,{[1,1]}); G:=ConjGR1(W,2); LegalC(2,G,[1,0,0,1]);
LegalC:=proc(k,G,w) local A,w1,lu,lu1:
if nops(w)<k then
 RETURN(FAIL):
fi:
A:=G[3]:
w1:=[op(nops(w)-k+1..nops(w),w)]:
if not member(w1,A) then
  RETURN(FAIL):
fi:

lu:=G[4][w1]:

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

end:


#Lang1(k,G,n): Given a type-3 grammar  with alphabet of length k, G=[St,En,A,CT] where St are the starters, En the enders,A, the alphabet, and  and n a positive integer n
#outputs all the prefixes of words of length n obeying this grammar. Try:
#W:=GW(10,{[1,1]}); G:=ConjGR1(W,2); Lang1(2,G,10);
Lang1:=proc(k,G,n)  local gu,gu1:
option remember:
if n<k then
 RETURN(FAIL):
fi:
if n=k then
 RETURN(G[1]):
fi:
gu:=Lang1(k,G,n-1):

{seq(op(LegalC(k,G,gu1)),gu1 in gu)}:

end:




#Lang(G,n): Given a type-3 grammar  with alphabet of length k, G=[k,[St,En,A,CT]] where St are the starters, En the enders,A, the alphabet, and  and n a positive integer n
#outputs all the  words of length n obeying this grammar. Try:
#W:=GW(10,{[1,1]}); G:=ConjGR(W,2); Lang(G,10);
Lang:=proc(G,n)  local lu,lu1,S,gu,k:
k:=G[1]:

lu:=Lang1(k,G[2],n):

S:=G[2][2]:

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

gu:={}:

for lu1 in lu do
if member([op(n-k+1..n,lu1)],S) then
 gu:=gu union {lu1}:
fi:
od:
gu:
end:






#ConjGR(W,K): A conjectured grammar for the set of words in {0,1} all of same length . w with memory<=K 
#it rertuns the memory-size, followed by the grammar ([k,G])
#Try:
#Try: W:=MGW(10,{[1,1]}); ConjGR(W,3); 
ConjGR:=proc(W,K) local n,k,w,G1:

if W={} then
 RETURN(FAIL):
fi:

if nops({seq(nops(w), w in W)})<>1 then
 RETURN(FAIL):
fi:
n:=nops(W[1]):

for k from 1 to K do

G1:=ConjGR1(W,k):
if Lang([k,G1],n)=W then
 RETURN([k,G1]):
fi:
od:

FAIL:

end:



#GGav(S,K,N1,N2): Given a set S of words in {0,1} and a positive integers K, and N, guesses a type-3 grammer for
#of memory <=K for the
#the set of binary words avoiding S as facors (consecutive substrings). It tests it up to words of length N. Try:
#GGav({[1,1]},2,8,12);
GGav:=proc(S,K,N1,N2) local s,G,W,n,K1:

if  N1>=N2  then
  RETURN(FAIL):
fi:

W:=GW(N1,S):

K1:=max(seq(nops(s), s in S)):

G:=ConjGR(W,max(K,K1)):

if G<>FAIL then

 for n from G[1]+1 to N2 do
  if Lang(G,n)  <> GW(n,S) then
   print(G, `did not work out `):
   RETURN(FAIL):
 fi:
od:
fi:
G:
end:




#GGmav(S,K,N1,N2): Given a set S of words in {0,1} and a positive integers K, and N, guesses a type-3 grammer for
#of memory <=K for the
#the set of MAXIMAL binary words avoiding S as facors (consecutive substrings). It tests it up to words of length N. Try:
#GGmav({[1,1]},2,8,12);
GGmav:=proc(S,K,N1,N2) local s,G,W,n,K1:

if  N1>=N2  then
  RETURN(FAIL):
fi:

W:=MGW(N1,S):

K1:=max(seq(nops(s), s in S)):

G:=ConjGR(W,max(K,K1)):

if G<>FAIL then

 for n from G[1]+1 to N2 do
  if Lang(G,n)  <> MGW(n,S) then
   print(G, `did not work out `):
   RETURN(FAIL):
 fi:
od:
fi:
G:
end:




#GFgr(G,x,t): Given a grammar finds the weight enumerator of the set of words generated by it with weight t^length*x^(#1). Try:
#G:=GGmav({[1,1]},2,10,14); GFgr(G,x,t);
GFgr:=proc(G,x,t) local G1,eq,var,X,w,a,lu1,w1,eq1,lu,var1,i:

G1:=G[2]:

var:={seq(X[a], a in G1[3])}:





eq:={}:

for w in G1[3] do

 if member(w,G1[1]) then
  eq1:=X[w]-t^nops(w)*x^convert(w,`+`):
 else
 eq1:=X[w]:
fi:


lu:=G1[4][w]:


 for lu1 in lu do
  w1:=[op(2..nops(w),w),lu1]:
  eq1:=eq1 - t*x^w[1]*X[w1]:
 od:
eq:=eq union {eq1}:
od:
  
var1:=solve(eq,var):

if var1=NULL then
  RETURN(FAIL):
fi:

lu:=normal(add(subs(var1,X[w]), w in G1[1])):
lu:

end:


#GFav(S,K,N1,N2,x,t): the generating function for words avoiding S. K, N1,N2 are guessing parameters. Try:
#GFav({[1,1]},3,10,14,x,t);
GFav:=proc(S,K,N1,N2,x,t) local G,i,S1,w,lu,kat:
option remember:
G:=GGav(S,K,N1,N2); 
if G=FAIL then
 RETURN(FAIL):
fi:

kat:={seq(nops(w),w in G[2][3])}:
if nops(kat)<>1 then
 RETURN(FAIL):
fi:
kat:=kat[1]-1:


lu:=0:

for i from 0 to kat do
 S1:=GW(i,S):
 lu:=lu+t^i*add(x^convert(w,`+`),w in S1):
od:

normal(lu+GFgr(G,x,t));

end:



#GFmav(S,K,N1,N2,x,t): the generating function for words avoiding S. K, N1,N2 are guessing parameters. Try:
#GFmav({[1,1]},3,10,14,x,t);
GFmav:=proc(S,K,N1,N2,x,t) local G,kat,w,i,lu,S1:
option remember:
G:=GGmav(S,K,N1,N2); 
if G=FAIL then
 RETURN(FAIL):
fi:


kat:={seq(nops(w),w in G[2][3])}:
if nops(kat)<>1 then
 RETURN(FAIL):
fi:
kat:=kat[1]-1:


lu:=0:

for i from 0 to kat do
 S1:=MGW(i,S):
 lu:=lu+t^i*add(x^convert(w,`+`),w in S1):
od:

normal(lu+GFgr(G,x,t));
end:





#CheckGFav(S,K,N1,N2,N3): checks GFav(S,K,N1,N2,x,t) for the coefficients up t^N3. Try
#CheckGFav({[1,1]},3,10,14,16);
CheckGFav:=proc(S,K,N1,N2,N3) local gu,x,t,ku,ku1,i,gu1,S1,w:
gu:=GFav(S,K,N1,N2,x,t):
gu1:=taylor(gu,t=0,N3+3):
for i from 0 to N3 do
ku:=expand(coeff(gu1,t,i)):
S1:=GW(i,S):
ku1:=expand(add(x^convert(w,`+`), w in S1)):
if ku<>ku1 then
 RETURN(false):
fi:
od:
true:

end:


#CheckGFmav(S,K,N1,N2,N3): checks GFmav(S,K,N1,N2,x,t) for the coefficients up t^N3. Try
#CheckGFmav({[1,1]},3,10,14,16);
CheckGFmav:=proc(S,K,N1,N2,N3) local gu,x,t,ku,ku1,i,gu1,S1,w:
gu:=GFmav(S,K,N1,N2,x,t):
gu1:=taylor(gu,t=0,N3+3):
for i from 0 to N3 do
ku:=expand(coeff(gu1,t,i)):
S1:=MGW(i,S):
ku1:=expand(add(x^convert(w,`+`), w in S1)):
if ku<>ku1 then
print(`i is`, i):
print(ku,ku1):
 RETURN(false):
fi:
od:
true:

end:


#RSA1(pi,S): Inputs a permutation pi and set of forbidden factors S, and a positive integer N, uses Random Sequential Adsorption to find the corresponding member of MGW(nops(pi),S). Try
#RSA1({[1,1]},[2,1,3,4]);
RSA1:=proc(S,pi) local N,L1,L2,a,i:

N:=nops(pi):

L1:=[0$N]:

for i from 1 to nops(pi) do
 a:=pi[i]:
 L2:=[op(1..a-1,L1),1,op(a+1..N,L1)]:
 if not IsBad(L2,S) then
  L1:=L2:
 fi:
od:
L1:
end:

#MGWrsa(k,S): MGW(k,S) via RSA, VERY slow. Only for checking. Try:
#MGWrsa(6,{[1,1]});
MGWrsa:=proc(k,S) local gu,pi:
gu:=permute(k):
{seq(RSA1(S,pi), pi in gu)}:
end:


#RSAex(k,S,z): The exact prob. generating function for words of length k avoding S in the variable z. Try:
#RSAex(7,{[1,1]},z);

RSAex:=proc(k,S,z) local gu,f,pi:
gu:=permute(k):
f:=add(z^convert(RSA1(S,pi),`+`), pi in gu):
f/subs(z=1,f):
end:

#RSAsimu(k,S,K1,K2): The approximate RSA statistics for maximal words in {0,1} according to the number of ones, simulationg K1 times. Try:
#RSAsimu(100,{[1,1]},K1,K2);
RSAsimu:=proc(k,S,K1,K2)  local f, z,i,Dist,lu:
f:=add(z^convert(RSA1(S,randperm(k)),`+`),i=1..K1)/K1:
f:=evalf(f/subs(z=1,f)):

Dist:=[]:

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

lu:=evalf(Alpha(f,z,K2),20):
lu:=[lu[1]/k,lu[2]/k,op(3..K2,lu)]:
lu,Dist:

end:

#UNI(k,S,z): The distribution of the number of ones according to the uniform distribution of S-avoiding maximal words. Try:
#UNI(100,{[1,1]},z);
UNI:=proc(k,S,z) local f,t:
f:=GFmav(S,1,10,14,z,t):
f:=expand(coeff(taylor(f,t=0,k+1),t,k)):
f:=evalf(f/subs(z=1,f)):
end:




#GFmavP(S,K,x,t): GFmav(S,3,12,15,z,t):
#GFmavP({[1,1]},x,t);
GFmavP:=proc(S,x,t):GFmav(S,6,12,15,x,t):end:


#GFmavApc(i,x,t): The precomputed values of GFmavP({[1$i]},x,t);  In other words the bi-variate generating function in x and t whose coefficient of x^a*t^b 
#is the number of maximal words in {0,1} of length a avoiding i consective 1's with b 1s. i must be from 1 to 7.
#Try:
#GFmavApc(7,x,t);
GFmavApc:=proc(i,x,t) local L:

if not (i>=1 and i<=7) then
 print(i, `should have been between 1 and 6`):
 RETURN(FAIL):
fi:

L:=
[1/(1-t),
-(t^2*x+t*x+1)/(t^3*x+t^2*x-1), -(t^5*x^3-t^3*x^2-t^2*x^2+t^2*x-t*x-1)/(t^6*x^3-t^4*x^2-t^3*x^2-t^2*x+1), 
-(t^9*x^6+t^7*x^5-t^5*x^4+t^5*x^3-2*t^4*x^3-t^3*x^3+t^3*x^2-t^2*x^2-t*x-1)/(t^10*x^6+t^8*x^5-t^6*x^4-2*t^5*x^3-t^4*x^3-t^3*x^2+1), 
-(t^14*x^10-t^11*x^8-2*t^9*x^7+t^9*x^6-2*t^8*x^6+t^6*x^5-t^6*x^4+2*t^5*x^4+t^4*x^4-2*t^4*x^3+t^3*x^3-t^3*x^2+t^2*x^2+t*x+1)/(t^15*x^10-t^12*x^8-2*t^10*x^7-2*t^9*x^6+t^7*x^5+2*t^6*x^4+t^5*x^4+t^4*x^3+t^3*x^2-1), 
-(t^20*x^15+t^17*x^13-2*t^14*x^11+t^14*x^10-3*t^13*x^10-2*t^11*x^9+t^11*x^8-2*t^10*x^8+t^8*x^7-2*
t^8*x^6+2*t^7*x^6-2*t^7*x^5+3*t^6*x^5+t^5*x^5-2*t^5*x^4+t^4*x^4-t^4*x^3+t^3*x^3+t^2*x^2+t*x+1)/(t^21*x^15+t^18*x^13-2*t^15*x^11-3*t^14*x^10-2*t^12*x^9-2*t^11*x^8+t^9*x^7+2*t^8*x^6+3*t^7*x^5+t^6*x^5+t^5*x^4+t^4*x^3-1),
-(t^27*x^21-t^23*x^18-3*t^20*x^16+t^20*x^15-3*t^19*x^15+2*t^16*x^13-t^16*x^12+3
*t^15*x^12+3*t^13*x^11-3*t^13*x^10+4*t^12*x^10-2*t^12*x^9+3*t^11*x^9-t^9*x^8+2*
t^9*x^7-2*t^8*x^7+2*t^8*x^6-3*t^7*x^6-t^6*x^6+3*t^6*x^5-t^5*x^5+2*t^5*x^4-t^4*x
^4+t^4*x^3-t^3*x^3-t^2*x^2-t*x-1)/(t^28*x^21-t^24*x^18-3*t^21*x^16-3*t^20*x^15+
2*t^17*x^13+3*t^16*x^12+3*t^14*x^11+4*t^13*x^10+3*t^12*x^9-t^10*x^8-2*t^9*x^7-3
*t^8*x^6-t^7*x^6-t^6*x^5-t^5*x^4-t^4*x^3+1)
]:

L[i]:
end:

#AsyAv(f,x,t): Given a rational generating  function f of x and t (t for length, x for number of 1s) 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*t-x^2),x,t);
AsyAv:=proc(f,x,t) local a,f0,f1,A1,A2,lu,i,aluf,si:

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

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


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(x=1,diff(f,x))):



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

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

-A2/(A1*a):

end:


#LimAvDen(S): The limit of the average density of maximal words in {0,1} avoiding the patterns in S. Try:
#LimAvDen({[1,1]}:
LimAvDen:=proc(S) local f,x,t:

if nops(S)=1 and nops(S[1])>=2 and nops(S[1])<=7 and {op(S[1])}={1} then
f:=GFmavApc(nops(S[1]),x,t):
else
f:=GFmavP(S,x,t):
fi:


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

evalf(AsyAv(f,x,t)):
end:


#elif nargs=1 and args[1]=DensStat then
#print(`DensStat(S,k,K): the density statistics for maximal  words in {0,1} avoiding S, up to the K-th moment. It also returns the actual distribution. `):
#print(`Try:`):
#print(`DensStat({[1,1]},1000,6);`):





#DensStat(S,s,K): the density statistics for maximal configurations of words of length s, avoiding S up to the K-th moment.
#Try:
#DensStat({[1,1]},1000,6);
DensStat:=proc(S,s,K) local f,x,t,lu,Dist,i:
f:=GFmavP(S,x,t):
f:=coeff(taylor(f,t=0,s+1),t,s):
f:=f/subs(x=1,f):

Dist:=[]:

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

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

end:



#PaperS(S,x,t,K1,K2,K3,K4): A paper about maximal S-avoiding maximal words in {0,1} comparing the uniform distribution of the density with
#RSA simulation K4 times. K1 is the number of terms for the OEIS. K2 is the length of the row taken both for the uniform distribution and
#the RSA simulation, K3 is the number of scaled moments to take. Try:
#PaperS({[1,1]},x,t,30,100,6,10000):

PaperS:=proc(S,x,t,K1,K2,K3,K4) local A,f,c,f1,lu,mu,t0,avden,i:
t0:=time():
f:=GFmavP(S,x,t):

print(`Investigating Maximal`,S, `avoiding  words in {0,1}`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Let `, A(c)(x), ` be the weight-enumerator of MAXIMAL words in {0,1} of length c, avoiding, as factors (consecutive subwords) in`, S):
print(``):
print(`according to the weight x^NumberOfOnes (or equivalently, x^SumOfEntries)`):
print(``):
print(`Then `):
print(``):
print(Sum(A(c)(x)*t^c,c=0..infinity)=f):
print(``):
print(`and in Maple notation `):
print(``):
lprint(f):
print(``):
f1:=subs(x=1,f):
print(`It  follows that the generating function for the total number is`):
print(``):
print(Sum(A(c)(1)*t^c,c=0..infinity)=f1):
print(``):
lu:=[seq(coeff(taylor(f1,t=0,K1+1),t,i),i=1..K1)]:
print(`For the sake of the OEIS, the first`, K1, `terms are `):
print(``):
print(lu):
print(``):
print(`-------------------------------------------------`):
print(``):




avden:=evalf(AsyAv(f,x,t)):

print(`First let us note that the EXACT limiting average density as the length of the word goes to infinity is`, avden):
print(``):

print(``):

print(`Now let's specialize and look at the statistical distribution of words of length`, K2 , `and look at the staistics of the densitiy`):
print(``):

print(``):
lu:=DensStat(S,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`, K2):
print(``):

mu:=RSAsimu(K2,S,K4,K3):
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(`and The ratio of the exact s.d (under the uniform distribution) and the s.d of RSA simulations is`):
print(``):
print(mu[1][2]/lu[1][2]):
print(``):
print(`-------------------------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `seconds to generate most of the time was spent by the RSA simulation`):
print(``):


end:








#PaperSpc(b,x,t,K1,K2,K3,K4): For b from 2 to 6, A paper about maximal [1$A] -avoiding maximal words in {0,1} comparing the uniform distribution of the density with
#RSA simulation K4 times. K1 is the number of terms for the OEIS. K2 is the length of the row taken both for the uniform distribution and
#the RSA simulation, K3 is the number of scaled moments to take. Try:
#PaperSpc(2,x,t,30,100,6,10000):

PaperSpc:=proc(b,x,t,K1,K2,K3,K4) local A,f,c,f1,lu,mu,t0,avden,i,S:
if not (b>=2 and b<=7) then
 print(`The first argument should be between 2 and 6`):
 RETURN(FAIL):
fi:

t0:=time():

f:=GFmavApc(b,x,t);

print(`Investigating Maximal words  in {0,1} avoiding `,b, `consecutive ones `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Let `, A(c)(x), ` be the weight-enumerator of MAXIMAL words in {0,1} of length c, avoiding`, b, `consecutive ones `):
print(``):
print(`according to the weight x^NumberOfOnes (or equivalently, x^SumOfEntries)`):

print(``):
print(`Then `):
print(``):
print(Sum(A(c)(x)*t^c,c=0..infinity)=f):
print(``):
print(`and in Maple notation `):
print(``):
lprint(f):
print(``):
f1:=subs(x=1,f):
print(`It  follows that the generating function for the total number is`):
print(``):
print(Sum(A(c)(1)*t^c,c=0..infinity)=f1):
print(``):
lu:=[seq(coeff(taylor(f1,t=0,K1+1),t,i),i=1..K1)]:
print(`For the sake of the OEIS, the first`, K1, `terms are `):
print(``):
print(lu):

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



S:={[1$b]}:


avden:=evalf(AsyAv(f,x,t)):

print(`First let us note that the EXACT limiting average density as the length of the word goes to infinity is`, avden):
print(``):

print(``):

print(`Now let's specialize and look at the statistical distribution of words of length`, K2 , `and look at the staistics of the densitiy`):
print(``):

print(``):
lu:=DensStat(S,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`, K2):
print(``):

mu:=RSAsimu(K2,S,K4,K3):
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(`and The ratio of the exact s.d (under the uniform distribution) and the s.d of RSA simulations is`):
print(``):
print(mu[1][2]/lu[1][2]):
print(``):
print(`-------------------------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `seconds to generate most of the time was spent by the RSA simulation`):
print(``):


end: