######################################################################
##LatinRectangles                                                    #
#Save this file as LatinRectangles.   To use it, stay in the         #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type "read LatinRectangles:"  <Enter>                     #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: 
 
print(`Created:  Jan. 10, 2007`):
print(` This is LatinRectangles`):
print(`to count latin rectangles with a fixed width`):
print(`It is one of the two Maple packages that accompany the paper `):
print(` "In How Many Ways Can n (Straight) Men and n (Straight) `):
print(` Women Get Married, 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: Derangements `):
 print( ` Oper `):
else
ezra(args):
fi:

end:



ezra:=proc()

if args=NULL then
 print(`The main procedures are: LatinRecs , Seq `):
 print(`  `):


elif nops([args])=1 and op(1,[args])=Derangements then
print(`Derangements(n): The set of all the dernagments of {1, ... , n}`):
print(`by brute-force. For example, try:`):
print(`Derangements(5);`):
print(``):

elif nops([args])=1 and op(1,[args])=LatinRecs then
print(`LatinRecs(n,k): all the reduced Latin rectangles with k rows`):
print(`For example: try: LatinRecs(4,3);`):


elif nops([args])=1 and op(1,[args])=Oper then   
print(`Oper(a,k): The operator that links the number of`):
print(`k by n-1 generalized Latin rectangles to  k by n Latin rectangles`):
print(`For example, try: Oper(a,3);`):

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 k by n Latin rectangles (divided by n!)`):
print(` for 1<=n<=N `):
print(`For example, try: Seq(3,5);`):

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

with(combinat):
#Derangements(n): The set of all the dernagments of {1, ... , n}
#by brute-force. For example, try:
#Derangements(5);
Derangements:=proc(n)
local gu,mu,p,i:
option remember:
mu:=permute(n):
gu:={}:
 for p in mu do
  if {seq(evalb(p[i]=i),i=1..n)}={false} then
   gu:=gu union {p}:
  fi:
 od:

gu:
end:

#LatinRecs(n,k): all the reduced Latin rectangles with k rows
#For example: try: LatinRecs(4,3);
LatinRecs:=proc(n,k) local j,i,gu,mu,lu,g,m:
option remember:
if k=1 then 
RETURN ({[[seq(i,i=1..n)]]}):
fi:

gu:=LatinRecs(n,k-1):
mu:=Derangements(n):
lu:={}:

 for g in gu do
  for m in mu do
   if {seq(seq(evalb(g[j][i]=m[i]),i=1..n),j=1..k-1)}={false} then
    lu:=lu union {[op(g),m]}:
   fi:
  od:
 od:

lu:

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): The operator that links the number of
#k-1 by n generalized Latin rectangles to  k by n Latin rectangles
#For example, try: Oper(a,3);
Oper:=proc(a,k)
local gu,mu,x,T,i,E,b,Degs,A,j,s:
mu:=powerset({seq(i,i=1..k)}) minus {{}}:
mu:=
[{seq(i,i=1..k)},op(convert(mu minus {{seq(i,i=1..k)}},list))]:


gu:=mul((1+add(x[i]*E[T minus {i}]/E[T],i in T))^b[T], T in mu):

for i from 1 to k do
gu:=taylor(gu,x[i],2):
gu:=coeff(gu,x[i],1):
od:

gu:=expand(subs(E[{}]=1,gu)):
gu:=subs({seq(b[mu[i]]=a[i],i=1..nops(mu)),
seq(E[mu[i]]=A[i],i=1..nops(mu))},gu):
gu:=expand(gu):

if type(gu,`+`) then
 Degs:={seq([seq(degree(op(i,gu),A[j]),j=1..nops(mu))],i=1..nops(gu))}:
else
 Degs:={[seq(degree(gu,A[j]),j=1..nops(mu))]}:
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,j,s:
option remember:
if min(op(vec))<0 then
  RETURN(0):
fi:

if vec=[0$nops(vec)] then
 RETURN(1):
fi:


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

end:





#Seq(k,N): The first N terms of the sequence:
#number of k by n Latin rectangles (divided by n!)
#Seq(3,5);
Seq:=proc(k,N) local a,i:
[seq(Rod(Oper(a,k),a,[i,0$(2^k-2)])/i!,i=1..N)]:
end: