######################################################################
##Bipartite: Save this file as Bipartite   To use it, stay in the    #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read Bipartite : <Enter>                           #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Dec. 29, 2006
 
print(`Created: Dec 29, 2006.`):
print(` This is Bipartite `):
print(`to automatically generate schemes, and then automatically`):
print(`use them, in order to count enumerating sequences`):
print(`for the number of regular bipartite graphs of degree k `):
print(`for any fixed inputted k`):
print(`It accompanies the paper `):
print(` "In How Many Ways Can n (Straight) Men and n (Straight) Women get Married,`):
print(` if Each Person Has Exactly k Spouses" `):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):
print(`and also available from his website`):
lprint(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
lprint(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg  .`):
 print(`For a list of the procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Mekadem, Rod, Oper, Seq, Sipur `):
else
ezra(args):
fi:

end:



ezra:=proc()

if args=NULL then
 print(`The main procedures are: OperR, OperRnice, SeqR, SipurA, SipurR `):


elif nops([args])=1 and op(1,[args])=Mekadem then
print(`Mekadem(P,A,deg): Given a polynomial P in A[1]. ..., A[k]`):
print(`extracts the coefficient of A[1]^deg[1]*...A[k]^deg[k]`):
print(`For example, try:`):
print(`Mekadem(3*A[1]*A[2]+A[1],A,[1,1]);`):



elif nops([args])=1 and op(1,[args])=Oper then
print(`Oper(a,k): Given  a symbol a, and a positive`):
print(`integer k, finds the recurrence operator with poly coeff.`):
print(`Oper(A[1], ..., A[k],a[1], ..., a[k])`):
print(`given in terms of a set {[poly,shif]}`):
print(`where [poly, shift] corresponds to`):
print(` poly*A[1]^shift[1]...A[k]^shift[k], `):
print(`such that the function`):
print(`F(n;a[1], ..., a[k]):=coeff of x1^k ... xn^k`):
print(`in e_k(x1, ..., xn)^a[k]*e_(k-1)(x1, ..., xn)^a[k-1]* ...`):
print(`*e_1(x1,...,xn)^a[1]`):
print(`safisfies F(n;a[1], ..., a[k])=Oper[F(n-1;a[1], ..., a[k])]`):
print(`For example, try:`):
print(`Oper(a,1);`):

elif nops([args])=1 and op(1,[args])=OperR then
print(`OperR(a,k,R): Given two symbols, A, and a, and a positive`):
print(`integer k, and a positive integer R,`):
print(`finds the recurrence operator with poly coeff.`):
print(`Oper(A[1], ..., A[k],a[1], ..., a[k])`):
print(`such that the function`):
print(`F_r(n;a[1], ..., a[k]):=coeff of x1^k ... xn^k`):
print(`in e_{k,r}(x1, ..., xn)^a[k]*e_{(k-1),r}(x1, ..., xn)^a[k-1]* ...`):
print(`*e_{1,r}(x1,...,xn)^a[1]`):
print(`where e_{i,r}(x1, ..., xn) is the weight-enumerator of`):
print(`of all sequences of integers between 0 and r that add`):
print(`up to i`):
print(`safisfies F_r(n;a[1], ..., a[k])=Oper[F_r(n-1;a[1], ..., a[k])]`):
print(`The representation of the operator is in the form as a set`):
print(`{[poly,shift]}`):
print(`For example, try:`):
print(`OperR(a,1,1);`):



elif nops([args])=1 and op(1,[args])=OperRnice then
print(`OperRnice(A,a,k,R): Given two symbols, A, and a, and a positive`):
print(`integer k, and a positive integer R,`):
print(`finds the recurrence operator with poly coeff.`):
print(`Oper(A[1], ..., A[k],a[1], ..., a[k])`):
print(`such that the function`):
print(`F_r(n;a[1], ..., a[k]):=coeff of x1^k ... xn^k`):
print(`in e_{k,r}(x1, ..., xn)^a[k]*e_{(k-1),r}(x1, ..., xn)^a[k-1]* ...`):
print(`*e_{1,r}(x1,...,xn)^a[1]`):
print(`where e_{i,r}(x1, ..., xn) is the weight-enumerator of`):
print(`of all sequences of integers between 0 and r that add`):
print(`up to i`):
print(`safisfies F_r(n;a[1], ..., a[k])=Oper[F_r(n-1;a[1], ..., a[k])]`):
print(`For example, try:`):
print(`OperRnice(A,a,1,1);`):

elif nops([args])=1 and op(1,[args])=Rod then
print(`Rod(Ope,a,k,n,vec): Given an operator `):
print(`Ope(a[1], ..., a[k],A[1], ..., A[k])`):
print(`and a vector vec, computes F(vec), where`):
print(`F(n;a[1], ..., a[k]):=Oper[F(n-1; a[1], ..., a[k])]`):
print(`For example, try:`):
print(`Rod(Oper(a,1),a,[5]);`):

elif nops([args])=1 and op(1,[args])=Seq then
print(`Seq(k,N): The first N terms of the sequence:`):
print(`number of (0,1) n by n matrices all whose rows and columns`):
print(` add up to k`):
print(`equivalently regular (n,n) bipartite graphs of degree k`):
print(`For example, try:`):
print(`Seq(3,10);`):

elif nops([args])=1 and op(1,[args])=SeqR then
print(`SeqR(k,N,R): The first N terms of the sequence:`):
print(`number of n by n matrices, with entries that `):
print(`nonnegative integers <=R`):
print(` all whose rows and columns add up to k.`):
print(`For example, try:`):
print(`SeqR(3,10,2);`):

elif nops([args])=1 and op(1,[args])=Sipur then
print(`Sipur(K,N): all the sequences for the number of`):
print(`n by n 0-1 matrices each of whose rows and columns`):
print(`add up to k, for k=1..K, and n=1..N`):
print(`For example, try:`):
print(`Sipur(4,20);`):

elif nops([args])=1 and op(1,[args])=SipurA then
print(`SipurA(K,N): all the sequences for the number of`):
print(`n by n matrices, whose entries are non-negative integers `):
print(`each of whose rows and columns`):
print(`add up to k, for k=1..K, and n=1..N`):
print(`For example, try:`):
print(`SipurA(4,20);`):

elif nops([args])=1 and op(1,[args])=SipurR then
print(`SipurR(K,N,R): all the sequences for the number of`):
print(`n by n matrices, whose entries are non-negative integers <=R,`):
print(`each of whose rows and columns`):
print(`add up to k, for k=R..K, and n=1..N`):
print(`For example, try:`):
print(`SipurR(4,20,2);`):

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


#Mekadem(P,A,deg): Given a polynomial P in A[1]. ..., A[k]
#extracts the coefficient of A[1]^deg[1]*...A[k]^deg[k]
#For example, try:
#Mekadem(3*A[1]*A[2]+A[1],A,[1,1]);
Mekadem:=proc(P,A,deg) local i,gu:
gu:=P:

for i from 1 to nops(deg) do
 gu:=coeff(gu,A[i],deg[i]):
od:
gu:
end:

#Oper(a,k): Given two symbols, A, and a, and a positive
#integer k, finds the recurrence operator with poly coeff.
#Oper(A[1], ..., A[k],a[1], ..., a[k])
#such that the function
#F(n;a[1], ..., a[k]):=coeff of x1^k ... xn^k
#in e_k(x1, ..., xn)^a[k]*e_(k-1)(x1, ..., xn)^a[k-1]* ...
#*e_1(x1,...,xn)^a[1]
#safisfies F(n;a[1], ..., a[k])=Oper[F(n-1;a[1], ..., a[k])]
#For example, try:
#Oper(a,1);
Oper:=proc(a,k) local gu,i,xn,Degs,s,j,A:
option remember:
gu:=mul((1+xn*A[i-1]/A[i])^a[i],i=1..k):
gu:=subs(A[0]=1,gu):
gu:=taylor(gu,xn=0,k+1):
gu:=expand(coeff(gu,xn,k)):

if type(gu,`+`)  then
 Degs:={seq([seq(degree(op(i,gu),A[j]),j=1..k)],i=1..nops(gu))}:
else
 Degs:={[seq(degree(gu,A[j]),j=1..k)]}:
fi:

{seq([Mekadem(gu,A,s),s], s in Degs)}:

end:



#Rod(Ope,a,k,n,vec): Given an operator 
#Ope(a[1], ..., a[k],A[1], ..., A[k])
#and a vector vec, computes F(vec), where
#F(n;a[1], ..., a[k]):=Oper[F(n-1; a[1], ..., a[k])]
#For example, try:
#Rod(Oper(a,1),a,[5]);
Rod:=proc(Ope,a,vec) local k,i,n,j,s:
option remember:
if min(op(vec))<0 then
  RETURN(0):
fi:
k:=nops(vec):
n:=add(i*vec[i],i=1..k):
if not type(n/k,integer) then
 RETURN(FAIL):
fi:
if n=0 then
 RETURN(1):
fi:

add(subs({seq(a[j]=vec[j],j=1..k)},s[1])*
Rod(Ope,a,[seq(vec[j]+s[2][j],j=1..k)]),s in Ope):

end:





#Seq(k,N): The first N terms of the sequence:
#number of (0,1) n by n matrices all whose rows and columns add up to k
#equivalently regular (n,n) bipartite graphs of degree k
#For example, try:
#Seq(3,10);
Seq:=proc(k,N) local A,a,i:
[seq(Rod(Oper(a,k),a,[0$(k-1),i]),i=1..N)]:
end:

#SeqR(k,N,R): The first N terms of the sequence:
#number of n by n matrices with entries from {0,1, ..., R}
#all of whose rows and columns add up to k
#For example, try:
#SeqR(3,10,1);
SeqR:=proc(k,N,R) local A,a,i:
[seq(Rod(OperR(a,k,R),a,[0$(k-1),i]),i=1..N)]:
end:


#Sipur(K,N): all the sequences for the number of
#n by n 0-1 matrices each of whose row and colum
#add to k, for k=1..K, and n=1..N
#For example, try:
#Sipur(4,20);
Sipur:=proc(K,N) local k:

for k from 1 to K do
print(`The sequence for the number of n by n matrices`):
print(`each of whose row and columns add up to `, k):
print(`for n between 1 and `, N, `is the following sequence`):
print(Seq(k,N)):
od:
end:

#SipurR(K,N,R): all the sequences for the number of
#n by n matrices whose entries are non-negative integers
#<=R, each of whose rows and columns
#add to k, for k=R..K, and n=1..N
#For example, try:
#SipurR(4,20,1);
SipurR:=proc(K,N,R) local k:

for k from R to K do
print(`The sequence for the number of n by n matrices,`):
 print(`whose entries are non-negative integers <= `, R , `and `):
print(`each of whose rows and columns add up to `, k):
print(`for n between 1 and `, N, `is the following sequence`):
print(SeqR(k,N,R)):
od:
end:




#OperR(a,k,R): Given two symbols, A, and a, and a positive
#integer k, and a positive integer R,
#finds the recurrence operator with poly coeff.
#Oper(A[1], ..., A[k],a[1], ..., a[k])
#such that the function
#F_r(n;a[1], ..., a[k]):=coeff of x1^k ... xn^k
#in e_{k,r}(x1, ..., xn)^a[k]*e_{(k-1),r}(x1, ..., xn)^a[k-1]* ...
#*e_{1,r}(x1,...,xn)^a[1]
#where e_{i,r}(x1, ..., xn) is the weight-enumerator of
#of all sequences of integers between 0 and r that add
#up to i
#safisfies F_r(n;a[1], ..., a[k])=Oper[F_r(n-1;a[1], ..., a[k])]
#For example, try:
#OperR(a,1,1);
OperR:=proc(a,k,R) local gu,i,xn,Degs,s,j,i1,A:
option remember:
gu:=mul((add(xn^i1*A[i-i1]/A[i],i1=0..R))^a[i],i=1..k):
gu:=subs({A[0]=1,seq(A[-i1]=0,i1=1..R-1)},gu):
gu:=taylor(gu,xn=0,k+1):
gu:=expand(coeff(gu,xn,k)):

if type(gu,`+`)  then
 Degs:={seq([seq(degree(op(i,gu),A[j]),j=1..k)],i=1..nops(gu))}:
else
 Degs:={[seq(degree(gu,A[j]),j=1..k)]}:
fi:

{seq([Mekadem(gu,A,s),s], s in Degs)}:

end:


#SipurA(K,N): all the sequences for the number of
#n by n matrices whose entries are non-negative integers
#each of whose rows and columns
#add to k, for k=1..K, and n=1..N
#For example, try:
#SipurA(4,20);
SipurA:=proc(K,N) local k:

for k from 1 to K do
print(`The sequence for the number of n by n matrices,`):
print(`with entries that are non-negative integers,`):
print(`each of whose rows and columns add up to `, k):
print(`for n between 1 and `, N, `is the following sequence`):
print(SeqR(k,N,k)):
od:
end:


#OperRnice(A,a,k,R): Given two symbols, A, and a, and a positive
#integer k, and a positive integer R,
#finds the recurrence operator with poly coeff.
#Oper(A[1], ..., A[k],a[1], ..., a[k])
#such that the function
#F_r(n;a[1], ..., a[k]):=coeff of x1^k ... xn^k
#in e_{k,r}(x1, ..., xn)^a[k]*e_{(k-1),r}(x1, ..., xn)^a[k-1]* ...
#*e_{1,r}(x1,...,xn)^a[1]
#where e_{i,r}(x1, ..., xn) is the weight-enumerator of
#of all sequences of integers between 0 and r that add
#up to i
#safisfies F_r(n;a[1], ..., a[k])=Oper[F_r(n-1;a[1], ..., a[k])]
#For example, try:
#OperRnice(A,a,1,1);
OperRnice:=proc(A,a,k,R) local gu,i,xn,i1:
option remember:
gu:=mul((add(xn^i1*A[i-i1]/A[i],i1=0..R))^a[i],i=1..k):
gu:=subs({A[0]=1,seq(A[-i1]=0,i1=1..R-1)},gu):
gu:=taylor(gu,xn=0,k+1):
gu:=expand(coeff(gu,xn,k)):

end: