######################################################################
## Hardin.txt Save this file as  Hardin.txt  to use it,              #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read  `Hardin.txt` <Enter>                        #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
# and Robert Dougherty-Bliss, Darmouth College                       #
# and Natalya Ter-Saakov      , Rutgers University ,                 #
##  DoronZeil at gmail dot com                                       # 
######################################################################


print(`First Written: Aug. 20,  2024: tested for Maple 2020 `):
print(` This Version  : Nov. 8, 2024 (adding A172555)`):
print():
print(`This is Hardin.txt, A Maple package`):
print(`accompanying Robert Dougherty-Bliss, Christoph Koutschan, Natalya Ter-Saakov  and Doron Zeilberger's  article: `):
print(` The Symbolic and Numeric Computational Challenges of Counting 0-1 balanced matrices`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/mamarim/mamarimhtml/hardin.html `):
print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/Hardin.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezra1();`):
print():
print(`-------------------`):




ezraS:=proc()
if args=NULL then

print(`The Story procedures are: ` ):
print(``):

else
ezra(args):

fi:


end:



ezra1:=proc()

if args=NULL then

print(`The SUPPORTING procedures are: BiWords, ekn, GoodWords,  GoodBmatrices, GoodBmatrices1, GoodBwords, GoodMatrices, IsGood, IsGood1, LF, LFb, Milim,  Nu1, Nu11,  SeqBstupid, Wt1, Wt1E, Wt2, Wt2E  `):
print(``):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` Hardin.txt: A Maple package for  counting 0-1 matrices dear to R. H. Hardin`):

print(`The MAIN procedures are:  A172555, GF1t, GF1xt, GF2xt, Opek4, SeqB, SeqBf`):


print(``):


elif nargs=1 and args[1]=A172555 then
print(`A172555(K): the first K terms of OEIS sequence A172555, the number of 8 by 2*nx8 binary arrays with row-sums 4 and column-sums n, using the CONJECTURED (non-rigorous but absolutely certain)`):
print(`9th-order linear recurrence found by Christoph Koutschan using the Kauers-Koutschan amazing method guessing with little date.  Try:`):
print(`A172555(150);`):

elif nargs=1 and args[1]=BiWords then
print(`BiWords(n,k): The set of all words in {0,1} of length n with k 1's. Try:`):
print(`BiWords(4,2);`):

elif nargs=1 and args[1]=CheckGF1xt then
print(`CheckGF1xt(A,H,V,n,m): checks GF1xt(A,H,V,n,x,t) for the first m coefficient of powers of t. Do`):
print(`CheckGF1xt({1,2},{[1,1]},{[2,2]},4,4);`):


elif nargs=1 and args[1]=CheckGF2xt then
print(`CheckGF2xt(H,V,n,m): checks GF2xt(H,V,n,x,t) for the first m coefficient of powers of t. Do`):
print(`CheckGF2xt({[0,1,0],[1,0,1]},({[0,1,0],[1,0,1]},4,4);`):


elif nargs=1 and args[1]=ekn then
print(`ekn(x,k,n): The elementary symmetric function e_k(x[1],...,x[n]). Try:`):
print(`ekn(x,3,6);`):

elif nargs=1 and args[1]=GF1t then
print(`GF1t(A,H,V,n,t): inputs a (finite) alphabet A, a set of forbidden horizontal factors H, a set of forbiddedn vertical factors V, a pos. integer n`):
print(`and variables x,t. Outputs the rational function in t such that the coeff. of t^m in its Maclaurin expansion`):
print(`equals  the cardinality of GoodMatrices(A,H,V,m,n). Try:`):
print(`GF1t({0,1},{[0,1,0],[1,0,1]}$2,3,t);`):

elif nargs=1 and args[1]=GF1xt then
print(`GF1xt(A,H,V,n,x,t): inputs a (finite) alphabet A, a set of forbidden horizontal factors H, a set of forbiddedn vertical factors V, a pos. integer n`):
print(`and variables x,t. Outputs the rational function in {seq(x[a],a in A)} and t such that the coeff. of t^m in its Maclaurin expansion`):
print(`equals  the weight enumerator of the set GoodMatrices(A,H,V,m,n) according to the weight Wt1(M,x). Try: `):
print(`GF1xt({0,1},{[0,1,0],[1,0,1]}$2,3,x,t);`):


elif nargs=1 and args[1]=GF2xt then
print(`GF2xt(H,V,n,x,t): inputs a set of forbidden horizontal factors H, a set of forbiddedn vertical factors V, a pos. integer n`):
print(`and variables x,t. Outputs the rational function in x[1],..., x[n] and t such that the coeff. of t^m in its Maclaurin expansion`):
print(`equals  the weight enumerator of the set GoodBmatrices1(H,V,m,n) according to the weight Wt2(M,x). Try:`):
print(`GF2xt({[0,1,0],[1,0,1]}$2,4,x,t);`):

elif nargs=1 and args[1]=GoodBmatrices then
print(`GoodBmatrices(H,V,m,n): all the m by n 0-1 matrices whose rows AND columns are balanced and that avoid H and colums avoid V. Try:`):
print(`GoodBmatrices({[0,1,0],[1,0,1]},{[0,1,0],[1,0,1]},6,6);`):

elif nargs=1 and args[1]=GoodBmatrices1 then
print(`GoodBmatrices1(H,V,m,n): all the m by n 0-1 matrices whose rows are balanced and that avoid H and colums avoid V`):
print(`GoodBmatrices1({[0,1,0],[1,0,1]}$2,4,4);`):

elif nargs=1 and args[1]=GoodBwords then
print(`GoodBwords(H,n): all the words of length n in the alphabet {0,1} that avoid H, the clever way. Try`):
print(`GoodBwords({[0,1,0],[1,0,1]},6);`):

elif nargs=1 and args[1]=GoodMatrices then
print(`GoodMatrices(A,H,V,m,n): all the m by n matrices whose rows avoid H and colums avoid V. Try`):
print(`GoodMatrices({0,1},{[1,0,1],[0,1,0]},{[1,0,1],[0,1,0]} ,4,4);`):

elif nargs=1 and args[1]=GoodWords then
print(`GoodWords(A,H,n): all the words of length n in the alphabet H that avoid H.  Try`):
print(`GoodWords({0,1},{[0,1,0],[1,0,1]},7);`):

elif nargs=1 and args[1]=IsGood then
print(`IsGood(w,H): Given a word w and a set of forbidden subwords,H,  finds out whether the word w contains any of them as factors. Try:`):
print(`IsGood([1,0,1,0,1],{[1,0,1]});`):

elif nargs=1 and args[1]=IsGood1 then
print(`IsGood1(w,h): Given a word w and a word h finds out whether h is not a factor. Try:`):
print(`IsGood1([1,0,1,0,1],[1,0,1]);`):

elif nargs=1 and args[1]=LF then
print(`LF(A,H,V,n,M): Given  matrix with nops(M) rows and n columns that is good finds all the members`):
print(`of GoodWords(A,H,n) that are good continuations (i.e. not introducing a violation). Try:`):
print(`LF({1,2},{[1,1]},{[2,2]},4,[[1,2,1,2]]);`):

elif nargs=1 and args[1]=LFb then
print(`LFb(H,V,n,M): Given  a 0-1 matrix with nops(M) rows and n columns that is ballanced and good finds all the members`):
print(`of GoodBwords(A,H,n) that are good continuations (i.e. not introducing a violation). Try:`):
print(`LFb({[1,0,1],[0,1,0]}$2,4,[[1,1,0,0],[0,0,1,1]]);`):

elif nargs=1 and args[1]=Milim then
print(`Milim(A,n): all the words of length n in the alphabet A. Try:`):
prin(`Milim({0,1},4);`):


elif nargs=1 and args[1]=Nu1 then
print(`Nu1(n,k,A,v): The number of n by k 0-1 matrices whose rows are from the set of words of length k,A, that end with a in A and add-up to v. Try:`):
print(`Nu1(4,2,{[0,0],[1,1]},[1,1]);`):

elif nargs=1 and args[1]=Nu11 then
print(`Nu11(n,k,A,a,v): The number of n by k 0-1 matrices whose rows are from the set of words of length k,A, that end with a in A and add-up to v. Try:`):
print(`Nu11(4,2,{[0,0],[1,1]},[1,1],[1,1]);`):

elif nargs=1 and args[1]=Opek4 then
print(`Opek4(n,N): The linear recurrence operator of order 9 that (conjecturally) anniliates OEIS sequence A17255, the number of 8 by 2*nx8 binary arrays with row-sums 4 and column-sums n. Try:`):
print(`Opek4(n,N);`):

elif nargs=1 and args[1]=SeqB then
print(`SeqB(k,N): For a fixed k, the first N terms of the sequence number of (2k)*(2n) balanced matrices, by brute force . Try:`):
print(`SeqB(2,10);`):

elif nargs=1 and args[1]=SeqBstupid then
print(`Seqstupid(k,N): For a fixed k, the first N terms of the sequence number of (2k)*(2n) balanced matrices, VERY STUPIDLY. ONLY FOR CHECKING. . Try:`):
print(`SeqBstupid(2,10);`):

elif nargs=1 and args[1]=SeqBf then
print(`SeqBf(H,V,n,M):   the first M terms, m=1,m-2, ..., of the sequence enumerating  2*n by 2*m balanced 0-1 matrices`):
print(`avoiding the horizontal restrictions H and the vertical restrictions V. Try:`):
print(`SeqBf({[0,1,0],[1,0,1]}$2,2,2); `):

elif nargs=1 and args[1]=Wt1 then
print(`Wt1(M,x): The weight, accordint to the simple weight of the matrix M using x as an index-variable. Try:`):
print(`Wt1([[1,2],[2,3]],x);`):

elif nargs=1 and args[1]=Wt1E then
print(`Wt1E(A,H,V,m,n): The sequence, starting at m=1 of the weight enumerators of the set GoodMatrices(A,H,V,m1,n) for m1=1..m. Try:`):
print(`Wt1E({0,1},{[0,1,0],[1,0,1]},{[0,1,0],[1,0,1]},4,6,x);`):


elif nargs=1 and args[1]=Wt2 then
print(`Wt2(M,x): The weight, accordint to the sophisticated weight of the 0-1  matrix M using x[j] to keep track of the excess of 1s to 0s. Try:`):
print(`Wt2([[0,1,1],[1,0,1]],x);`):

elif nargs=1 and args[1]=Wt2E then
print(`Wt2E(H,V,m,n): The sequence, starting at m=1 of the weight enumerators according to Wt2(M,x), of the set GoodMatrices({0,1},H,V,m1,n) for m1=1..m. Only applicable to 0-1 matrices.Try:`):
print(`Wt2E({[0,1,0],[1,0,1]},{[0,1,0],[1,0,1]},4,6,x);`):

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

fi:


end:

ez:=proc(): 
rint(`GF1xt(A,H,V,n,x,t);`):
print(`CheckGF1xt(A,H,V,n,m)`):
print(`GF1xt({1,2},{[1,1]},{[2,2]},2,x,t);`):
print(`CheckGF1xt({1,2},{[1,1]},{[2,2]},4,4);`):
print(`GF2xt(H,V,n,x,t)`):
print(`GF2xt({[0,1,0],[1,0,1]}$2,4,x,t);`):
print(`CheckGF2xt({[0,1,0],[1,0,1]}$2,4,4);`):
end:

#IsGood1(w,h): Given a word w and a word h finds out whether h is not a factor. Try:
#IsGood1([1,0,1,0,1],[1,0,1]);
IsGood1:=proc(w,h) local i:
if nops(w)<nops(h) then
 RETURN(true):
fi:

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

#IsGood(w,H): Given a word w and a set of forbidden subwords,H,  finds out whether the word w contains any of them as factors. Try:
#IsGood([1,0,1,0,1],{[1,0,1]});
IsGood:=proc(w,H) local h:

for h in H do
 if not IsGood1(w,h) then
  RETURN(false):
 fi:
od:
true:
end:

#Milim(A,n): all the words of length n in the alphabet A. Try:
#Milim({0,1},4);
Milim:=proc(A,n) local gu,a,gu1:
option remember:
if n=0 then
 RETURN({[]}):
fi:
gu:=Milim(A,n-1):
{seq(seq([op(gu1),a], a in A),gu1 in gu)}:
end:

#GoodWordsS(A,H,n): all the words of length n in the alphabet H that avoid H, the stupid way. Try
#GoodWordsS({0,1},{[0,1,0],[1,0,1]},7);
GoodWordsS:=proc(A,H,n) local mu,mu1,gu:
mu:=Milim(A,n):
gu:={}:

for mu1 in mu do
 if IsGood(mu1,H) then
  gu:=gu union {mu1}:
 fi:
od:
gu:
end:


#GoodWords(A,H,n): all the words of length n in the alphabet H that avoid H, the clever way. Try
#GoodWords({0,1},{[0,1,0],[1,0,1]},7);
GoodWords:=proc(A,H,n) local ga, mu,mu1, h,a,hal,gu:
if H={} then
 RETURN(Milim(A,n)):
fi:
ga:=max(seq(nops(h),h in H))-1:
if n<=ga then
 RETURN(GoodWordsS(A,H,n)):
fi:

mu:=GoodWords(A,H,n-1):

gu:={}:

for mu1 in mu do
 for a in A do
 hal:=[op(mu1),a]:
 if IsGood([op(n-1-ga..n,hal)],H) then
  gu:=gu union {hal}:
 fi:
 od:
od:

gu:
end:



 

#GoodMatrices(A,H,V,m,n): all the m by n matrices whose rows avoid H and colums avoid V
GoodMatrices:=proc(A,H,V,m,n) local ga,v,lu,gu1,lu1,gu,i1,j1:
if m=0 then
 RETURN({[]}):
 fi:

ga:=max(seq(nops(v),v in V))-1:

lu:=GoodWords(A,H,n):



if m<=ga then
 gu:=GoodMatrices(A,H,V,m-1,n):
 RETURN({seq(seq([op(gu1),lu1],gu1 in gu),lu1 in lu)}):
fi:

 gu:=GoodMatrices(A,H,V,m-1,n):
gu:={seq(seq([op(gu1),lu1],gu1 in gu),lu1 in lu)}:

lu:={}:

for gu1 in gu do
 if {seq(IsGood([seq(gu1[i1][j1],i1=1..nops(gu1))],V),j1=1..n)}={true} then
  lu:=lu union {gu1}:
 fi:
od:
lu:
end:


#Wt1(M,x): The weight, accordint to the simple weight of the matrix M using x as an index-variable. Try:
#Wt1([[1,2],[2,3]],x);
Wt1:=proc(M,x) local i,j:
mul(mul(x[M[i][j]],j=1..nops(M[i])),i=1..nops(M)):
end:

#Wt1E(A,H,V,m,n): The sequence, starting at m=1 of the weight enumerators of the set GoodMatrices(A,H,V,m1,n) for m1=1..m. Try:
#Wt1E({0,1},{[0,1,0],[1,0,1]},{[0,1,0],[1,0,1]},4,6,x);
Wt1E:=proc(A,H,V,m,n,x) local gu,m1,lu,lu1,gu1:
gu:=[]:

for m1 from 1 to m do
 lu:=GoodMatrices(A,H,V,m1,n):
 gu1:=add(Wt1(lu1,x),lu1 in lu):
 gu:=[op(gu),gu1]:
od:
gu:
end:


#LF(A,H,V,n,M): Given  matrix with nops(M) rows and n columns that is good finds all the members
#of GoodWords(A,H,n) that are good continuations (i.e. not introducing a violation). Try:
#LF({1,2},{[1,1]},{[2,2]},4,[[1,2,1,2]]);
LF:=proc(A,H,V,n,M) local gu,lu,M1,lu1:
if not member(M,GoodMatrices(A,H,V,nops(M),n) ) then
  RETURN(FAIL):
fi:

lu:=GoodWords(A,H,n):

gu:={}:

for lu1 in lu do
M1:=[op(M),lu1]:
if member(M1,GoodMatrices(A,H,V,nops(M1),n)) then
 gu:=gu union {lu1}:
fi:
od:
gu:
end:

#GF1xt(A,H,V,n,x,t): inputs a (finite) alphabet A, a set of forbidden horizontal factors H, a set of forbiddedn vertical factors V, a pos. integer n
#and variables x,t. Outputs the rational function in {seq(x[a],a in A)} and t such that the coeff. of t^m in its Maclaurin expansion
#equals  the weight enumerator of the set GoodMatrices(A,H,V,m,n) according to the weight Wt1(M,x).
GF1xt:=proc(A,H,V,n,x,t) local ga,v,eq,var,X,gu,gu1,eq1,lu,lu1,w,lu1N,F,var1,F0,i1:
ga:=max(seq(nops(v),v in V))-1:
gu:=GoodMatrices(A,H,V,ga,n):
var:={seq(X[gu1],gu1 in gu)}:
eq:={}:

for gu1 in gu do
 eq1:=X[gu1]-Wt1(gu1,x)*t^nops(gu1):

 lu:=LF(A,H,V,n,gu1):

w:=Wt1([gu1[1]],x):
 for lu1 in lu do
  lu1N:=[op(2..nops(gu1),gu1),lu1]:
   eq1:=eq1-w*t*X[lu1N]:
 od:
eq:=eq union {eq1}:
od:


F0:=Wt1E(A,H,V,ga-1,n,x):


F0:=1+add(F0[i1]*t^i1,i1=1..nops(F0)):


F:=add(X[gu1],gu1 in gu):
var1:=solve(eq,var):

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


F:=subs(var1,F):




F:=normal(F+F0):

end:


#GF1t(A,H,V,n,t): inputs a (finite) alphabet A, a set of forbidden horizontal factors H, a set of forbiddedn vertical factors V, a pos. integer n
#and variables x,t. Outputs the rational function in t such that the coeff. of t^m in its Maclaurin expansion
#equals  the cardinality of GoodMatrices(A,H,V,m,n). Try:
#GF1t({[0,1,0],[1,0,1]}$2,3,t);
GF1t:=proc(A,H,V,n,t) local f,x,a:
f:=GF1xt(A,H,V,n,x,t):
normal(subs({seq(x[a]=1,a in A)},f)):
end:


#CheckGF1xt(A,H,V,n,m): checks GF1xt(A,H,V,n,x,t) for the first m coefficient of powers of t. Do
#CheckGF1xt({1,2},{[1,1]},{[2,2]},4,4);
CheckGF1xt:=proc(A,H,V,n,m) local f,gu,x,i,t:
f:=GF1xt(A,H,V,n,x,t):
if f=FAIL then
 RETURN(FAIL):
fi:
gu:=Wt1E(A,H,V,m,n,x):
f:=taylor(f,t=0,m+3):
evalb([seq(expand(coeff(f,t,i)),i=1..m)]=gu):
end:




###start more refined weight


#GoodBwords(H,n): all the balanced words of length n in the alphabet {0,1} that avoid H, the clever way. Try
#GoodBwords({[0,1,0],[1,0,1]},7);
GoodBwords:=proc(H,n) local mu,mu1,gu:
mu:=GoodWords({0,1},H,n):


gu:={}:

for mu1 in mu do
if  convert(subs(0=-1,mu1),`+`)=0 then
 gu:=gu union {mu1}:
fi:
od:

gu:
end:

#Wt2(M,x): The weight, according to location for of the 0-1 matrix M using x as an index-variable but taking care of location. 
#The power of x[j] is the excess of 1's to 0's in column j Try:
#Wt2([[1,2],[2,3]],x);
Wt2:=proc(M,x) local i,j,gu,m,n:
m:=nops(M):
n:=nops(M[1]):
gu:=mul(mul(x[M[i][j],j],j=1..n),i=1..m):
subs({seq(x[0,j]=1/x[j],j=1..n), seq(x[1,j]=x[j],j=1..n)},gu):
end:


#Wt2E(H,V,m,n): The sequence, starting at m=1 of the refined weight (for 0-1 matrices) enumerators of the set GoodMatrices({0,1},H,V,m1,n)  for m1=1..m. Try:
#Wt2E({0,1},{[0,1,0],[1,0,1]},{[0,1,0],[1,0,1]},4,6,x);
Wt2E:=proc(H,V,m,n,x) local gu,m1,lu,lu1,gu1:
gu:=[]:

for m1 from 1 to m do
 lu:=GoodMatrices({0,1},H,V,m1,n):
 gu1:=add(Wt2(lu1,x),lu1 in lu):
 gu:=[op(gu),gu1]:
od:
gu:
end:



#CheckGF2xt(H,V,n,m): checks GF2xt(H,V,n,x,t) for the first m coefficient of powers of t. Do
#CheckGF2xt({[1,0,1],[0,1,0]}$2,4,4);
CheckGF2xt:=proc(H,V,n,m) local f,gu,x,i,t:
f:=GF2xt(H,V,n,x,t):
if f=FAIL then
 RETURN(FAIL):
fi:
gu:=Wt2E(H,V,m,n,x):
f:=taylor(f,t=0,m+3):
evalb([seq(expand(coeff(f,t,i)),i=1..m)]=gu):
end:



#GF1xt(A,H,V,n,x,t): inputs a (finite) alphabet A, a set of forbidden horizontal factors H, a set of forbiddedn vertical factors V, a pos. integer n
#and variables x,t. Outputs the rational function in {seq(x[a],a in A)} and t such that the coeff. of t^m in its Maclaurin expansion
#equals  the weight enumerator of the set GoodMatrices(A,H,V,m,n) according to the weight Wt1(M,x).
GF1xt:=proc(A,H,V,n,x,t) local ga,v,eq,var,X,gu,gu1,eq1,lu,lu1,w,lu1N,F,var1,F0,i1:
ga:=max(seq(nops(v),v in V))-1:
gu:=GoodMatrices(A,H,V,ga,n):
var:={seq(X[gu1],gu1 in gu)}:
eq:={}:

for gu1 in gu do
 eq1:=X[gu1]-Wt1(gu1,x)*t^nops(gu1):

 lu:=LF(A,H,V,n,gu1):

w:=Wt1([gu1[1]],x):
 for lu1 in lu do
  lu1N:=[op(2..nops(gu1),gu1),lu1]:
   eq1:=eq1-w*t*X[lu1N]:
 od:
eq:=eq union {eq1}:
od:


F0:=Wt1E(A,H,V,ga-1,n,x):


F0:=1+add(F0[i1]*t^i1,i1=1..nops(F0)):


F:=add(X[gu1],gu1 in gu):
var1:=solve(eq,var):

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


F:=subs(var1,F):




F:=normal(F+F0):

end:


#GoodBmatrices1(H,V,m,n): all the m by n 0-1 matrices whose rows are balanced and that avoid H and colums avoid V
GoodBmatrices1:=proc(H,V,m,n) local ga,v,lu,gu1,lu1,gu,i1,j1:
if m=0 then
 RETURN({[]}):
 fi:

ga:=max(seq(nops(v),v in V))-1:

lu:=GoodBwords(H,n):



if m<=ga then
 gu:=GoodBmatrices1(H,V,m-1,n):
 RETURN({seq(seq([op(gu1),lu1],gu1 in gu),lu1 in lu)}):
fi:

 gu:=GoodBmatrices1(H,V,m-1,n):
gu:={seq(seq([op(gu1),lu1],gu1 in gu),lu1 in lu)}:

lu:={}:

for gu1 in gu do
 if {seq(IsGood([seq(gu1[i1][j1],i1=1..nops(gu1))],V),j1=1..n)}={true} then
  lu:=lu union {gu1}:
 fi:
od:
lu:
end:



#LFb(H,V,n,M): Given  a 0-1 matrix with nops(M) rows and n columns that is ballanced and good finds all the members
#of GoodBwords(A,H,n) that are good continuations (i.e. not introducing a violation). Try:
#LFb({[1,0,1],[0,1,0]}$2,[[1,0,1,0],[1,1,0,0]]);
LFb:=proc(H,V,n,M) local gu,lu,M1,lu1:
if not member(M,GoodBmatrices1(H,V,nops(M),n) ) then
  RETURN(FAIL):
fi:

lu:=GoodBwords(H,n):

gu:={}:

for lu1 in lu do
M1:=[op(M),lu1]:
if member(M1,GoodBmatrices1(H,V,nops(M1),n)) then
 gu:=gu union {lu1}:
fi:
od:
gu:
end:


#GF2xt(H,V,n,x,t): inputs a set of forbidden horizontal factors H, a set of forbiddedn vertical factors V, a pos. integer n
#and variables x,t. Outputs the rational function in {seq(x[a],a in A)} and t such that the coeff. of t^m in its Maclaurin expansion
#equals  the weight enumerator of the set GoodBmatrices({0,1},H,V,m,n) according to the weight Wt2(M,x).
GF2xt:=proc(H,V,n,x,t) local ga,v,eq,var,X,gu,gu1,eq1,lu,lu1,w,lu1N,F,var1,F0,i1:
ga:=max(seq(nops(v),v in V))-1:
gu:=GoodBmatrices1(H,V,ga,n):


var:={seq(X[gu1],gu1 in gu)}:
eq:={}:

for gu1 in gu do
 eq1:=X[gu1]-Wt2(gu1,x)*t^nops(gu1):

 lu:=LFb(H,V,n,gu1):

w:=Wt2([gu1[1]],x):
 for lu1 in lu do
  lu1N:=[op(2..nops(gu1),gu1),lu1]:
   eq1:=eq1-w*t*X[lu1N]:
 od:
eq:=eq union {eq1}:
od:


F0:=Wt2E(H,V,ga-1,n,x):


F0:=1+add(F0[i1]*t^i1,i1=1..nops(F0)):


F:=add(X[gu1],gu1 in gu):
var1:=solve(eq,var):

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


F:=subs(var1,F):




F:=normal(F+F0):

end:




#GoodBmatrices(H,V,m,n): all the m by n 0-1 matrices whose rows AND columns are balanced and that avoid H and colums avoid V. Try:
#GoodBmatrices({[0,1,0],[1,0,1]},{[0,1,0],[1,0,1]},6,6);
GoodBmatrices:=proc(H,V,m,n) local lu,lu1,gu,i1,j1:

if not m mod 2=0 then
 print(`m `,m, `should be even`):
 RETURN(FAIL):
fi:

lu:=GoodBmatrices1(H,V,m,n):

gu:={}:

for lu1 in lu do
  if {seq(add(lu1[i1][j1],i1=1..m),j1=1..n)}={m/2} then
   gu:=gu union {lu1}:
  fi:
od:

gu:

end:




#ekn(x,k,n): The elementary symmetric function e_k(x[1],...,x[n]). Try:
#ekn(x,3,6);
ekn:=proc(x,k,n)  local i:
option remember:
if k=1 then
 RETURN(add(x[i],i=1..n)):
fi:

if n=1 then
if k=1 then
 RETURN(x[1]):
else
 RETURN(0):
fi:
fi:

expand(ekn(x,k,n-1) +x[n]*ekn(x,k-1,n-1)):
end:


#SeqBold(k,N): For a fixed k, the first N terms of the sequence number of (2k)*(2n) balanced matrices, by brute force . Try:
#SeqBold(2,10);
SeqBold:=proc(k,N) local x,f,gu,n1,i1,f1:

f:=ekn(x,k,2*k):

gu:=[]:

for n1 from 1 to N do
 f1:=expand(f^(2*n1)):
  for i1 from 1 to 2*k do
   f1:=coeff(f1,x[i1],n1):
  od:
gu:=[op(gu),f1]:
od:
gu:
end:



#SeqBf(H,V,n,M):   the first M terms, m=1,m-2, ..., of the sequence enumerating  2*n by 2*m balanced 0-1 matrices
#avoiding the horizontal restrictions H and the vertical restrictions V. Try:
#SeqB({[0,1,0],[1,0,1]}$2,2,2)
SeqBf:=proc(H,V,n,M) local gu,x,t,L,i,gu1,m:

 gu:=GF2xt(H,V,2*n,x,t):
 gu:=taylor(gu,t=0,2*M+2):

L:=[]:

for m from 1 to M do
 gu1:=expand(coeff(gu,t,2*m)):
 for i from 1 to 2*n do
  gu1:=coeff(gu1,x[i],0):
 od:
L:=[op(L),gu1]:
od:
L:
end:






#BiWords(n,k): The set of all words in {0,1} of length n with k 1's. Try:
#BiWords(4,2);
BiWords:=proc(n,k) local Gu1,Gu2,w:
option remember:

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

Gu1:=BiWords(n-1,k):
Gu2:=BiWords(n-1,k-1):

{seq([op(w),0],w in Gu1), seq([op(w),1],w in Gu2)}:
end:




#Nu11(n,k,A,a,v): The number of n by k 0-1 matrices whose rows are from the set of words of length k,A, that end with a in A and add-up to v. Try:
#Nu11(4,2,{[0,0],[1,1]},[1,1],[1,1]);
Nu11:=proc(n,k,A,a,v) local w,b:
option remember:
if {seq(nops(w), w in A)}<>{k} then
  RETURN(FAIL):
fi:
if not member(a,A) then
 RETURN(FAIL):
fi:

if min(op(v))<0 then
  RETURN(0):
fi:

if n=1 then
if a=v then
  RETURN(1):
else
 RETURN(0):
fi:
fi:

add(Nu11(n-1,k,A,b,v-a), b in A):
end:


#Nu1(n,k,A,v): The number of n by k 0-1 matrices whose rows are from the set of words of length k, that add-up to v. Try:
#Nu11(4,2,{[0,0],[1,1]},[1,1]);
Nu1:=proc(n,k,A,v) local a:
option remember:
add(Nu11(n,k,A,a,v), a in A):
end:





#SeqBstupid(k,N): For a fixed k, the first N terms of the sequence number of (2k)*(2n) balanced matrices, cleverly . Try:
#SeqBstupid(2,10);
SeqBstupid:=proc(k,N) local A,n1:
A:=BiWords(2*k,k):
[seq(Nu1(2*n1,2*k,A,[n1$(2*k)]),n1=1..N)]:
end:





#SeqB(k,N): For a fixed k, the first N terms of the sequence number of (2k)*(2n) balanced matrices, by brute force . Try:
#SeqB(2,10);
SeqB:=proc(k,N) local x,f,gu,n1,i1,f1:

f:=ekn(x,k,2*k):

gu:=[]:

for n1 from 1 to N do
 f1:=f^(2*n1):
  for i1 from 1 to 2*k do
   f1:=coeff(f1,x[i1],n1):
  od:
gu:=[op(gu),f1]:
od:
gu:
end:

#Opek4(n,N): The linear recurrence operator of order 9 that (conjecturally) anniliates OEIS sequence A17255, the number of 8 by 2*nx8 binary arrays with row-sums 4 and column-sums n. Try:
#Opek4(n,N);
Opek4:=proc(n,N):
-496125000000000*(1+n)*(2+n)*(3+n)*(4+n)*(1+2*n)*(3+2*n)*(5+2*n)*(7+2*n)*(9+2*n)*(11+2*n)*(13+2*n)*(15+2*n)*(17+2*n)*(1708369933613326663680*n^23+255401305075192336220160*n^22+18222343765549559379394560*n^21+825585244086079515474788352*n^20+26659931622963281030074073088*n^
19+652880304667928117936189079552*n^18+12595230124506862351125236154368*n^17+196275127066044457020180222377984*n^16+2513375241520987501607113962110976*n^15+26762020178485189208591017270616064*n^14+238840293875155345925965469462368256*n^13+
1795411794455612692492065984729558016*n^12+11394951871120481737104432307024736064*n^11+61059117549298299205982739568911311776*n^10+275581033499618888334332524246956934944*n^9+1042470173007586769022365064935495981808*n^8+3278923215358322490867617390292630707992*n^7+
8474441209752165916584686593489897696448*n^6+17691002481605703461451617209981375383622*n^5+29091257191392529216620930498704807839965*n^4+36275761013903607492425291639403051165825*n^3+32238337642618415900882675854119580117125*n^2+18191320641499286998161262039631309285625*n+
4898365292987617609676333150712805443750)+22500000000*(2+n)*(3+n)*(4+n)*(3+2*n)*(5+2*n)*(7+2*n)*(9+2*n)*(11+2*n)*(13+2*n)*(15+2*n)*(17+2*n)*(23286790565083255752622080*n^25+3597809142305363013780111360*n^24+265947489714600097696123453440*n^23+
12518192548307278073543987822592*n^22+421288284055214512389176327405568*n^21+10789754377512053953992942194786304*n^20+218550306104887732632392482955198464*n^19+3591824497932180221162044781642711040*n^18+48755070909494577497350222255706128384*n^17+
553512625909260100450743990756884094976*n^16+5302626919716634381775701695712989605888*n^15+43125712935482569224575468596160023036928*n^14+298876158689185025610727511075068857960640*n^13+1768141174828344914341029658546946681910048*n^12+
8928252680488018747786285958461370805846752*n^11+38404771595398636290659238721615901656915856*n^10+140168597190672045590905274764441695077758696*n^9+431362321021059782746945619728502858226275560*n^8+1109202646434225274017994655405725656710433226*n^7+
2352868537001625347533816777095357631728572409*n^6+4043777096617207712823952138388524484230272507*n^5+5487372816107522596170768453441872959294963752*n^4+5656099613769357129030370130386094712783664980*n^3+4159530249102375554089857103375055810568010425*n^2+
1943329630343148480909354041811724772295661875*n+433286032074174796514012692442911970512301250)*N-500000000*(3+n)*(4+n)*(5+2*n)*(7+2*n)*(9+2*n)*(11+2*n)*(13+2*n)*(15+2*n)*(17+2*n)*(85995925718227637596323840*n^27+13888342003493763471806300160*n^26+
1076269485254494305234242764800*n^25+53280341772456945587622915342336*n^24+1892459266928007427920327198375936*n^23+51351665225985281457443899703820288*n^22+1106721861154992634290250493443702784*n^21+19444413125915414233017836516063707136*n^20+
283643240234436472823145722241080328192*n^19+3480994039135509457081902214193218764800*n^18+36287507156652026494960774227357829586944*n^17+323545650337044214302381647478402081642496*n^16+2479227822233408888078938902293854590661248*n^15+
16375908894496278175222297537711170935687232*n^14+93373658678186737039422157129020872741680384*n^13+459570685209165389574283108635360379321134816*n^12+1949509842293969790211135621062918803814263504*n^11+7105528477831712091828556470089814945735575984*n^10+
22143490514872634047879888007965152008273004064*n^9+58591193986328462305260891961389685019179357006*n^8+130360381994295938164307578921642833634214698460*n^7+240666222054644916081900418235523333451282193761*n^6+361954942514678749910376677047346405222558047524*n^5+
432027359098806536257782836660287623106109244210*n^4+393607625714951350311287525698807066566307053570*n^3+257053796012665600195609689721569091535999123025*n^2+107133298007006278241600505542568004815482341500*n+21402278524083561049599995803378759509036037500)*N^2-10000000*(
4+n)*(7+2*n)*(9+2*n)*(11+2*n)*(13+2*n)*(15+2*n)*(17+2*n)*(262258701988716229448171520*n^29+44715108689076117120913244160*n^28+3668703205864601929224636334080*n^27+192875914340856511052486402899968*n^26+7299544405839804767945733455216640*n^25+
211807114254879654312799797423636480*n^24+4900435371804469588713322633546956800*n^23+92821159761359016162491962263237165056*n^22+1466559751138970915319903980843200249856*n^21+19594015457226928079958020961794785755136*n^20+223623054440090975709124262676072394285056*n^19+
2196577001480509479875537257949520853456896*n^18+18672488450019054334765274518614045142822016*n^17+137897951093808390633995994025266417406248896*n^16+886901827250605121602012818616206467769276416*n^15+4973562265572502613909642566448385205683677536*n^14+
24317818877254235589073270234669578447532088016*n^13+103541908946332621246221974372798849349573244896*n^12+382967237971613498217870145290648946638457027616*n^11+1225656913541855645334956119857888333067130486786*n^10+3375444680620227963760389109023155136928624029176*n^9+
7938948888326572364136808136097106639831340160991*n^8+15785169763932060794766418870128526458484267666566*n^7+26172342121537160549989099114590339772726948186023*n^6+35514410632347855704573407331783063377673799067026*n^5+38410586974831070242819681787714379159155490925191*n^4
+31837904680271210221948188114852600093151355119716*n^3+18988756022307766893492591781408519605558366436485*n^2+7253518484549736587612215117433437584422260009500*n+1332648850866969781650466919612784681375563155500)*N^3-200000*(9+2*n)*(11+2*n)*(13+2*n)*(15+2*n)*(17+2*n)*(
2752626430863875106794373120*n^31+499601697201793331883178721280*n^30+43760685557381616964928305889280*n^29+2463662764563037891897659834236928*n^28+100172543110915776477044356085710848*n^27+3133708065014176560444885365460566016*n^26+78458713100336395643436286934656483328*n
^25+1614678564438710649355389663423015223296*n^24+27838810426560219418274102574281606479872*n^23+407772451917604785112177265828780604809216*n^22+5128132419841019732660052631221880725041152*n^21+55813549021064580192335298397172514642076672*n^20+
528902673633061763083170896015839644831117248*n^19+4383343466590264290822298600432055635327072480*n^18+31871331993802595079401376182774994900036603168*n^17+203721842475199970097720484523940707459515318096*n^16+1145897049212681583233066271149240949087360964904*n^15+
5671828971752077998173985719588380635143018164624*n^14+24679707274024235050155282928657622031264677832934*n^13+94214190035021759772679786312997848535182806258395*n^12+314541885094928078004567092587367396698766026037145*n^11+
914295033294423593837720059861275029666434916746038*n^10+2299963957703050805992469855535590820529662412475475*n^9+4967241787341303361104239989115851737677844082561036*n^8+9113713182130918778452680618018592805171822279926993*n^7+
14007994953398262662741651126482561733145911183074596*n^6+17696908746509679109765547801181303794758429777099413*n^5+17892371633195428996233844105126203219539932540054827*n^4+13917292904523702142362651131270964535410107730813120*n^3+
7817783228239131602382248249375498756951195436352500*n^2+2822402219978909537358814388144095747887345747516000*n+491705576267742740668687287902657299818290786400000)*N^4+800*(5+n)*(11+2*n)*(13+2*n)*(15+2*n)*(17+2*n)*(56065272522732619155363594240*n^31+
10344042780444168234164583137280*n^30+921226848300012768572606597038080*n^29+52744134225275882827255197714087936*n^28+2181445155830546677118085139662372864*n^27+69429940301784979794801162034232426496*n^26+1768925770193478095017003703849234202624*n^25+
37052676081110943650811464058654878269440*n^24+650326375309449917376842634107300940759040*n^23+9698955348543166512998220785217964553461760*n^22+124213497622663285088419390705876287482941440*n^21+1376970202417390494232795771788518547126246400*n^20+
13292463860007607593208451100294945878480686400*n^19+112239574000557864589264267154562308846116975200*n^18+831598721939357671416809565093067186723506869600*n^17+5417299212553518174854989710194108315456570832080*n^16+31058259334229216052671851027084367288628997462520*n^15+
156708150832012538977730084612418422896820689140680*n^14+695171561790705271432667410289469050383367122897370*n^13+2705791143956678704367533084442447792246048394535295*n^12+9211292875584371663876242081207330367848737532399515*n^11+
27304015934951291724048210068369056759396519477232840*n^10+70046887236907125631073472501072787912942197611459725*n^9+154289359055268421127765681677713019117402568584682754*n^8+288728559381544710851895828349520883168621563599468871*n^7+
452647942398049248581796568334584127447401547628798994*n^6+583289669106687506809820620444850739191618647669969911*n^5+601540481498485553214734040327408384774529947106377345*n^4+477271248581648477923382004154390794580437322315837800*n^3+
273468255129560862533004913229282331143525191327495500*n^2+100704769427745716161263054919066655324841615547140000*n+17895137437841140825745636970614851539889297579200000)*N^5+16*(5+n)*(6+n)*(13+2*n)*(15+2*n)*(17+2*n)*(7196990105667417528871157760*n^31+
1349435644812640786663342080000*n^30+122146385542362293023755542200320*n^29+7108628882028533911636400494608384*n^28+298880256264990204333131598557872128*n^27+9671250943921781791490854649945653248*n^26+250534806684411701692101354292614529024*n^25+
5336275369548235127826091118843285078016*n^24+95245760149039526908108812377093666406400*n^23+1444668153545436161610510184937244953231360*n^22+18817924664449648911496374543342178118082560*n^21+212185578704098379021134878827072234499696640*n^20+
2083585726653990457353034777520937354453416320*n^19+17897408841432147319372293289907265200417433280*n^18+134901688163022700165742302488415707129848353920*n^17+894055273173932773313508894468870595244813032800*n^16+5214978398983584983250940291779994993246327757360*n^15+
26771579280716628972321368662429849544260377953520*n^14+120834723425746121124702161449860259193924694482840*n^13+478540668648075209571342423630852651466177261720090*n^12+1657580151215781745901566538110074425990294007336620*n^11+
4999325011911335893749740867402580411054410366038315*n^10+13049775203003340786829611016747499182332072623901530*n^9+29246561157980405097807474860990372918003597935691351*n^8+55685906229755149692841467113716129565716446343194062*n^7+
88822211704470699203584286253746111365539085757839027*n^6+116449117759209746565804497340843550493082593721721056*n^5+122177048227950589386315422528598393503797118250506469*n^4+98614633946508581199652924403692591875974591975228100*n^3+
57479009553418393424253587273904850677419208241477500*n^2+21530297609932645266588657702813781082132306946360000*n+3891351686564358063191469660858623417943135630960000)*N^6-8*(5+n)*(6+n)*(7+n)*(15+2*n)*(17+2*n)*(2189545992374989025118781440*n^31+
417108511547435409285127864320*n^30+38355025667357192636849392189440*n^29+2267352840872690548944664082251776*n^28+96819917282297302158404596911833088*n^27+3181432791970598950450994456494080000*n^26+83678722732943532016515326383169208320*n^25+
1809348448652326558487085457399340728320*n^24+32778774774772269859484446602783226036224*n^23+504545483047281698935787036893132292243456*n^22+6668179021745548175927636208108316153470976*n^21+76272675186176555901547808994622892461365248*n^20+
759613369550247298718090568627478404675950720*n^19+6616165277921750021160075927920368254562900416*n^18+50555896761990089580031308858960358231388868736*n^17+339591563273384807107953904589714251542515682528*n^16+2007149189763090610325749317167246003608416682640*n^15+
10438283927376161615399593888019433066696220910976*n^14+47716126430148518169345370060744640569257601151736*n^13+191336615086316865305790791064634504637947931046538*n^12+670879320830714945184414928946372469398092743854060*n^11+
2047642948925248640605968356994916729776088338765421*n^10+5407530505638185971250710984259046766744118670515526*n^9+12257500170546340504788701326630559718835841945743717*n^8+23598250507244583321677670140114040743655066122829602*n^7+
38048560370258750244494491040831413790654466500584421*n^6+50409037427826176939887158425999078009232343504912416*n^5+53430567815334907072517070636823241262928904257674123*n^4+43555313223924190572522423021252330390861305831113476*n^3+
25631843445854552902279750997858412209618156492254740*n^2+9690860917565823091556634053697700420973479967361600*n+1767362308993618469321480623903583635034472036144000)*N^7+2*(5+n)*(6+n)*(7+n)*(8+n)^3*(17+2*n)*(22398438199604325887508480*n^29+3975722780429767845032755200*n^
28+339690492840187044804710891520*n^27+18603409376580638976672584957952*n^26+733649225995577952763276690980864*n^25+22189400779795700463916879647080448*n^24+535286044118169199520991392421380096*n^23+10574945171435588058138730228603224064*n^22+
174318785633245263331852062649517522944*n^21+2430604803850992277552370633077150474240*n^20+28959222262428697503733095780966269448192*n^19+297049515675305722778680833370339212555264*n^18+2637722622903340676154712351306394415591104*n^17+
20354542960562825184809545430884477789295200*n^16+136831937745091631528274766540034373584674912*n^15+802267598042943030377706156870124986015824464*n^14+4102485049949791899845918522125402535644559784*n^13+18274253945561423033814500713106918892323753280*n^12+
70732011003348632853947112905205964760739880802*n^11+236964494430635662243068683861444375313214538779*n^10+683334862462281505196293598108083457702810418619*n^9+1683375853820331731022446217534133320427688135815*n^8+3506783911711734968687217777553957587813718457667*n^7+
6093528657551925134728980155822720469989194433337*n^6+8668037255522269976531937101616670215747314238105*n^5+9830568813989012954695019483223113374730200271817*n^4+8546797943251860211583607063375016167113743061811*n^3+5348157926480993451925976008113224780089144936140*n^2+
2143977249971319532946730180933642849474255905100*n+413488176150319412586890844611086224835553244000)*N^8-9*(5+n)*(6+n)*(7+n)*(8+n)^3*(9+n)^7*(1708369933613326663680*n^23+216108796602085822955520*n^22+13035732647099499628462080*n^21+498888203329478996652982272*n^20+
13596728537953448000079003648*n^19+280777858606528849440941015040*n^18+4563573214041399099165613752320*n^17+59861219294009027529267916767232*n^16+644654295874091892639003131854848*n^15+5767423429843304321241914825891840*n^14+43208030699189337528136840512245760*n^13+
272403057743763525537209617959222272*n^12+1448589264911577756626446150131138368*n^11+6497689096131722570692369446389199840*n^10+24525693851840842091943891464123460000*n^9+77514784387741037883769443039027850672*n^8+203509237081451490395762679019778104728*n^7+
438608431094715348469097797970402022040*n^6+762800151857017231327762036260844513790*n^5+1043978045152125077047565766566495579927*n^4+1082413567759839362750214677564636470353*n^3+799050142492970208675152589943906208220*n^2+374168579040629551760577258256728027300*n+
83528294022396251739836927583932052000)*N^9
:
end:




#From FindRec.txt
Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:

#SeqFromRec(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

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

gu:=Ini:

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

gu:

end:

#End From FindRec.txt

#A172555(K): the first K terms of OEIS sequence A172555, the number of 8 by 2*nx8 binary arrays with row-sums 4 and column-sums n, using the CONJECTURED (non-rigorous but absolutely certain)
#9th-order linear recurrence found by Christoph Koutschan using the Kauers-Koutschan amazing method guessing with little date.  Try:
#A172555(150);
A172555:=proc(K) local Ini,ope,n,N:
ope:=Opek4(n,N):
Ini:=[70, 44730, 60871300, 116963796250, 273957842462220, 729833528645228700, 2128370421607383353400, 6637298299496411164868250, 21796558191507153941744357500]:
SeqFromRec(ope,n,N,Ini,K):
end: