######################################################################
##HOOKER: Save this file as HOOKER   To use it, stay in the          #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read HOOKER : <Enter>                              #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: July 2010
with(combinat): 
print(`Created: July 2010`):
print(` This is HOOKER `):
print(`a Maple package to study the number of`):
print(`Young-tableaux in bounded hooks (and sum of the squares of f_lambda)`):
print(`It accompanies the paper `):
print(` "Refined Asymptotics and Linear Recurrences for the Numbers of `):
print(`Young Tableaux bounded in the hook H(k,l) for k+l<=5 "`):
print(`by Shalosh B. Ekhad and Amitai Regev,`):
print(`published in the Personal Journal of Ekhad and Zeilberger.`):
print(`http://www.math.rutgers.edu/~zeilberg/pj.html `):
lprint(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
lprint(``):
 print(`The most current version of this  package `):
 print(` is  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/tokhniot/HOOKER  .`):
 print(`For a list of the procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):




####Begin AsyRec
with(numtheory):
with(SolveTools):

ezra1AsyRec:=proc()
if args=NULL then

print(` `):
print(`The supporting procedures are: Asy1, Asy1special `):
print(` Atom, CODV, CODV1, Finda, FindExpP, FindExpP1, FindExpP1g`):
print(` FindKk, Ksect,  `):
print(` NakedStirling, NewOpe, Nor, OneStepG, OpeBin,`):
print(`OpePer, OpePerG,OpePerN `):
print(` SeqFromRec `):
print():

else
 ezra(args):
fi:
end:
 
ezraAsyRec:=proc()
if args=NULL then

print(` AsyRec: A Maple package for finding the asymptotics of`):
print(` solutions of (homog.) linear recurrence equations with polynomial`):
print(` coefficents using (a variant) of Birkhoff-Trjitzinsky method `):
print():
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures:  `):
print(`Asy, AsyBin, AsyC, AsyF, AsyFC, AsyPerm, SipurBin, SipurPerm `):
print():


elif nargs=1 and args[1]=Asy then
print(`Asy(ope,n,N,M): the asymptotic expansion of solutions `):
print(`to ope(n,N)f(n)=0,where ope(n,N) is  a recurrence operator`):
print(`up to the M's term`):
print(`For example, try:`):
print(`Asy((n+2)^3*N^2-(34*n^3+153*n^2+231*n+117)*N+(n+1)^3,n,N,3);`):
print(`Asy((n+2)^2*N^2-(7*n^2+21*n+16)*N-8*(n+1)^2,n,N,3);`):
print(`Asy(N^2-N-(n+1),n,N,3);`):
print(`Asy(N^3-N^2-(n+2)*N-(n+2)*(n+1),n,N,3);`):

elif nargs=1 and args[1]=Asy1 then
print(`Asy1(ope,n,N,M): the asymptotic expansion of solutions `):
print(`to ope(n,N)f(n)=0,where ope(n,N) is  a recurrence operator`):
print(`up to the M's term, followed by the [K,k]`):
print(`such that the asymptotics is for a(Kn)/n!^k:`):
print(`For example, try:`):
print(`Asy1((n+2)^3*N^2-(34*n^3+153*n^2+231*n+117)*N+(n+1)^3,n,N,3);`):
print(`Asy1((n+2)^2*N^2-(7*n^2+21*n+16)*N-8*(n+1)^2,n,N,3);`):
print(`Asy1(N^2-N-(n+1),n,N,3);`):
print(`Asy1(N^3-N^2-(n+2)*N-(n+2)*(n+1),n,N,3);`):

elif nargs=1 and args[1]=Asy1special then
print(`Asy1special(ope,n,N,K,x): the asymptotic expansion of solutions `):
print(`to ope(n,N)f(n)=0,given as a list`):
print(`[pu,lu,expansion,r] where it is`):
print(`exp(pu)*lu^n*expansion(x), where x=1/n^(1/r), and r is`):
print(`a positive integer. It also returns [K,k] (see Asy1)`):
print(`where ope(n,N) is  a recurrence operator`):
print(`up to the K's term`):
print(`For example, try: `):
print(`Asy1special((n+2)^2*N^2-(7*n^2+21*n+16)*N-8*(n+1)^2,n,N,3,x);`):
print(`Asy1special(N^2-N-(n+1),n,N,3,x);`):

elif nargs=1 and args[1]=AsyBin then
print(`AsyBin(F,n,k,M,L): the asymptotics in n, up to order M,`):
print(`of a binomial coefficient`):
print(`sum Sum(F(n,k),k=0..n) (assuming F is supported between 0 and n`):
print(`where F is a hypergeometric term in n and k, `):
print(`and L is the number of terms in the sequence`):
print(`for estimating the constant in front`):
print(`For example, try:`):
print(`AsyBin(binomial(n,k)*binomial(n+k,k),n,k,5,1000);`):

elif nargs=1 and args[1]=AsyC then
print(`AsyC(ope,n,N,M,Ini,K): the asymptotic expansion of solutions `):
print(`to ope(n,N)f(n)=0, with the given initial conditions`):
print(`where ope(n,N) is  a recurrence operator`):
print(`up to the M's term`):
print(`and complete with an empirically derived constant in front`):
print(`using K terms `):
print(`For example, try:`):
print(`AsyC((n+2)^3*N^2-(34*n^3+153*n^2+231*n+117)*N+(n+1)^3,n,N,5,[5,73],1000);`):
print(`AsyC((n+2)^2*N^2-(7*n^2+21*n+16)*N-8*(n+1)^2,n,N,5,[2,10],1000);`):
print(`AsyC(N^2-N-(n+1),n,N,5,[1,2],1000);`):

elif nargs=1 and args[1]=AsyF then
print(`AsyF(ope,n,N,M): the asymptotic expansion of solutions `):
print(`to ope(n,N)f(n)=0,where ope(n,N) is  a recurrence operator`):
print(`up to the M's term, in terms of factorials`):
print(`For example, try:`):
print(`AsyF((n+2)^3*N^2-(34*n^3+153*n^2+231*n+117)*N+(n+1)^3,n,N,3);`):
print(`AsyF((n+2)^2*N^2-(7*n^2+21*n+16)*N-8*(n+1)^2,n,N,3);`):
print(`AsyF(N^2-N-(n+1),n,N,3);`):
print(`AsyF(N^3-N^2-(n+2)*N-(n+2)*(n+1),n,N,3);`):

elif nargs=1 and args[1]=AsyFC then
print(`AsyFC(ope,n,N,M,Ini,K): the asymptotic expansion of solutions `):
print(`to ope(n,N)f(n)=0, with the given initial conditions`):
print(`where ope(n,N) is  a recurrence operator`):
print(`up to the M's term, in terms of the factorial function`):
print(`and complete with an empirically derived constant in front`):
print(`using K terms `):
print(`For example, try:`):
print(`AsyFC((n+2)^3*N^2-(34*n^3+153*n^2+231*n+117)*N+(n+1)^3,n,N,5,[5,73],1000);`):
print(`AsyFC((n+2)^2*N^2-(7*n^2+21*n+16)*N-8*(n+1)^2,n,N,5,[2,10],1000);`):
print(`AsyFC(N^2-N-(n+1),n,N,5,[1,2],1000);`):

elif nargs=1 and args[1]=AsyPerm then
print(`AsyPerm(n,r,M,L): the asymptotics for the number of permutations`):
print(`whose r-th power is the identity permutation. `):
print(`M is the desired order, and L is the number of terms in the`):
print(`sequence used to estimate the constant`):
print(`For example, try AsyPerm(n,2,5,1000);`):

elif nargs=1 and args[1]=Atom then
print(`Atom(s,r,i,x,K): Expanding (n+i)^(s/r)-n^(s/r) in terms of`):
print(`x=n^(-1/r) up to K terms`):
print(`For example try:`):
print(`Atom(2,3,2,x,3);`):

elif nargs=1 and args[1]=CODV then
print(`CODV(ope,n,N,k,K): Given a linear recurrence operator ope(n,N)`):
print(`annihilating a[n], say, outputs the operator`):
print(`annihilating b[n]:=a[n*K]/n!^k . `):
print(`For example, try: CODV(N-(n+1),n,N,1,1):`):

elif nargs=1 and args[1]=CODV1 then
print(`CODV1(ope,n,N,k): Given a linear recurrence operator ope(n,N)`):
print(`annihilating a[n], say, outputs the operator`):
print(`annihilating b[n]:=a[n]/n!^k`):
print(`For example, try: CODV1(N-(n+1),n,N,1):`):


elif nargs=1 and args[1]=Finda then
print(`Finda(ope,N,x,r): finds the first power x^a in the`):
print(`asymptotic solution of ope(N,n)f(n)=0, where x=1/n^(1/r)`):
print(`For example, try:`):
print(`Finda((1+x)-(1+3*x)*N,N,x,1);`):


elif nargs=1 and args[1]=FindExpP then
print(`FindExpP1(ope,n,N): finds the exponential part of the asymptotics`):
print(`for the solutions of ope(n,N)a(n)=0 if it is normalized`):
print(`such that the leading asymp. is 1^n`):
print(`for example, try:`):
print(`FindExpP((n+1)*N^2-2*(n+5)*N+n+3,n,N);`):

elif nargs=1 and args[1]=FindExpP1 then
print(`FindExpP1(ope,n,N,r,x): finds the exponential part of the asymptotics`):
print(`in terms of x=n^(1/r)`):
print(`for the solutions of ope(n,N)a(n)=0 if it is of type r.`):
print(`for example, try:`):
print(`FindExpP1((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x);`):

elif nargs=1 and args[1]=FindExpP1g then
print(`FindExpP1g(ope,n,N,r,x,ds): finds the exponential part of the`):
print(` asymptotics`):
print(`in terms of x=n^(1/r) as a poly. of degree ds in x`):
print(`for the solutions of ope(n,N)a(n)=0 if it is of type r`):
print(`For example, try:`):
print(`FindExpP1g((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x,1);`):

elif nargs=1 and args[1]=FindKk then
print(`FindKk(ope,n,N): Given a linear recurrence operator with polynomial`):
print(`coefficients, ope(n,N), finds the integers K and k such`):
print(`that if a(n) is a solution of ope(n,N)a(n)=0, then`):
print(`b(n):=a(K*n)/n!^k is annihilated by a standard operator`):
print(`the output is the pair [k,K] and the transformed operator`):
print(`For exanple, try: FindKk(N^2-n,n,N);`):


elif nargs=1 and args[1]=Ksect then
print(`Ksect(ope,n,N,k,r): Given a linear recurrence operator ope(n,N)`):
print(`annihilating a sequence a[n], and pos. integer k, and`):
print(`non-neg. integer r, outputs the one, of the same`):
print(`order, that annihilates  a[k*n+r]`):
print(`For example try: Ksect(N^2-N-(n+1),n,N,2,0);`):

elif nargs=1 and args[1]=NakedStirling then 
print(`NakedStirling(n,K): the asymptotic expansion of`):
print(`n!/((n/e)^n*sqrt(2*Pi*n). For example, try:`):
print(`NakedStirling(n,5);`):


elif nargs=1 and args[1]=NewOpe then
print(`NewOpe(ope,n,N,K,x): the exponential part + the transformed operator`):
print(`(up to degree K asymp. in the coefficients), in terms of x=1/n^(1/r)`):
print(`For example, try:`):
print(`NewOpe((n+1)*N^2-2*(n+5)*N+n+3,n,N,4,x);`):

elif nargs=1 and args[1]=Nor then
print(`Nor(ope,n,N): the Normalizer of the linear recurrence`):
print(`operator with polynomial coefficients ope(n,N)`):
print(`followed by its exponential growth constant`):
print(`For example, try:`):
print(`Nor(N^2-3*N+2,n,N);`):

elif nargs=1 and args[1]=OneStepG then
print(`OneStepG(ope,N,x,r,Cu): Given the asymptotic expansion of`):
print(`a solution to ope(n,N) expressed in terms of x=1/n^(1/r)`):
print(`finds the next one`):
print(`finds one more term`):
print(`OneStepG((1+x)-(1+3*x)*N,N,x,1,x^2);`):

elif nargs=1 and args[1]=OpeBin then
print(`OpeBin(n,N,r): the operator for the sum of the r-th power of the`):
print(`binomial coefficients followed by their initial conditions`):
print(`For example, try:`):
print(`OpeBin(n,N,3);`):

elif nargs=1 and args[1]=OpePer then
print(`OpePer(N,n,r): the linear recurrenec operator that`):
print(`annihilates a(n,r):=The number of permutations`):
print(`with cycles of length<=r. For example, try:`):
print(`OpePer(N,n,2);`):

elif nargs=1 and args[1]=OpePerG then
print(`OpePerG(N,n,S): the linear recurrenec operator that`):
print(`annihilates a(n,S):=The number of permutations`):
print(`with cycles of lengths in S. It also outputs`):
print(`the initial conditions `):
print(`For example, try:`):
print(`OpePerG(N,n,{1,2});`):

elif nargs=1 and args[1]=OpePerN then
print(`OpePerN(N,n,r): the normalized linear recurrenec operator that`):
print(`annihilates a(n,r):=The number of permutations`):
print(`with cycles of length<=r. For example, try:`):
print(`OpePerN(N,n,2);`):

elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends it to the first K values`):
print(`For example, try:`):
print(`SeqFromRec(N-2,n,N,[1],10);`):
print(`SeqFromRec((n+2)^3*N^2-(34*n^3+153*n^2+231*n+117)*N+(n+1)^3,n,N,[5,73],20);`):
print(`SeqFromRec((n+2)^2*N^2-(7*n^2+21*n+16)*N-8*(n+1)^2,n,N,[2,10],20);`):
print(`SeqFromRec(N^2-N-(n+1),n,N,[1,2],20);`):

elif nargs=1 and args[1]=SipurBin then
print(`SipurBin(R,n,M,L): The story for sums of powers of binomial`):
print(`coefficients. For example, try:`):
print(`SipurBin(4,n,5,1000);`):


elif nargs=1 and args[1]=SipurPerm then
print(`SipurPerm(R,n,M,L): The story for the asymptotics for`):
print(`the number of premutations pi  of [1,n] such that pi^r=Identity`):
print(`for r from 2 to R (r=2 is involutions).`):
print(`M is the desired order and L is the number of terms in the`):
print(`squence for estimating the constant`):
print(`(for some reason, it only works for R<=6).`):
print(`For example, try:`):
print(`SipurPerm(6,n,5,1000);`):

else

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

fi:


end:




#MonoN(ope,n,N,k): Given a linear recurrence operator with
#poly coeffs. ope(n,N), and a pos. integer k, outputs the expression of
#N^k as a linear combination of 1, N, ..., N^(ORDER-1)
#For example, try
#MonoN(N-n-1,n,N,3);
MonoN:=proc(ope,n,N,k) local ORDER,coe0,i,lu1,lu2:
ORDER:=degree(ope,N):

if k<ORDER then
 RETURN(N^k):
fi:

if k=ORDER then
coe0:=coeff(ope,N,ORDER):
RETURN(add(normal(-coeff(ope,N,i)/coe0)*N^i,i=0..ORDER-1)):
fi:

lu1:=expand(MonoN(ope,n,N,k-1)):
lu2:=expand(MonoN(ope,n,N,ORDER)):
lu1:=expand(subs(n=n+1,lu1)*N):
lu1:=subs(N^ORDER=lu2,lu1):
add(normal(coeff(lu1,N,i))*N^i,i=0..ORDER-1):

end:


#Ksect(ope,n,N,k,r): Given a linear recurrence operator ope(n,N)
#annihilating a sequence a[n], and pos. integer k, and
#non-neg. integer r, outputs the one, of the same
#order, that annihilates  a[k*n+r]
#For example try: Ksect(N^2-N-(n+1),n,N,2,0);
Ksect:=proc(ope,n,N,K,r) local ORDER,eq,var,lu,mu,c,Ope,var1,i,vu,v:
 
ORDER:=degree(ope,N):

for i from 0 to ORDER do
 lu[i]:=MonoN(ope,n,N,r+K*i):
 lu[i]:=subs(n=n*K,lu[i]):
od:

Ope:=add(c[i]*N^i,i=0..ORDER):

mu:=expand(add(c[i]*lu[i],i=0..ORDER)):

eq:={seq(coeff(mu,N,i),i=0..ORDER)}:
var:={seq(c[i],i=1..ORDER)}:

var1:=solve(eq,var):


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



vu:={}:

for i from 1 to nops(var1) do
 if op(1,var1[i])=op(2,var1[i]) then
    vu:=vu union {op(1,var1[i])}:
 fi:
od:

Ope:=subs(var1,Ope):

for v in vu do
Ope:=subs(v=0,Ope):
od:



Ope:=subs(c[0]=1,numer(normal(Ope))):



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
#For example, try:
#SeqFromRec(N-2,n,N,[1],10);
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1,kap:
ope1:=Yafe(ope,N)[2]:



kap:={solve(denom(normal(ope1)),n)}:

kap:=evalf(max(op(kap))):


L:=degree(ope1,N):

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

if nops(Ini)<=kap then
RETURN(FAIL):
fi:
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:
 



#HafelOper(ope,n,N,L): applies the operator ope(n,N) to
#the sequence L, for example, try:
#HafelOper(N-(n+1),n,N,[1,2,6,24,120]);
HafelOper:=proc(ope,n,N,L) local i,n1:
[seq(add(subs(n=n1,coeff(ope,N,i))*L[n1+i],i=0..degree(ope,N)),
n1=1..nops(L)-degree(ope,N))]:
end:


#TestKsect(ope,n,N,K,r,M,Ini): tests procedure Ksect
#with initial conditions Ini up to M terms
TestKsect:=proc(ope,n,N,K,r,M,Ini)  local gu,Ope,n1:
gu:=SeqFromRec(ope,n,N,Ini,M*K+r):
Ope:=Ksect(ope,n,N,K,r):
gu:=[seq(gu[K*n1+r],n1=1..M-1)]:
evalb({op(HafelOper(Ope,n,N,gu))}={0}):
end:

rf:=proc(a,n) local i:mul(a+i,i=0..n-1):end:

#CODV1(ope,n,N,k): Given a linear recurrence operator ope(n,N)
#annihilating a[n], say, outputs the operator
#annihilating b[n]:=a[n]/n!^k
#For example, try: CODV1(N-(n+1),n,N,1):
CODV1:=proc(ope,n,N,k) local i:
add(coeff(ope,N,i)*(rf(n+1,i))^k*N^i,i=0..degree(ope,N)):
end:


#CODV(ope,n,N,k,K): Given a linear recurrence operator ope(n,N)
#annihilating a[n], say, outputs the operator
#annihilating b[n]:=a[n*K]/n!^k
#For example, try: CODV(N-(n+1),n,N,1,1):
CODV:=proc(ope,n,N,k,K) local Ope:
Ope:=Ksect(ope,n,N,K,0):
CODV1(Ope,n,N,k):
end:



#FindKk(ope,n,N): Given a linear recurrence operator with polynomial
#coefficients, ope(n,N), finds the integers K and k such
#that if a(n) is a solution of ope(n,N)a(n)=0, then
#b(n):=a(K*n)/n!^k is annihilated by a standard operator
#the output is the pair [k,K] and the transformed operator
#For exanple, try: FindKk(N^2-n,n,N);
FindKk:=proc(ope,n,N) local ope1,Ld,deg,sp,yakhas,K,k,pu:


pu:=Aope1(ope,n,N):

if {solve(pu,N)}<>{} and {solve(pu,N)}<>{0} then
 RETURN([1,0],ope):
fi:

ope1:=expand(ope):
Ld:=ldegree(ope,N):
ope1:=expand(subs(n=n-Ld,ope1)/N^Ld ):

deg:=degree(ope,N):

sp:=degree(coeff(ope,N,0),n)-degree(coeff(ope,N,deg),n):

yakhas:=sp/deg:

k:=numer(yakhas):
K:=denom(yakhas):
[K,k],CODV(ope,n,N,k,K):
end:


#Aope1(ope,n,N): the asymptotics of the difference operator with poly. coeffs. ope(n,N)
Aope1:=proc(ope,n,N)
local gu:
gu:=expand(numer(normal(ope))):
gu:=coeff(gu,n,degree(gu,n)):
gu:=expand(gu/coeff(gu,N,degree(gu,N))):
factor(gu):
end:

C3:=proc(a,b,c) option remember:

if [a,b,c]=[0,0,0] then
  RETURN(1):
fi:
if a<0 or b<0 or c<0 then
   RETURN(0):
fi:
C3(a-1,b,c)+C3(a,b-1,c)+C3(a,b,c-1)-4*C3(a-1,b-1,c-1):
end:

#Aluf1(pit): is it positive dominant?
Aluf1:=proc(pit)
local aluf,shi,i,makom:


shi:=evalf(abs(pit[1])):
aluf:={evalf(evalc(pit[1]))}:
makom:={1}:

for i from 2 to nops(pit) do

 if evalf(abs(pit[i]))>evalf(shi) then
   shi:=abs(evalf(evalc(pit[i]))):
  aluf:={pit[i]}:
  makom:={i}:
 elif evalf(abs(pit[i]))=shi then
   aluf:=aluf union {evalf(evalc(pit[i]))}:
   makom:=makom union {i}:
 fi:

od:

aluf,shi,makom:

end:



#OneStepAS1(ope1,n,N,alpha,f,S1): Given the partial asymptotic
#expansion of solutions of ope1(n,N) with exponent alpha
#extends it by one term
OneStepAS1:=proc(ope1,n,N,alpha,f,S1)
local x,f1,L,F,A,mu,i,A1:

f1:=subs(n=1/x,f):
L:=degree(f1,x):
F:=f1+A*x^(L+S1):
mu:=add(
subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha*subs(x=x/(1+i*x),F)
,i=ldegree(ope1,N)..degree(ope1,N)):
mu:=normal(mu):
mu:=taylor(mu,x=0,L+11):
mu:=simplify(mu):
for i from L+1 to L+9 while coeff(mu,x,i)=0 do od:
if i=L+10 then
 RETURN(FAIL):
fi:
mu:=coeff(mu,x,i):
A1:=subs(Linear({mu},{A}),A):
if A1=0 then
  RETURN(FAIL):
fi:


subs({x=1/n,A=A1},F):

end:


#OneStepA(ope1,n,N,alpha,f): Given the partial asymptotic
#expansion of solutions of ope1(n,N) with exponent alpha
#extends it by one term
OneStepA:=proc(ope1,n,N,alpha,f)
local S1,pu:
for S1 from 1 to 5 while OneStepAS1(ope1,n,N,alpha,f,S1)=FAIL do od:
pu:=OneStepAS1(ope1,n,N,alpha,f,S1):
if pu=FAIL then
 RETURN(f):
else
 RETURN(pu):
fi:
end:



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:






#Nor(ope,n,N): the Normalizer of the linear recurrence
#operator with polynomial coefficients ope(n,N)
#followed by its exponential growth constant
#For example, try:
#Nor(N^2-3*N+2,n,N);
Nor:=proc(ope,n,N)
local gu,pit,lu,makom,ope1:

gu:=Aope1(ope,n,N):
if degree(gu,N)=0 then
   print(`Not of exponential growth, case not implemented`):
   RETURN(FAIL):
fi:

if not type(expand(normal(ope/gu)),`+`) then
 pit:=[solve(gu,N)]:
 makom:=Aluf1(pit)[3]:

lu:=pit[makom[1]]:
ope1:=Yafe(subs(N=lu*N,ope),N)[2]:
RETURN(ope1,lu):

fi:

pit:=[solve(gu,N)]:
makom:=Aluf1(pit)[3]:

if nops(makom)<>1 then
lu:=evalf(evalc(pit[makom[1]])):
 if coeff(lu,I,1)>(1/10)^(Digits-2) then
 print(`Dominant roots are complex`):
 RETURN(FAIL):
  elif pit[makom[1]]+pit[makom[2]]=0 then
   print(`Dominant real roots are negatives of each other`):
   RETURN(FAIL):
 fi:
fi:


lu:=evalf(evalc(pit[makom[1]])):

if coeff(lu,I,1)>(1/10)^(Digits-2) then
 print(`Dominant root is complex`):
 RETURN(FAIL):
fi:

lu:=pit[makom[1]]:

ope1:=Yafe(subs(N=lu*N,ope),N)[2]:
ope1,lu:
end:









#Atom(s,r,i,x,K): Expanding (n+i)^(s/r)-n^(s/r) in terms of
#x=n^(-1/r) up to K terms
#For example try:
#Atom(2,3,2,x,3);
Atom:=proc(s,r,i,x,K) local gu,y,i1:
gu:=(1+i*y)^(s/r)-1:
gu:=taylor(gu,y=0,K+1):
gu:=add(coeff(gu,y,i1)*y^i1,i1=0..K):
expand(subs(y=x^r,gu)/x^s):
end:




#FindExpP1(ope,n,N,r,x): finds the exponential part of the asymptotics
#in terms of x=n^(1/r)
#for the solutions of ope(n,N)a(n)=0 if it is of type r
#For example, try:
#FindExpP1((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x);
FindExpP1old:=proc(ope,n,N,r,x) local c,eq,var,s,sof,gu,ope1,deg,i,i1,v:
sof:=add(c[s]*x^s,s=1..r-1):
var:={seq(c[i],i=1..r-1)}:

deg:=degree(ope,n):
ope1:=expand(ope/n^deg):
ope1:=expand(subs(n=1/x^r,ope1)):

gu:=0:

for i from 0 to degree(ope1,N) do
 gu:=gu+coeff(ope1,N,i)*exp(add(c[s]*Atom(s,r,i,x,2),s=1..r-1) ):
od:

gu:=taylor(gu,x=0,2*r-1):



eq:={seq(coeff(gu,x,i1),i1=r..2*r-2)}:


var:=[solve(eq,var)][1]:
for v  in var do
 if op(1,v)=op(2,v) then
  RETURN(FAIL):
 fi:
od:

subs(var,sof):

end:







#Findr(ope,n,N): Given an operator ope(n,N)
#such that the leading terms in n is a multiple of
#(N-1), finds the largest r such that it is
# a multiple of (N-1)^r. For example, try:
#Findr((N-1)^4,n,N);
Findr:=proc(ope,n,N) local ope1,r:
ope1:=coeff(expand(ope),n,degree(ope,n)):
if expand(subs(N=1,ope1))<>0 then
  RETURN(FAIL):
fi:

for r from 1 while expand(subs(N=1,diff(ope1,N$r)))=0 do od:
r:
end:


FindExpP:=proc(ope,n,N)  local r,x,lu:
 r:=Findr(ope,n,N):


  if r=FAIL then
    RETURN(FAIL):
  fi:
 lu:=FindExpP1(ope,n,N,r,x) :
  if lu=FAIL then
    RETURN(FAIL):
   fi:
  subs(x=n^(1/r),lu):
end:




#NewOpe(ope,n,N,K): the exponential part + the transformed operator
#(up to degree K asymp. in the coefficients)
#For example, try:
#NewOpe((n+1)*N^2-2*(n+5)*N+n+3,n,N,4);
NewOpe:=proc(ope,n,N,K,x)  local r,lu,lu1,ope1,deg,s,i,gu,mu:
 r:=Findr(ope,n,N):
  if r=FAIL then
    RETURN(FAIL):
  fi:
 lu:=FindExpP1(ope,n,N,r,x) :
 lu1:=FindExpP(ope,n,N) :

  if lu=FAIL then
    RETURN(FAIL):
   fi:

deg:=degree(ope,n):
ope1:=expand(ope/n^deg):
ope1:=subs(n=1/x^r,ope1):
ope1:=expand(ope1):


gu:=0:

for i from 0 to degree(ope1,N) do
  mu:=0:
 for s from 1 to r-1 do
     mu:=mu+coeff(lu,x,s)*Atom(s,r,i,x,K+3):
  od:
     
 gu:=gu+coeff(ope1,N,i)*exp(mu)*N^i:
od:

gu:=taylor(gu,x=0,K+3):
lu1,add(coeff(gu,x,i)*x^i,i=0..K):
end:



#Bdok(ope,n,N,K): the exponential part + the transformed operator
#(up to degree K asymp. in the coefficients)
#For example, try:
#Bdok((n+1)*N^2-2*(n+5)*N+n+3,n,N,4);
Bdok:=proc(ope,n,N,K)  local r,x,lu,lu1,ope1,deg,s,i,gu,mu:
 r:=Findr(ope,n,N):
  if r=FAIL then
    RETURN(FAIL):
  fi:
 lu:=FindExpP1(ope,n,N,r,x) :
 lu1:=FindExpP(ope,n,N) :
  if lu=FAIL then
    RETURN(FAIL):
   fi:

deg:=degree(ope,n):
ope1:=expand(ope/n^deg):
ope1:=subs(n=1/x^r,ope1):
ope1:=expand(ope1):


gu:=0:

for i from 0 to degree(ope1,N) do
  mu:=0:
 for s from 1 to r-1 do
     mu:=mu+coeff(lu,x,s)*Atom(s,r,i,x,K+3):
  od:
     
 gu:=gu+coeff(ope1,N,i)*exp(mu)*N^i:
od:

gu:=taylor(gu,x=0,K+3):
lu1,add(coeff(gu,x,i)*x^i,i=0..K):
end:




OpePerN:=proc(N,n,r) local i,j,ope:
ope:=1-add(N^(-i)*mul(n-j,j=1..i-1),i=1..r):
ope:=expand(N^r*subs(n=n+r,ope)):
ope:=add(factor(coeff(ope,N,i))*N^i,i=0..r):
ope:=expand(FindKk(ope,n,N)[2]):
ope:=add(factor(coeff(ope,N,i))*N^i,i=0..r):
ope:=numer(Nor(ope,n,N)[1]):
ope:=add(factor(coeff(ope,N,i))*N^i,i=0..r):
end:

OpePer:=proc(N,n,r) local i,j,ope:
ope:=1-add(N^(-i)*mul(n-j,j=1..i-1),i=1..r):
ope:=expand(N^r*subs(n=n+r,ope)):
end:

OpePerG:=proc(N,n,S) local i,j,ope,r,gu,t:
ope:=1-add(N^(-i)*mul(n-j,j=1..i-1),i in S):
r:=max(op(S)):
ope:=expand(N^r*subs(n=n+r,ope)):

gu:=exp(add(t^i/i,i in S)):
gu:=taylor(gu,t=0,degree(ope,N)+2):

ope,[seq(coeff(gu,t,i)*i!,i=1..degree(ope,N))]:
end:


#Finda(ope,N,x,r): finds the first power x^a in the
#asymptotic solution of ope(N,n)f(n)=0, where x=1/n^(1/r)
#For example, try:
#Finda((1+x)-(1+3*x)*N,N,x,1);
Finda:=proc(ope,N,x,r) local gu,i,a,a1:

gu:=add(coeff(ope,N,i)*(1+i*x^r)^(a/r),i=0..degree(ope,N)):

gu:=taylor(gu,x=0,5*r+5):
gu:=expand(add(coeff(gu,x,i)*x^i,i=0..5*r+4)):

for i from 0 to 5*r+4 while expand(simplify(coeff(gu,x,i)))=0 do od:

if i=5*r+5 then
  RETURN(FAIL):
fi:


a1:=[solve(coeff(gu,x,i),a)][1]:


if a1=NULL then
 RETURN(FAIL):
elif not type(a1,integer) then
 RETURN(FAIL):

else
 RETURN(-a1):
fi:
end:






#OneStepG(ope,N,x,r,Cu): Given a partial
#asympt. expansion, in terms of x=1/n^(1/r),
#for solutions of the recurrence equation
#ope(n,N)a(n)=0, where ope is normalized of type
#r (i.e. its leading coeff. is a multiple of (N-1)^r)
#finds the next term
#OneStepG((1+x)-(1+3*x)*N,N,x,1,x^2);
OneStepG:=proc(ope,N,x,r,Cu) local gu,i,a,c,c1,Cu1,K,ka:
K:=degree(Cu,x):
Cu1:=Cu+c*x^(K+1):
gu:=add(coeff(ope,N,i)*
add(coeff(Cu1,x,a)*x^a*(1+i*x^r)^(-a/r),a=ldegree(Cu1,x)..degree(Cu1,x)),
i=0..degree(ope,N)):
gu:=expand(gu):

gu:=taylor(gu,x=0,5*r+8+K):
gu:=
simplify(expand(add(coeff(gu,x,i)*x^i,i=0..5*r+7+K))):

gu:=pashet(gu,x,5*r+7+K,c):
ka:=expand(simplify(coeff(gu,x,ldegree(gu,x)))):
ka:=simplify(ka):

c1:=[solve(simplify(coeff(gu,x,ldegree(gu,x))),c)][1]:


if c1=NULL then
 RETURN(FAIL):
else
 RETURN(subs(c=c1,Cu1)):
fi:
end:



#Asy1(ope,n,N,K): the asymptotic expansion of solutions 
#to ope(n,N)f(n)=0,where ope(n,N) is  a recurrence operator
#up to the K's term
Asy1:=proc(ope,n,N,K)
local gu,lu,alpha,mu,ope1,ku,i,f,x,ka,vu,ope1A,r,Kk,pu,a:


gu:=FindKk(ope,n,N):
Kk:=gu[1]:

ope1:=gu[2]:

vu:=Nor(ope1,n,N):

if vu=FAIL then
  RETURN(FAIL):
fi:

ope1:=vu[1]:
lu:=vu[2]:

ope1A:=Aope1(ope1,n,N):



for r from 1 while subs(N=1,diff(ope1A,N$r))=0 do od:


if r=1 then

###r=1 case

mu:=add(
subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha,i=ldegree(ope1,N)..degree(ope1,N)):
mu:=normal(mu):
ku:=factor(coeff(taylor(mu,x=0,2),x,1)):


ka:=simplify([solve(ku,alpha)]):


if coeff(ka[1],I,1)<>0 then
  RETURN(FAIL):
fi:

alpha:=max(op(ka)):


if normal(subs(N=1,ope1))=0 then
 RETURN(lu^n*n^alpha,Kk):
fi:

f:=1:
for i from 1 to K do
f:=OneStepA(ope1,n,N,alpha,f):
od:

if f=FAIL then
 RETURN(FAIL):
fi:

RETURN(lu^n*n^alpha*f,Kk):

#end r=1 case

else



ope1:=numer(ope1):


if expand(diff(ope1,n))=0 then
 ka:=n^(r-1)*lu^n:
 RETURN(ka,Kk):
fi:

 pu:=FindExpP(ope1,n,N):

if pu=FAIL then
  RETURN(FAIL):
fi:

ope1:=NewOpe(ope1,n,N,K+5,x)[2]:

a:=Finda(ope1,N,x,r):

if a=FAIL then
 RETURN(FAIL):
fi:
ka:=x^a:

for i from a to K do
 ka:=OneStepG(ope1,N,x,r,ka):
od:

ka:=exp(pu)*subs(x=1/n^(1/r),ka)*lu^n:

RETURN(ka,Kk):
fi:
end:



pashet:=proc(gu,x,K,c) local gu1,i,pu:

Digits:=300:

gu1:=0:

for i from 0 to K do
 pu:=evalf(coeff(gu,x,i)):
if not (abs(coeff(pu,c,1))<10^(-20)  and abs(coeff(pu,c,0))<10^(-20) )
then
  gu1:=gu1+coeff(gu,x,i)*x^i:
fi:
od:
gu1:
end:


  
 
#NakedStirling(n,K): the asymptotic expansion of
#n!/((n/e)^n*sqrt(2*Pi*n). For example, try:
#NakedStirling(n,5);
NakedStirling:=proc(n,K) local gu,i:
gu:=n!/(n/exp(1))^n/sqrt(2*Pi*n):
gu:=expand(asympt(gu,n,K+1)):

(n/exp(1))^n*sqrt(2*Pi*n), add(coeff(op(i,gu),n,-(i-1))/n^(i-1),i=1..K+1):
end:


#AsyF(ope,n,N,M): the asymptotic expansion of solutions 
#to ope(n,N)f(n)=0,where ope(n,N) is  a recurrence operator
#up to the M's term, in terms of the factorial function
AsyF:=proc(ope,n,N,M) local K,k,gu:
gu:=Asy1(ope,n,N,M):

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

K:=gu[2][1]:
k:=gu[2][2]:
gu:=gu[1]:
(n/K)!^k*subs(n=n/K,gu):
end:



#AsyFC(ope,n,N,M,Ini,K): the asymptotic expansion of solutions 
#to ope(n,N)f(n)=0, with the given initial conditions
#where ope(n,N) is  a recurrence operator
#up to the M's term, in terms of the factorial function
#and complete with an empirically derived constant in front
#using K terms
AsyFC:=proc(ope,n,N,M,Ini,K) local gu,L,i,mu,er,C,D1:
Digits:=100:
gu:=AsyF(ope,n,N,M):
L:=SeqFromRec(ope,n,N,evalf(Ini),K):
mu:=[seq(evalf(L[i]/subs(n=i,gu)),i=K-10..K)]:
er:=mu[nops(mu)]-mu[nops(mu)-1]:
D1:=-trunc(log(abs(er))/log(10))-3:
if D1<2 then
 print(`can't determine the constant`):
RETURN(gu):
fi:

C:=evalf(mu[nops(mu)],D1):
C,gu:
end:





###Asy

#Asy1special(ope,n,N,K,x): the asymptotic expansion of solutions 
#to ope(n,N)f(n)=0,given as a list
#[pu,lu,expansion,r] where it is
#exp(pu)*lu^n*expansion(x), where x=1/n^(1/r), and r is
#a positive integer. It also returns [K,k] (see Asy1)
#where ope(n,N) is  a recurrence operator
#up to the K's term
Asy1special:=proc(ope,n,N,K,x)
local gu,lu,alpha,mu,ope1,ku,i,f,ka,vu,ope1A,r,Kk,pu,a:


gu:=FindKk(ope,n,N):



Kk:=gu[1]:

ope1:=gu[2]:

vu:=Nor(ope1,n,N):


if vu=FAIL then
  RETURN(FAIL):
fi:

ope1:=vu[1]:
lu:=vu[2]:

ope1A:=Aope1(ope1,n,N):



for r from 1 while subs(N=1,diff(ope1A,N$r))=0 do od:


if degree(ope1,n)=0 then
RETURN([0,lu,x^(-r+1),1,1],Kk):
fi:


if r=1 then

###r=1 case

mu:=add(
subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha,i=ldegree(ope1,N)..degree(ope1,N)):
mu:=normal(mu):
ku:=factor(coeff(taylor(mu,x=0,2),x,1)):


ka:=simplify([solve(ku,alpha)]):


if coeff(ka[1],I,1)<>0 then
  RETURN(FAIL):
fi:

alpha:=max(op(ka)):


if normal(subs(N=1,ope1))=0 then
RETURN([0,lu,x^(-alpha),1,1],Kk):
fi:

f:=1:
for i from 1 to K do
f:=OneStepA(ope1,n,N,alpha,f):
od:

if f=FAIL then
 RETURN(FAIL):
fi:

f:=subs(n=1/x,f):
RETURN([0,lu,x^(-alpha),f,1],Kk):
#end r=1 case

else



ope1:=numer(ope1):


if expand(diff(ope1,n))=0 then
 ka:=n^(r-1)*lu^n:
 RETURN(ka,Kk):
fi:




 pu:=FindExpP(ope1,n,N):



if pu=FAIL then
  RETURN(FAIL):
fi:


ope1:=NewOpe(ope1,n,N,K+5,x)[2]:

a:=Finda(ope1,N,x,r):

if a=FAIL then
 RETURN(FAIL):
fi:
ka:=x^a:

for i from a to K do
 ka:=OneStepG(ope1,N,x,r,ka):
od:


RETURN([pu,lu,1,ka,r],Kk):

fi:
end:

#IsMispar(a): is a (generalized) numeric type
IsMispar:=proc(a) :
if type(a,numeric) then
 true:
elif type(a,`^`) and type(op(1,a),numeric) and type(op(2,a),numeric)  then
 true:
elif a=Pi  then
true:
elif type(a,`^`) and op(1,a)=Pi and type(op(2,a),numeric)  then
 true:
elif type(a,function) and type(op(1,a),numeric) then
 true:
else
 false:
fi:
end:

#GetRidOfConst(P): gets rid of the constant in front
GetRidOfConst:=proc(P) local i,P1:
if not type(P,`*`) then
 RETURN(P):
fi:

P1:=1:

for i from 1 to nops(P) do
 if not IsMispar(op(i,P)) then
  P1:=P1*op(i,P):
 fi:
od:
P1:
end:

#Asy(ope,n,N,M): the asymptotic expansion of solutions 
#to ope(n,N)f(n)=0,where ope(n,N) is  a recurrence operator
#up to the M's term
Asy:=proc(ope,n,N,M) local K,k,gu,vu,vu1,vu2,pu,lu,ka1,ka2,r,x,i,vu2a:
gu:=Asy1special(ope,n,N,M,x):
if gu=FAIL then
 RETURN(FAIL):
fi:

K:=gu[2][1]:
k:=gu[2][2]:


gu:=gu[1]:
pu:=subs(n=n/K,gu[1]):
lu:=gu[2]:
ka1:=gu[3]:
ka2:=gu[4]:
r:=gu[5]:
ka2:=subs(x=K^(1/r)*x,ka2):
vu:=NakedStirling(n,M+1):
vu1:=vu[1]:
vu2:=vu[2]:


vu1:=subs(n=n/K,vu1)^k:
vu2:=subs(n=n/K,vu2)^k:
vu2:=expand(subs(n=1/x^r,vu2)):

vu2:=expand(ka2*vu2):


vu2:=add(simplify(coeff(vu2,x,i))*x^i,i=0..M*r):
vu2:=subs(x=1/n^(1/r),vu2):


ka1:=subs(x=1/n^(1/r),ka1):
ka1:=subs(n=n/K,ka1):
vu2a:=op(1,vu2):
vu2:=expand(normal(vu2/vu2a)):
gu:=simplify(vu2a*vu1*lu^(n/K))*exp(pu)*ka1*vu2:

GetRidOfConst(gu):

end:






#AsyC(ope,n,N,M,Ini,K): the asymptotic expansion of solutions 
#to ope(n,N)f(n)=0, with the given initial conditions
#where ope(n,N) is  a recurrence operator
#up to the M's term,
#and complete with an empirically derived constant in front
#using K terms
AsyC:=proc(ope,n,N,M,Ini,K) local gu,L,i,mu,er,C,D1:
Digits:=100:
gu:=Asy(ope,n,N,M):
if gu=FAIL then
 RETURN(FAIL):
fi:
L:=SeqFromRec(ope,n,N,evalf(Ini),K):

mu:=[seq(evalf(L[i]/subs(n=i,gu)),i=K-10..K)]:
er:=mu[nops(mu)]-mu[nops(mu)-1]:


if er=0 then
 D1:=Digits:
else
D1:=-trunc(log(abs(er))/log(10))-3:
fi:

if D1<2 then
 print(`can't determine the constant`):
RETURN(gu):
fi:

C:=evalf(mu[nops(mu)],D1):
C,gu:
end:








#FindExpP1g(ope,n,N,r,x,ds): finds the exponential part of the asymptotics
#in terms of x=n^(1/r) as a poly. of degree ds in x
#for the solutions of ope(n,N)a(n)=0 if it is of type r
#For example, try:
#FindExpP1g((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x,1);
FindExpP1g:=proc(ope,n,N,r,x,ds) local c,eq,var,s,sof,gu,ope1,
deg,i,i1,v,ka,i11:
sof:=add(c[s]*x^s,s=1..ds):
var:={seq(c[i],i=1..ds)}:

deg:=degree(ope,n):
ope1:=expand(ope/n^deg):
ope1:=expand(subs(n=1/x^r,ope1)):


gu:=0:

for i from 0 to degree(ope1,N) do
 gu:=gu+coeff(ope1,N,i)*exp(add(c[s]*Atom(s,r,i,x,2),s=1..ds) ):
od:


gu:=taylor(gu,x=0,3*r+10):


for i1 from 0 to 3*r+8 while coeff(gu,x,i1)=0 do od:

eq:={seq(coeff(gu,x,i11),i11=i1..i1+ds-1)}:
eq:={seq(coeff(gu,x,i11),i11=i1..i1+ds)}:

lprint(eq,var):
ka:=[solve(eq,var)]:


if ka=[] then
 RETURN(FAIL):
fi:


var:=ka[1]:
for v  in var do
 if op(1,v)=op(2,v) then
  RETURN(FAIL):
 fi:
od:

subs(var,sof):

end:


#FindExpP1(ope,n,N,r,x): finds the exponential part of the asymptotics
#in terms of x=n^(1/r)
#for the solutions of ope(n,N)a(n)=0 if it is of type r
#For example, try:
#FindExpP1((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x);
FindExpP1:=proc(ope,n,N,r,x) local s,gu,nosaf:

for s from 1 to r-1 do
 for nosaf from 0 to 1 do
 gu:=FindExpP1gExtra(ope,n,N,r,x,s,nosaf):
 if gu<>FAIL then
  RETURN(gu):
 fi:
od:
od:

FAIL:
end:





#FindExpP1gExtra(ope,n,N,r,x,ds,nosaf)
#: finds the exponential part of the asymptotics
#in terms of x=n^(1/r) as a poly. of degree ds in x
#for the solutions of ope(n,N)a(n)=0 if it is of type r
#For example, try:
#FindExpP1gExtra((n+1)*N^2-2*(n+5)*N+n+3,n,N,2,x,1);
FindExpP1gExtra:=proc(ope,n,N,r,x,ds,nosaf) local c,eq,var,s,sof,gu,ope1,
deg,i,i1,v,ka,i11,ku,varf:
sof:=add(c[s]*x^s,s=1..ds):
var:={seq(c[i],i=1..ds)}:

deg:=degree(ope,n):
ope1:=expand(ope/n^deg):
ope1:=expand(subs(n=1/x^r,ope1)):


gu:=0:

for i from 0 to degree(ope1,N) do
 gu:=gu+coeff(ope1,N,i)*exp(add(c[s]*Atom(s,r,i,x,2),s=1..ds) ):
od:


gu:=taylor(gu,x=0,3*r+10):


for i1 from 0 to 3*r+8 while coeff(gu,x,i1)=0 do od:

eq:={seq(coeff(gu,x,i11),i11=i1..i1+ds-1+nosaf)}:


ka:=[solve(eq,var)]:


if ka=[] then
 RETURN(FAIL):
fi:


var:=ka[1]:
for v  in var do
 if op(1,v)=op(2,v) then
  RETURN(FAIL):
 fi:
od:

#varf:=evalf(var):

#ku:=[seq(subs(varf,c[s]),s=1..ds)]:

#print(`ku is`, ku):
#add(ku[s]*x^s,s=1..ds):

subs(var,sof):

end:


 
 

#OpeBin(n,N,r): the operator for the sum of the r-th power of the
#binomial coefficients followed by their initial conditions
#For example, try:
#OpeBin(n,N,3);
OpeBin:=proc(n,N,r) local ope,k,n1:
ope:=zeil(binomial(n,k)^r/2^(n*r),k,n,N)[1]:

ope,[seq(add((n1!/k!/(n1-k)!)^r/2^r,k=0..n1),n1=1..degree(ope,N))]:
end:



#AsyBin(F,n,k,M,L): the asymptotics in n, up to order M,
#of a binomial coefficient
#sum Sum(F(n,k),k=0..n) (assuming F is supported between 0 and n
#where F is a hypergeometric term in n and k, 
#and L is the number of terms in the sequence
#for estimating the constant in front
#For example, try:
#AsyBin(binomial(n,k)*binomial(n+k,k),n,k,5,1000);
AsyBin:=proc(F,n,k,M,L) local k1,n1,ope,Ini,N:
ope:=zeil(F,k,n,N):
if ope=FAIL then
RETURN(FAIL):
fi:
ope:=ope[1]:
Ini:=[seq(add(subs({n=n1,k=k1},F),k1=0..n1),n1=1..degree(ope,N))]:

AsyC(ope,n,N,M,Ini,L):

end:


#SipurBin(R,n,M,L): The story for sums of powers of binomial
#coefficients. For example, try:
#SipurBin(4,n,5,1000);
SipurBin:=proc(R,n,M,L) local r,k,lu,t0:
print(`These are the asymptotics of the sum of the powers of Pascal's`):
print(`triangles for powers between 2 to `, R):
t0:=time():
for r from 2 to R do
print(`The asymtotics of`):
print(Sum(binomial(n,k)^r,k=0..n)):
print(`up to order`, M, `is equal to `):
lu:=AsyBin(binomial(n,k)^r,n,k,M,L):

print(lu[1]*lu[2]):

od:

print(`Everything is rigorous, but the constants in front are`):
print(`non-rigorous (yet fairly reliable) estimates`):

print(`This took`, time()-t0, `seconds of CPU time`):

end:



#AsyPerm(n,r,M,L): the asymptotics for the number of permutations
#whose r-th power is the identity permutation. 
#M is the desired order, and L is the number of terms in the
#sequence used to estimate the constant
#For example, try AsyPerm(n,2,5,1000);
#
AsyPerm:=proc(n,r,M,L) local N,ope:
ope:=OpePerG(N,n,divisors(r)):
AsyC(ope[1],n,N,M,ope[2],L):
end:


#SipurPerm(R,n,M,L): The story for the asymptotics for
#the number of premutations pi  of [1,n] such that pi^r=Identity
#for r from 2 to R (r=2 is involutions).
#M is the desired order and L is the number of terms in the
#squence for estimating the constant
#(for some reason, it only works for R<=6).
#For example, try:
#SipurPerm(6,n,5,1000);
SipurPerm:=proc(R,n,M,L) local r,lu,t0,pi, IdentityPerm:
print(`These are the asymptotics for the number of permutations`):
print(`whose r-th power is the identity permutation, for r between 2`):
print(` and `, R):
t0:=time():
for r from 2 to R do
print(`The asymtotics of the number of permutations pi of [1,n] such that`):
print(pi^r=IdentityPerm):
print(`up to order`, M, `is equal to `):
lu:=AsyPerm(n,r,M,L):

print(lu[1]*lu[2]):
print(`which in floating-point is`):
print(evalf(lu[1]*lu[2])):
od:

print(`Everything is rigorous, but the constants in front are`):
print(`non-rigorous (yet fairly reliable) estimates`):

print(`This took`, time()-t0, `seconds of CPU time`):

end:

####End AsyRec

 
##begin Findrec
ezraFindrec:=proc()
if args=NULL then

print(` FindRec: A Maple package for empirically guessing partial recurrence`):
print(`equations satisfied by Discrete Functions of TWO Variables`):
print():
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures:  `):
print(`  findrec, Findrec, FindrecF`):
print():

elif nargs=1 and args[1]=findrec then
print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`):
print(`the sequence f of degree DEGREE and order ORDER.`):
print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`):

elif nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N): Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `):
print(`of maximum DEGREE+ORDER<=MaxC`):
print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N);`):

elif nargs=1 and args[1]=FindrecF then
print(`FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with`):
print(`poly coffs. .g. try FindrecF(i->i!,n,N);`):


elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends it to the first K values`):
print(`For example, try:`):
print(`SeqFromRec(N-n-1,n,N,[1],10);`):


fi:
 
end:


###Findrec 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrecVerbose:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
print(`There is some slack, there are `, nops(ope)):
print(ope):
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:


if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then
 RETURN(FAIL):
fi:


if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 


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:


#Findrec(f,n,N): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N);
Findrec:=proc(f,n,N)
local DEGREE, ORDER,ope,L,MaxC,lu:

for MaxC from 1 do
lu:=max(seq((2+d1)*(1+MaxC-d1)+6,d1=0..MaxC)):
if lu>nops(f) then
break:
fi:
od:

MaxC:=MaxC-1:

for L from 0 to MaxC do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL then
     if SeqFromRec(ope,n,N,[op(1..ORDER,f)],nops(f))=f then
       RETURN(ope):
      fi:
     fi:
 od:
od:
FAIL:

end:





 
#SeqFromRecOld(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
SeqFromRecOld:=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 Findrec


with(linalg):

#findrecEx(f,DEGREE,ORDER,m1): Explores whether thre
#is a good chance that there is a recurrence of degree DEGREE
#and order ORDER, using the prime m1
#For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003);
findrecEx:=proc(f,DEGREE,ORDER,m1)
local ope,var,eq,i,j,n0,eq1,a,A1,
D1,E1,Eq,Var,f1,n,N:
option remember:
f1:=f mod m1:
if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1:
  od:
 
   eq:= eq union {eq1}:
od:


Eq:= convert(eq,list):
Var:= convert(var,list):


D1:=nops(Var):
E1:=nops(Eq):
if E1<D1 then
  RETURN(true):
fi:

A1:=matrix(D1,D1):

for i from 1 to D1-1 do
for j from 1 to D1 do
  A1[i,j]:=coeff(Eq[i],Var[j]):
od:
od:

for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1],Var[j]):
od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:

if E1-D1>=1 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+1],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

if E1-D1>=2 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+2],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

true:

end:

##End Findrec



with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are:  Check21a, CdsaBs3 `):
 print(` CdsaBselberg, GW,hafokh, InfoSlatex, InfoTs3, `):
 print(` MamarKatanTs3, partition1, InfoTs3, IsHkl `):
 print(` ParBox, ParH `):
else
ezra(args):
fi:

end:



ezra:=proc()

if args=NULL then
 print(` aklz, Cdsa,  dSA, InfoS, InfoT,`):
 print(` MamarKatan,MamarKatanT, SeqS, SeqT, Sklz, Tdsa, YF`):


elif args=aklz then
print(`akl1(k,l,z): The constant in front of the asympt. expansion`):
print(`for S^(z)(k,l)(n) in the Berele-Regev paper`):
print(`"Asymptotics for Young tableaux in the (k,l) hook" .`):
print(`For example, try:`);
print(`aklz(2,2,1);`):

elif args=Cdsa then
print(`Cdsa(d,s,a): the constant in front of Regev's asymptotic formula`):
print(`for the (d-sum) Tdsa. For example, try:`);
print(` Cdsa(2,2,1); `):


elif args=CdsaBselberg then
print(`CdsaBselberg(d,s,a): CdsaB via Selberg. For example, try:`):
print(`CdsaBselberg(2,3,1);`):

elif args=CdsaBs3 then
print(`CdsaBs3(d,a,K): the secdond constant in front of Regev's`):
print(` asymptotic formula `):
print(`for the (d-sum) when s=3. For example, try:`):
print(`CdsaBs3(2,1,K);`):

elif args=Check21a then
print(`Check21a(m1,s,d,b): checks proposition 2.1 for d,s, and`):
print(`m=s^2*m1^2, and a vector of integers b, fixed and small of`):
print(`length s. For example, try:`):
print(`Check21a(5,3,2,[1,0,-1]);`):

elif args=Check21b then
print(`Check21b(m1,s,d,b): checks proposition 2.1, second form`):
print(` for d,s, and`):
print(`m=s^2*m1^2, and a vector of integers b, fixed and small of`):
print(`length s. For example, try:`):
print(`Check21a(5,3,2,[1,0,-1]);`):

elif args=dSA then
print(`dSA(d,s,a,m):The Amitai Regev d-sum asymptotics`):
print(`for T_{d,ds}^(a)(dm) in Amitai Regev's paper`):
print(`"Asymptotics of Young tableaux if the strip, the d-sums"`):
print(` For example try:  `):
print(`dSA(2,3,1,m);`):

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

elif args=InfoS then
print(`InfoS(k,l,z,n0,M,K,n,N): the linear recurrence operator `):
print(`ope(n,N), phrased in terms of the discrete variable n`):
print(`and the forward shift operator N (Nx(n):=x(n+1)) `):
print(`annihilating`):
print(`the sequence S(k,l,z)(n):=The sum of (f^lam)^z over all`):
print(`partions n bounded in the Berele-Regev hook H(k,l),`):
print(`followed by the asympotics to order M (if available)`):
print(`followed by the first K terms of that sequence`):
print(`using n0 values of the sequence.`):
print(`[Note: The output is semi-rigorous. The recurrence`):
print(`is obtained by pure guessing. We know from general nonsense`):
print(`that such a recurrence exists, and it is extremely unlikely`):
print(`that the recurrence found is a different one, but`):
print(`to prove it rigorously you need the multi-variable Zeilberger`):
print(`algorithm (or the Wegscheider Sister Celine algorithm)].`):
print(`For example, try:`):
print(`InfoS(2,1,1,40,3,100,n,N);`):

elif args=InfoSlatex then
print(`InfoSlatex(k,l,z,n0,M,K,n,N): A latex version of InfoS`):
print(`InfoSlatex(2,1,1,40,3,100,n,N);`):

elif args=InfoT then
print(`InfoT(d,s,a,n0,M,K,n,N): the linear recurrence operator `):
print(`ope(n,N), phrased in terms of the discrete variable n`):
print(`and the forward shift operator N (Nx(n):=x(n+1)) `):
print(`annihilating`):
print(`the sequence T(d,s,a)(n):=The sum of (f^(d*mu))^a over all`):
print(`partions n bounded with at most s rows,`):
print(`followed by the first K terms of that sequence`):
print(`using n0 values of the sequence.`):
print(`followed by the asympotics to order M (if available)`);
print(`followed by the leading asymptotics predicted by`):
print(`Amitai Regev in his paper`):
print(`"Asymptotics for Young tableaux in the strip, the d-sums"`):
print(`[Note: The output (except for the last item)`):
print(`is semi-rigorous. The recurrence`):
print(`is obtained by pure guessing. We know from general nonsense`):
print(`that such a recurrence exists, and it is extremely unlikely`):
print(`that the recurrence found is a different one, but`):
print(`to prove it rigorously you need the multi-variable Zeilberger`):
print(`algorithm (or the Wegscheider Sister Celine algorithm)].`):
print(`For example, try:`):
print(`InfoT(2,1,1,40,3,100,n,N);`):


elif args=InfoTs3 then
print(`InfoTs3(d,a,n0,M,K,n,N): like `):
print(`InfoT(d,s,a,n0,M,K,n,N): with s=3 without using the `):
print(`Selberg evaluation `):
print(`For example, try:`):
print(`InfoTs3(2,1,40,3,100,n,N);`):

elif args=MamarKatan then
print(`MamarKatan(k,l,z,n0,M,K,n,A): Verbose version of `):
print(`InfoS(k,l,z,n0,M,K,n,N)`):
print(`A is a symbol for the sequence studied A(n) `):
print(`For example, try:`):
print(`MamarKatan(2,1,1,40,3,100,n,A);`):

elif args=MamarKatanT then
print(`MamarKatanT(d,s,a,n0,M,K,n,A): Verbose version of`):
print(`InfoT(d,s,a,n0,M,K,n,N)`):
print(`A is a symbol for the sequence studied A(n) `):
print(`For example, try:`):
print(`MamarKatanT(2,1,1,40,3,100,n,A);`):

elif args=MamarKatanTs3 then
print(`MamarKatanTs3(d,a,n0,M,K,n,A): Like`):
print(`MamarKatanT(d,s,a,n0,M,K,n,A) but only for s=3`):
print(`not using Selberg's evaluation`):
print(`For example, try:`):
print(`MamarKatanTs3(2,1,40,10,300,n,A);`):

elif args=IsHkl then
print(`IsHkl(L,k,l): is the shape L in H(k,l)? For example, try:`):
print(`IsHkl([4,3],2,2);`):


elif args=ParBox then
print(`ParBox(m,n): all the partitions bounded in an m by n box`):
print(`For example, try:`):
print(`ParBox(3,4);`):

elif args=ParH then
print(`ParH(n,k,l): all the partitions of n bounded in the hook H(k,l)`):
print(`For example, try: ParH(7,2,2);`):

elif args=SeqS then
print(`SeqS(k,l,z,n0): the first n0 terms (starting at 1)`):
print(`of the sequence for the sums of f_lam^z`):
print(`for lam a partition of n bounded in the (k,l) hook`):
print(`For example, try:`):
print(`SeqS(1,1,1,20);`):

elif args=SeqT then
print(`SeqT(d,s,a,m0): The first m0 terms of the sequence of`):
print(`sums of the a-th power of f_(d*mu)`):
print(`for mu partions of m with at most  s parts .`):
print(`For example, try:`):
print(`SeqT(1,2,1,10);`):


elif args=Sklz then
print(`Sklz(k,l,z,n): The sum of f_lam^z`):
print(`for lam a partiton of n bounded in the (k,l) hook`):
print(`For example, try:`):
print(`Sklz(1,1,1,4);`):


elif args=Tdsa then
print(`Tdsa(d,s,a,m): The sum of the a-th power of f_(d*mu)`):
print(`for mu partions of m with at most  s parts .`):
print(`For example, try:`):
print(`Tdsa(1,1,1,4);`):

elif args=YF then
print(`YF(b): The number of Young tableaux of shape b using`):
print(`the Young-Frobenius-MacMahon formula`):
print(`For example, try:`):
print(`Yf([2,2]);`):

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



#GWold(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]);
GWold:=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:
 
 

hafokh:=proc(L) local i: [seq(L[nops(L)+1-i],i=1..nops(L))]: end:

with(combinat):

partition1:=proc(n) local gu,g:
gu:=partition(n):
{seq(hafokh(g), g in gu)}:
end:

#IsHkl(L,k,l): is the shape L in H(k,l)? For example, try:
#IsHkl([4,3],2,2);
IsHkl:=proc(L,k,l) :
if nops(L)<=k or L[k+1]<=l then
 true:
else
 false:
fi:

end:



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

if min(op(b))<0 then
   RETURN(0):
fi:
 if b=[0$k] then
      1:
 elif not {seq(evalb(b[i]-b[i+1]>=0),i=1..k-1)}={true} then
     0:
 else
    add(GW([op(1..i-1,b),b[i]-1,op(i+1..k,b)]),i=1..k):
 fi:
end:
 

#SklzSlow(k,l,z,n): The sum of the z-th power of f_lam^z
#for lam bounded in the (k,l) hook
#For example, try:
#SklzSlow(1,1,1,4);
SklzSlow:=proc(k,l,z,n) local gu,g,mu:
gu:=partition1(n):
mu:=0:

for g in gu do
 if IsHkl(g,k,l)  then
  mu:=mu+YF(g)^z:
 fi:
od:
mu:
end:

SeqSslow:=proc(k,l,z,n0) local n:
[seq(SklzSlow(k,l,z,n), n=1..n0)]:
end:


#Sklz(k,l,z,n): The sum of the z-th power of f_lam^z
#for lam bounded in the (k,l) hook
#For example, try:
#Sklz(1,1,1,4);
Sklz:=proc(k,l,z,n) local gu,g:
gu:=ParH(n,k,l):
add(YF(g)^z, g in gu):
end:

SeqS:=proc(k,l,z,n0) local n:
[seq(Sklz(k,l,z,n), n=1..n0)]:
end:




#YF(b): The number of Young tableaux of shape b using
#the Young-Frobenius-MacMahon formula
#For example, try:
#Yf([2,2]);
YF:=proc(b) local i,j,k:
k:=nops(b):
convert(b,`+`)!/mul((b[i]+k-i)!,i=1..k)*
mul(mul(b[i]-b[j]+j-i,j=i+1..k),i=1..k):
end:


#ParBox(m,n): all the partitions bounded in an m by n box
#For example, try:
#ParBox(3,4);
ParBox:=proc(m,n) local gu,mu,mu1,i:
option remember:
gu:={[]}:

if n=1 then
 RETURN({[],seq([i],i=1..m)}):
fi:

if m=1 then
 RETURN({[],seq([1$i],i=1..n)}):
fi:

gu:=ParBox(m-1,n):
mu:=ParBox(m,n-1):
gu union {seq([m,op(mu1)],mu1 in mu)}:

end:

#Par1(n,m,k):partitions of n into at most k parts
#with largest part exactly m
Par1:=proc(n,m,k) local gu,mu,mu1,m1:

if n<0 then
 RETURN({}):
fi:

if k=0 then
  RETURN({}):
fi:

if k=1 then
 if n=m then
     RETURN({[m]}):
 else
  RETURN({}):
 fi:
fi:

gu:={}:
for m1 from 1 to m do
mu:=Par1(n-m,m1,k-1):
gu:=gu union {seq([m,op(mu1)],mu1 in mu)}:
od:
gu:
end:

#Par(n,k): partitions of n into at most k parts
Par:=proc(n,k) local m1,k1:
{seq(seq(op(Par1(n,m1,k1)),m1=1..n),k1=1..k)}:
end:

#Conj1(L): the conjugate partition of L
Conj1:=proc(L) local i:
hafokh([seq(i$(L[i]-L[i+1]),i=1..nops(L)-1),nops(L)$L[nops(L)]]):
end:


ParHslow:=proc(n,k,l) local gu,L,mu:
gu:=partition1(n):
mu:={}:

for L in gu do
if IsHkl(L,k,l) then
 mu:=mu union {L}:
fi:
od:
mu:
end:


#ParH(n,k,l): all the partitions of n bounded in the hook H(k,l)
#For example, try: ParH(7,2,2);
ParH:=proc(n,k,l) local gu,m1,mu,mu1,ru,ru1:
gu:=Par(n,k):

for m1 from 1 to n-1 do
mu:=Par(n-m1,k) minus Par(n-m1,k-1):

 for mu1 in mu do
  ru:=Par(m1,min(mu1[k],l)):
  gu:=gu union {seq([op(mu1),op(Conj1(ru1))],ru1 in ru)}:
 od:
od:
gu:
end:



#akl1(k,l): The constant in front of the asympt. expansion
#for S^(1)(k,l)(n) in the Berele-Regev paper
#"Asymptotics for Young tableaux in the (k,l) hook"
#For example, try:
#akl1(2,2);
akl1:=proc(k,l) local i:
(1/2)^(k*l-k-l)/Pi^((k+l)/2)*(k+l)^((k*(k-1)+l*(l-1))/4)/k!/l!*
mul(GAMMA(1+i/2),i=1..k)*mul(GAMMA(1+i/2),i=1..l):
end:



#aklz(k,l,z): The constant in front of the asympt. expansion
#for S^(z)(k,l)(n) in the Berele-Regev paper
#"Asymptotics for Young tableaux in the (k,l) hook"
#For example, try:
#aklz(2,2,1);
aklz:=proc(k,l,z) local i,lu,gu:

lu:=(1/sqrt(2*Pi))^(k+l-1)/2^(k*l)*(k+l)^((k^2+l^2)/2):

gu:=lu^z:

lu:=(z/2*(k*(k-1)+l*(l-1))+k+l)/2:
gu:=gu*(1/(k+l))^lu:

gu:=gu/k!/l!:

gu:=gu*sqrt(z/2/Pi):
gu:=gu*(sqrt(2*Pi)^(k+l)):
lu:=-1/2*(z/2*(k*(k-1)+l*(l-1))+k+l):
gu:=gu*z^lu:
gu:=gu/GAMMA(1+z/2)^(k+l):
gu:=gu*mul(GAMMA(1+z/2*i),i=1..k)*mul(GAMMA(1+z/2*i),i=1..l):
gu:

end:

#InfoS(k,l,z,n0,M,K,n,N): the linear recurrence operator 
#ope(n,N), phrased in terms of the discrete variable n
#and the forward shift operator N (Nx(n):=x(n+1)) 
#annihilating
#the sequence S(k,l,z)(n):=The sum of (f^lam)^z over all
#partions n bounded in the Berele-Regev hook H(k,l),
#followed by the first K terms of that sequence
#using n0 values of the sequence.
#followed by the asympotics to order M (if available)
#followed by the leading asymptotics predicted by
#Allan Berele and Amitai Regev in their paper
#"Asymptotics for Young tableaux in the (k,l) hook"
#[Note: The output (except for the last item)
#is semi-rigorous. The recurrence
#is obtained by pure guessing. We know from general nonsense
#that such a recurrence exists, and it is extremely unlikely
#that the recurrence found is a different one, but
#to prove it rigorously you need the multi-variable Zeilberger
#algorithm (or the Wegscheider Sister Celine algorithm)].
#For example, try:
#InfoS(2,1,1,40,3,100,n,N);
InfoS:=proc(k,l,z,n0,M,K,n,N) local gu,ope,lu1,lu2,BR:
gu:=SeqS(k,l,z,n0):

BR:=aklz(k,l,z)*(1/n)^gklz(k,l,z)*(k+l)^(z*n):



ope:=Findrec(gu,n,N):

if ope=FAIL then
 RETURN([gu,BR]):
fi:




lu1:=SeqFromRec(ope,n,N,gu,K):

if lu1=FAIL then
 RETURN([gu,BR,ope]):
fi:

lu2:=Asy(ope,n,N,M):


if lu2=FAIL then
 RETURN([gu,BR,ope,lu1]):
else
lu2:=lu2*aklz(k,l,z):
RETURN([gu,BR,ope,lu1,lu2]):
fi:

end:





#gklz(k,l,z): g(k,l,z) in the Berele-Regev paper in Theorem 1.1. of
#"Asymptotics for Young tableaux in the (k,l) hook"
gklz:=proc(k,l,z)
(z/2*(k*(k+1)+l*(l+1)-2)-(k+l-1))/2:
end:



#MamarKatan(k,l,z,n0,M,K,n,A): Verbose version of InfoS
#For example, try:
#MamarKatan(2,1,1,40,3,100,n,A);
MamarKatan:=proc(k,l,z,n0,M,K,n,A) local gu,N,Lambda,H,i,ope,ord1,BR,opee,
amal,gu5,i1,gue:


print(`On the Berele-Regev Sequence for`, H(k,l)):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let 's define the sequence`, A(n) ):
print(`to be the sum of the number of Standard Young tableaux raised`):
print(`to the`, z, `-th power , over all shapes with n cells`):
print(`that are bounded in the`, H(k,l), `hook , i.e. such that`):
print(Lambda[k+1]<=l):

gu:=InfoS(k,l,z,n0,M,K,n,N):

print(`The first`, n0, `terms of this sequence are`):
print(gu[1]):

BR:=gu[2]:
print(`The Berele-Regev asympototics is`):
print(BR):

print(`the sequence of ratios of the terms of the sequence`):
print(`to the asymptotics are`):
print([seq(evalf(gu[1][i]/subs(n=i,BR)),i=1..nops(gu[1]))]):
print(`this should tend to 1 `):

if nops(gu)=2 then
print(`We were unable to guess a linear recurrence satisfied by the sequence`):
print(`even though one definitely exists,`):
print(` you need to make the 4th argument`, n0, ` larger. `):
RETURN(` Sorry `):
fi:

ope:=gu[3]:
ord1:=degree(ope,N):
ope:=expand(subs(n=n-ord1,ope)/N^(ord1)):

ope:=Yafe(ope,N)[2]:

print(`The sequence`, A(n), `satisfies the following order-`,  ord1):
print(`linear recurrence equation`):
print(coeff(ope,N,0)*A(n)=add(-coeff(ope,N,-i)*A(n-i),i=1..-ldegree(ope,N))):



if nops(gu)=3 then
print(`We were unable to use this recurrence to get more terms`):
fi:


print(`Thanks to the recurrence`):
print(` we can now easily get many more terms of the sequence`):
print(`For example, the first `, K, `terms are `):
print(gu[4]):

print(`The ratio of the last term to the Berele-Regev asymptotics`):
print(`is: `):
print(evalf(gu[4][nops(gu[4])]/subs(n=nops(gu[4]),BR))):


if nops(gu)=4 then
print(`We were unable to find the finer asymptotics due to the`):
print(`fact that this recurrence has dominant roots that are equal but`):
print(`opposite, one could try to do the even entries separately, `):
print(`and guess a recurrence for them`):
print(`Let's do it`):
gue:=[seq(gu[1][2*i1],i1=1..nops(gu[1])/2)]:


opee:=Findrec(gue,n,N):

if opee=FAIL then
print(`There are not enough terms, make n0, currently, `, n0, `larger `):
RETURN(` Sorry `):
fi:

amal:=Asy(opee,n,N,M):

if amal=FAIL then
print(`There are not enough terms, make n0, currently, `, n0, `larger `):
RETURN(` Sorry `):
fi:

gu5:=aklz(k,l,z)/2^gklz(k,l,z)*subs(n=n/2,amal):


print(`Using the Poincare-Birkhoff-Trjitzinsky method, and using`):
print(`the Berele-Regev constant`):
print(`one can deduce the following asymptotics expression to order`, M):
print(gu5):
print(`The ratio of the last term to this more refined asymptotics`):
print(`is: `):
print(evalf(gu[4][nops(gu[4])]/subs(n=nops(gu[4]),gu5))):
RETURN(` Thank you `):
fi:

print(`Using the Poincare-Birkhoff-Trjitzinsky method, and using`):
print(`the Berele-Regev constant`):
print(`one can deduce the following asymptotics expression to order`, M):
print(gu[5]):
print(`The ratio of the last term to this more refined asymptotics`):
print(`is: `):
print(evalf(gu[4][nops(gu[4])]/subs(n=nops(gu[4]),gu[5]))):


end:








#Tdsa(d,s,a,m): The sum of the z-th power of f_(d*mu)^a
#for mu partions of m with at most  s parts 
#For example, try:
#Tdsa(1,1,1,4);
Tdsa:=proc(d,s,a,m) local gu,g,mu,g1,i:
option remember:
gu:=ParH(m,s,0):
mu:=0:
for g in gu do
g1:=[seq(g[i]$d,i=1..nops(g))]:
mu:=mu+YF(g1)^a:
od:
mu:
end:

SeqT:=proc(d,s,a,m0) local m:
option remember:
[seq(Tdsa(d,s,a,m), m=1..m0)]:
end:


#InfoT(d,s,a,n0,M,K,n,N): the linear recurrence operator 
#ope(n,N), phrased in terms of the discrete variable n
#and the forward shift operator N (Nx(n):=x(n+1)) 
#annihilating
#the sequence T(d,s,a)(n):=The sum of (f^(d*mu))^a over all
#partions n bounded with at most s rows,
#followed by the first K terms of that sequence
#using n0 values of the sequence.
#followed by the asympotics to order M (if available)
#followed by the leading asymptotics predicted by
#Amitai Regev in his paper
#"Asymptotics for Young tableaux in the strip, the d-sums"
#[Note: The output (except for the last item)
#is semi-rigorous. The recurrence
#is obtained by pure guessing. We know from general nonsense
#that such a recurrence exists, and it is extremely unlikely
#that the recurrence found is a different one, but
#to prove it rigorously you need the multi-variable Zeilberger
#algorithm (or the Wegscheider Sister Celine algorithm)].
#For example, try:
#InfoT(2,1,1,40,3,100,n,N);
InfoT:=proc(d,s,a,n0,M,K,n,N) local gu,ope,lu1,lu2,BR:
gu:=SeqT(d,s,a,n0):

BR:=Cdsa(d,s,a)*(n)^(-a*(d^2*s^2+d^2*s-2)/4+(s-1)/2)*(d*s)^(a*d*n):



ope:=Findrec(gu,n,N):

if ope=FAIL then
 RETURN([gu,BR]):
fi:




lu1:=SeqFromRec(ope,n,N,gu,K):

if lu1=FAIL then
 RETURN([gu,BR,ope]):
fi:

lu2:=Asy(ope,n,N,M):


if lu2=FAIL then
 RETURN([gu,BR,ope,lu1]):
else
lu2:=lu2*Cdsa(d,s,a):
RETURN([gu,BR,ope,lu1,lu2]):
fi:

end:




#Cdsa(d,s,a): the constant in front of Regev's asymptotic formula
#for the (d-sum) Tdsa. For example, try:
#Cdsa(2,2,1);
Cdsa:=proc(d,s,a) local gu,i,j:
gu:=(1/sqrt(2*Pi))^(d*s-1)*sqrt(d)*s^(d^2*s^2/2)*
mul(i!,i=2..d-1)^s:
gu:=gu^a:

gu:=gu*(d/s)^((s-1)*(d^2*a*s+2)/4)*d/sqrt(s)*sqrt(a/2/Pi)/s!*
(2*Pi)^(s/2)*(d^2*a)^(-s/2-d^2*a*s*(s-1)/4)*
GAMMA(1+d^2*a/2)^(-s)*
mul(GAMMA(1+d^2*a*j/2),j=1..s):
end:


#MamarKatanT(d,s,a,n0,M,K,n,A): Verbose version of InfoT
#For example, try:
#MamarKatanT(d,s,a,40,3,100,n,A);
MamarKatanT:=proc(d,s,a,n0,M,K,n,A) local gu,N,i,ope,ord1,BR,amitai:



print(`On the Regev Sequence for d-sums for d=`,d, `s= `, s, `a= `, a):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let 's define the sequence`, A(n) ):
print(`to be the sum of the number of Standard Young tableaux raised`):
print(`to the`, a, `-th power , over all magnifications by  a factor`):
print(`of `, d , ` of all shapes with at most `, s, `rows `):
gu:=InfoT(d,s,a,n0,M,K,n,N):

print(`The first`, n0, `terms of this sequence are`):
print(gu[1]):

BR:=gu[2]:
print(`The Amitai Regev asympototics is`):
print(BR):

print(`the sequence of ratios of the terms of the sequence`):
print(`to the asymptotics are`):
print([seq(evalf(gu[1][i]/subs(n=i,BR)),i=1..nops(gu[1]))]):
print(`this should tend to 1 `):

if nops(gu)=2 then
print(`We were unable to guess a linear recurrence satisfied by the sequence`):
print(`even though one definitely exists,`):
print(` you need to make the 4th argument`, n0, ` larger. `):
RETURN(` Sorry `):
fi:

ope:=gu[3]:
ord1:=degree(ope,N):
ope:=expand(subs(n=n-ord1,ope)/N^(ord1)):

ope:=Yafe(ope,N)[2]:

print(`The sequence`, A(n), `satisfies the following order-`,  ord1):
print(`linear recurrence equation`):
print(coeff(ope,N,0)*A(n)=add(-coeff(ope,N,-i)*A(n-i),i=1..-ldegree(ope,N))):



if nops(gu)=3 then
print(`We were unable to use this recurrence to get more terms`):
fi:


print(`Thanks to the recurrence`):
print(` we can now easily get many more terms of the sequence`):
print(`For example, the first  100 terms are `):
print([op(1..100,gu[4])]):

print(`The ratio of the last term to the Amitai Regev asymptotics`):
print(`is: `):
print(evalf(gu[4][nops(gu[4])]/subs(n=nops(gu[4]),BR))):


if nops(gu)=4 then
print(`We were unable to find the finer asymptotics due to the`):
print(`fact that this recurrence has dominant roots that are equal but`):
print(`opposite, one could try to do the even entries separately, `):
print(`and guess a recurrence for them`):
RETURN(`Thank You `):

fi:

print(`Using the Poincare-Birkhoff-Trjitzinsky method, and using`):
print(`the Berele-Regev constant`):
print(`one can deduce the following asymptotics expression to order`, M):
print(gu[5]):
print(`The ratio of the last term to this more refined asymptotics`):
print(`is: `):
amitai:=evalf(gu[4][nops(gu[4])]/subs(n=nops(gu[4]),gu[5])):
print(amitai):
amitai:


end:





#dSA(d,s,a,m):The Amitai Regev d-sum asymptotics
#for T_{d,ds}^(a)(dm) in Amitai Regev's paper
#"Asymptotics of Young tableaux if the strip, the d-sums"
#For example try: 
#dSA(2,3,1,m);
dSA:=proc(d,s,a,n):
Cdsa(d,s,a)*(n)^(-a*(d^2*s^2+d^2*s-2)/4+(s-1)/2)*(d*s)^(a*d*n):
end:

#Check21a(m1,s,d,b): checks proposition 2.1 for d,s, and
#m=s^2*m1^2, and a vector of integers b, fixed and small of
#length s. For example, try:
#Check21a(5,3,2,[1,1,1]);
Check21a:=proc(m1,s,d,b) local m,L,yemin,smol,Ds,i,j:
if nops(b)<>s then
print(b , `should have `, s, `parts `):
 RETURN(FAIL):
fi:
m:=s^2*m1^2:
L:=[seq(m/s+b[i]*sqrt(m),i=1..s)]:

if {seq(evalb(L[i]-L[i+1]>=0),i=1..s-1)}<>{true} then
 print(L, `should be a partition `):
 RETURN(FAIL):
fi:

L:=[seq(L[i]$d,i=1..nops(L))]:

smol:=YF(L):

Ds:=mul(mul(b[i]-b[j],j=i+1..s),i=1..s):
yemin:=(1/sqrt(2*Pi))^(d*s-1)
*d^((d^2*s^2+d^2*s)/4)
*s^(d^2*s^2/2)
*mul(i!,i=2..d-1)^s*
(1/sqrt(d*m))^((d^2*s^2+d^2*s-2)/2)*
(d*s)^(d*m)*
Ds^(d^2)*
exp(-(d*s/2)*add(b[i]^2,i=1..s)):


print(`the right side is`, evalf(yemin)):
print(`the ratio is` , evalf(yemin/smol)):
evalf(yemin/smol):

end:



#Check21b(m1,s,d,b): checks proposition 2.1 for d,s, and
#m=s^2*m1^2, and a vector of integers b, fixed and small of
#length s. For example, try:
#Check21b(5,3,2,[1,1,1]);
Check21b:=proc(m1,s,d,b) local m,L,yemin,smol,Ds,i,j:

if nops(b)<>s then
print(b , `should have `, s, `parts `):
 RETURN(FAIL):
fi:

if {seq(evalb(b[i]-b[i+1]>=0),i=1..nops(b)-1)}<>{true} then
 print(b, `should be weakly increasing `):
 RETURN(FAIL):
fi:

if convert(b,`+`)<>0 then
 print(b, `should add-up to zero `):
 RETURN(FAIL):
fi:

m:=s^2*m1^2:
L:=[seq(m/s+b[i]*sqrt(m),i=1..s)]:

if {seq(evalb(L[i]-L[i+1]>=0),i=1..s-1)}<>{true} and nops(L)>1 then
 print(L, `should be a partition `):
 RETURN(FAIL):
fi:

L:=[seq(L[i]$d,i=1..nops(L))]:

smol:=YF(L):


Ds:=mul(mul(b[i]-b[j],j=i+1..s),i=1..s):


yemin:=(1/sqrt(2*Pi))^(d*s-1)
*sqrt(d)
*s^(d^2*s^2/2)
*mul(i!,i=2..d-1)^s*
(1/sqrt(m))^((d^2*s^2+d^2*s-2)/2)*
(d*s)^(d*m)*
Ds^(d^2)*
exp(-(d*s/2)*add(b[i]^2,i=1..s)):


print(`the right side is`, evalf(yemin)):
print(`the ratio is` , evalf(yemin/smol)):
evalf(yemin/smol):

end:



#CdsaA(d,s,a): the first
#constant in front of Regev's asymptotic formula
#for the (d-sum) Tdsa. For example, try:
#CdsaA(2,2,1);
CdsaA:=proc(d,s,a) local gu,i,j:
gu:=(1/sqrt(2*Pi))^(d*s-1)*sqrt(d)*s^(d^2*s^2/2)*
mul(i!,i=2..d-1)^s:
gu:=gu^a:
end:


#CdsaBs3(d,a,K): the secdond constant in front of Regev's asymptotic formula
#for the (d-sum) when s=3. For example, try:
#CdsaBs3(2,1);
CdsaBs3:=proc(d,a,K) local gu,x1,x2,x3,Ds,s:
option remember:
s:=3:
Ds:=(x1-x2)*(x1-x3)*(x2-x3):
gu:=Ds^(d^2)*exp(-d*s/2*(x1^2+x2^2+x3^2)):
gu:=gu^a:
gu:=subs(x3=-x1-x2,gu):
evalf(int(int(gu,x2=-x1/2..x1),x1=0..K)):
end:



#CdsaBselberg(d,s,a): CdsaB via Selberg. For example, try:
#CdsaBselberg(2,3,1);
CdsaBselberg:=proc(s,d,a) local j:
(d/s)^((s-1)*(a*s+2)/4)*d/sqrt(s)*sqrt(a/2/Pi)/s!*
(2*Pi)^(s/2)*(d^2*a)^(-s/2-d^2*a*s*(s-1)/4)*
GAMMA(1+d^2*a/2)^(-s)*
mul(GAMMA(1+d^2*a*j/2),j=1..s):

end:



#InfoTs3(d,a,n0,M,K,n,N): Like InfoT but only
#for s=3, and without using the Selberg evaluation
#For example, try:
#InfoTs3(2,1,40,3,100,n,N);
InfoTs3:=proc(d,a,n0,M,K,n,N) local gu,ope,lu1,lu2,BR,s:
s:=3:
gu:=SeqT(d,s,a,n0):

BR:=CdsaA(d,s,a)*CdsaBs3(d,a,20)*
(n)^(-a*(d^2*s^2+d^2*s-2)/4+(s-1)/2)*(d*s)^(a*d*n):



ope:=Findrec(gu,n,N):

if ope=FAIL then
 RETURN([gu,BR]):
fi:




lu1:=SeqFromRec(ope,n,N,gu,K):

if lu1=FAIL then
 RETURN([gu,BR,ope]):
fi:

lu2:=Asy(ope,n,N,M):


if lu2=FAIL then
 RETURN([gu,BR,ope,lu1]):
else
lu2:=lu2*CdsaA(d,s,a)*CdsaBs3(d,a,20):
RETURN([gu,BR,ope,lu1,lu2]):
fi:

end:





#MamarKatanTs3(d,a,n0,M,K,n,A): Like
#MamarKatanT(d,s,a,n0,M,K,n,A) but only for s=3
#not using Selberg's evaluation
#For example, try:
#MamarKatanTs3(d,a,40,3,100,n,A);
MamarKatanTs3:=proc(d,a,n0,M,K,n,A) local gu,N,i,ope,ord1,BR,amitai,s:
s:=3:



print(`On the Regev Sequence for d-sums for d=`,d, `s= `, s, `a= `, a):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let 's define the sequence`, A(n) ):
print(`to be the sum of the number of Standard Young tableaux raised`):
print(`to the`, a, `-th power , over all magnifications by  a factor`):
print(`of `, d , ` of all shapes with at most `, s, `rows `):
gu:=InfoTs3(d,a,n0,M,K,n,N):

print(`The first`, n0, `terms of this sequence are`):
print(gu[1]):

BR:=gu[2]:
print(`The Amitai Regev asympototics is`):
print(BR):

print(`the sequence of ratios of the terms of the sequence`):
print(`to the asymptotics are`):
print([seq(evalf(gu[1][i]/subs(n=i,BR)),i=1..nops(gu[1]))]):
print(`this should tend to 1 `):

if nops(gu)=2 then
print(`We were unable to guess a linear recurrence satisfied by the sequence`):
print(`even though one definitely exists,`):
print(` you need to make the 4th argument`, n0, ` larger. `):
RETURN(` Sorry `):
fi:

ope:=gu[3]:
ord1:=degree(ope,N):
ope:=expand(subs(n=n-ord1,ope)/N^(ord1)):

ope:=Yafe(ope,N)[2]:

print(`The sequence`, A(n), `satisfies the following order-`,  ord1):
print(`linear recurrence equation`):
print(coeff(ope,N,0)*A(n)=add(-coeff(ope,N,-i)*A(n-i),i=1..-ldegree(ope,N))):



if nops(gu)=3 then
print(`We were unable to use this recurrence to get more terms`):
fi:


print(`Thanks to the recurrence`):
print(` we can now easily get many more terms of the sequence`):
print(`For example, the first  100 terms are `):
print([op(1..100,gu[4])]):

print(`The ratio of the last term to the Amitai Regev asymptotics`):
print(`is: `):
print(evalf(gu[4][nops(gu[4])]/subs(n=nops(gu[4]),BR))):


if nops(gu)=4 then
print(`We were unable to find the finer asymptotics due to the`):
print(`fact that this recurrence has dominant roots that are equal but`):
print(`opposite, one could try to do the even entries separately, `):
print(`and guess a recurrence for them`):
RETURN(`Thank You `):

fi:

print(`Using the Poincare-Birkhoff-Trjitzinsky method, and using`):
print(`the Berele-Regev constant`):
print(`one can deduce the following asymptotics expression to order`, M):
print(gu[5]):
print(`The ratio of the last term to this more refined asymptotics`):
print(`is: `):
amitai:=evalf(gu[4][nops(gu[4])]/subs(n=nops(gu[4]),gu[5])):
print(amitai):
amitai:


end:


#InfoSlatex(k,l,z,n0,M,K,n,N): A latex version of InfoS
#InfoSlatex(2,1,1,40,3,100,n,N);
InfoSlatex:=proc(k,l,z,n0,M,K,n,N) local gu,ope,lu1,lu2,BR:
gu:=SeqS(k,l,z,n0):

BR:=aklz(k,l,z)*(1/n)^gklz(k,l,z)*(k+l)^(z*n):



ope:=Findrec(gu,n,N):

if ope=FAIL then
print(gu):
print(latex(BR))
fi:




lu1:=SeqFromRec(ope,n,N,gu,K):

if lu1=FAIL then
print(gu):
print(latex(BR)):
print(latex(ope)):
fi:

lu2:=Asy(ope,n,N,M):


if lu2=FAIL then
print(gu):
print(latex(BR)):
print(latex(ope)):
else
lu2:=lu2*aklz(k,l,z):
print(gu):
print(latex(ope)):
print(latex(lu2)):
fi:

end: