######################################################################
##AMITAI: Save this file as AMITAI   To use it, stay in the          #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read AMITAI : <Enter>                              #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Nov. 2006
with(combinat): 
print(`Created: Nov. 2006`):
print(` This is AMITAI `):
print(`a Maple package to study the number of`):
print(`Young-tableaux of bounded height.`):
print(`It accompanies the paper `):
print(` Proof of a Conjecture of Amitai Regev about Three-Rowed`):
print(`Young Tableaux (and much more!) `):
print(`by Shalosh B. Ekhad and Doron Zeilberger,`):
print(`also available from the authors' websites`):
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: ApplyOper, Axkn, Bxkn, CF1,`):
 print(` GenOper, GuessOper1, GW1, HafelOper, nuYTwith, Pars1, Pars`):
 print( ` ParsWith, , RegevStyle3YTwith1, RGF, RGFt, SYTkn, YTwith `):
else
ezra(args):
fi:

end:



ezra:=proc()

if args=NULL then
 print(`The main procedures are: AR, ARa, ARaS, ARaf, CF,`):
 print(`  GuessRegev3D, GuessRegev4D, GuessRegev4Dnice, `):
 print(`GuessAndProveRegev3D, `):
 print(` GuessAndProveRegev3Dterse`):
 print(`GW, GWf `):
 print(` ProveRegev, R , R1, RegevSummand, RegevStyle3YTwith`):
 print(` Sipur3, Sipur3with, Sipur4, `):
 print(`  TnuYTodd, TnuYToddS, TnuYTwith  `):



elif nops([args])=1 and op(1,[args])=ApplyOper then
print(`ApplyOper(S,ope,n,N): applies the operator ope(n,1/N) to the`):
print(`sequence S as soon as it makes sense`):
print(`For example, try: ApplyOper([a,b,c],1-n/N,n,N);`):

elif nops([args])=1 and op(1,[args])=AR then
print(`AR(n,k): total number of Young tableaux with at most k parts`):
print(`with n cells. For example, try: AR(8,2);`):

elif nops([args])=1 and op(1,[args])=ARa then
print(`ARa(a,n): total number of Young tableaux with at most nops(a) parts`):
print(`with n cells and such that the cells labelled 1, ... sum(a) `):
print(`form shape a`):
print(`For example, try: ARa([0,0],8);`):

elif nops([args])=1 and op(1,[args])=ARaf then
print(`ARaf(a,n): total number of Young tableaux with at most nops(a) parts`):
print(`with n cells and such that the cells labelled 1, ... sum(a) `):
print(`form shape a`):
print(`using redunant generating functions`):
print(`For example, try: ARaf([0,0],8);`):

elif nops([args])=1 and op(1,[args])=ARaS then
print(`ARaS(a,N): the sequence, for n=1..N for the `):
print(`total number of Young tableaux with at most nops(a) parts`):
print(`with n cells and such that the cells labelled 1, ... sum(a) `):
print(`form shape a`):
print(`For example, try: ARaS([0,0],8);`):

elif nops([args])=1 and op(1,[args])=Axkn then
print(`Axkn(x,k,n): The generating function for all partitions`):
print(`of n, for example, try:Axkn(x,2,5);`):

elif nops([args])=1 and op(1,[args])=Bxkn then
print(`Bxkn(x,k,n): The generating function for all partitions`):
print(`of n, times (1-x[1]/x[k)*(1-x[2]/x[k-1])*...`):
print(`where all negative exponets are discarded`):
print(`For example, try:Bxkn(x,2,20);`):


elif nops([args])=1 and op(1,[args])=CF then 
print(`CF(a,A): The closed-form expression `):
print(`(where k=nops(a)) for the number of good walks`):
print(`that start at a`):
print(`For example, try: CF([0,0],A);`):

elif nops([args])=1 and op(1,[args])=CF1 then 
print(`CF1(a,A): The closed-form expression (divided by the`):
print(`multinomial coefficient (A[1]+...+A[k])!/A[1]!/.../A[k]!`):
print(`(where k=nops(a)) for the number of good walks`):
print(`that start at a`):
print(`For example, try: CF1([0,0],A);`):



elif nops([args])=1 and op(1,[args])=GuessAndProveRegev3D then 
print(` GuessAndProveRegev3D(n,k,a): inputs symbols `):
print(` n and k (for phrasing the output) and a lattice `):
print(` point a=[a1,a2,a3] with a1>=a2>=a3>=0 integers `):
print(` conjectures an expression (in terms of the Motzkin `):
print(` numbers M(n)) for the number of walks of length `):
print(` n-(a1+a2+a3) starting at [a1,a2,a3], with `):
print(`unit positive steps ([1,0,0],[0,1,0],[0,0,1]) `):
print(`always staying in x>=y>=z, and proceeds to prove`):
print(`it completely rigorously. For example,for `):
print(` the original Regev case try: `):
print(` GuessAndProveRegev3D(n,k,[2,1,0]); `):


elif nops([args])=1 and op(1,[args])=GuessAndProveRegev3Dterse then 
print(` GuessAndProveRegev3Dterse(n,a,M): The inputs is a symbol `):
print(` n  (for phrasing the output) and a lattice `):
print(` point a=[a1,a2,a3] with a1>=a2>=a3>=0 integers `):
print(` and a symbol M for denoting the Motzkin numbers.`):
print(` The output is`):
print(` an expression (in terms of the Motzkin `):
print(` numbers M(n)) for the number of walks of length `):
print(` n-(a1+a2+a3) starting at [a1,a2,a3], with `):
print(`unit positive steps ([1,0,0],[0,1,0],[0,0,1]) `):
print(`always staying in x>=y>=z, and proceeds to prove`):
print(`it completely rigorously, but internally. `):
print(`It only outputs the expression in terms of the Motzkin numbers`):
print(`No news is good news. If it fails to prove it, it will tell you so`):
print(` For example,for `):
print(` the original Regev case try: `):
print(` GuessAndProveRegev3Dterse(n,[2,1,0],M); `):

elif nops([args])=1 and op(1,[args])=GenOper then 
print(`GenOper(n,N,Ord,Deg,a): a generic recurrence operator with`):
print(`polynomial coefficients with symbol n and shift operator N`):
print(`of order Ord and and deggree Deg, followed by the set of`):
print(`coefficients. For example try: GenOper(n,N,1,1,a);`):


elif nops([args])=1 and op(1,[args])=GuessOper1 then
print(`GuessOper1(S1,S2,Ord,Deg,n,N): Given two sequences S1 and S2`):
print(`guessesan operator P such that S1=P(n,N)S1 for example try:`):
print(`GuessOper1(ARaS([0,0,0],20),ARaS([2,1,0],20),3,0,n,N);`):


elif nops([args])=1 and op(1,[args])=GuessRegev3D then
print(`GuessRegev3D(a,N): Given a lattice point a=[a1,a2,a3]`):
print(`with a1>=a2>=a3 guesses a recurrence operaor (let's call`):
print(` it P(N)`):
print(`with constant coefficients (where N is the shift in n)`):
print(`such that the number of walks with n-sum(a) steps starting`):
print(`at a is P(n)M(n) where M(n) is the n-th Motzkin number`):
print(`For example, try: GuessRegev3D([2,1,0],N);`):


elif nops([args])=1 and op(1,[args])=GuessRegev4D then
print(`GuessRegev4D(a,n,N): Given a lattice point a=[a1,a2,a3,a4]`):
print(`with a1>=a2>=a3>=a4 guesses a recurrence operaor (let's call`):
print(` it P(n,N)`):
print(`with coefficients that are polynomials of degree 1 in n`):
print(`(where N is the shift in n)`):
print(`such that the number of walks with n-sum(a) steps starting`):
print(`at a is P(n,N)M(n) where M(n) is the number of`):
print(`n-step ballot-walks that start at [0,0,0,0]`):
print(`For example, try: GuessRegev4D([3,2,1,0],n,N);`):

elif nops([args])=1 and op(1,[args])=GuessRegev4Dnice then
print(`GuessRegev4Dnice(a,n,M): Given a lattice point a=[a1,a2,a3,a4]`):
print(`with a1>=a2>=a3>=a4 guesses an expression`):
print(`for the number of walks with n-sum(a) steps starting`):
print(`at a `):
print(` in terms of M(n): the analogous quantity starting at [0,0,0,0],`):
print(`For example, try: GuessRegev4Dnice([3,2,1,0],n,M);`):

elif nops([args])=1 and op(1,[args])=GW then
print(`GW(a,b): the number of walks from a to b with unit steps`):
print(`always staying in x_1>=x_2>=...>=x_k (where k=nops(a))`):
print(`For example try: GW([0,0],[4,4]);`):

elif nops([args])=1 and op(1,[args])=GW1 then
print(`GW1(a,b): the number of walks from a to b with unit steps`):
print(`always staying in x_1>=x_2>=...>=x_k (where k=nops(a))`):
print(`Using Closed-Form `):
print(`For example try: GW([0,0],[4,4]);`):

elif nops([args])=1 and op(1,[args])=GWf then
print(`GWf(a,b): the number of walks from a to b with unit steps`):
print(`always staying in x_1>=x_2>=...>=x_k (where k=nops(a))`):
print(`Using the redundant generating function`):
print(`For example try: GWf([0,0],[4,4]);`):

elif nops([args])=1 and op(1,[args])=HafelOper then 
print(`HafelOper(Oper,x,k,a): applies the recurrence operator Oper`):
print(`(phrased  as a monomial in x[1], ..., x[k]) to the`):
print(`multinomial coefficient (a[1]+...+a[k])!/a[1]!/.../a[k]!`):
print(`and divides by the latter`):
print(`For example, try: HafelOper(1-x[2]/x[1],x,2,a);`):


elif nops([args])=1 and op(1,[args])=nuYTwith then
print(`nuYTwith(shape,cell,m): the number of Young-tableaux of shape`):
print(`shape with the entry in cell cell being m.`):
print(`For example, try: nuYTwith([2,2],[2,2],4); `):

elif nops([args])=1 and op(1,[args])=Pars then
print(`Pars(n,k): all partitions of n with at most k parts.`):
print(`For example, try: Pars(7,3);`):

elif nops([args])=1 and op(1,[args])=Pars1 then
print(`Pars1(K,n,k): all partitions of n with at most k parts.`):
print(`with largest part K`):
print(`For example, try: Pars1(4,7,3);`):

elif nops([args])=1 and op(1,[args])=ParsWith then
print(`ParsWith(n,k,cell): all partitions of n with at most k parts.`):
print(`that contain the cell cell=[i,j]`):
print(`For example, try: ParsWith(7,3,[2,2]);`):

elif nops([args])=1 and op(1,[args])=ProveRegevOriginal then
print(`ProveRegevOriginal(n,k): proves the original`):
print(`Regev conjecture verbosely, for paths of length`):
print(` n that start at [2,1,0]`):
print(`using the symbols n and k`):

elif nops([args])=1 and op(1,[args])=RegevStyle3YTwith  then
print(`RegevStyle3YTwith(cell,m,M,n): the explicit expression`):
print(`in terms of the Motzkin number sequence M(n) for`):
print(`the expression for the total number of`):
print(`of 3-rowed Young-tableaux with n cells whose`):
print(`entry in cell is m. For example, try:`):
print(`RegevStyle3YTwith([2,1],3,M,n); `):

elif nops([args])=1 and op(1,[args])=RegevStyle3YTwith1  then
print(`RegevStyle3YTwith1(cell,m,N): the explicit operator`):
print(`applied to the Motzkin number sequence M(n) for`):
print(`the expression for the total number of`):
print(`of 3-rowed Young-tableaux with n cells whose`):
print(`entry in cell is m. For example, try:`):
print(`RegevStyle3YTwith1([1,1],1,N); `):


elif nops([args])=1 and op(1,[args])=RGF then
print(`RGF(a,x): The redundant generating function (in x[1],x[2], ..)`):
print(`for good walks`):
print(`that start at a`):
print(`For example, try: RGF([0,0],x);`):

elif nops([args])=1 and op(1,[args])=RGFt then
print(`RGFt(a,x,t): The redundant generating function (in x[1],x[2], ..)`):
print(`for the total number of good walks`):
print(`that start at a`):
print(`For example, try: RGFt([0,0],x,t)`):


elif nops([args])=1 and op(1,[args])=R then
print(`R(a,A): The expression `):
print(`in A[1],A[2],A[3]`):
print(`such that the number of good walks from`):
print(`a to [A[1],A[2],A[3]] can be written as`):
print(`R(a;A[1],A[2],A[3])-R(a;A[1]+1,A[2],A[3]-1):`):
print(`For example, try: R([0,0,0],A);`):

elif nops([args])=1 and op(1,[args])=R1 then
print(`R1(a,k,n): the two summands for the parts`):
print(`k<=n/3 and n/3<=k<=n/2 for the Regev-style`):
print(`single-sum for the the number of good walks`):
print(`of length n that start at a, for example`):
print(`try: R1([0,0,0],k,n);`):


elif nops([args])=1 and op(1,[args])=RegevSummand then
print(`RegevSummand(a,n,k): The summand for the Regev-style`):
print(`single-binomial coefficient sum for good walks that`):
print(`start at a. For example, try:`):
print(`RegevSummand([0,0,0],n,k):`):

elif nops([args])=1 and op(1,[args])=Sipur3 then
print(`Sipur3(n,M,L,K): All the expressions for the number of lattice  paths`):
print(`of length n-(a1+a2) that start at [a1,a2,0] and stay`):
print(`in x>=y>=z for L>=a1>=a2>0. `):
print(`The program rigorously(!) proves each statement`):
print(`It also prints out the first K terms`):
print(`For example, try:`):
print(`Sipur3(n,M,3,20);`):



elif nops([args])=1 and op(1,[args])=Sipur3with then
print(`Sipur3with(n,M,L,K): All the expressions for the number of `):
print(`<=3-rowed Young-tableaux with n cells, in terms of`):
print(`the Motzkin numbers M, whose [i,j] cell has m`):
print(`for m<=L and i*j<=m`):
print(`it all also prints out the first K terms`):
print(`For example, try:`):
print(`Sipur3with(n,M,3,20);`):

elif nops([args])=1 and op(1,[args])=Sipur4 then
print(`Sipur4(n,M,L,K): All the expressions for the number of lattice  paths`):
print(`of length n-(a1+a2+a3) that start at [a1,a2,a3,0] and stay`):
print(`in x>=y>=z>=w for L>=a1>=a2>=a3>0. `):
print(`The program semi-rigorously proves each statement`):
print(`It also prints out the first K terms`):
print(`For example, try:`):
print(`Sipur4(n,M,3,20);`):

elif nops([args])=1 and op(1,[args])=SYTkn then
print(`SYTkn(k,n): the set of Standard Young Tableaux of <=k rows`):
print(`with n cells. For example, try`):
print(`SYTkn(3,5);`):

elif nops([args])=1 and op(1,[args])=TnuYTodd then
print(`TnuYTodd(k,n,cell): the total number of <=k-rowed Young-tableaux`):
print(`of n cells such that cell has an odd entry`):
print(`For example, try: TnuYTodd(3,10,[1,2]);`):

elif nops([args])=1 and op(1,[args])=TnuYToddS then
print(`TnuYToddS(k,n,cell,m): the sequence, for n=1..N`):
print(`for the total number of <=k-rowed Young-tableaux`):
print(`of n cells such that cell cell has an odd entry in it.`):
print(`For example, try: TnuYToddS(3,10,[1,2],3);`):

elif nops([args])=1 and op(1,[args])=TnuYTwith then
print(`TnuYTwith(k,n,cell,m): the total number of <=k-rowed Young-tableaux`):
print(`of n cells such that cell cell has m in it.`):
print(`For example, try: TnuYTwith(3,10,[1,2],3);`):


elif nops([args])=1 and op(1,[args])=TnuYTwithS then
print(`TnuYTwithS(k,n,cell,m): the sequence, for n=m..N+m`):
print(`for the total number of <=k-rowed Young-tableaux`):
print(`of n cells such that cell cell has m in it.`):
print(`For example, try: TnuYTwithS(3,10,[1,2],3);`):


elif nops([args])=1 and op(1,[args])=YTwith then
print(`YTwith(k,n,cell,m): the set of <=k-rowed Young-tableaux`):
print(`of n cells such that cell cell has m in it.`):
print(`for checking purposes)`):
print(`For example, try: YTwith(3,10,[1,2],3);`):

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



#GW(a,b): the number of walks from a to b with unit steps
#always staying in x_1>=x_2>=...>=x_k (where k=nops(a))
#For example try: GW([0,0],[4,4]);
GW:=proc(a,b) local k,i:
option remember:
k:=nops(a):

if nops(b)<>k then
 ERROR(`Bad input`):
fi:

 if not {seq(evalb(b[i]>=a[i]),i=1..k)}={true} then
      0:
 elif b=a then
      1:
 elif not {seq(evalb(b[i]-b[i+1]>=0),i=1..k-1)}={true} then
     0:
 else
    add(GW(a,[op(1..i-1,b),b[i]-1,op(i+1..k,b)]),i=1..k):
 fi:
end:
 
 
 
#Pars1(K,n,k): all partitions of n with at most k parts.
#and largest part K
#For example, try: Pars1(3,7,3);
Pars1:=proc(K,n,k) local n1,K1,mu,gu,i:
option remember:
  if k=1 then
    if n=K then
      RETURN({[n]}):
    else
      RETURN({}):
    fi:
 fi:
if n<K then
 RETURN({}):
elif n=K then
  RETURN({[K,0$(k-1)]}):
else
   n1:=n-K:
    gu:={}:
    for K1 from 0 to K do
     mu:=Pars1(K1,n1,k-1):
     gu:=gu union {seq([K,op(mu[i])],i=1..nops(mu))}:
    od:
 RETURN(gu):
fi:

end:

#Pars(n,k): all partitions of n with at most k parts.
#For example, try: Pars(7,3);
Pars:=proc(n,k) local K:
option remember:
{seq(op(Pars1(K,n,k)),K=0..n)}:

end:

#AR(n,k): total number of Young tableaux with at most k parts
#with n cells. For example, try: AR(8,2);
AR:=proc(n,k) local gu,b:
gu:=Pars(n,k):
add(GW([0$k],b), b in gu):
end:

#ARa(a,n): total number of Young tableaux with at most k=nops(a) parts
#such that the the cells labelled 1, 2, ..., sum(a) has shape a
#and with n cells. For example, try: ARa([0,0],8);
ARa:=proc(a,n) local gu,b,k:
k:=nops(a):
gu:=Pars(n,k):
add(GW(a,b), b in gu):
end:

#Ref(x,i): The reflection of a vector x with respect to
#x[i]-x[i+1]=-1
Ref:=proc(x,i) local k:
k:=nops(x):

if not (i>=1 and i<=k-1) then
 ERROR(`Bad input`):
fi:

[op(1..i-1,x),x[i+1]-1,x[i]+1,op(i+2..k,x)]:
end:


#Images(a): all the images of the vector a under reflections
#w.r.t. to the hyperplanes x[i]-x[i+1]=-1
#Try: Images([a,b,c]);
Images:=proc(a) local gu,gu1:
gu:={[1,a]}:
gu1:=OneStepImages(gu):

while gu<>gu1 do
gu:=gu1:
gu1:=OneStepImages(gu):
od:
gu1:
end:



OneStepImages1:=proc(a) local k,i:
k:=nops(a):
{seq(Ref(a,i) ,i=1..k-1)}:
end:

#OneStepImages(S): all the reflections of elements of S
#w.r.t. to ONE of the hyperplanes x[i]-x[i+1]=-1 (at a time)
#Try: OneStepImages({[a,b,c]});
OneStepImages:=proc(S) local gu,w,mu,j,a,s:
gu:={}:

for s in S do
 w:=s[1]:
 a:=s[2]:
 mu:=OneStepImages1(a):
 gu:=gu union {seq([-w,mu[j]],j=1..nops(mu))}:
od:

gu union S:
end:


#RGF(a,x): The redundant generating function (in x[1],x[2], ..)
#for good walks
#that start at a
#For example, try: RGF([0,0],x);
RGF:=proc(a,x) local gu,k,i,j:
 k:=nops(a):
 gu:=Images(a):
add(gu[i][1]*mul(x[j]^gu[i][2][j],j=1..k),i=1..nops(gu))/
(1-add(x[j],j=1..k)):
end:

#RGF1(a,x): The numerator of the
#redundant generating function (in x[1],x[2], ..)
#for good walks
#that start at a
#For example, try: RGF1([0,0],x);
RGF1:=proc(a,x) local gu,k,i,j:
 k:=nops(a):
 gu:=Images(a):
add(gu[i][1]*mul(x[j]^gu[i][2][j],j=1..k),i=1..nops(gu)):
end:


#GWf(a,b): the number of walks from a to b with unit steps
#always staying in x_1>=x_2>=...>=x_k (where k=nops(a))
#using constant terms and generating functions
#For example try: GWf([0,0],[4,4]);
GWf:=proc(a,b) local gu,x,i,k:
k:=nops(a):
gu:=RGF(a,x):

for i from 1 to k do
gu:=coeff(series(gu,x[i]=0,b[i]+abs(degree(numer(gu),x[i]))+2),x[i],b[i]):
gu:=normal(gu):
od:
gu:

end:

#RGFt(a,x,t): The redundant generating function (in x[1],x[2], ..)
#for the total number of good walks
#that start at a
#For example, try: RGFt([0,0],x,t);
RGFt:=proc(a,x,t) local gu,i,j,k:
 k:=nops(a):
gu:=RGF(a,x):
gu/mul(1-mul(t/x[j],j=1..i),i=1..k):
end:

#CT1(F,x,k): The constant term of any rational function of
#form Laurent(x[1], ...,x[k])/(1-x[1]-...-x[k])
CT1:=proc(F,x,k) local mu,gu,mu1,c,i,i1:
mu:=expand(normal(F*(1-add(x[i],i=1..k)))):
gu:=0:


for i from 1 to nops(mu) do
mu1:=op(i,mu):
if {seq(evalb(degree(mu1,x[i1])<=0),i1=1..k)}={true} then
 c:=normal(mu1/mul(x[i1]^degree(mu1,x[i1]),i1=1..k)):
  if not type(c,integer) then
    RETURN(FAIL):
   fi:
 gu:=gu+c*add(-degree(mu1,x[i1]),i1=1..k)!/
mul((-degree(mu1,x[i1]))!,i1=1..k):
fi:
od:
gu:
end:



#ARf(n,k): total number of Young tableaux with at most k parts
#with n cells. Using redudant generating function
#For example, try: ARf(8,2);
ARf:=proc(n,k) local gu,x,t:
gu:=RGFt([0$k],x,t):
gu:=taylor(gu,t=0,n+1):
gu:=normal(coeff(gu,t,n)):
CT1(gu,x,k):
end:


#ARaf(a,n): total number of Young tableaux with at most k parts
#with n cells. 
#such that the the cells labelled 1, 2, ..., sum(a) has shape a
#Using redudant generating function
#For example, try: ARaf([0,0],8);
ARaf:=proc(a,n) local gu,x,t,k:
k:=nops(a):
gu:=RGFt(a,x,t):
gu:=taylor(gu,t=0,n+1):
gu:=normal(coeff(gu,t,n)):
CT1(gu,x,k):
end:

#Sym(f,x): the symmetrizer of f(x) (x is a list of variables)
#For example, try: Sym(1/(1-x-2*y),[x,y]):
Sym:=proc(f,x) local k,gu,mu,pi,i:
k:=nops(x):
mu:=permute(k):

gu:=0:
for pi in mu do
gu:=gu+subs({seq(x[i]=x[pi[i]],i=1..k)},f):
od:
gu:=normal(gu/k!):
end:



#RGFtSym(a,x,t): The Symmetrized
#redundant generating function (in x[1],x[2], ..)
#for the total number of good walks
#that start at a
#For example, try: RGFtSym([0,0],x,t);
RGFtSym:=proc(a,x,t) local i:
Sym(RGFt(a,x,t),[seq(x[i],i=1..nops(a))]):
end:

#Axkn(x,k,n): The generating function for all partitions
#of n, for example, try:Axkn(x,2,5);
Axkn:=proc(x,k,n) local gu,i,j,t:
gu:=1/mul(1-mul(t*x[j],j=1..i),i=1..k):
gu:=taylor(gu,t=0,n+1):
expand(coeff(gu,t,n)):
end:


F1abc:=proc(a1,a2,a3):
(1-a2/(a1+1)-a3/(a1+1)+a2*a3/(a1+1)/(a1+2))*
(a1+a2+a3)!/a1!/a2!/a3!:
end:

Gabc:=proc(a,b,c): (a-b+1)*(a-c+2)*(b-c+1)*
(a+b+c)!/(a+2)!/(b+1)!/c!:
end:

GabcOld:=proc(a,b,c): (a-b+1)*(a-c+2)*(b-c+1)/(a+1)/(a+2)/(b+1)*
(a+b+c)!/a!/b!/c!:
end:

Regev:=proc(n) local gu,a1,m:
gu:=0:
for a1 from 0 to trunc(n/3) do
m:=n-a1:
 if m mod 2=0 then
  gu:=gu+F1abc(a1,m/2,m/2):
 else
 gu:=gu+F1abc(a1,(m+1)/2,(m-1)/2):
 fi:
od:

end:

  
AR3:=proc(n) local gu,a,b,c,m:
gu:=0:

for a from trunc(n/3) to n do
  m:=n-a:
   for b from trunc(m/2) to a do
   c:=n-a-b:
 if a>=0 and b>=0 and c>=0 and a>=b and b>=c then  
    gu:=gu+Gabc(a,b,c):
 fi:
 od:
od:

gu:
end:


AR3new:=proc(n) local gu,a:
gu:=0:

for a from trunc(n/3) to n do
 print(ifactor(AR3new1(n,a))):
  gu:=gu+AR3new1(n,a):
od:

gu:
end:


AR3new1:=proc(n,a) local b,m,gu,c:
gu:=0:
 m:=n-a:
   for b from trunc(m/2) to a do
   c:=n-a-b:
 if a>=0 and b>=0 and c>=0 and a>=b and b>=c then  
    gu:=gu+Gabc(a,b,c):
 fi:
 od:
gu:
end:


ARaS:=proc(a,N) local n:[seq(ARa(a,n),n=1..N)]:end:


Bxkn:=proc(x,k,n) local gu,i,j:
  gu:=Axkn(x,k,n):
  gu:=expand(gu*mul(1-x[i]/x[k+1-i],i=1..trunc(k/2))):
  
 for i from 1 to k do
  gu:=add(coeff(gu,x[i],j)*x[i]^j,j=0..degree(gu,x[i])):
 od:

gu:=expand(gu):

gu:=sort(gu):

end:

CleanUp:=proc(g,x,k) local i,j,gu:
gu:=g:
 for i from 1 to k do
  gu:=add(coeff(gu,x[i],j)*x[i]^j,j=0..degree(gu,x[i])):
 od:
gu:=expand(gu):
collect(gu,x[1]):
end:

#HafelMono(Mono,x,k,a): applies the shift-monomial Mono
#(phrased  as a monomial in x[1], ..., x[k]) to the
#multinomial coefficient (a[1]+...+a[k])!/a[1]!/.../a[k]!
#and divides by the latter
#For example, try: HafelMono(x[2]/x[1],x,2,a)
HafelMono:=proc(Mono,x,k,a) local F,G,i:
F:=add(a[i],i=1..k)!/mul(a[i]!,i=1..k):

G:=subs({seq(a[i]=a[i]-degree(Mono,x[i]),i=1..k)},F):
normal(simplify(expand(G/F))):
end:

#Mekadem(Mono,x,k): the coeff. of Mono
Mekadem:=proc(Mono,x,k) local i:
normal(Mono/mul(x[i]^degree(Mono,x[i]),i=1..k)):
end:


  
  
 
#HafelOper(Oper,x,k,a): applies the recurrence operator Oper
#(phrased  as a monomial in x[1], ..., x[k]) to the
#multinomial coefficient (a[1]+...+a[k])!/a[1]!/.../a[k]!
#and divides by the latter
#For example, try: HafelOper(1-x[2]/x[1],x,2,a)
HafelOper:=proc(Oper,x,k,a) local i,mu,gu:
mu:=expand(Oper):

if not type(mu,`+`) then
 RETURN(HafelMono(mu,x,k,a)):
fi:

gu:=0:
for i from 1 to nops(mu) do
 gu:=gu+Mekadem(op(i,mu),x,k)*HafelMono(op(i,mu),x,k,a):
od:

factor(normal(gu)):

end:



#CF1(a,A): The closed-form expression (divided by the
#multinomial coefficient (A[1]+...+A[k])!/A[1]!/.../A[k]!
#(where k=nops(a)) for the number of good walks
#that start at a
#For example, try: RGF([0,0],A);
CF1:=proc(a,A) local gu,k,x:
k:=nops(a):
gu:=RGF1(a,x):
HafelOper(gu,x,k,A):
end:




#CF(a,A): The closed-form expression 
#(where k=nops(a)) for the number of good walks
#that start at a
#For example, try: CF([0,0],A);
CF:=proc(a,A) local i,k:
k:=nops(a):
CF1(a,A)*add(A[i],i=1..k)!/mul(A[i]!,i=1..k):
end:




#GW1(a,b): the number of walks from a to b with unit steps
#always staying in x_1>=x_2>=...>=x_k (where k=nops(a))
#using CF
#For example try: GW1([0,0],[4,4]);
GW1:=proc(a,b) local gu,A,i,k:
gu:=CF(a,A):
k:=nops(a):
eval(subs({seq(A[i]=b[i],i=1..k)},gu)):
end:



#R(a,A): The expression 
#in A[1],A[2],A[3]
#such that the number of good walks from
#a to [A[1],A[2],A[3]] can be written as
#R(a;A[1],A[2],A[3])-R(a;A[1]+1,A[2],A[3]-1):
#For example, try: R([0,0,0],A);
R:=proc(a,A)  local gu,x:
if nops(a)<>3 then
 ERROR(`Bad input`):
fi:
gu:=normal(RGF1(a,x)/(1-x[3]/x[1])):
normal(HafelOper(gu,x,3,A)):
end:



#R1(a,k,n): the two summands for the parts
#k<=n/3 and n/3<=k<=n/2 for the Regev-style
#single-sum for the the number of good walks
#of length n that start at a, for example
#try: R1([0,0,0],k,n);
R1:=proc(a,k,n)  local A,gu:
gu:=R(a,A);
factor([subs({A[1]=n-2*k,A[2]=k,A[3]=k},gu),
subs({A[1]=k,A[2]=k,A[3]=n-2*k},gu)]):
end:

#RegevSummand(a,n,k): The summand for the Regev-style
#single-binomial coefficient sum for good walks that
#start at a. For example, try:
#RegevSummand([0,0,0],n,k):
RegevSummand:=proc(a,n,k) local gu:
gu:=R1(a,k,n):

if gu[1]<>gu[2] then
 RETURN(FAIL):
fi:
gu[1]*n!/k!^2/(n-2*k)!:
end:


ProveRegevOriginal:=proc(n,k) local gu,mu,gu1:
print(`Theorem: Let M(n) be the number of walks of length n`):
print(`from (0,0,0) using unit positive steps and staying in`):
print(`the region x>=y>=z>=0. Let A(n) be the number of such walks`):
print(`with n-3 steps, that start at (2,1,0). Then `):
print(`A(n)=M(n-1)-M(n-3) . `):
print():
print(`Proof:`):
print(`The summand for M(n),`):
print(`let's call it F(n,k), is`):
gu:=RegevSummand([0,0,0],n,k):
print(gu):

print(`The summand for A(n), `):
print(`let's call it G(n,k), is`):
mu:=RegevSummand([2,1,0],n,k):
print(mu):
print(`To prove Amitai Regev's conjecture we must show that`):
print(`Sum(F(n-1,k)-F(n-3,k),k=0..n)=Sum(G(n,k),k=0..n).`):
print(`Hence the summand of the difference,`):
print(`F(n-1,k)-F(n-3,k)-G(n,k), is `):
gu:=normal(simplify(expand(subs(n=n-1,gu)/gu))-
simplify(expand(subs(n=n-3,gu)/gu))
-
simplify(expand(mu/gu))
)*gu:
print(gu):
print(`But this is a telescoping sum, the summand being H(n,k+1)-H(n,k)`):
print(`where H(n,k) is`):
gu1:=sum(gu,k):
if simplify(subs(k=k+1,gu1)-gu1-gu)<>0 then
 RETURN(FAIL):
fi:
print(gu1):
print(`Check! (or believe me). This concludes the proof of Amitai Regev's`):
print(`conjecture. QED.`):
end:

#GenOper(n,N,Ord,Deg,a): a generic recurrence operator with
#polynomial coefficients with symbol n and shift operator N
#of order Ord and and deggree Deg, followed by the set of
#coefficients. For example try: GenOper(n,N,1,1,a);
GenOper:=proc(n,N,Ord,Deg,a) local ope,i,j,var:
ope:=0:
var:={}:

for i from 0 to Ord do
for j from 0 to Deg do
 var:=var union {a[i,j]}:
 ope:=ope+a[i,j]*n^j*N^(-i):
od:
od:

ope,var:

end:


#ApplyOper(S,ope,n,N): applies the operator ope(n,1/N) to the
#sequence S as soon as it makes sense
ApplyOper:=proc(S,ope,n,N) local gu,L,kha,n1,i:
gu:=[]:

L:=-ldegree(ope,N):

for n1 from L+1 to nops(S) do
 kha:=add(S[n1-i]*subs(n=n1,coeff(ope,N,-i)),i=0..L):
 gu:=[op(gu),kha]:
od:
gu:
end:



#GuessOper1(S1,S2,Ord,Deg,n,N): Given two sequences S1 and S2
#guessesan operator P such that S2=P(n,N)S1 for example try:
#GuessOper1(ARaS([0,0,0],20),ARaS([2,1,0],20),3,0,n,N);
GuessOper1:=proc(S1,S2,Ord,Deg,n,N) local ope, var,S1t,S2t,L,var1,a,i,eq:

if nops(S1)<>nops(S2) then
 ERROR(`Bad input`):
fi:

ope:=GenOper(n,N,Ord,Deg,a):
var:=ope[2]:
ope:=ope[1]:
L:=-ldegree(ope,N):
if nops(S1)-L-1-nops(var)<9 then
  ERROR(`Insufficient data`):
fi:

S1t:=ApplyOper(S1,ope,n,N):

S2t:=[op(L+1..nops(S2),S2)]:


eq:={seq(S1t[i]-S2t[i],i=3..nops(S1t))}:

var1:=solve(eq,var):

if var1=NULL then
 RETURN(FAIL):
else
 RETURN(subs(var1,ope)):
fi:
end:



#GuessRegev3D(a,N): Given a lattice point a=[a1,a2,a3]
#with a1>=a2>=a3 guesses a recurrence operaor (let's call
# it P(N)
#with constant coefficients (where N is the shift in n)
#such that the number of walks with n-sum(a) steps starting
#at a is P(n)M(n) where M(n) is the n-th Motzkin number
#For example, try: GuessRegev3D([2,1,0],N);
GuessRegev3D:=proc(a,N) local n,L,orekh,gu:
L:=convert(a,`+`):

orekh:=11+2*L:
gu:=GuessOper1(ARaS([0,0,0],orekh),ARaS(a,orekh),L,0,n,N) :

if gu=FAIL then
 RETURN(FAIL):
fi:
sort(gu):
end:



#GuessAndProveRegev3D(n,k,a): inputs symbols
#n and k (for phrasing the output) and a lattice
#point a=[a1,a2,a3] with a1>=a2>=a3>=0 integers
#conjectures an expression (in terms of the Motzkin
#numbers M(n)) for the number of walks of length
#n-(a1+a2+a3) starting at [a1,a2,a3], with
#unit positive steps ([1,0,0],[0,1,0],[0,0,1])
#always staying in x>=y>=z, and proceeds to prove
#it completely rigorously. For example,for
#the original Regev case try:
#GuessAndProveRegev3D(n,k,[2,1,0]);
GuessAndProveRegev3D:=proc(n,k,a) local gu,mu,gu1,ope,N,i:

ope:=GuessRegev3D(a,N):
if ope=FAIL then
 RETURN(FAIL):
fi:
print(`Theorem: Let M(n) be the number of walks of length n`):
print(`from  [0,0,0], using unit positive steps and staying in`):
print(`the region x>=y>=z>=0. These are the famous Motzkin numbers,`):
print(`as first proved by Amitai Regev in 1981.`):
print(`Let A(n) be the number of such walks`):
print(`with `, n-convert(a,`+`),  `steps, that start at`,a,  ` Then `):
print(A(n)=add(coeff(ope,N,-i)*M(n-i),i=-degree(ope,N)..-ldegree(ope,N))):
print():
print(`Proof:`):
print(` Thanks to Regev's seminal result (that can be reproved `):
print(` with our approach), the summand for M(n),`):
print(`let's call it F(n,k), is`):
gu:=RegevSummand([0,0,0],n,k):
print(gu):

print(`The summand for A(n), `):
print(`let's call it G(n,k), is`):
mu:=RegevSummand(a,n,k):
print(mu):
print(`To prove the Theorem, we must show that`):
print(`the sum (w.r.t. k), whose summand is`):

print(G(n,k)-add(coeff(ope,N,-i)*F(n-i,k),i=-degree(ope,N)..-ldegree(ope,N))):

print(`is identically zero`):
gu:=
normal(
add(
coeff(ope,N,-i)*
normal(
simplify(expand(subs(n=n-i,gu)/gu))
      ),
i=-degree(ope,N)..-ldegree(ope,N))
-expand(simplify(mu/gu)) )*gu:
#HERE
print(`But this is a telescoping sum, the summand being H(n,k+1)-H(n,k)`):
print(`where H(n,k) is`):
gu1:=sum(gu,k):
if simplify(subs(k=k+1,gu1)-gu1-gu)<>0 then
 RETURN(FAIL):
fi:
print(gu1):
print(`Check! (or believe me). This concludes the proof of our theorem`):
print(` that is an analog of Amitai Regev's`):
print(` original conjecture (whose starting point was [2,1,0]). QED.`):
end:



#GuessAndProveRegev3Dterse(n,k,a,M): inputs symbols
#n (for phrasing the output) and a lattice
#point a=[a1,a2,a3] with a1>=a2>=a3>=0 integers
#conjectures an expression (in terms of the Motzkin
#numbers M(n)) for the number of walks of length
#n-(a1+a2+a3) starting at [a1,a2,a3], with
#unit positive steps ([1,0,0],[0,1,0],[0,0,1])
#always staying in x>=y>=z, and proceeds to prove
#it completely rigorously, but the output is terse. 
#the output is simply the expression in the Motzkin
#numbers M(n), where M is part of the input indicationg
#the symbol to be used, and if it does not complain
#it means that it proved it. 
#the original Regev case try:
#GuessAndProveRegev3Dterse(n,k,[2,1,0],M);
GuessAndProveRegev3Dterse:=proc(n,a,M) local gu,mu,gu1,ope,N,i,k:

ope:=GuessRegev3D(a,N):
if ope=FAIL then
 RETURN(FAIL):
fi:
gu:=RegevSummand([0,0,0],n,k):

mu:=RegevSummand(a,n,k):

gu:=
normal(
add(
coeff(ope,N,-i)*
normal(
simplify(expand(subs(n=n-i,gu)/gu))
      ),
i=-degree(ope,N)..-ldegree(ope,N))
-expand(simplify(mu/gu)) )*gu:
gu1:=sum(gu,k):
if simplify(subs(k=k+1,gu1)-gu1-gu)<>0 then
 RETURN(FAIL):
fi:

RETURN(add(coeff(ope,N,-i)*M(n-i),i=-degree(ope,N)..-ldegree(ope,N))):
end:


#Sipur3(n,M,L,K): All the expressions for the number of lattice  paths
#of length n-(a1+a2) that start at [a1,a2,0] and stay
#in x>=y>=z for L>=a1>=a2>. 
#it all also prints out the first K terms
#For example, try:
#Sipur3(n,M,3,20);
Sipur3:=proc(n,M,L,K) local a1, a2:
print(`Let `, M(n), ` be the Motzkin numbers `):
for a1 from 0 to L do
for a2 from 0 to a1 do
print(`The number of`,n-a1-a2, `-step 3D ballot-lattice paths starting at`, [a1,a2,0] , ` equals (and (rigorously!) proved)`):
print(GuessAndProveRegev3Dterse(n,[a1,a2,0],M)):
print(`The first`, K-a1-a2, `terms are `):
print(op(a1+a2..K,ARaS([a1,a2,0],K))):
od:
od:
end:





#GuessRegev4D(a,n,N): Given a lattice point a=[a1,a2,a3,a4]
#with a1>=a2>=a3>=a4 guesses a recurrence operaor (let's call
# it P(n,N)
#with  coefficients that are polynomials of degree 1 in n
#(where N is the shift in n)
#such that the number of walks with n-sum(a) steps starting
#at a is P(n,N)M(n) where M4(n) is this quantity when it
#starts at [0,0,0,0]
#For example, try: GuessRegev4D([3,2,1,0],n,N);
GuessRegev4D:=proc(a,n,N) local L,orekh,gu:
L:=convert(a,`+`):

orekh:=20+4*L:
gu:=GuessOper1(ARaS([0,0,0,0],orekh),ARaS(a,orekh),L,1,n,N) :

if gu=FAIL then
 RETURN(FAIL):
fi:
gu:
end:



#GuessRegev4Dnice(a,n,M): Given a lattice point a=[a1,a2,a3,a4]
#with a1>=a2>=a3>=a4 guesses an expression
#for the number of walks with n-sum(a) steps starting
#in terms of M(n), which is this quantity when it
#starts at [0,0,0,0]
#For example, try: GuessRegev4D([3,2,1,0],n,N);
GuessRegev4Dnice:=proc(a,n,M) local ope,N,i:
ope:=GuessRegev4D(a,n,N) :

RETURN(add(factor(coeff(ope,N,-i))*M(n-i),i=-degree(ope,N)..-ldegree(ope,N))):

end:



#Sipur4(n,M,L,K): All the expressions for the number of lattice  paths
#of length n-(a1+a2+a3) that start at [a1,a2,a3,0] and stay
#in x>=y>=z>=w for L>=a1>=a2>a3. For example, try:
#it also prints out the first K terms
#Sipur4(n,M,3,20);
Sipur4:=proc(n,M,L,K) local a1, a2,a3:
print(`Let `, M(n), ` be the Motzkin numbers `):
for a1 from 0 to L do
for a2 from 0 to a1 do
for a3 from 0 to a2 do
print(`The number of`,n-a1-a2-a3, `-step 3D ballot-lattice paths starting at`, [a1,a2,a3,0] , ` equals (and is provable, but we only proved it semi-rigorously)`):
print(GuessRegev4Dnice([a1,a2,a3,0],n,M)):
print(`The first`, K-a1-a2-a3, `terms are `):
print(op(a1+a2+a3..K,ARaS([a1,a2,a3,0],K))):

od:
od:
od:

end:




#ParsWith(n,k,cell): all partitions of n with at most k parts.
#that contain the cell cell=[i,j] as a corner-cell
#For example, try: ParsWith(7,3,[2,2]);
ParsWith:=proc(n,k,cell) local gu,s,i,j,mu:
i:=cell[1]: j:=cell[2]:

mu:=Pars(n,k):
gu:={}:
for s in mu do
if i<=k and j=s[i] then
  gu:=gu union {s}:
fi:
od:

gu:

end:



#nuYTwith(shape,cell,m): the number of Young-tableaux of shape
#shape with the entry in cell cell being m.
#For example, try: nuYTwith([2,2],[2,2],4) 
nuYTwith:=proc(shape,cell,m)
local gu,mu,i,j,k,sh,sh1:
k:=nops(shape):
i:=cell[1]: j:=cell[2]:

if m<i*j then
  RETURN(0):
fi:

mu:=ParsWith(m,k,cell):

gu:=0:
for sh in mu do
 sh1:=[op(1..i-1,sh),sh[i]-1,op(i+1..k,sh)]:
 gu:=gu+GW([0$k],sh1)*GW(sh,shape):
od:
gu:
end:





#TnuYTwith(k,n,cell,m): the total number of <=k-rowed Young-tableaux
#of n cells such that cell cell has m in it.
#For example, try: TnuYTwith(3,10,[1,2],3);
TnuYTwith:=proc(k,n,cell,m) local gu,mu,i,sh,sh1:

mu:=ParsWith(m,k,cell):
i:=cell[1]:
gu:=0:
for sh in mu do
 sh1:=[op(1..i-1,sh),sh[i]-1,op(i+1..k,sh)]:
 gu:=gu+GW([0$k],sh1)*ARa(sh,n):
od:
gu:

end:

#TnuYTwithS(k,N,cell,m): the sequence 
#for n=m..N+m, for the total number of <=k-rowed Young-tableaux
#of n cells such that cell cell has m in it.
#For example, try: TnuYTwithS(3,10,[1,2],3);
TnuYTwithS:=proc(k,N,cell,m) local n:

[seq(TnuYTwith(k,n,cell,m),n=m..N+m)]:

end:


#TnuYTodd(k,n,cell): the total number of <=k-rowed Young-tableaux
#of n cells such that cell has an odd entry
#For example, try: TnuYTodd(3,10,[1,2]);
TnuYTodd:=proc(k,n,cell) local m:
add(TnuYTwith(k,n,cell,2*m+1),m=0..trunc(n/2)+1):
end:

#TnuYToddS(k,N,cell): the sequence 
#for n=1..N, for the total number of <=k-rowed Young-tableaux
#of n cells such that cell cell has an odd entry  in it.
#For example, try: TnuYToddS(3,10,[1,2]);
TnuYToddS:=proc(k,N,cell) local n:

[seq(TnuYTodd(k,n,cell),n=1..N)]:

end:



#GuessAndProveRegev3DVeryTerse(a,N): inputs 
#a symbol for the shift operator N and a
#point a=[a1,a2,a3] with a1>=a2>=a3>=0 integers
#conjectures an expression (in terms of the Motzkin
#numbers M(n)) for the number of walks of length
#n-(a1+a2+a3) starting at [a1,a2,a3], with
#unit positive steps ([1,0,0],[0,1,0],[0,0,1])
#always staying in x>=y>=z, and proceeds to prove
#it completely rigorously, but the output
#just the operator (constant-coeff. in N)
#it means that it proved it. 
#the original Regev case try:
#GuessAndProveRegev3DVeryTerse([2,1,0],N);
GuessAndProveRegev3DVeryTerse:=proc(a,N) local gu,mu,gu1,ope,i,n,k:

ope:=GuessRegev3D(a,N):
if ope=FAIL then
 RETURN(FAIL):
fi:
gu:=RegevSummand([0,0,0],n,k):

mu:=RegevSummand(a,n,k):

gu:=
normal(
add(
coeff(ope,N,-i)*
normal(
simplify(expand(subs(n=n-i,gu)/gu))
      ),
i=-degree(ope,N)..-ldegree(ope,N))
-expand(simplify(mu/gu)) )*gu:
gu1:=sum(gu,k):
if simplify(subs(k=k+1,gu1)-gu1-gu)<>0 then
 RETURN(FAIL):
fi:

ope:
end:



#RegevStyle3YTwith(cell,m,N): the explicit operator
#applied to the Motzkin number sequence M(n) for
#the expression for the total number of
#of 3-rowed Young-tableaux with n cells whose
#entry in cell is m. For example, try:
#RegevStyle3YTwith([1,1],1,N):

RegevStyle3YTwith1:=proc(cell,m,N) local gu,mu,i,sh,sh1,ope:
gu:=0:
if cell[1]>3 then
 ERROR(`Bad input`):
fi:
mu:=ParsWith(m,3,cell):
i:=cell[1]:
gu:=0:

for sh in mu do
 sh1:=[op(1..i-1,sh),sh[i]-1,op(i+1..3,sh)]:
 ope:=GuessAndProveRegev3DVeryTerse(sh,N):
 
 if ope=FAIL then
  print(`Something went wrong`):
 fi:

 gu:=gu+GW([0$3],sh1)*ope:
od:
gu:

gu:
end:


#RegevStyle3YTwith(cell,m,M,n): gives 
#a (rigorously proved!) expression
#in terms of the Motzkin numbers M(n) 
#(M and n are inputted symbols)
#for the number of n-celled three-rowed Young tableaux
#whose entry at cell cell is m, for example, try:
#RegevStyle3YTwith([1,2],3,M,n);
RegevStyle3YTwith:=proc(cell,m,M,n) local N,ope,i:
ope:=RegevStyle3YTwith1(cell,m,N):
add(coeff(ope,N,-i)*M(n-i),i=-degree(ope,N)..-ldegree(ope,N)):
end:


#Sipur3with(n,M,L,K): All the expressions for the number of 
#<=3-rowed Young-tableaux with n cells, in terms of
#the Motzkin numbers M, whose [i,j] cell has m
#for m<=L and i*j<=m
#it all also prints out the first K terms
#For example, try:
#Sipur3with(n,M,3,20);
Sipur3with:=proc(n,M,L,K) local m,Y,i,j,lu:
print(`Let `, M(n), ` be the Motzkin numbers `):
print(`The following are (rigorously proved!) explicit expressions`):
print(`in terms of the Motzkin numbers M(n), for the number of 3-rowed`):
print(` n-celled `):
print(`Young-tableaux whose [i,j]-entry has m in it, for m between`):
print(`1 and `, L, `and the cell [i,j], with i<=3 (of course) and`):
print(` i*j<=m `):

for m from 1 to L do
for i from 1 to 3 do
for j from 1 to trunc(m/i) do
lu:=RegevStyle3YTwith([i,j],m,M,n):
if lu<>0 then
print(`The number of three-rowed n-celled Young tableaux Y such that`):
print(Y[i,j]=m, ` equals `):
print(lu):
print(`The first`, K, `non-zero terms are `):
print(TnuYTwithS(3,K,[i,j],m)):
fi:
od:
od:
od:



end:





Children:=proc(Lam)
local i, temp,Lam1,k:

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

temp:={}:

for i from 1 to nops(Lam)-1 do
 if Lam[i]>Lam[i+1] then
   Lam1:=[op(1..i-1,Lam),Lam[i]-1,op(i+1..nops(Lam),Lam)]:
    temp:=temp union {[Lam1,i]}:
 fi:
od:


k:=nops(Lam):

if Lam[k]>1 then
 Lam1:=[op(1..k-1,Lam),Lam[k]-1]:
    temp:=temp union {[Lam1,k]}:
 else

 Lam1:=[op(1..k-1,Lam)]:
    temp:=temp union {[Lam1,k]}:
fi:

temp:

end:

SYT:=proc(Lam)
local ch,i,n,j,temp,temp1:
option remember:

n:=convert(Lam,`+`):

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

temp:={}:
ch:=Children(Lam):

for i from 1 to nops(ch) do
temp1:=SYT(ch[i][1]):
temp:=temp union {seq(Ap(temp1[j],ch[i][2],n),j=1..nops(temp1))}:
od:
temp:
end:

Ap:=proc(Y,i,n) :
if i<=nops(Y) then
[op(1..i-1,Y),[op(Y[i]),n],op(i+1..nops(Y),Y)]:
else
[op(Y),[n]]:
fi:
end:

Shrink:=proc(lam) local i:
for i from 1 to nops(lam) do
 if lam[i]=0 then
  RETURN([op(1..i-1,lam)]):
 fi:
od:
lam:
end:

#SYTkn(k,n): the set of Standard Young Tableaux of <=k rows
#with n cells
SYTkn:=proc(k,n)
local gu,mu,lam:
option remember:
mu:=Pars(n,k):

gu:={}:
for lam in mu do
 
 gu:=gu union SYT(Shrink(lam)):
od:
gu:
end:

 


#YTwith(k,n,cell,m): the set of <=k-rowed Young-tableaux
#of n cells such that cell cell has m in it.
#for checking purposes)
#For example, try: YTwith(3,10,[1,2],3);
YTwith:=proc(k,n,cell,m) local gu,mu,Y,i,j:

i:=cell[1]: j:=cell[2]:

mu:=SYTkn(k,n):

gu:={}:

for Y in mu do
 if i<=nops(Y) and j<=nops(Y[i]) and Y[i][j]=m then
     gu:=gu union {Y}:
 fi:
od:

gu:
end: