######################################################################
##TANGLED: Save this file as TANGLED To use it, stay in the          #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read TANGLED : <Enter>                             #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Dec. 10, 2007
 
print(`Created: Dec. 10, 2007`):
print(` This is TANGLED, `):
print(`to enumerate the number of k-noncrossing tangled-diagrams`):
print(`It accompanies the paper `):
print(`Efficient Counting of k-noncrossing tangled diagrams`):
print(`by William Y.C. Chen, Jing Qin,`):
print(`Christian M. Reidys and Doron Zeilberger`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 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 different helps type ezra();, for help with`):
 print(``):

with(combinat):


ezra:=proc(): 

if nargs=0 then
print(`The MAIN procedure is: Seq, SeqD, SeqDt, Sipur `):

elif nops([args])=1 and op(1,[args])=Seq then
print(`Seq(S,F,k,L):`):
print(`the sequence whose i^th term (for i=1..L), is the number`):
print(`of simple random walks`):
print(`staying in x1>=x2>= ...>=xk>=0 `):
print(`from  the origin back to the origin using 2i steps`):
print(`For example, try:`):
print(`Seq(2,20);`):

elif nops([args])=1 and op(1,[args])=SeqD then
print(`SeqD(k,L): the sequence whose i^th term is D3kn(k,i) `):
print(`For i=1..L. For example, try:`):
print(`SeqD(2,20);`):

elif nops([args])=1 and op(1,[args])=SeqDt then
print(`SeqDt(k,L): the sequence whose i^th term is D3tkn(k,i) `):
print(`For i=1..L. For example, try:`):
print(`SeqDt(2,20);`):

elif nops([args])=1 and op(1,[args])=Sipur then
print(`Sipur(k,n,N,a,L,L1): The story about (k+1)-non-crossing `):
print(`tangled diagrams. For example, try:`):
print(`Sipur(2,n,N,a,100,1000);`):

else
 print(`there is no help for `, args ):
fi:
end:



#####Begin Walks in x1>=x2>=...xk>=0
#JatN(a,t,N): The first N terms of Sum(t^(2*n+a)/n!/(n+a)!, n=0..N);
#For example, try: JatN(1,t,10);
JatN:=proc(a,t,N) local n,gu:
gu:=0:
for n from 0 while 2*n+a<=N do
 gu:=gu+t^(2*n+a)/n!/(n+a)!:
od:

gu:

end:

#NuW3(S,F,n): The number of walks in x1>=x2>=...>=xk>=0 with
#unit steps in all directions, from pt. S to pt. F
#having n steps. For example, try:
#NuW3([0,0],[5,3],5);
NuW3:=proc(S,F,n) local k,i:
option remember:
k:=nops(S):
if nops(F)<>k then
 RETURN(FAIL):
fi:

if n=0 then 
 if F=S then
    1:
 else
   0:
 fi:

elif {seq(evalb(F[i]>=F[i+1]),i=1..k-1),evalb(F[k]>=0)}<>{true} then
    0:
else
  add(NuW3(S,[op(1..i-1,F),F[i]-1,op(i+1..k,F)],n-1), i=1..k)+
  add(NuW3(S,[op(1..i-1,F),F[i]+1,op(i+1..k,F)],n-1), i=1..k):
fi:

end:



#NuW3CT(S,F,n): The number of walks in Z^k with
#unit steps in all directions, staying in type 3 wedge
#from pt. S to pt. F
#having n steps. 
#Using the Constant-Term formula
#For example, try:
#NuW3CT([0,0],[5,3],5);
NuW3CT:=proc(S,F,n) local k,i,gu,t,lu:
k:=nops(S):
gu:=add(t[i]+1/t[i],i=1..k):
gu:=expand(gu^n):
lu:=Phi3G(t,S):
gu:=expand(gu*lu):
for i from 1 to k do
 gu:=coeff(gu,t[i],F[i]):
od:
gu:
end:


#GF3(S,F,J): The formal power series
#(as a polynomial in J[a]'s, whose coeff. of t^n/n!
#is the number of simple random walks
#staying in x[1]>=...>=x[k]>=0 starting at S ending 
#in F. It also outputs the largest index of J
#For example, try:
#GF3([0,0],[5,3],J);
GF3:=proc(S,F,J) local k,i,gu,t,lu,gadol:
k:=nops(S):
lu:=Phi3G(t,S):
lu:=expand(lu/mul(t[i]^F[i],i=1..k)):
gadol:=max(seq(max(degree(lu,t[i]),-ldegree(lu,t[i])),i=1..k)):
PolyToJ(lu,t,k,J),gadol:
end:

#CheckNuW3CT(S,N,T): checks that NuW3 and NuW3CT coincide
#for starting point S and all points in Pars(N,k) and time t<=T
#For example, try: CheckNuW3CT([3,2],7,10);
CheckNuW3CT:=proc(S,N,T) local k,gu,F,t:
k:=nops(S):
gu:=Pars(N,k) minus Pars(N,k-1):

[seq(seq(NuW3(S,gu[j],t)-NuW3CT(S,gu[j],t), j=1..nops(gu)), 
t=1..T)]:
end:




#GF3(S,F,J): The formal power series
#(as a polynomial in J[a]'s, whose coeff. of t^n/n!
#is the number of simple random walks
#staying in x[1]>=...>=x[k]>=0 starting at S ending 
#in F. It also outputs the largest index of J
#For example, try:
#GF3([0,0],[5,3],J);
GF3:=proc(S,F,J) local k,i,gu,t,lu,gadol:
k:=nops(S):
lu:=Phi3G(t,S):
lu:=expand(lu/mul(t[i]^F[i],i=1..k)):
gadol:=max(seq(max(degree(lu,t[i]),-ldegree(lu,t[i])),i=1..k)):
PolyToJ(lu,t,k,J),gadol:
end:

#Seq(S,F,k,L): Like Seqslow but using (exp.) generating function
#For example, try: Seq([0,0],[0,0],2,20);
Seqg:=proc(S,F,k,L) local gu,J,gadol,t,i: 
if nops(S)<>k or nops(F)<>k then
 RETURN(FAIL):
fi:

gu:=GF3(S,F,J):
gadol:=gu[2]:
gu:=gu[1]:

for i from 0 to gadol do
 gu:=subs(J[i]=JatN(i,t,L),gu):
od:

gu:=taylor(gu,t=0,L+1):
[seq(coeff(gu,t,i)*i!,i=1..L)]:
end:

Seq:=proc(k,L) local gu,i: 
gu:=Seqg([0$k],[0$k],k,2*L):
[seq(gu[2*i],i=1..L)]:
end:




#f3kn(k,n): The number of ways of going from the
#origin back to the origin in n steps, staying
#in x1>=x2>=...>=xk>=0 and using unit steps in
#all directions. For example, try:
#f3kn(2,10);
f3kn:=proc(k,n): NuW3([0$k],[0$k],n) :end:


#g3kn(k,n): The number of ways of going from the
#origin back to the origin in n steps, staying
#in x1>=x2>=...>=xk>=0 and using unit steps in
#all directions plus the lazy option. For example, try:
#g3kn(2,10);
g3kn:=proc(k,n) local i:  add(binomial(n,i)*f3kn(k,n-i),i=0..n) :end:

#SeqW3(k,L): the sequence whose i^th term is the number
#of non-lazy simple random walks from the origin back to the origin with
#2*i steps staying at x1>=x2>=...>=xk>=0. For example, try:
#SeqW3(2,20);
SeqW3:=proc(k,L) local i: [seq(f3kn(k,2*i),i=1..L)]:end:

#Seqslow(S,F,k,L):
#the sequence whose i^th term (for i=1..L), is the number
#of simple random walks
#staying in x1>=x2>= ...>=xk>=0 
#from  S to F`):
#For example, try:
#Seqslow([0,0],[1,2], 2,20);`):
Seqslow:=proc(S,F,k,L) local i: [seq(NuW3(S,F,i),i=1..L)]:end:


#NuLW3(S,F,n): The number of walks in x1>=x2>=...>=xk>=0 with
#unit steps in all directions, from pt. S to pt. F
#as well as the lazy (do-nothing) step
#having n steps. For example, try:
#NuLW3([0,0],[5,3],5);
NuLW3:=proc(S,F,n) local k,i:
option remember:
k:=nops(S):
if nops(F)<>k then
 RETURN(FAIL):
fi:

if n=0 then 
 if F=S then
    1:
 else
   0:
 fi:

elif {seq(evalb(F[i]>=F[i+1]),i=1..k-1),evalb(F[k]>=0)}<>{true} then
    0:
else
  add(NuLW3(S,[op(1..i-1,F),F[i]-1,op(i+1..k,F)],n-1), i=1..k)+
  add(NuLW3(S,[op(1..i-1,F),F[i]+1,op(i+1..k,F)],n-1), i=1..k)+
  NuLW3(S,F,n-1):
:
fi:

end:

#SeqLW3Old(k,L): the sequence whose i^th term is the number
#of non-lazy simple random walks from the origin back to the origin with
#i steps staying at x1>=x2>=...>=xk>=0. For example, try:
#SeqLW3Old(2,20);
SeqLW3Old:=proc(k,L) local i: [seq(NuLW3([0$k],[0$k],i),i=1..L)]:end:


#SeqLW3(k,L): the sequence whose i^th term is the number
#of non-lazy simple random walks from the origin back to the origin with
#i steps staying at x1>=x2>=...>=xk>=0. For example, try:
#SeqLW3(2,20);
SeqLW3:=proc(k,L) local i: [seq(g3kn(k,i),i=1..L)]:end:


#D3tkn(k,n): The number of (k+1)-noncrossing tangled diagrams over [n]
#For example, try: D3tkn(2,10);
D3tkn:=proc(k,n) local i: add(binomial(n,i)*f3kn(k,2*n-i),i=0..n):end:


#SeqDt(k,L): The sequence of D3tkn(k,i) for i=1..L. For example, try:
#SeqDt(2,20);
SeqDtSlow:=proc(k,L) local i: [seq(D3tkn(k,i),i=1..L)]:end:

#SeqDt(k,L): The sequence of D3tkn(k,i) for i=1..L. For example, try:
#SeqDt(2,20);
SeqDt:=proc(k,L) local i,gu: 
gu:=Seqg([0$k],[0$k],k,2*L):
[seq(add(binomial(n,i)*gu[2*n-i],i=0..n),n=1..L)]:
end:


#D3kn(k,n): The number of (k+1)-tangled diagrams over [n]
#For example, try: D3kn(2,10);
D3kn:=proc(k,n) local i: add(binomial(n,i)*D3tkn(k,n-i),i=0..n):end:

#SeqD(k,L): The sequence of D3kn(k,i) for i=1..L. For example, try:
#SeqD(2,20);
SeqDSlow:=proc(k,L) local i: [seq(D3kn(k,i),i=1..L)]:end:

#SeqD(k,L): The sequence of D3kn(k,i) for i=1..L. For example, try:
#SeqD(2,20);
SeqD:=proc(k,L) local i,gu: 
gu:=SeqDt(k,L):
[seq(1+add(binomial(n,i)*gu[i],i=1..n),n=1..L)]:
end:

#Phi3G(t,S): The anti-symmetrizerof t^S  for the group of reflections
#w.r.t. x[i]-x[i+1]=-1 and x[k]=-1 . For example, try:
#Phi3G(t,[0,0]);
Phi3G:=proc(t,S) local S1,x,i,j,k:
k:=nops(S):
S1:={seq([seq(x[j], j=1..i-1),x[i+1]-1,x[i]+1, seq(x[j],j=i+2..k)],i=1..k-1),
[seq(x[j],j=1..k-1),-2-x[k]]}:
ASx0(S1,k,t,x,S):
end:

#AS(S,k,t,x): The anti-symmetrizer of the monomial 
#t[1]^x[1]*...*t[k]^x[k] under the Coxeter group
#generated by the set of reflections in S. For example, try:
#AS({[x[2],x[1],x[3]], [x[1],x[3],x[2]]} ,3,t,x);
AS:=proc(S,k,t,x) local gu,d,i,m,mu:
gu:=0:
for d from 0 while Trans1(S,k,x,d)<>{} do
mu:=Trans1(S,k,x,d): 
gu:=gu+(-1)^d*add(mul(t[i]^m[i],i=1..k), m in mu):
od:
gu:
end:


#ASx0(S,k,t,x,x0): The anti-symmetrizer of the monomial 
#t[1]^x0[1]*...*t[k]^x0[k] under the Coxeter group
#generated by the set of reflections in S, where
#x0 is a specific numerical point in Z^k.
# For example, try:
#ASx0({[x[2],x[1],x[3]], [x[1],x[3],x[2]]} ,3,t,x,[2,1,0]);
ASx0:=proc(S,k,t,x,x0) local gu,i:
gu:=AS(S,k,t,x):
factor(subs({seq(x[i]=x0[i],i=1..k)},gu)):
end:

#Trans1(S,k,x,d): Given a set S of affine-linear transformations
#in x[1], ..., x[k] (x is a symbol and k is a positive integer)
#and a non-negative integer d, finds all transformations whose
#(minimal) length is d. For example, try:
#Trans1({[x[2],x[1]]} ,2,x,1);
Trans1:=proc(S,k,x,d) local gu,i,mu,lu,s,kha,m:
option remember:
if d=0 then
  RETURN({[seq(x[i],i=1..k)]}):
fi:

mu:=Trans1(S,k,x,d-1):

lu:={seq(op(Trans1(S,k,x,i)),i=0..d-1)}:

gu:={}:

for m in mu do
 for s in S do
  kha:=subs({seq(x[i]=s[i],i=1..k)},m):
  if not member(kha,lu) then
    gu:=gu union {kha}:
  fi:
 od:
od:

gu:
end:

#MonoToJ(M,t,k,J): Given a monomial M in t[1], ..., t[k]
#t[1]^a1*...*t[k]^ak, and a symbol J  outputs J[a1]*...J[ak]
#For example, try:
#MonoToJ(t[1]^2*t[2],t,2,J);
MonoToJ:=proc(M,t,k,J) local lu,i,d,gu:
gu:=M:
lu:=1:
for i from 1 to k do
 d:=degree(M,t[i]):
 lu:=lu*J[abs(d)]:
 gu:=normal(gu/t[i]^d):
od:

if not type(gu,integer) then
 RETURN(FAIL):
fi:

lu*gu:
end:

#PolyToJ(P,t,k,J): Given a polynomial P in t[1], ..., t[k]
#and a symbol J  outputs the umbral image t^a->J[a]
#For example, try:
#PolyToJ(t[1]+t[2],t,2,J);
PolyToJ:=proc(P,t,k,J) local i,P1:
P1:=expand(P):

if not type(P1,`+`) then
 RETURN(MonoToJ(P1,t,k,J)):
else
add(MonoToJ(op(i,P1),t,k,J),i=1..nops(P1)) :
fi:

end:






###GuessHolo without its ezra
with(SolveTools):


#####One variable rational-function geussing

#Bdok1(DS1,P,x): checks that P in the single variable x  indeed fits the data set DS1
Bdok1:=proc(DaS,P,x) local i:

for i from 1 to nops(DaS) do
if subs(x=DaS[i][1],denom(P))<>0 and
normal(subs(x=DaS[i][1],P)-DaS[i][2])<>0 then

  RETURN(false):
fi:
od:
true:
end:

#GenRat1(x,dxt,dxb,a): a generic (monic top)
#rational function of degree of top dxt 
#and degree of bottom dxb
#with coeffs.parametrized by a[i] followed by the set
#of coeffs.
GenRat1:=proc(x,dxt,dxb,a) local i:

(add(a[i]*x^i,i=0..dxt))/
add(a[i+dxt+1]*x^i,i=0..dxb)
,{seq(a[i],i=0..dxt+dxb+1)}:
end:


#GR1d1d2(S1,x,d1,d2): guesses a rational function in x of degree d1/d2
#that fits the data in DS1
GR1d1d2:=proc(S1,x,d1,d2)
local eq,var,P,a,var1,i:
P:=GenRat1(x,d1,d2,a):
var:=P[2]:
P:=P[1]:

if nops(S1)<=d1+d2+3 then
# print(`Insufficient data`):
 RETURN(FAIL):
fi:

eq:={seq(numer(normal(subs(x=S1[i][1],P)-S1[i][2])),i=1..d1+d2+3)}:

var1:=Linear(eq,var):
if expand(subs(var1,numer(P)))=0 then
  RETURN(FAIL):
fi:
P:=subs(var1,P):
if Bdok1(S1,P,x) then 
   RETURN(factor(normal(P))):
 else
 RETURN(FAIL):
fi:

end:


#GR1(S1,x): guesses a polynomial from the data-set S1
GR1:=proc(S1,x) local d1,d2,gu,L,i1:
if {seq(S1[i1][2],i1=1..nops(S1))}={0} then
   RETURN(0):
fi:

for L from 0 to nops(S1)-4 do
 for d1 from 0 to L do
   d2:=L-d1:
 gu:=GR1d1d2(S1,x,d1,d2):
  if gu<>FAIL then
    RETURN(gu):
  fi:
 od:
od:

FAIL:
end:

#GR1start(S1,x,D1,D2): guesses a polynomial from the data-set S1
#starting with degree D1/D2
GR1start:=proc(S1,x,D1,D2) local d1,d2,gu,L,i1:
if {seq(S1[i1][2],i1=1..nops(S1))}={0} then
   RETURN(0):
fi:
 for d1 from D1 to nops(S1)-3-D2 do
  for d2 from D2 to nops(S1)-3-d1 do
 gu:=GR1d1d2(S1,x,d1,d2):
  if gu<>FAIL then
    RETURN(gu):
  fi:
 od:
od:

FAIL:
end:





#GR1Pol(S1,N,x): Guesses a poly in N rational in x for
#the data set S1 of the form {[i,POL(N)]}:
GR1Pol:=proc(S1,N,x)
local gu,S1a,Deg,lu,i,i1,d1,d2,lDeg:
if {seq(S1[i1][2],i1=1..nops(S1))}={0} then
   RETURN(0):
fi:

Deg:=max(seq(degree(S1[i1][2],N),i1=1..nops(S1))):
lDeg:=min(seq(ldegree(S1[i1][2],N),i1=1..nops(S1))):
gu:=0:

S1a:= {seq([S1[i1][1],coeff(S1[i1][2],N,lDeg)],i1=1..nops(S1))}:
 lu:=GR1(S1a,x):
  if lu=FAIL  then
   RETURN(FAIL):
  else
  gu:=gu+factor(normal(lu)*N^lDeg):
  fi:


d1:=degree(numer(lu),x):
d2:=degree(denom(lu),x):
for i from lDeg+1 to Deg do
S1a:= {seq([S1[i1][1],coeff(S1[i1][2],N,i)],i1=1..nops(S1))}:
 lu:=GR1start(S1a,x,d1,d2):
  if lu=FAIL then
   RETURN(FAIL):
  else
  gu:=gu+factor(normal(lu))*N^i:
  fi:
 od:

gu:
end:


#GR1PolRat(S1,N,x,M): Guesses a poly in N rational in x for
#the data set S1 of the form {[i,POL(N,Rat(M))]}:
GR1PolRat:=proc(S1,N,x,M)
local gu,S1a,Deg,lu,i,i1:
option remember:
if {seq(S1[i1][2],i1=1..nops(S1))}={0} then
   RETURN(0):
fi:

Deg:=max(seq(degree(S1[i1][2],N),i1=1..nops(S1))):
gu:=0:
for i from 0 to Deg do
S1a:= {seq([S1[i1][1],coeff(S1[i1][2],N,i)],i1=1..nops(S1))}:
 lu:=GR1Rat(S1a,M,x):
  if lu=FAIL then
   RETURN(FAIL):
  else
  gu:=gu+factor(normal(lu))*N^i:
  fi:
 od:
gu:
end:



####End one-variable rational function guessing




###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,MaxC): 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,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:

for L from 0 to MaxC do
for ORDER from 1 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
       RETURN(ope):
      fi:
 od:
od:
FAIL:

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
SeqFromRec:=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:

#findrecK(f,DEGREE,ORDER): Is there a unique recurrence if degree DEGREE and order ORDER
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrecK([seq(i,i=1..10)],0,2);
findrecK:=proc(f,DEGREE,ORDER)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,n,N,a:
option remember:
if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:

#if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then
# RETURN(FAIL):
#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(false):
fi:

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

if nops(ope)=1 then
 RETURN(true):
else
 RETURN(false):
fi:

end:
 
 
#end Findrec



####Begin Asymptotics
#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:
for S1 from 1 while OneStepAS1(ope1,n,N,alpha,f,S1)=FAIL do od:
OneStepAS1(ope1,n,N,alpha,f,S1):
end:

#Asy(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
Asy:=proc(ope,n,N,K)
local gu,pit,lu,makom,alpha,mu,ope1,ku,i,f,x,ka:

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]:

     RETURN(pit[makom[1]]^n*expand((nops(makom)-1)!*binomial(nops(makom)-1+n,nops(makom)-1))):
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(`Not only that but Dominant roots are complex`):
 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]:

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,nops(makom)+1),x,nops(makom))):

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

alpha:=max(op(ka)):

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

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

lu^n*n^alpha*f:
end:

###EndAsymptotics



#####Coeffs. of reciprocal of a polynomial


#PolToSteps(POL,xList): Inputs a polynomial POL in
#the set of variables xList and returns the 
#Weighted Steps for the scheme followed by the reciprocal
#of the constant term
PolToSteps:=proc(POL,xList) local k,i,POL1,mu,Steps,ku,ku1,i1,pu:
option remember:
POL1:=expand(POL):
k:=nops(xList):
if not type(POL,`+`) then
 print(`POL must have at least two terms`):
 RETURN(FAIL):
fi:

pu:=POL1:
for i from 1 to nops(xList) do
pu:=coeff(pu,xList[i],0):
od:

if pu=0 then
 print(`The polynomial should have a non-zero constant term`):
 RETURN(FAIL):
fi:

POL1:=expand(POL1/pu):

Steps:={}:

for i from 1 to nops(POL1) do
 mu:=op(i,POL1):

 ku:=[seq(degree(mu,xList[i1]),i1=1..nops(xList))]:
 ku1:=normal(mu/convert([seq(xList[i1]^ku[i1],i1=1..nops(xList))],`*`)):
 Steps:=Steps union {[ku,-ku1]}:
od:
Steps:=Steps minus {[[0$k],-1]}:
Steps,1/pu:
end:


#Coeff1(Steps,a): 
#The number of weighter-lattice walks in the k-dim lattice with pos. steps
#in the set of steps Steps ending at the point a
#For example, try: Coeff1({[[1,0],1],[[0,1],1]},[1,1]);
Coeff1:=proc(Steps,a) local i,j,k,Prev:
option remember:
k:=nops(a):
if a=[0$k] then
  RETURN(1):
fi:
if min(op(a))<0 then
  RETURN(0):
fi:

Prev:=
{seq([[seq(a[j]-Steps[i][1][j],j=1..k)],Steps[i][2]],i=1..nops(Steps))}:
expand(add(Prev[i][2]*Coeff1(Steps,Prev[i][1]),i=1..nops(Prev))):

end:



#Coeff1S(Steps,a): 
#The number of weighted-lattice walks in the k-dim lattice with pos. steps
#in the set of symmetric steps Steps ending at the point a
#For example, try: Coeff1({[[1,0],1],[[0,1],1]},[1,1]);
Coeff1S:=proc(Steps,a) local i,j,k,Prev:
option remember:
k:=nops(a):
if a=[0$k] then
  RETURN(1):
fi:
if min(op(a))<0 then
  RETURN(0):
fi:

Prev:=
{seq([[seq(a[j]-Steps[i][
1][j],j=1..k)],Steps[i][2]],i=1..nops(Steps))}:
expand(add(Prev[i][2]*Coeff1S(Steps,sort(Prev[i][1])),i=1..nops(Prev))):

end:


#CoeffP(POL,xList,aList): The coeff. of xList^aList in POL
#(a polynomial in the variables of xList)
CoeffP:=proc(POL,xList,aList)
local Steps,coe:
Steps:=PolToSteps(POL,xList):
coe:=Steps[2]:
Steps:=Steps[1]:
coe*Coeff1(Steps,aList): 
end:

#CoeffPs(POL,xList,aList): The coeff. of xList^aList in 1/POL
#(a Symmetric polynomial in the variables of xList)
CoeffPs:=proc(POL,xList,aList)
local Steps,coe:
Steps:=PolToSteps(POL,xList):
coe:=Steps[2]:
Steps:=Steps[1]:
coe*Coeff1S(Steps,aList): 
end:

#####

#LaW(Steps,a): The number of lattice walks in the k-dim lattice with pos. steps
#in the set of steps Steps ending at the point a
#For example, try: LaW({[1,0],[0,1]},[1,1]);
LaW:=proc(Steps,a) local i,j,k,Prev:
option remember:
k:=nops(a):
if a=[0$k] then
  RETURN(1):
fi:
if min(op(a))<0 then
  RETURN(0):
fi:

Prev:={seq([seq(a[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}:
add(LaW(Steps,Prev[i]),i=1..nops(Prev)):

end:

#LaWp(Steps,a): The probability of lattice walks in the k-dim lattice with pos. steps
#in the set of steps Steps with attached probabilites
#to get to the point a
#For example, try: LaWp({[[1,0],1/2],[[0,1],1/2]},[1,1]);
LaWp:=proc(Steps,a) local i,j,k,Prev:
option remember:
k:=nops(a):
if a=[0$k] then
  RETURN(1):
fi:
if min(op(a))<0 then
  RETURN(0):
fi:

Prev:={seq([[seq(a[j]-Steps[i][1][j],j=1..k)],Steps[i][2]],i=1..nops(Steps))}:
add(Prev[i][2]*LaWp(Steps,Prev[i][1]),i=1..nops(Prev)):

end:

#BaW(Steps,a): The number of ballot-lattice walks in the k-dim lattice with pos. steps
#in the set of steps Steps ending at the point a
#For example, try: BaW({[1,0],[0,1]},[1,1]);
BaW:=proc(Steps,a) local i,j,k,Prev:
option remember:
k:=nops(a):
if a=[0$k] then
  RETURN(1):
fi:
if min(op(a))<0 then
  RETURN(0):
fi:

if min(seq(a[i]-a[i+1],i=1..nops(a)-1))<0 then
 RETURN(0):
fi:

Prev:={seq([seq(a[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}:
add(BaW(Steps,Prev[i]),i=1..nops(Prev)):

end:


#BaWp(Steps,a): The probability of ballot-lattice walks in the k-dim lattice with pos. steps
#in the set of steps Steps, with assoc. prob. ending at the point a
#For example, try: BaWp({[[1,0],1/2],[[0,1],1/2]},[1,1]);
BaWp:=proc(Steps,a) local i,j,k,Prev:
option remember:
k:=nops(a):
if a=[0$k] then
  RETURN(1):
fi:
if min(op(a))<0 then
  RETURN(0):
fi:

if min(seq(a[i]-a[i+1],i=1..nops(a)-1))<0 then
 RETURN(0):
fi:

Prev:={seq([[seq(a[j]-Steps[i][1][j],j=1..k)],Steps[i][2]],i=1..nops(Steps))}:
add(Prev[i][2]*BaWp(Steps,Prev[i][1]),i=1..nops(Prev)):

end:


#BaWg(Steps,a,g): The number of generalized ballot-lattice walks 
#(i.e. that stay in a[i]-a[i+1]>=-g
#in the k-dim lattice with pos. steps
#in the set of steps Steps ending at the point a
#For example, try: BaWg({[1,0],[0,1]},[1,1],2);
BaWg:=proc(Steps,a,g) local i,j,k,Prev:
option remember:
k:=nops(a):
if a=[0$k] then
  RETURN(1):
fi:
if min(op(a))<0 then
  RETURN(0):
fi:

if min(seq(a[i]-a[i+1]+g,i=1..nops(a)-1))<0 then
 RETURN(0):
fi:

Prev:={seq([seq(a[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}:
add(BaWg(Steps,Prev[i],g),i=1..nops(Prev)):

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
SeqFromRec:=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:





#FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with
#poly coffs.
#e.g. try FindrecF(i->i!,n,N);
FindrecF:=proc(f,n,N)
local DEGREE, ORDER,ope,L,i,t0:
option remeber:
t0:=time():

for L from 0 do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:

 ope:=findrec([seq(f(i),i=1..(2+DEGREE)*(1+ORDER)+4)],DEGREE,ORDER,n,N):
     if ope<>FAIL then
       RETURN(ope):
      fi:
 od:
od:
FAIL:

end:


#GH2MPol(POL,x,y,m,n,M,Patience,s1): The recurrence in the m-direction  of
#the coeff. of x^m*y^n of 1/POL(x,y)
#For example, try: GH2MPol(1-x-y,x,y,m,n,M,10);
GH2MPol:=proc(P,x,y,m,n,M,Patience,s1) local i,j:
if expand(subs({x=y,y=x},P)-P)=0 then
GH2M((i,j)->CoeffPs(P,[x,y],[i,j]),m,n,M,Patience,s1):
else
GH2M((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,Patience,s1):
fi:
end:


GH2:=proc(F,m,n,M,N,Patience,s1) local i,j: GH2M(F,m,n,M,Patience,s1),GH2M((i,j)->F(j,i),n,m,N,Patience,s1):end:

GH2dir:=proc(F,m,n,M,N,MaxC) local i,j: GH2Mdir(F,m,n,M),GH2Mdir((i,j)->F(j,i),n,m,N):end:


GH2PolDir:=proc(P,x,y,m,n,M,N) local ope:
if expand(subs({x=y,y=x},P)-P)=0 then
ope:=GH2Mdir((i,j)->CoeffPs(P,[x,y],[i,j]),m,n,M):
RETURN(ope,subs({m=n,n=m,M=N},ope)):
else
GH2dir((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,N):
fi:
end:

#GH2c(F,m,n,M,N,Patience,s1): Canonical Form of GH2
GH2c:=proc(F,m,n,M,N,Patience,s1) local lu,lu1,lu2: 
lu:=GH2(F,m,n,M,N,Patience,s1):
lu1:=numer(normal(lu[1])):lu2:=numer(normal(lu[2])):
Yafe1(lu1,M,N),Yafe1(lu2,M,N):end:


GH2Polc:=proc(P,x,y,m,n,M,N)local lu:
lu:=GH2c((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,N):
Yafe1(lu[1],M,N),Yafe1(lu[2],M,N):end:




GH2DPol:=proc(P,x,y,m,n,M,N):GH2D((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,N):end:


#EvalOpeF2(F,ope,m,n,M,N,m0,n0): evaluates ope(m,n,M,N)F(m0,n0)
#For example try: 
#EvalOpeF2((i,j)->binomial(i+j,j),M*N-M-N,m,n,M,N,3,5);
EvalOpeF2:=proc(F,ope,m,n,M,N,m0,n0)
local gu,i,j:
gu:=0:
for i from 0 to degree(ope,M) do
for j from 0 to degree(ope,N) do
 gu:=gu+subs({m=m0,n=n0},coeff(coeff(ope,M,i),N,j))*F(m0+i,n0+j):
od:
od:
gu:
end:

#CanSmn(OpeM,OpeN,m,n,M,N,i0,j0,F): Given a holonomic representaion
#of a discrete function (let's call it F(m,n)), in terms
#of the two operators OpeM(m,n,M) and OpeN(m,n,N), expresses
#F(m+i0,n+j0) in terms of F(m+a,n+b), with 0<=a<degree(opeM,M)
#and 0<=b<=degree(opeN,N). It first goes in the m direction
#and then in the n direction. For example, try:
#CanSmn(M-1,N-1,m,n,M,N,1,1,F);
CanSmn:=proc(OpeM,OpeN,m,n,M,N,i0,j0,F)
local gu,Ordm,Ordn,i1,gu1,j1:
option remember:
Ordm:=degree(OpeM,M):
Ordn:=degree(OpeN,N):


if coeff(OpeM,M,Ordm)<>1 then
 print(OpeM, `should be monic in `, M):
 RETURN(FAIL):
fi:

if coeff(OpeN,N,Ordn)<>1 then
 print(OpeN, `should be monic in `, N):
 RETURN(FAIL):
fi:


if i0<Ordm and j0<Ordn then
 RETURN(F(m+i0,n+j0)):
fi:

if j0<Ordn then
 if i0=Ordm then
  RETURN(add(subs(n=n+j0,normal(-coeff(OpeM,M,i1)))*F(m+i1,n+j0),i1=0..Ordm-1)):
 fi:

gu:=subs(m=m+1,CanSmn(OpeM,OpeN,m,n,M,N,i0-1,j0,F)):
gu:=expand(subs(F(m+Ordm,n+j0)=CanSmn(OpeM,OpeN,m,n,M,N,Ordm,j0,F),gu)):
gu:=add(factor(normal(coeff(gu,F(m+i1,n+j0),1)))*F(m+i1,n+j0),i1=0..Ordm-1):
RETURN(gu):
fi:

#Now j0>=Ordn
if i0<Ordm then


 if j0=Ordn then
  RETURN(add(subs(m=m+i0,normal(-coeff(OpeN,N,i1)))*F(m+i0,n+i1),i1=0..Ordn-1)):
 fi:

gu:=subs(n=n+1,CanSmn(OpeM,OpeN,m,n,M,N,i0,j0-1,F)):
gu:=expand(subs(F(m+i0,n+Ordn)=CanSmn(OpeM,OpeN,m,n,M,N,i0,Ordn,F),gu)):
gu:=add(factor(normal(coeff(gu,F(m+i0,n+i1),1)))*F(m+i0,n+i1),i1=0..Ordn-1):
RETURN(gu):
fi:

#Now i0>=Ordm and j0>=Ordn

gu:=subs(m=m+1,CanSmn(OpeM,OpeN,m,n,M,N,i0-1,j0,F)):

for j1 from 0 to Ordn-1 do
gu:=expand(subs(F(m+Ordm,n+j1)=CanSmn(OpeM,OpeN,m,n,M,N,Ordm,j1,F),gu)):
od:
gu1:=0:

for i1 from 0 to Ordm-1 do
for j1 from 0 to Ordn-1 do
gu1:=gu1+factor(coeff(gu,F(m+i1,n+j1),1))*F(m+i1,n+j1):
od:
od:
RETURN(gu1):

end:





#CanSmnO(OpeM,OpeN,m,n,M,N,i0,j0):
#CanSmnO(M-1,N-1,m,n,M,N,1,1);
CanSmnO:=proc(OpeM,OpeN,m,n,M,N,i0,j0)
local F,i1,j1,gu,gu1,Ordm,Ordn:

gu:=CanSmn(OpeM,OpeN,m,n,M,N,i0,j0,F):

Ordm:=degree(OpeM,M):
Ordn:=degree(OpeN,N):

gu1:=0:

for i1 from 0 to Ordm-1 do
for j1 from 0 to Ordn-1 do
gu1:=gu1+factor(coeff(gu,F(m+i1,n+j1),1))*M^i1*N^j1:
od:
od:
gu1:=M^i0*N^j0-gu1:
RETURN(gu1):
end:



#CheckCanSmnO(P,x,y,i0,j0): Checks CanSmnO on the reciprocals of P
#For example:
#CheckCanSmnO(1-x-y-x*y,x,y,3,4);
CheckCanSmnO:=proc(P,x,y,i0,j0)
local lu,m,n,M,N,Ope,i,j,m0,n0:

lu:=GH2Pol(P,x,y,m,n,M,N,5):
Ope:=CanSmnO(lu[1],lu[2],m,n,M,N,i0,j0):

evalb(
{seq(seq(EvalOpeF2((i,j)->CoeffP(P,[x,y],[i,j]),Ope,m,n,M,N,m0,n0),m0=0..40),n0=0..40)}
={0})
:
end:




#CanSmnO1(OpeM,OpeN,m,n,M,N,i0,j0):
#CanSmnO1(M-1,N-1,m,n,M,N,1,1);
CanSmnO1:=proc(OpeM,OpeN,m,n,M,N,i0,j0)
local F,i1,j1,gu,gu1,Ordm,Ordn:

gu:=CanSmn(OpeM,OpeN,m,n,M,N,i0,j0,F):

Ordm:=degree(OpeM,M):
Ordn:=degree(OpeN,N):

gu1:=0:

for i1 from 0 to Ordm-1 do
for j1 from 0 to Ordn-1 do
gu1:=gu1+factor(coeff(gu,F(m+i1,n+j1),1))*M^i1*N^j1:
od:
od:
gu1:=gu1:
RETURN(gu1):
end:






#CanSnm(OpeM,OpeN,m,n,M,N,i0,j0,F): Given a holonomic representaion
#of a discrete function (let's call it F(m,n)), in terms
#of the two operators OpeM(m,n,M) and OpeN(m,n,N), expresses
#F(m+i0,n+j0) in terms of F(m+a,n+b), with 0<=a<degree(opeM,M)
#and 0<=b<=degree(opeN,N). It first goes in the n direction
#and then in the m direction. For example, try:
#CanSnn(M-1,N-1,m,n,M,N,1,1,F);
CanSnm:=proc(OpeM,OpeN,m,n,M,N,i0,j0,F)
local gu,Ordm,Ordn,i1,gu1,j1:
option remember:
Ordm:=degree(OpeM,M):
Ordn:=degree(OpeN,N):


if coeff(OpeM,M,Ordm)<>1 then
 print(OpeM, `should be monic in `, M):
 RETURN(FAIL):
fi:

if coeff(OpeN,N,Ordn)<>1 then
 print(OpeN, `should be monic in `, N):
 RETURN(FAIL):
fi:


if i0<Ordm and j0<Ordn then
 RETURN(F(m+i0,n+j0)):
fi:

if j0<Ordn then
 if i0=Ordm then
  RETURN(add(subs(n=n+j0,normal(-coeff(OpeM,M,i1)))*F(m+i1,n+j0),i1=0..Ordm-1)):
 fi:

gu:=subs(m=m+1,CanSnm(OpeM,OpeN,m,n,M,N,i0-1,j0,F)):
gu:=expand(subs(F(m+Ordm,n+j0)=CanSnm(OpeM,OpeN,m,n,M,N,Ordm,j0,F),gu)):
gu:=add(factor(normal(coeff(gu,F(m+i1,n+j0),1)))*F(m+i1,n+j0),i1=0..Ordm-1):
RETURN(gu):
fi:

#Now j0>=Ordn
if i0<Ordm then


 if j0=Ordn then
  RETURN(add(subs(m=m+i0,normal(-coeff(OpeN,N,i1)))*F(m+i0,n+i1),i1=0..Ordn-1)):
 fi:

gu:=subs(n=n+1,CanSnm(OpeM,OpeN,m,n,M,N,i0,j0-1,F)):
gu:=expand(subs(F(m+i0,n+Ordn)=CanSnm(OpeM,OpeN,m,n,M,N,i0,Ordn,F),gu)):
gu:=add(factor(normal(coeff(gu,F(m+i0,n+i1),1)))*F(m+i0,n+i1),i1=0..Ordn-1):
RETURN(gu):
fi:

#Now i0>=Ordm and j0>=Ordn

gu:=subs(n=n+1,CanSmn(OpeM,OpeN,m,n,M,N,i0,j0-1,F)):


for i1 from 0 to Ordm-1 do
gu:=expand(subs(F(m+i1,n+Ordn)=CanSnm(OpeM,OpeN,m,n,M,N,i1,Ordn,F),gu)):
od:


gu1:=0:

for i1 from 0 to Ordm-1 do
for j1 from 0 to Ordn-1 do
gu1:=gu1+factor(coeff(gu,F(m+i1,n+j1),1))*F(m+i1,n+j1):
od:
od:
RETURN(gu1):

end:


#CanSnmO(OpeM,OpeN,m,n,M,N,i0,j0):
#CanSnmO(M-1,N-1,m,n,M,N,1,1);
CanSnmO:=proc(OpeM,OpeN,m,n,M,N,i0,j0)
local F,i1,j1,gu,gu1,Ordm,Ordn:

gu:=CanSnm(OpeM,OpeN,m,n,M,N,i0,j0,F):

Ordm:=degree(OpeM,M):
Ordn:=degree(OpeN,N):

gu1:=0:

for i1 from 0 to Ordm-1 do
for j1 from 0 to Ordn-1 do
gu1:=gu1+factor(coeff(gu,F(m+i1,n+j1),1))*M^i1*N^j1:
od:
od:
gu1:=M^i0*N^j0-gu1:
RETURN(gu1):
end:


#CheckCanSnmO(P,x,y,i0,j0): Checks CanSmnO on the reciprocals of P
#For example:
#CheckCanSnmO(1-x-y-x*y,x,y,3,4);
CheckCanSnmO:=proc(P,x,y,i0,j0)
local lu,m,n,M,N,Ope,i,j,m0,n0:

lu:=GH2Pol(P,x,y,m,n,M,N,5):
Ope:=CanSnmO(lu[1],lu[2],m,n,M,N,i0,j0):

evalb(
{seq(seq(EvalOpeF2((i,j)->CoeffP(P,[x,y],[i,j]),Ope,m,n,M,N,m0,n0),m0=0..40),n0=0..40)}
={0})
:
end:

#Yafe1(Ope,M,N): Makes Ope nice
Yafe1:=proc(Ope,M,N)
local d1,d2,Ope1,i1,i2,gu,Ld1,Ld2:
#Ope1:=numer(normal(Ope)):
Ope1:=Ope:
d1:=degree(Ope1,M):
Ld1:=ldegree(Ope1,M):
d2:=degree(Ope1,N):
Ld2:=ldegree(Ope1,N):
gu:=0:
for i1 from Ld1 to d1 do
for i2 from Ld2 to d2 do
 gu:=gu+factor(coeff(coeff(Ope1,M,i1),N,i2))*M^i1*N^i2:
od:
od:

gu:

end:


#Mult(Ope1,Ope2,m,n,M,N): Ope1(m,n,M,N)*Ope2(m,n,M,N)
Mult:=proc(Ope1,Ope2,m,n,M,N) local d1,d2,i1,i2,gu,Ld1,Ld2:
d1:=degree(Ope1,M):
Ld1:=ldegree(Ope1,M):
d2:=degree(Ope1,N):
Ld2:=ldegree(Ope1,N):
gu:=0:
for i1 from Ld1 to d1 do
for i2 from Ld2 to d2 do
 gu:=expand(gu+factor(coeff(coeff(Ope1,M,i1),N,i2))*M^i1*N^i2
  *subs({m=m+i1,n=n+i2},Ope2) ):
od:
od:
Yafe1(gu,M,N):
end:

#Comm(Ope1,Ope2,m,n,M,N): The commutator of Ope1, Ope2
#i.e. Ope1*Ope2-Ope2*Ope1
Comm:=proc(Ope1,Ope2,m,n,M,N)
Yafe1(Mult(Ope1,Ope2,m,n,M,N)-Mult(Ope2,Ope1,m,n,M,N),M,N):
end:


#CommSeq(OPE0,Ope1,m,n,M,N): Given a partial recurrence operator
#with constant coefficients, OPE0, and a partial recurrence
#operator with polynomial coefficients, finds the
#sequence of commuators with OPE0: [Ope1,[OPE0,Ope1],[OPE0,%], ...]
#and checks that the last one is a multiple of OPE0. Otherwise
#it returns FAIL
CommSeq:=proc(OPE0,Ope1,m,n,M,N)
local gu,lu,gu1:
lu:=Ope1:
gu:=[lu]:
while lu<>0 do
lu:=Comm(OPE0,lu,m,n,M,N):
gu:=[op(gu),lu]:
od:
gu:=[op(1..nops(gu)-1,gu)]:
gu:=[op(1..nops(gu)-1,gu),factor(gu[nops(gu)])]:
gu1:=denom(ormal(gu[nops(gu)]/OPE0));
if gu1<>1 then
 FAIL:
else
 gu:
fi:
end:



Diag1:=proc(OpeM,OpeN,m,n,M,N,k0)
local gu,a,var,eq,eqTry,i1,j1,OpeD,var1T,i,var1:

OpeD:=M^k0*N^k0:
var:={}:

 gu:=CanSmnO1(OpeM,OpeN,m,n,M,N,k0,k0):
for i from 0 to k0-1 do
 gu:=gu+a[i]*CanSmnO1(OpeM,OpeN,m,n,M,N,i,i):
 var:= var union {a[i]}:
 OpeD:=OpeD+a[i]*M^i*N^i:
od:
gu:=expand(gu):
 eq:={seq(seq(coeff(coeff(gu,M,i1),N,j1),j1=0..degree(gu,N)),
i1=0..degree(gu,M))}:

eqTry:=subs({m=11/5,n=17/19},eq):

var1T:=Linear(eqTry,var):

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

var1:=Linear(eq,var):

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

end:

#DiagOper(OpeM,OpeN,m,n,M,N):The operator in the MN direction
#deduced from OpeM(M,m,n) and OpeN(N,m,n)
#For example, try: DiagOper(M-1,N-1,m,n,M,N);
DiagOper:=proc(OpeM,OpeN,m,n,M,N) local k0,gu:
for k0 from 1 do
gu:=Diag1(OpeM,OpeN,m,n,M,N,k0):
 if gu<>FAIL then
   RETURN(gu):
fi:
od:
end:




#GH2MPolPar(POL,x,y,m,n,M,par,Patience,s1): The recurrence in the m-direction  of
#the coeff. of x^m*y^n of 1/POL(x,y)
#where POL also depends (in a rational way)
#on P
#For example, try: GH2MPolPar(1-x-y-a*x*y,x,y,m,n,M,a,10,11);
GH2MPolPar:=proc(P,x,y,m,n,M,par,Patience,s1) local S1,i,Ope,Ope1,i1,j1,m0,n0,N:
S1:={}:
for i from 2 to 5 do
 S1:=S1 union {[i,GH2MPol(subs(par=i,P),x,y,m,n,M,Patience,s1)]}:
od:

for i from 6 do
Ope:= GR1Pol(S1,M,par):
if Ope<>FAIL then

  Ope1:=numer(normal(Ope)):

 if {seq(seq(normal(EvalOpeF2((i1,j1)->CoeffP(P,[x,y],[i1,j1]),
  Ope1,m,n,M,N,m0,n0)),m0=0..10),n0=0..10)}<>{0} then
print( {seq(seq(normal(EvalOpeF2((i1,j1)->CoeffP(P,[x,y],[i1,j1]),
  Ope1,m,n,M,N,m0,n0)),m0=0..10),n0=0..10)}):

  print(Ope, ` failed `):
 else
  RETURN(Ope):
 fi:
fi:
 S1:=S1 union {[i,GH2MPol(subs(par=i,P),x,y,m,n,M,Patience,s1)]}:
od:

end:

#GH2PolPar(P,x,y,m,n,M,N,par,Patience,s1): Like GH2Pol, but P depends
#on a parameter par
GH2PolPar:=proc(P,x,y,m,n,M,N,par,Patience,s1) local ope:
if expand(subs({x=y,y=x},P)-P)=0 then
ope:=GH2MPolPar(P,x,y,m,n,M,par,Patience,s1):
RETURN(ope,subs({m=n,n=m,M=N},ope)):
else
RETURN(GH2MPolPar(P,x,y,m,n,M,par,Patience,s1),
subs({m=n,n=m,M=N},GH2MPolPar(subs({x=y,y=x},P),x,y,m,n,M,par,Patience,s1)))
:
fi:
end:



#GuessDiag(F,n,N,a): Given a function F of two discrete variables, guesses
#a recurrence for  the diag. F(n+a,n), where N is the shift in N
#For example, try: GuessDiag((i,j)->i!+j!,n,N,0);
GuessDiag:=proc(F,n,N,a) local i:
FindrecF(i->F(i+a,i),n,N):
end:

#GuessDiagPol(P,x,y,n,N,a): Given a polynomial of two  variables x and y, 
#guesses a recurrence for  the coeffs. of x^n*y^(n+a)
#, where N is the shift in N
#For example, try: GuessDiagPol(1-x-y+2*x*y,x,y,n,N,0);
GuessDiagPol:=proc(P,x,y,n,N,a) local i,j:
GuessDiag((i,j)->CoeffP(P,[x,y],[i,j]),n,N,a):
end:


#GuessDiagPar(F,n,N,a,p): Given a function of two variables (m,n) 
#(featuring a paramter p), guesss a recurrence
#in the diagonal (n+a,n)
#For example, try: GuessDiagPar((i,j)->p*i!+j!,n,N,1,p);
GuessDiagPar:=proc(F,n,N,a,p)
local i,ope,h,ope1,ope2,ope3,DEG,ORD,gu1,gu2,gu3,S1,
Ope,j:

ope[1]:=FindrecF(h->subs(p=1,F(h+a,h)),n,N):
ope[2]:=FindrecF(h->subs(p=2,F(h+a,h)),n,N):
ope[3]:=FindrecF(h->subs(p=3,F(h+a,h)),n,N):
ope1:=numer(normal(ope[1])):
ope2:=numer(normal(ope[2])):
ope3:=numer(normal(ope[3])):

for i from 4 do

ope1:=numer(normal(ope[i-1])):
ope2:=numer(normal(ope[i-2])):
ope3:=numer(normal(ope[i-3])):
  if nops({degree(ope1,n),degree(ope2,n),degree(ope2,n)})=1 and
   nops({degree(ope1,N),degree(ope2,N),degree(ope2,N)})=1  then
   DEG:=degree(ope3,n):
   ORD:=degree(ope3,N):

gu1:=[seq(subs(p=i-1,F(h+a,h)),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
gu2:=[seq(subs(p=i-2,F(h+a,h)),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
gu3:=[seq(subs(p=i-3,F(h+a,h)),h=1..(1+DEG)*(1+ORD)+5+ORD)]:

if findrecK(gu1,DEG,ORD) and findrecK(gu2,DEG,ORD) and
    findrecK(gu3,DEG,ORD) then
      break:

  fi:
fi:

ope[i]:=FindrecF(h->subs(p=i,F(h+a,h)),n,N):

od:

S1:={}:

for j from i-2 to i+6 do

  gu1:=[seq(subs(p=j,F(h+a,h)),h=1..(1+DEG)*(1+ORD)+5+ORD)]:

   if findrecK(gu1,DEG,ORD) then
    ope1:=findrec(gu1,DEG,ORD,n,N):
    S1:=S1 union {[j,ope1]}:

   fi:

od:

for j from i+7 do
Ope:= GR1Pol(S1,N,p):
 if Ope<>FAIL then
  RETURN(Ope):
 fi:

  gu1:=[seq(subs(p=j,F(h+a,h)),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
   if findrecK(gu1,DEG,ORD) then
    ope1:=findrec(gu1,DEG,ORD,n,N):
    S1:=S1 union {[j,ope1]}:
   fi:

od:
end:




#GuessDiagPolPar(P,x,y,n,N,a,p): Given a polynomial of two  variables x and y, 
#also depending (rationally) on a parameter p
#guesses a recurrence for  the coeffs. of x^n*y^(n+a)
#, where N is the shift in N
#For example, try: GuessDiagPolPar(1-x-y+p*x*y,x,y,n,N,0,p);
GuessDiagPolPar:=proc(P,x,y,n,N,a,p):
GuessDiagPar((i,j)->CoeffP(P,[x,y],[i,j]),n,N,a,p):
end:




#GR1Rat(S1,N,x): Guesses a rational function in N rational in x for
#the data set S1 of the form {[i,RAT(N)]}:
GR1Rat:=proc(S1,N,x)
local Deg1,Deg2,i,i1,S1Top,S1Bot,R1,c,Top,Bot:

if {seq(S1[i1][2],i1=1..nops(S1))}={0} then
   RETURN(0):
fi:

Deg1:=max(seq(degree(numer(S1[i1][2]),N),i1=1..nops(S1))):

Deg2:=max(seq(degree(denom(S1[i1][2]),N),i1=1..nops(S1))):

S1Top:={}:
S1Bot:={}:

for i from 1 to nops(S1) do
R1:=S1[i][2]:
 if degree(numer(R1),N)=Deg1 and degree(denom(R1),N)=Deg2  then
 c:=coeff(numer(R1),N,Deg1):
 S1Top:=S1Top union {[S1[i][1],numer(R1)/c]}:
 S1Bot:=S1Bot union {[S1[i][1],denom(R1)/c]}:
 fi:
od:
Top:=GR1Pol(S1Top,N,x):
Bot:=GR1Pol(S1Bot,N,x):

 if Top<>FAIL and Bot<>FAIL then
  RETURN(Top/Bot):
 else
  RETURN(FAIL):
 fi:

end:


#####Direct Way for GH2M: GH2Mdirect

#GenPol2(x,y,d,a): a generic Polynomial
#of dgree d in (x,y) with coeffs. indexed by a
#followed by the coefficient GenPol2(x,y,1,a)
#For example, try
GenPol2:=proc(x,y,d,a) local i,k:
convert([seq(seq(a[i,k-i]*x^i*y^(k-i),i=0..k),k=0..d)],`+`),
{seq(seq(a[i,k-i],i=0..k),k=0..d)}:
end:

#GenOper2(m,n,M,deg,ORD,a): a generic operator of degree deg in (m,n)
#and order ORD in M
GenOper2:=proc(m,n,M,deg,ORD,a) local i,var,gu,lu:
gu:=0:
var:={}:
for i from 0 to ORD do
lu:=GenPol2(m,n,deg,a[i]):
gu:=gu+lu[1]*M^i:
 var:=var union lu[2]:
od:
gu,var:
end:





#GH2Mdir1(F,m,n,M,deg,ORD), guesses an operator ope(m,n,M) of degree deg in (m,n)
#and order ORD in M annihilating the function F
#for example, try
#GH2Mdir1((i,j)->(i+j!,m,n,M,2,2);
GH2Mdir1:=proc(F,m,n,M,deg,ORD) local ope,var,eq,a,m0,n0,k,N,var1,ind1,i:
ope:=GenOper2(m,n,M,deg,ORD,a):
var:=ope[2]:
ope:=ope[1]:

eq:={}:
for k from 5 while nops(eq)<=nops(var)+5 do
for m0 from 0 to k do
 n0:=k-m0:
 eq:=eq union {EvalOpeF2(F,ope,m,n,M,N,m0+2,n0+4)}:
od:
od:

var1:=Linear(eq,var):

ind1:={}:

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

if nops(ind1)=0 then
 RETURN(FAIL):
fi:
if nops(ind1)>1 then
 RETURN(FAIL):
fi:

ope:=subs(ind1[1]=1,subs(var1,ope)):
if {seq(seq(EvalOpeF2(F,ope,m,n,M,N,m0,n0),m0=10..20),n0=10..20)}<>{0} then
 RETURN(FAIL):
else
RETURN(Yafe1(ope,M,N)):
fi:

end:



#GH2Mdir(F,m,n,M), guesses an operator ope(m,n,M) of degree deg in (m,n)
#and order ORD in M annihilating the function F
#for example, try
#GH2Mdir((i,j)->(i+j!,m,n,M);
GH2Mdir:=proc(F,m,n,M) local ope,DEG,ORD,k,LT,i,ope1,ope2,ope3,ope4,j:

ope1:=numer(normal(FindrecF(i->F(i,20),n,N,MaxC))):
ope2:=numer(normal(FindrecF(i->F(i,21),n,N,MaxC))):
ope3:=numer(normal(FindrecF(i->F(i,22),n,N,MaxC))):

for j from 23 do
if nops({degree(ope1,N),degree(ope2,N),degree(ope3,N)})=1 and
  nops({degree(ope1,n),degree(ope2,n),degree(ope3,n)})=1 then
   DEG:=degree(ope3,n):
   ORD:=degree(ope3,N):
 break:
fi:
ope4:=numer(normal(FindrecF(i->F(i,j),n,N,MaxC))):
ope1:=ope2:
ope2:=ope3:
ope3:=ope4:
od:

 ope:=GH2Mdir1(F,m,n,M,DEG,ORD):
  if ope<>FAIL then
 LT:=coeff(ope,M,degree(ope,M)):
 ope:=add(normal(coeff(ope,M,i)/LT)*M^i,i=0..degree(ope,M)):
    RETURN(ope):
  fi:

FAIL:
end:


#GH2MPolDir(POL,x,y,m,n,M,MaxC): The recurrence in the m-direction  of
#the coeff. of x^m*y^n of 1/POL(x,y) the direct way
#For example, try: GH2MPolDir(1-x-y,x,y,m,n,M,10);
GH2MPolDir:=proc(P,x,y,m,n,M,MaxC) local i,j:
if expand(subs({x=y,y=x},P)-P)=0 then
GH2Mdir((i,j)->CoeffPs(P,[x,y],[i,j]),m,n,M,MaxC):
else
GH2Mdir((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,MaxC):
fi:
end:



#GH2D(F,m,n,M,N): Given a function of two variables (m,n) guesss a recurrence
#in the diagonal direction (MN), where M and N are the shift-operator in the m-direction 
#and n-directions 
#For example, try: GH2D((i,j)->i!+j!,m,n,M,N);
GH2D:=proc(F,m,n,M,N)
local i,ope,h,ope1,ope2,ope3,DEG,ORD,gu1,gu2,gu3,S1,
Ope,j,K,x,Ope1,m0,n0,i1:

ope[1]:=FindrecF(h->F(h+1,h),n,K):
ope[2]:=FindrecF(h->F(h+2,h),n,K):
ope[3]:=FindrecF(h->F(h+3,h),n,K):
ope1:=numer(normal(ope[1])):
ope2:=numer(normal(ope[2])):
ope3:=numer(normal(ope[3])):

for i from 4 do

ope1:=numer(normal(ope[i-1])):
ope2:=numer(normal(ope[i-2])):
ope3:=numer(normal(ope[i-3])):
  if nops({degree(ope1,n),degree(ope2,n),degree(ope2,n)})=1 and
   nops({degree(ope1,K),degree(ope2,K),degree(ope2,K)})=1  then
   DEG:=degree(ope3,n):
   ORD:=degree(ope3,K):

gu1:=[seq(F(h+i-1,h),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
gu2:=[seq(F(h+i-2,h),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
gu3:=[seq(F(h+i-3,h),h=1..(1+DEG)*(1+ORD)+5+ORD)]:

if findrecK(gu1,DEG,ORD) and findrecK(gu2,DEG,ORD) and
    findrecK(gu3,DEG,ORD) then
      break:

  fi:
fi:

ope[i]:=FindrecF(h->F(h+i,h),n,K):

od:

S1:={}:

for j from i-2 to i+6 do

  gu1:=[seq(F(h+j,h),h=1..(1+DEG)*(1+ORD)+5+ORD)]:

   if findrecK(gu1,DEG,ORD) then
    ope1:=findrec(gu1,DEG,ORD,n,K):
    S1:=S1 union {[j,ope1]}:

   fi:

od:


for j from i+7 do
Ope:= GR1PolRat(S1,K,x,n):
 if Ope<>FAIL then

 Ope:=subs({x=m-n},Ope):
Ope:=add(factor(normal(coeff(Ope,K,i1)))*K^i1,i1=0..degree(Ope,K)):
 Ope:=subs({K=M*N},Ope):

 Ope1:=numer(normal(Ope)):
 if {seq(seq(normal(EvalOpeF2(F,Ope1,m,n,M,N,m0,n0)),m0=0..40),n0=0..40)}<>{0} then
  print(Ope, ` failed `):
 else
  RETURN(Ope):
 fi:
 fi:

  gu1:=[seq(F(h,j),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
   if findrecK(gu1,DEG,ORD) then
    ope1:=findrec(gu1,DEG,ORD,n,K):
    S1:=S1 union {[j,ope1]}:
   fi:

od:
end:


#GH2M(F,m,n,M,Patience,s1): 
#Given a function of two variables (m,n) guesses a recurrence
#in the m direction, where M is the shift-operator in the m-direction 
#It does it the fast way (using GR1Rat instead of GR1)
#s1 is where it starts exploring
#For example, try: GH2M((i,j)->i!+j!,m,n,M,10);
GH2M:=proc(F,m,n,M,Patience,s1)
local i,ope,h,ope1,ope2,ope3,ope4,DEG,ORD,gu1,gu2,gu3,gu4,S1,Ope,j,
m0,n0,N,Ope1:
ope[1]:=FindrecF(h->F(h,s1+1),m,M):
ope[2]:=FindrecF(h->F(h,s1+2),m,M):
ope[3]:=FindrecF(h->F(h,s1+3),m,M):
ope[4]:=FindrecF(h->F(h,s1+4),m,M):
ope1:=numer(normal(ope[1])):
ope2:=numer(normal(ope[2])):
ope3:=numer(normal(ope[3])):
ope4:=numer(normal(ope[4])):



for i from 5 to Patience do

ope1:=numer(normal(ope[i-1])):
ope2:=numer(normal(ope[i-2])):
ope3:=numer(normal(ope[i-3])):
ope4:=numer(normal(ope[i-4])):


  if nops({degree(ope1,m),degree(ope2,m),degree(ope3,m),degree(ope4,m)})=1 and
   nops({degree(ope1,M),degree(ope2,M),degree(ope3,M),degree(ope4,M)})=1  then
   DEG:=degree(ope4,m):
   ORD:=degree(ope4,M):

gu1:=[seq(F(h,s1+i-1),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
gu2:=[seq(F(h,s1+i-2),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
gu3:=[seq(F(h,s1+i-3),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
gu4:=[seq(F(h,s1+i-4),h=1..(1+DEG)*(1+ORD)+5+ORD)]:


if findrecK(gu1,DEG,ORD) and findrecK(gu2,DEG,ORD) and
    findrecK(gu3,DEG,ORD)  and  findrecK(gu4,DEG,ORD)  then
      break:
  fi:
fi:
ope[i]:=FindrecF(h->F(h,i+s1),m,M):
od:

if i-1=Patience then
 print(`Couldn't do it with patience`, Patience, `and Starting place`, s1):
 RETURN(FAIL):
fi:

S1:={}:

for j from i-3 to i+6 do
  gu1:=[seq(F(h,j+s1),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
   if findrecK(gu1,DEG,ORD)  then
    ope1:=findrec(gu1,DEG,ORD,m,M):
    S1:=S1 union {[j+s1,ope1]}:
   fi:
od:

if nops( {seq(S1[i][2],i=1..nops(S1))})=1 then
   RETURN(S1[1][2]):
fi:

for j from i+7 do
Ope:= GR1PolRat(S1,M,n,m):
 if Ope<>FAIL then
  Ope1:=numer(normal(Ope)):

 if {seq(seq(normal(EvalOpeF2(F,Ope1,m,n,M,N,m0,n0)),m0=0..40),n0=0..40)}<>{0} then
  print(Ope, ` failed `):
  RETURN(FAIL):
 else

   RETURN(Ope):
 fi:
fi:
  gu1:=[seq(F(h,j+s1),h=1..(1+DEG)*(1+ORD)+5+ORD)]:

   if findrecK(gu1,DEG,ORD) then
    ope1:=findrec(gu1,DEG,ORD,m,M):
    S1:=S1 union {[j+s1,ope1]}:
   fi:

od:
end:


GH2Pol:=proc(P,x,y,m,n,M,N,Patience,s1) local ope:
if expand(subs({x=y,y=x},P)-P)=0 then
ope:=GH2M((i,j)->CoeffPs(P,[x,y],[i,j]),m,n,M,Patience,s1):
RETURN(ope,subs({m=n,n=m,M=N},ope)):
else
GH2M((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,Patience,s1),
GH2M((i,j)->CoeffP(P,[x,y],[j,i]),m,n,N,Patience,s1)
fi:
end:


#AF(F,n,k,m0,n0): Given a function F(n,k) finds
#sum(F(n0+m0,k),k=0..m0)
AF:=proc(F,n,k,m0,n0)
option remember:

if m0<0 then
0:
else
AF(F,n,k,m0-1,n0+1)+eval(subs({k=m0,n=m0+n0},F)):
fi:
end:

#HypS(F,n,k,m,M,Patience,s1): The Operator Oper(M,m,n) annihilationg
#sum(F(n+m,k),k=0..m)
#For example, try: HypS(binomial(n,k),n,k,m,M,10,10);
HypS:=proc(F,n,k,m,M,Patience,s1):
GH2M((i,j)->AF(F,n,k,i,j),m,n,M,Patience,s1):
end:






 
Zinn:=proc(resh)
local s1,s2,n:
n:=nops(resh)-2:
s1:=sn(resh,n):
s2:=sn(resh,n-1):
evalf(2*(s1+s2)/(s1-s2)^2),
evalf(sqrt(op(n+1,resh)/op(n-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))):
end:
 
 
sn:=proc(resh,n):
-1/log(op(n+1,resh)*op(n-1,resh)/op(n,resh)^2):
#evalf(%):
end:
 

#CanOpe(Ope,m,n,M,N): Given a partial  linear-recurrence operator with polynomial coefficients
#divides it by an appropriate shiftp operaor to find an equivalent polynomial
#in the form of a polynomial in (m,n,1/M,1/N)
#For example, try: CanOpe(M*N-M-N,m,n,M,N);
CanOpe:=proc(Ope,m,n,M,N) local gu,i,j,d1,d2,Ope1,i1,j1:
Ope1:=numer(Ope):
d1:=degree(Ope1,M):
d2:=degree(Ope1,N):
Ope1:=expand(subs({m=m-d1,n=n-d2},Ope1)/M^d1/N^d2):

gu:=0:

for i1 from ldegree(Ope1,M) to degree(Ope1,M) do
 for j1 from ldegree(Ope1,N) to degree(Ope1,N) do
   gu:=gu+factor(coeff(coeff(Ope1,M,i1),N,j1))*M^i1*N^j1:
 od:
od:
gu:
end:

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



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:
 

#FindrecR(f,n,N,MaxC): 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,2);
FindrecR:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:

for L from 0 to MaxC do
for ORDER from 1 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:=findrecR([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL then
       RETURN(ope):
      fi:
 od:
od:
FAIL:

end:




#Tp(f,x,p): The transformation T_p(f(x))
Tp:=proc(f,x,p)
local i:
normal(
subs({seq(x[i]=p*x[i]/(1-(1-p)*x[i]),i=1..nops(x))},f)/
mul(1-(1-p)*x[i],i=1..nops(x))):
end:

#Sp(f,x,p): The transformation S_p(f(x))
Sp:=proc(f,x,p)
local i:
normal(
subs({seq(x[i]=x[i]/(p+(1-p)*x[i]),i=1..nops(x))},f)/
mul(1+(1-p)*x[i]/p,i=1..nops(x))):
end:






#ApplyOper(Seq1,ope,n,N,i0): Given a sequence Seq1, applies the operator 
#ope(n,N) at Seq1[i]
ApplyOper:=proc(Seq1,ope,n,N,i0) local i:
add(Seq1[i0+i]*subs(n=i0,coeff(ope,N,i)),i=ldegree(ope,N)..degree(ope,N)):
end:


ez:=proc():print(`Mult1(ope1,ope2,n,N), Div1(ope1,ope2,n,N)`):end:

#Div1(ope1,ope2,n,N): Given two operators ope1(n,N) and ope2(n,N)
#outputs operators Q(n,N) and R(n,N) such that ope1=Q*ope2+R
#and the order of R is less than the order of ope2
Div1:=proc(ope1,ope2,n,N) local d1,d2,i,ope1a,gu:
d1:=degree(ope1,N):
d2:=degree(ope2,N):

if d1<d2 then
 RETURN(0,ope1):
fi:

ope1a:=expand(ope1- coeff(ope1,N,d1)/subs(n=n+d1-d2,coeff(ope2,N,d2))
                  *N^(d1-d2)*subs(n=n+d1-d2,ope2) ):
ope1a:=add(normal(coeff(ope1a,N,i))*N^i,i=0..degree(ope1a,N)):
ope1a:

gu:=Div1(ope1a,ope2,n,N):

normal(coeff(ope1,N,d1)/subs(n=n+d1-d2,coeff(ope2,N,d2)))*N^(d1-d2)+gu[1],gu[2]:
end:




#Mult1(ope1,ope2,n,N): ope1(n,N)*ope2(n,N)
Mult1:=proc(ope1,ope2,n,N) local d1,i1,gu:
d1:=degree(ope1,N):

gu:=0:
for i1 from 0 to d1 do
 gu:=expand(gu+coeff(ope1,N,i1)*N^i1*subs(n=n+i1,ope2)):
od:

add(normal(coeff(gu,N,i1))*N^i1,i1=0..degree(gu,N)):
end:



###-----------------3D tableaux
ezY3:=proc(): print(`IsPar(L), IsGood(L), Shrink(L), Yeladim(L), Y3(L)`): end:

IsPar:=proc(L) local i:
if member(0,{op(L)}) then
 RETURN(false):
fi:
for i from 1 to nops(L)-1 do
 if L[i]-L[i+1]<0 then
   RETURN(false):
 fi:
od:
true:
end:

#IsGood(L): Given a list of partitions (each of them a list) outputs whether it is good
IsGood:=proc(L) local i,j,sh:
for i from 1 to nops(L) do
 if not IsPar(L[i]) then
  RETURN(false):
 fi:
od:

sh:=[seq(nops(L[i]),i=1..nops(L))]:
if not IsPar(sh) then
  RETURN(false):
fi:

for i from 1 to nops(L)-1 do
for j from 1 to nops(L[i]) do
 if j<=nops(L[i+1]) and L[i][j]-L[i+1][j]<0 then
     RETURN(false):
 fi:
od:
od:
true:
end:


Shrink1:=proc(L):

if L[nops(L)]=0 then
 RETURN([op(1..nops(L)-1,L)]):
fi:

L:
end:

Shrink:=proc(L) local i,L1:

L1:=[seq(Shrink1(L[i]),i=1..nops(L))]:

if L[nops(L1)]=[] or L[nops(L1)]=[0] then
L1:=[op(1..nops(L1)-1,L1)]:
fi:
L1:
end:






#Yeladim(L): all the legal children of the PP L
Yeladim:=proc(L) local gu,L1,i,j:


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

gu:={}:
for i from 1 to nops(L) do
for j from 1 to nops(L[i]) do
 L1:=[op(1..i-1,L),[op(1..j-1,L[i]),L[i][j]-1,op(j+1..nops(L[i]),L[i])],
op(i+1..nops(L),L)]:

 L1:=Shrink(L1):
  if IsGood(L1) then
    gu:=gu union {L1}:
  fi:
od:
od:
gu:
end:



#Y3(L):The number of 3-D YT of shape L
Y3:=proc(L) local gu,i:
option remember:


if L=[] then
 RETURN(1):
fi:

if not IsGood(L) then
 RETURN(0):
fi:

gu:=Yeladim(L):

add(Y3(gu[i]) ,i=1..nops(gu)):
end:



##------------end 3D tableaux

###end GuessHolo without its ezra

#Sipur(k,n,N,a,L,L1): The story about (k+1)-non-crossing 
#tangled diagrams. For example, try:
#Sipur(2,n,N,a,100,1000);
Sipur:=proc(k,n,N,a,L,L1)
local gu,ope,C,i,lu,gu1:
print(`The sequence counting the number`):
print(`of ` , k+1, ` - non-crossing tangled diagrams  of [n] `):
print(`(allowing isolated points) (let's call it a(n))`):
print(`from 1 to` , L, ` is: ` ):
gu:=SeqD(k,L):
print(gu):
C:=trunc(2*(sqrt(L-4)-3/2))-1:
ope:=Findrec(gu,n,N,C):

if ope<>FAIL then
 print(`This sequence is P-finite, and is annihilated by the linear`):
 print(`recurrence operator (N is the shift operator in n)`):
 lprint(ope):
 print(`In everyday language, this means that a(n), satisfies the recurrence`):
 print(add(coeff(ope,N,i)*a(n+i),i=0..degree(ope,N)) =0  ):
 lu:=Asy(ope,n,N,3):

  if lu<>FAIL then
   print(`The asymptotics in n, is a constant times`):
   print(lu):
  fi:
 print(`Using the recurrence we can crank out many more terms`):
 print(`of the sequence`):

gu1:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],L1):

if gu<>[op(1..nops(gu),gu1)] then
 ERROR(`Something terrible happended `):
fi:

print(gu1):

fi:

end: