#This is a short Maple package accompanying Manuel Kauers and Doron Zeilberger's article
#Counting Standard Young Tableaux With Restricted Runs 
#https://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/cyt.html
#The main procedures re SeqG(N) and SeqH(N) producing the two sequences discussed at the end of the paper

Help:=proc(): print(`G(a,b,c), SeqG(N), H(a,b,c), SeqH(N), Zinn(List) `): end:

###start G sequence 
#G(a,b,c): the number of walks from [0,0,0] to [a,b,c] with unit positive steps in the cubic lattice staying in x>=y>=z
#AND such that each run has length at least 2
G:=proc(a,b,c) option remember:
G1(a,b,c)+G2(a,b,c)+G3(a,b,c):
end:

G1:=proc(a,b,c) local a1:
option remember:

if a<0 or b<0 or c<0 or a<b or b<c then
0:
elif b=0 and c=0 and a>1 then
1:
else
add(G2(a1,b,c)+G3(a1,b,c) ,a1=b..a-2):
fi:
end:


G2:=proc(a,b,c) local b1:
option remember:

if a<0 or b<0 or c<0 or a<b or b<c then
0:
else
add(G1(a,b1,c)+G3(a,b1,c) ,b1=c..b-2):
fi:
end:


G3:=proc(a,b,c) local c1:
option remember:

if a<0 or b<0 or c<0 or a<b or b<c then
0:
else
add(G1(a,b,c1)+G2(a,b,c1) , c1=0..c-2):
fi:
end:

SeqG:=proc(N) local i:

[seq(G(i,i,i),i=1..N)]:

end:

######################End G seqeunce



###start H sequence 
#H(a,b,c): the number of walks from [0,0,0] to [a,b,c] with unit positive steps in the cubic lattice staying in x>=y>=z
#AND such that each run is ODD (i.e. NO even runs)
H:=proc(a,b,c) option remember:
H1(a,b,c)+H2(a,b,c)+H3(a,b,c):
end:

H1:=proc(a,b,c) local r1:
option remember:

if a<0 or b<0 or c<0 or a<b or b<c then
0:
elif b=0 and c=0 and a mod 2 =0 then
0:
elif b=0 and c=0 and a mod 2 =1 then
1:
else
add(H2(a-2*r1-1,b,c)+H3(a-2*r1-1,b,c) , r1=0..trunc((a-b)/2) ):
fi:
end:



H2:=proc(a,b,c) local r1:
option remember:

if a<0 or b<0 or c<0 or a<b or b<c then
0:
else
add(H1(a,b-2*r1-1,c)+H3(a,b-2*r1-1,c) , r1=0..trunc((b-c)/2) ):
fi:
end:



H3:=proc(a,b,c) local r1:
option remember:

if a<0 or b<0 or c<0 or a<b or b<c then
0:
else
add(H1(a,b,c-2*r1-1)+H2(a,b,c-2*r1-1) , r1=0..trunc(c/2) ):
fi:
end:



SeqH:=proc(N) local i:

[seq(H(i,i,i),i=1..N)]:

end:

######################End H seqeunce


sn:=proc(resh,n1):
-1/log(op(n1+1,resh)*op(n1-1,resh)/op(n1,resh)^2):
end:
 
#Zinn(resh): Zinn-Justin's method to estimate
#the C,mu, and theta such that
#resh[i] is appx. Const*mu^i*i^theta
#For example, try:
#Zinn([seq(5*i*2^i,i=1..30)]);
Zinn:=proc(resh)
local s1,s2,theta,mu,n1,i:
if nops({seq(sign(resh[i]),i=1..nops(resh))})<>1 then
 RETURN(FAIL):
fi:

n1:=nops(resh)-1:
s1:=sn(resh,n1):
s2:=sn(resh,n1-1):
theta:=evalf(2*(s1+s2)/(s1-s2)^2):
mu:=evalf(sqrt(op(n1+1,resh)/op(n1-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))):
[theta,mu]:
end: