######################################################################
## DyckClever.txt Save this file as  DyckClever.txt  to use it,      #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read DyckClever.txt <Enter>                       #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  DoronZeil at gmail dot com                                       # 
######################################################################

with(plots):
with(combinat):

print(`First Written May 2020: tested for Maple 18 `):
print(`Version   May 22, 2020  `):
print():
print(`This is DyckClever.txt, One of the  Maple packages`):
print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's  article: `):
print(`Automatic Counting of Restricted Dyck Paths via (Numeric and Symbolic) Dynamic Programming `):
print(`Available from:`):
print(` http://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/dyck.html .`):
print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/Dyck.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):



print():
print(`For a list of the  STORY functions,`):
print(` type "ezraStory();". For specific help type "ezra(procedure_name);" `):
print():


print():
print(`For a list of the  Checking functions,`):
print(` type "ezraCheck();". For specific help type "ezra(procedure_name);" `):
print():



 
print():
print(`For a list of the SUPPORTING functions,`):
print(` type "ezra1();". For specific help type "ezra(procedure_name);" `):
print():

 

print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
#print(`For a list of the supporting functions type: ezra1();`):
print():
 




ezraCheck:=proc()
if args=NULL then

print(` The Checking procedures are :  CheckfAB, CheckfABr, CheckfCD, CheckfCDr `):

else
ezra(args):
fi:

end:

ezraStory:=proc()
if args=NULL then

print(` The Story procedures are : `):

else
ezra(args):
fi:

end:

ezra1:=proc()
if args=NULL then

print(` The supporting procedures are : algtorec,  CheckOpe, ChildAB, ChildABr, ChildCD, OperToRecEq, SeqFromRec `):

else
ezra(args):

fi:


end:


ezra:=proc()
if args=NULL then

print(` DyckClever.txt:  A Maple package for automatic enumeration of families of Dyck paths obeying certain kinds of restrictions  `):
print(`but NOT using the KISS method, but rather symbolic Dynamic Programing `):

print(`The MAIN procedures are: fAB, fABcat, fABr, fCD, fCDr, InfoAB, InfoABr `):
print(` `):



elif nargs=1 and args[1]=algtorec then
print(`algtorec(F,P,x,n,N): Inputs a polynomial F in the two variables P and x, where P is the ordinary generating function`):
print(`of a sequence, let's call it a(n), that satifies the algebraic equation. To get a recurrence operator`):
print(`for the Catalan numbers when they are defined by P=1+x*P^2, try`):
print(`algtorec(P-1-x*P^2,P,x,n,N);`):

elif nargs=1 and args[1]=CheckfAB then
print(`CheckfAB(A,B,K): Checks procedure fAB(A,B,P,x) up to K terms. Try:`):
print(`CheckfAB({1},{1},50);`):


elif nargs=1 and args[1]=CheckfABr then
print(`CheckfABr(A,B,r,K): Checks procedure fABr(A,B,r,P,x) up to K terms. Try:`):
print(`CheckfABr({2*r+1},{3*r+1},r,50);`):


elif nargs=1 and args[1]=CheckfCD then
print(`CheckfCD(C,D,K): Checks procedure fCD(C,D,P,x) up to K terms. Try:`):
print(`CheckfCD({1},{1},50);`):


elif nargs=1 and args[1]=CheckfCDr then
print(`CheckfCDr(C,D,r,K): Checks procedure fCDr(C,D,r,P,x) up to K terms. Try:`):
print(`CheckfCDr({2*r+1},{},r,50);`):

elif nargs=1 and args[1]=CheckOpe then
print(`CheckOpe(L,ope,n,N):plugs in the sequence to the operator ope, in terms of n and N. Try: `):
print(`CheckOpe([1$10],N-1,n,N);`):

elif nargs=1 and args[1]=ChildAB then
print(`ChildAB(A1,B1,x,Z): inputs a pair of sets of non-negative integers A1,B1, and symbols x and Z`):
print(`outputs their "children" A2,B2, and the resulting equation. Try:`):
print(`ChildAB({1},{1},x,Z);`):



elif nargs=1 and args[1]=ChildABr then
print(`ChildABr(A1,B1,r,x,Z): inputs a pair of sets of affine-linear expressions A1,B1, in the variable r, and symbols x and Z`):
print(`outputs their "children" A2,B2, and the resulting equation. Try:`):
print(`ChildABr({3*r+1},{5*r+1,7*r+1},x,Z);`):


elif nargs=1 and args[1]=ChildCD then
print(`ChildCD(STATE,Z,x): Inputs a state STATE=[C,C1,D,D1,IR]: where IR is 0 ot 1, outputs`):
print(`the equation satisfied by Z[STATE] and the set of children states that feature in the`):
print(`formula. Try:`):
print(`ChildCD([{1},{1},{1},{1},0],Z,x);`):

elif nargs=1 and args[1]=fAB then
print(`fAB(A,B,x,P): inputs sets of positive integers A, and B, and variables x and P`):
print(` and outputs a polynomial in x and P , let's call it F(x,P) such that`):
print(`if you plug-in for P the  ordinary generating function for`):
print(`the sequence "number of Dyck path of semi-length n with no peak-heights in A `):
print(`and no valley-heights in B `):
print(`you would get 0. `):
print(`Try: `):
print(`fAB({1},{1},x,P); `):



elif nargs=1 and args[1]=fABr then
print(`fAB(A,B,r,x,P): inputs sets of affine linear expressions A, and B in thevariable r, and variables x and P`):
print(` and outputs a polynomial in x and P , let's call it F(x,P) such that`):
print(`if you plug-in for P the  ordinary generating function for`):
print(`the sequence "number of Dyck path of semi-length n with no peak-heights in the range of A `):
print(`and no valley-heights in the range of B `):
print(`you would get 0. `):
print(`Try: `):
print(`fABr({2*r+1},{3*r+1},r,x,P); `):


elif nargs=1 and args[1]=fABcat then
print(`fABcat(A,B,x,C): inputs sets of positive integers A, and B, and variables x and C`):
print(` and outputs a polynomial in x and C , where C is the Catalan generating function`):
print(`Try: `):
print(`fABcat({1},{1},x,C); `):

elif nargs=1 and args[1]=fCD then
print(`fCD(C,D,x,P): Inputs finite sets of positive integers C and D and symbols x and P. Outputs `):
print(`The algebraic equation satisfied by the ordinary generating function, P(x), of the`):
print(`sequence enumerating Dyck paths with no upward runs in C and no downward runs in D. Try:`):
print(`fCD({1},{1},x,P);`):


elif nargs=1 and args[1]=fCDr then
print(`fCDr(C,D,r,x,P): Inputs finite sets of arithmetical progressions of the form ar+b, with a and b positive integers`):
print( `C and D and symbols r, x and P. Outputs `):
print(`The algebraic equation satisfied by the ordinary generating function, P(x), of the`):
print(`sequence enumerating Dyck paths with no upward runs in the range of C and no downward runs in the range of D. Try:`):
print(`fCDr({2*r+2},{2*r+2},r,x,P);`):




elif nargs=1 and args[1]=InfoAB then
print(`InfoAB(A,B,x,P,C,a,n,M): inputs two sets of positive integers A and B, symbols x, P , C, a, and n, and a positive  integer`):
print(`M outputs information about the integer sequence, let's call it a(n)`):
print(`that enumerates Dyck paths of semi-length n none of whose peaks have heights in the set A`):
print(`and none of of whose valleys have heights in the set B. `):
print(`(i) The first entry is the algebraic equation satisfied by the ordinary generating function P(x)=Sum(a(n)*x^n,n=0..infinity)`):
print(` (ii) The second entry is the expression in terms of C=(1-sqrt(1-4x))/(2*x), the Catalan generating function `):
print(`(iii) the third entry is the recurrence and initial conditions satisfied by a(n)`):
print(`(iv):The fourth entry is the list of first 50 terms from n=1 to n=50 (for the OEIS)`):
print(` (v) the fifth entry is the  value of a(M), just for fun. Try: `):
print(` InfoAB({1},{1},x,P,C,a,n,1000); `):



elif nargs=1 and args[1]=InfoABr then
print(`InfoABr(A,B,r,x,P,a,n,M): Inputs two sets of affine-linear expressions in the variable r, `):
print(`A and B, symbols x and P, a, and n, and a positive  integer`):
print(`M. Outputs information about the integer sequence, let's call it a(n)`):
print(`that enumerates Dyck paths of semi-length n none of whose peaks have heights in the range of the set A`):
print(`and none of of whose valleys have heights in the range of the set B. `):
print(`(i) The first entry is the algebraic equation satisfied by the ordinary generating function P(x)=Sum(a(n)*x^n,n=0..infinity)`):
print(`(ii) the second entry is the recurrence and initial conditions satisfied by a(n)`):
print(`(iii):The third entry is the list of first 50 terms from n=1 to n=50 (for the OEIS)`):
print(` (iv) the fourthh entry is the  value of a(M), just for fun. Try: `):
print(` InfoABr({2*r+3},{},r,x,P,a,n,1000); `):


elif nargs=1 and args[1]=OperToRecEq then
print(`OperToRecEq(ope,n,N,a): inputs a linear recurrence operator ope in n and the forward shift operator N,`):
print(`and a symbol a, outputs the linear recurrence equation that asserts that the sequence is annihilated by`):
print(`the operator, in terms of decreasing arguments of n. Try`):
print(`OperToRecEq(N-(n+1),n,N,a);`):

elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,Ini,K): The first K terms of the sequence given by the recurrence ope. Try:`):
print(`SeqFromRec(N^2-N-1,n,N,[1,1],20);`):


else
 print(`There is no such thing as`, args):

fi:


end:

#ChildAB(A1,B1,x,Z): inputs a pair of sets of non-negative integers A1,B1, and symbols x and Z
#outputs their "children" A2,B2, and the resulting equation. Try:
#ChildAB({1},{1},x,Z);
ChildAB:=proc(A1,B1,x,Z) local A2,B2,eq,m:

if member(0,A1) then
A2:=A1 minus {0}:
B2:=B1:
eq:=Z[A1,B1]-Z[A2,B1]+1:
else
A2:={seq(m-1, m in A1)}:
B2:={seq(m-1, m in B1)} minus {-1}:

if member(0,B1) then
  eq:=Z[A1,B1]-1-x*Z[A2,B2]:
else
 eq:=Z[A1,B1]-1-x*Z[A2,B2]*Z[A1,B1]:
fi:

fi:
[A2,B2,eq]:
end:



#fAB(A,B,x,P): inputs sets of positive integers A, and B, and variables x and P
# and outputs a polynomial in x and P , let's call it F(x,P) such that
#if you plug-in for P the  ordinary generating function for
#the sequence "number of Dyck path of semi-length n with no peak-heights in A 
#and no valley-heights in B
#you would get 0.
#Try:
#fAB({1},{1},x,P);
fAB:=proc(A,B,x,P) local eq,var,var1,Z,A1,B1,mu,gu:
A1:=A:
B1:=B:
var:={}:
var1:={Z[A1,B1]}:
eq:={}:

while var1<>var do

var:=var1:
mu:=ChildAB(A1,B1,x,Z):

var1:=var1 union {Z[mu[1],mu[2]]}:

eq:=eq union {mu[3]}:
A1:=mu[1]:
B1:=mu[2]:
od:

var1:=var minus {Z[A,B]}:

gu:=eliminate(eq,var1)[2]:

if nops(gu)<>1 then
 RETURN(FAIL):
fi:

gu:=gu[1]:

gu:=subs(Z[A,B]=P,gu):

add(factor(coeff(gu,P,i))*P^i,i=0..degree(gu,P)):

end:




#fABcat(A,B,x,C): inputs sets of positive integers A, and B, and variables x and P
# and outputs a polynomial in x and C , where C is the Catalan generating function
fABcat:=proc(A,B,x,C) local eq,var,var1,Z,A1,B1,mu,gu:
if A={} and B={} then
 RETURN(C):
fi:

A1:=A:
B1:=B:
var:={}:
var1:={Z[A1,B1]}:
eq:={}:

while var1<>var do

var:=var1:
mu:=ChildAB(A1,B1,x,Z):

var1:=var1 union {Z[mu[1],mu[2]]}:

eq:=eq union {mu[3]}:
A1:=mu[1]:
B1:=mu[2]:
od:


eq:=eq minus {Z[{},{}]-1-x*Z[{},{}]^2}:

var:=var minus {Z[{},{}]}:

eq:=subs(Z[{},{}]=C,eq):


var:=solve(eq,var):

gu:=subs(var,Z[A,B]):

collect(numer(gu),C)/collect(denom(gu),C):

end:


#ChildABr(A1,B1,r,x,Z)
ChildABr:=proc(A1,B1,r,x,Z) local A2,B2,eq,m,a1,b1:


if member(0,{seq(coeff(a1,r,0),a1 in A1)}) then
  A2:={}:
 for a1 in A1 do
    if coeff(a1,r,0)=0 then
     A2:=A2 union {subs(r=r+1,a1)}:
     else
      A2:=A2 union {a1}:
    fi:
  od:

B2:=B1:
eq:=Z[A1,B1]-Z[A2,B2]+1:
else

A2:={seq(m-1, m in A1)}:

if member(0,{seq(coeff(b1,r,0),b1 in B1)}) then
  B2:={}:
 for b1 in B1 do
    if coeff(b1,r,0)=0 then
     B2:=B2 union {subs(r=r+1,b1)}:
     else
      B2:=B2 union {b1}:
    fi:
  od:

 B2:={seq(m-1, m in B2)}:
eq:=Z[A1,B1]-1-x*Z[A2,B2]:
else

 B2:={seq(m-1, m in B1)}:

 eq:=Z[A1,B1]-1-x*Z[A2,B2]*Z[A1,B1]:
fi:

fi:
[A2,B2,eq]:
end:



#fABr(A,B,r,x,P): inputs sets of linear affine expressions A, and B in the variable r, and variables x and P
# and outputs a polynomial in x and P , let's call it F(x,P) such that
#if you plug-in for P the  ordinary generating function for
#the sequence "number of Dyck path of semi-length n with no peak-heights in A 
#and no valley-heights in B
#you would get 0.
#Try:
#fABr({2*r+1},{},r,x,P);
fABr:=proc(A,B,r,x,P) local eq,var,var1,Z,A1,B1,mu,gu:
if not ( type(A,set)  and type(B,set) and type(r,symbol)) then
 print(`Bad input`):
 RETURN(FAIL):
fi:


A1:=A:
B1:=B:
var:={}:
var1:={Z[A1,B1]}:
eq:={}:

while var1<>var do

var:=var1:
mu:=ChildABr(A1,B1,r,x,Z):

var1:=var1 union {Z[mu[1],mu[2]]}:

eq:=eq union {mu[3]}:
A1:=mu[1]:
B1:=mu[2]:
od:

var1:=var minus {Z[A,B]}:

gu:=eliminate(eq,var1)[2]:

if nops(gu)<>1 then
 RETURN(FAIL):
fi:

gu:=gu[1]:

gu:=subs(Z[A,B]=P,gu):

add(factor(coeff(gu,P,i))*P^i,i=0..degree(gu,P)):

end:



###Start FROM SCHUZENBERGER.txt



#OperToRecEq(ope,n,N,a): inputs a linear recurrence operator ope in n and the forward shift operator N,
#and a symbol a, outputs the linear recurrence equation that asserts that the sequence is annihilated by
#the operator, in terms of decreasing arguments of n. Try
#OperToRecEq(N-(n+1),n,N,a);

OperToRecEq:=proc(ope,n,N,a) local ope1,ORD,i:

ORD:=degree(ope,N):

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

add(factor(coeff(ope1,N,-i))*a(n-i),i=0..ORD)=0:

end:

 
simp:=proc(pa,deg,bitui,P,x)
local mone,mekh,gu:
gu:=normal(bitui):
mone:=numer(gu):
mekh:=denom(gu):
mone:=simp1(pa,deg,mone,P,x):
mekh:=simp1(pa,deg,mekh,P,x):
expand(mone)/expand(mekh):
end:
 
simp1:=proc(pa,deg,bitui,P,x)
local i,gu:
 
gu:=expand(bitui):
 
for i from deg to 3*deg+1 do
 gu:=subs(P(x)^i=pa[i],gu):
od:
 
gu:
 
end:

#MyGcd(L): inputs a list of polynomials and finds their gcd. Try
#MyGcd([n,n^2,n^3]);
MyGcd:=proc(L) local gu:
if L=[] then
 RETURN(FAIL):
elif nops(L)=1 then
 RETURN(L[1]):
else
 gu:=MyGcd([op(1..nops(L)-1,L)]):
 RETURN(gcd(gu,L[nops(L)])):
fi:
end:
 
pashet:=proc(p,n,N)
local i,gu1,gu,p1,ra:
p1:=normal(p):
gu1:=denom(p1):
ra:=degree(gu1,N):
p1:=subs(n=n+ra,numer(p1)):
p1:=expand(p1):
 
gu:=0:
for i from 0 to degree(p1,N) do
gu:=gu+factor(coeff(p1,N,i))*N^i:
od:
RETURN(gu):
end:
 
 
 
difftorec:=proc(ope1,D,x,n,N)
local i,j,gu,ord,ope,r,mu:
 
ope:=expand(ope1):
gu:=0:
ord:=degree(ope,D):
 
for i from 0 to  ord do
 mu:=coeff(ope,D,i):
 
for j from 0 to degree(mu,x) do
  gu:=gu+coeff(mu,x,j)*product(n+r-j,r=1..i)*N^(i-j):
od:
od:
 
 
pashet(gu,n,N):
end:
 
algtodiff:=proc(F,P,x,ORDER)
local i,gu,a,lu,pa,deg,eq1,lu1,y,var,eq,ope,pip:
deg:=degree(F,P):
pa:=array(deg..3*deg+1):
eq1:=subs(P^deg=y,F):
pa[deg]:=solve(eq1=0,y):
pa[deg]:=subs(P=P(x),pa[deg]):
 
for i from 1 to 2*deg+1 do
   lu:=pa[deg+i-1]*P(x):
   lu:=expand(lu):
   lu:=subs(P(x)^deg=pa[deg],lu):
   pa[deg+i]:=lu:
od:
 
 
gu:=a[0]*P(x):
lu:=subs(P=P(x),F):
lu1:=diff(lu,x):
lu1:=subs(diff(P(x),x)=y,lu1):
lu1:=solve(lu1=0,y):
lu1:=simp(pa,deg,lu1,P,x):
 lu1:=normal(lu1):
 
 
gu:=gu+a[1]*lu1:
gu:=normal(gu):
 gu:=simp(pa,deg,gu,P,x):
 
var:={a[0],a[1]}:
ope:=a[0]+a[1]*D:
lu:=lu1:
 for i from 2 to ORDER do
 lu:=diff(lu,x):
  lu:=subs(diff(P(x),x)=lu1,lu):
  lu:=normal(lu):
 lu:=simp(pa,deg,lu,P,x):
 gu:=gu+a[i]*lu:
 gu:=normal(gu):
 gu:=simp(pa,deg,gu,P,x):
 var:=var union {a[i]}:
 
 ope:=ope+a[i]*D^i:
od:
 
gu:=normal(gu):
 
 
gu:=numer(gu):
 
gu:=expand(gu):
eq:={}:
 
for i from 0 to deg-1 do
 eq:=eq union {coeff( gu,P(x),i)=0}:
od:
 
 
var:=solve(eq,var):

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

ope:=subs(var,ope):

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

ope:=normal(ope):
ope:=numer(ope):
ope:=expand(ope):
pip:=coeff(ope,D,degree(ope,D)):
 
pip:=coeff(pip,x,degree(pip,x)):
normal(ope/pip):
end:
 

#algtorec(F,P,x,n,N): Inputs a polynomial F in the two variables P and x, where P is the ordinary generating function
#of a sequence, let's call it a(n), that satifies the algebraic equation. To get a recurrence operator
#for the Catalan numbers when they are defined by P=1+x*P^2, try
#algtorec(P-1-x*P^2,P,x,n,N);
algtorec:=proc(F,P,x,n,N) local ope,d,i, mu:

if not type(F,`+`) then
  RETURN(FAIL):
fi:

for d from 1 to degree(F,P) do
ope:=algtodiff(F,P,x,d):
 if ope<>FAIL then
   ope:=difftorec(ope,D,x,n,N):
  mu:=MyGcd([seq(coeff(ope,N,i),i=0..degree(ope,N))]):
  ope:=normal(ope/mu):
  ope:=add(factor(coeff(ope,N,i))*N^i,i=0..degree(ope,N)):
 RETURN(ope):   
 fi:
od:
FAIL:
end:

###End FROM SCHUZENBERGER.txt

###From Dyck.txt


#IsRange1(f,k,n): inputs a linear expression f in k and a non-negative integer n
#decides whether f(k0)=n for some non-negative integer k0.
#Try:
#IsRange1(3*k+1,k,7);

IsRange1:=proc(f,k,n) local k0:
k0:=solve(f=n,k):

if k0>=0 and type(k0,integer) then
 true:
else
 false:
fi:
end:


#IsRange(F,k,n): inputs a set of linear expressions F={f} in k and a non-negative integer n
#decides whether f(k0)=n for some non-negative integer k0 and some f in F
#Try:
#IsRange({3*k+1,3*k+2},k,9);

IsRange:=proc(F,k,n) local f:

for f in F do
 if IsRange1(f,k,n) then
  RETURN(true):
 fi:
od:
false:
end:




#NRVr(m,n,A,B,C,D,r):  Inputs positive integers m>=n and sets A and B of linear expressions in the variable r
#Outputs
#The number of walks from the origin to [m,n] ending with a vertical step
#where no H to  V corner [i,j] can be such that i-j belongs to the range of A and
#no V to H corner [i,j] can be such that  i-j belongs to the range of B. Try:
#
#NRVr(10,10,{2*r+3},{},{2*r+1},{2*r+5},r);
NRVr:=proc(m,n,A,B,C,D,r) local gu,k:
option remember:

if (m<0 or n<0 or m<n) then
 RETURN(0):
fi:

if m=0 and n=0 then
 RETURN(1):
fi:

if IsRange(B,r,m-n) then
  RETURN(0):
fi:

gu:=0:

for k from 1 to n do
 if not IsRange(D,r,k) then
   gu:=gu+NRHr(m,n-k,A,B,C,D,r):
 fi:
od:

gu:

end:





#NRHr(m,n,A,B,C,D,r):  Inputs positive integers m>=n and sets A and B of linear expressions in the variable r
#Outputs
#The number of walks from the origin to [m,n] ending with a horizontal step
#where no H to  V corner [i,j] can be such that i-j belongs to the range of A and
#no V to H corner [i,j] can be such that  i-j belongs to the range of B. Try:
#
#NRHr(10,10,{2*r+3},{},{},{}, r);
NRHr:=proc(m,n,A,B,C,D,r) local gu,k:
option remember:

if (m<0 or n<0 or m<n) then
 RETURN(0):
fi:

if m=0 and n=0 then
 RETURN(0):
fi:


if IsRange(A,r,m-n) then
  RETURN(0):
fi:

gu:=0:

for k from 1 to m-n do
 if not IsRange(C,r,k) then
  gu:=gu+NRVr(m-k,n,A,B,C,D,r):
 fi:

od:

gu:

end:

#NRV(m,n,A,B,C,D,): The number of walks from the origin to [m,n] ending with a vertical step
#where no H to  V corner [i,j] can be such that i-j belongs to A and
#no V to H corner [i,j] can be such that  i-j belongs to B. 
#No upward run belongs to C and No downward run in D
#Try:
#NRV(10,10,{3,5,7,9},{},{1},{4});
NRV:=proc(m,n,A,B,C,D) local gu,k:
option remember:

if (m<0 or n<0 or m<n) then
 RETURN(0):
fi:

if m=0 and n=0 then
 RETURN(1):
fi:

if member(m-n,B) then
  RETURN(0):
fi:

gu:=0:

for k from 1 to n do
 if not member(k,D) then
   gu:=gu+NRH(m,n-k,A,B,C,D):
 fi:
od:

gu:

end:



#NRH(m,n,A,B,C,D): The number of walks from the origin to [m,n] ending with a horizontal step
#where no H to  V corner [i,j] can be such that i-j belongs to A and
#no V to H corner [i,j] can be such that  i-j belongs to B. 
#No upward run belongs to C and No downward run in D. 
NRH:=proc(m,n,A,B,C,D) local gu,k:
option remember:

if (m<0 or n<0 or m<n) then
 RETURN(0):
fi:

if m=0 and n=0 then
 RETURN(0):
fi:


if member(m-n,A) then
  RETURN(0):
fi:

gu:=0:

for k from 1 to m-n do
 if not member(k, C) then
   gu:=gu+NRV(m-k,n,A,B,C,D):
  fi:
od:

gu:

end:


#SeqABCD(A,B,C,D,K): Given four finite sets of non-negative integers A, B, C,D, and a positive integer K outputs
#the first K terms (starting at n=1) of the sequence of Dyck paths  (but not listing the zeros)
#where no H to  V corner [i,j] can be such that i-j belongs to A and
#no V to H corner [i,j] can be such that  i-j belongs to  B. Try:
#No upward run belongs to C and No downward run in D. Try:
#SeqABCD({1},{1},{1},{1},20);
SeqABCD:=proc(A,B,C,D,K) local i,L:
L:=[seq(NRV(i,i,A,B,C,D),i=0..K)]:

end:


#SeqABCDr(A,B,C,D,r,K): Given two finite sets of linear expressions A, B in the variable r, and a positive integer K outputs
#the first K terms (starting at n=1) of the sequence of Dyck paths (but not listing the zeros)
#where no H to  V corner [i,j] can be such that i-j belongs to the range of A and
#no V to H corner [i,j] can be such that  i-j belongs to  the range of B. 
#
#The output is a pair [Sequene,start] where the start is the first index where it is non-zero
#All the terms  with n<start are zero.
#Try:
#SeqABCDr({2*r+3},{},r,20);
SeqABCDr:=proc(A,B,C,D,r,K) local i,L:

L:=[seq(NRVr(i,i,A,B,C,D,r),i=0..K)]:

end:

#CheckOpe(L,ope,n,N):plugs in the sequence to the operator ope, in terms of n and N
CheckOpe:=proc(L,ope,n,N) local n1,i1:
{seq(add(subs(n=n1,coeff(ope,N,i1))*L[n1+i1],i1=0..degree(ope,N)),n1=1..nops(L)-degree(ope,N))}:
end:

###End From Dyck.txt


#SeqFromRec(ope,n,N,Ini,K): The first K terms of the sequence given by the recurrence ope. Try:
#SeqFromRec(N^2-N-1,n,N,[1,1],20);
SeqFromRec:=proc(ope,n,N,Ini,K) local L,ope1,n1,i1,gu,kha:

L:=degree(ope,N):


ope1:=expand(subs(n=n-L,ope)/N^L):
ope1:=ope1/coeff(ope1,N,0):
ope1:=1-ope1:


gu:=Ini:


for n1 from nops(Ini)+1 to K do
 kha:=add(subs(n=n1,coeff(ope1,N,-i1))*gu[n1-i1] ,i1=1..L  ):
 gu:=[op(gu),kha]:
od:

gu:

end:


#InfoAB(A,B,x,P,C,a,n,M): inputs two sets of positive integers A and B, symbols x and P, a, and n, and a positive  integer
#M outputs information about the integer sequence, let's call it a(n)
#that enumerates Dyck paths of semi-length n none of whose peaks have heights in the set A
#and none of of whose valleys have heights in the set B. 
#(i) The first entry is the algebraic equation satisfied by the ordinary generating function P(x)=Sum(a(n)*x^n,n=0..infinity)
#(ii) The second entry is the expression in terms of C=(1-sqrt(1-4x))/(2*x), the Catalan generating function
#(iii) the third entry is the recurrence and initial conditions satisfied by a(n)
#(iv):The fourth entry are the first 50 terms from n=1 to n=50 (for the OEIS)
#(v) the fifth entry is the  value of a(M). Try:
#InfoAB({1},{1},x,P,C,a,n,1000);

InfoAB:=proc(A,B,x,P,C,a,n,M) local EQ1,EQ2,Ini,L1,L2,ope,N,i1:

EQ1:=fAB(A,B,x,P):

EQ2:=fABcat(A,B,x,C):


ope:=algtorec(EQ1,P,x,n,N):

Ini:=SeqABCD(A,B,{},{},degree(ope,N)):

Ini:=[op(2..nops(Ini),Ini)]:

L1:=SeqFromRec(ope,n,N,Ini,20+nops(Ini)):
L2:=SeqABCD(A,B,{},{},20+nops(Ini)):
L2:=[op(2..nops(L2),L2)]:
if L1<>L2 then
 print(`Something terrible happened`):
print(L1):
print(L2):
 RETURN(FAIL):
fi:


[EQ1,EQ2,{OperToRecEq(ope,n,N,a), seq(a(i1)=Ini[i1],i1=1..nops(Ini))}, SeqFromRec(ope,n,N,Ini,50), SeqFromRec(ope,n,N,Ini,M)[M]]:

end:




#InfoABr(A,B,r,x,P,a,n,M): inputs two sets of affine linear expressions A and B in the variable r, 
#symbols x and P, a, and n, and a positive  integer
#M outputs information about the integer sequence, let's call it a(n)
#that enumerates Dyck paths of semi-length n none of whose peaks have heights in the range of the set A
#and none of of whose valleys have heights in the set range of the set B. 
#(i) The first entry is the algebraic equation satisfied by the ordinary generating function P(x)=Sum(a(n)*x^n,n=0..infinity)
#(ii) the second entry is the recurrence and initial conditions satisfied by a(n)
#(iii):The third entry are the first 50 terms from n=1 to n=50 (for the OEIS)
#(v) the fourth entry is the  value of a(M). Try:
#InfoABr({2*r+1},{},r,x,P,a,n,1000);

InfoABr:=proc(A,B,r,x,P,a,n,M) local EQ1,Ini,L1,L2,ope,N,i1:

EQ1:=fABr(A,B,r,x,P):

ope:=algtorec(EQ1,P,x,n,N):

Ini:=SeqABCDr(A,B,{},{},r,degree(ope,N)):

Ini:=[op(2..nops(Ini),Ini)]:

L1:=SeqFromRec(ope,n,N,Ini,20+nops(Ini)):
L2:=SeqABCDr(A,B,{},{},r,20+nops(Ini)):
L2:=[op(2..nops(L2),L2)]:
if L1<>L2 then
 print(`Something terrible happened`):
print(L1):
print(L2):
 RETURN(FAIL):
fi:

if ope<>FAIL then
[EQ1,{OperToRecEq(ope,n,N,a), seq(a(i1)=Ini[i1],i1=1..nops(Ini))}, SeqFromRec(ope,n,N,Ini,50), SeqFromRec(ope,n,N,Ini,M)[M]]:
else
[EQ1,FAIL, SeqFromRec(ope,n,N,Ini,50), SeqFromRec(ope,n,N,Ini,M)[M]]:
fi:

end:


#ChildCD(STATE,Z,x): Inputs a state STATE=[C,C1,D,D1,IRI]: outputs
#the equation satisfied by Z[STATE] and the set of children states that feature in the
#formula. Try:
#ChildCD([{1},{1},{1},{1},0],Z,x);
ChildCD:=proc(STATE,Z,x)
local C,C1,D,D1,eq,C2,D2,c,d,IR:

if nops(STATE)<>5 then
 RETURN(FAIL):
fi:

C:=STATE[1]: C1:=STATE[2]: D:=STATE[3]:D1:=STATE[4]: IR:=STATE[5]:

if IR=0 then
 eq:=Z[C,C1,D,D1,0]-1-Z[C,C1,D,D1,1]-Z[C,C1,D,D,1]*Z[C,C,D,D,0]*Z[C,C,D,D1,1]:
 RETURN([eq,{[C,C1,D,D1,1],  [C,C1,D,D,1],[C,C,D,D,0],[C,C,D,D1,1]}]):
fi:

if IR=1 then

if not member(1,C1) and not member(1,D1) then
 C2:={seq(c-1,c in C1)}:
 D2:={seq(d-1,d in D1)}:
eq:=Z[C,C1,D,D1,1]-x*Z[C,C2,D,D2,0]:
 RETURN([eq,{[C,C2,D,D2,0]}]):
fi:

if  member(1,C1) and not member(1,D1) then
eq:=Z[C,C1,D,D1,1]-Z[C,C1 minus {1},D,D1,1] +x:
 RETURN([eq,{[C,C1 minus {1},D,D1,1]} ]):
fi:


if  not member(1,C1) and member(1,D1) then
eq:=Z[C,C1,D,D1,1]-Z[C,C1 ,D,D1 minus {1},1] +x:
 RETURN([eq,{[C,C1 ,D,D1 minus {1},1]}]):
fi:

if   member(1,C1) and member(1,D1) then
eq:=Z[C,C1,D,D1,1]-Z[C,C1 minus {1} ,D,D1 minus {1},1] +x:
 RETURN([eq,{[C,C1 minus {1} ,D,D1 minus {1},1]}]):

fi:

fi:



end:








#fCD(C,D,x,P): The algebraic equation satisfied by the ordinary generating function of the
#sequence enumerating Dyck paths with no upward runs in C and no downward runs in D. Try:
#fCD({1},{1},x,P);
fCD:=proc(C,D,x,P) local eq,gu,STILLTODO, ALREADYDONE,Z,a,STATE,var,var1,lu,gu1,i,i1:

ALREADYDONE:={}:

STILLTODO:={[C,C,D,D,0]}:

eq:={}:

while STILLTODO<>{} do
 STATE:=STILLTODO[1]:

ALREADYDONE:=ALREADYDONE union {STATE}:

 gu:=ChildCD(STATE,Z,x):
 eq:=eq union {gu[1]}:

 STILLTODO:=(STILLTODO union gu[2]) minus ALREADYDONE:
od:


var:={seq(Z[op(a)], a in ALREADYDONE)}:

var1:=var minus {Z[C,C,D,D,0]}:

var:=[op(var1),Z[C,C,D,D,0]]:


gu:=subs(Z[C,C,D,D,0]=P,Groebner[Basis](eq,plex(op(var)))[1]);

gu:=factor(gu):

if not type(gu,`*`) then
 RETURN(gu):
fi:

lu:=SeqABCD({},{},C,D,30):
lu:=add(lu[i]*x^(i-1),i=1..nops(lu)):

for i from 1 to nops(gu) do
gu1:=op(i,gu):
if {seq(coeff(taylor(subs(P=lu,gu1),x=0,31),x,i1),i1=0..30)}={0} then
gu1:=add(factor(coeff(gu1,P,i))*P^i,i=0..degree(gu1,P)):
 RETURN(gu1):
fi:
od:

FAIL:
end:


#CheckfAB(A,B,K): Checks procedure fAB(A,B,P,x) up to K terms. Try:
#CheckfAB({1},{1},50);
CheckfAB:=proc(A,B,K) local P,x,eq,lu,mu,i:
eq:=fAB(A,B,x,P):

lu:=SeqABCD(A,B,{},{},K):

mu:=add(lu[i]*x^(i-1),i=1..nops(lu)):

eq:=taylor(subs(P=mu,eq),x=0,K+1):

evalb({seq(coeff(eq,x,i),i=0..K)}={0}):

end:


#CheckfABr(A,B,K): Checks procedure fAB(A,B,P,x) up to K terms. Try:
#CheckfABr({2*r+1},{2*r+3},50);
CheckfABr:=proc(A,B,r,K) local P,x,eq,lu,mu,i:
eq:=fABr(A,B,r,x,P):

lu:=SeqABCDr(A,B,{},{},r,K):

mu:=add(lu[i]*x^(i-1),i=1..nops(lu)):

eq:=taylor(subs(P=mu,eq),x=0,K+1):

evalb({seq(coeff(eq,x,i),i=0..K)}={0}):

end:




#CheckfCD(C,D,K): Checks procedure fCD(C,D,P,x) up to K terms. Try:
#CheckfCD({1},{1},50);
CheckfCD:=proc(C,D,K) local P,x,eq,lu,mu,i:
eq:=fCD(C,D,x,P):

lu:=SeqABCD({},{},C,D,K):

mu:=add(lu[i]*x^(i-1),i=1..nops(lu)):

eq:=taylor(subs(P=mu,eq),x=0,K+1):

evalb({seq(coeff(eq,x,i),i=0..K)}={0}):

end:

#ChildCDr(STATE,r,Z,x): Inputs a state STATE=[C,C1,D,D1,IR], where C, C1, D, D1 are sets of arithmetical progresions.Outputs
#the equation satisfied by Z[STATE] and the set of children states that feature in the
#formula. Try:
#ChildCDr([{2*r+1},{2*r},{4*r+4},{4*r+3}],Z,x);
ChildCDr:=proc(STATE,r,Z,x)
local C,C1,D,D1,eq,C2,D2,c,d,IR,c1,d1:

if nops(STATE)<>5 then
 RETURN(FAIL):
fi:

C:=STATE[1]: C1:=STATE[2]: D:=STATE[3]:D1:=STATE[4]: IR:=STATE[5]:

if IR=0 then
 eq:=Z[C,C1,D,D1,0]-1-Z[C,C1,D,D1,1]-Z[C,C1,D,D,1]*Z[C,C,D,D,0]*Z[C,C,D,D1,1]:
 RETURN([eq,{[C,C1,D,D1,1],  [C,C1,D,D,1],[C,C,D,D,0],[C,C,D,D1,1]}]):
fi:

if IR=1 then

if not member(1, {seq(coeff(c1,r,0),c1 in C1)}) and not  member(1, {seq(coeff(d1,r,0),d1 in D1)})  then

 C2:={seq(c-1,c in C1)}:
 D2:={seq(d-1,d in D1)}:

eq:=Z[C,C1,D,D1,1]-x*Z[C,C2,D,D2,0]:
 RETURN([eq,{[C,C2,D,D2,0]}]):
fi:



if  member(1, {seq(coeff(c1,r,0),c1 in C1)}) and not  member(1, {seq(coeff(d1,r,0),d1 in D1)})  then
C2:={}:

 for c1 in C1 do
  if coeff(c1,r,0)=1 then
   C2:=C2 union {subs(r=r+1,c1)}:
   else
    C2:=C2 union {c1}:
 fi:

 od:

eq:=Z[C,C1,D,D1,1]-Z[C,C2,D,D1,1] +x:
 RETURN([eq,{[C,C2,D,D1,1]} ]):
fi:





if  member(1, {seq(coeff(d1,r,0),d1 in D1)}) and not  member(1, {seq(coeff(c1,r,0),c1 in C1)})  then
D2:={}:

 for d1 in D1 do
  if coeff(d1,r,0)=1 then
   D2:=D2 union {subs(r=r+1,d1)}:
   else
    D2:=D2 union {d1}:
 fi:

 od:

eq:=Z[C,C1,D,D1,1]-Z[C,C1,D,D2,1] +x:
 RETURN([eq,{[C,C1,D,D2,1]} ]):
fi:


if  member(1, {seq(coeff(d1,r,0),d1 in D1)}) and  member(1, {seq(coeff(c1,r,0),c1 in C1)})  then


C2:={}:

 for c1 in C1 do
  if coeff(c1,r,0)=1 then
   C2:=C2 union {subs(r=r+1,c1)}:
   else
    C2:=C2 union {c1}:
 fi:

 od:


D2:={}:

 for d1 in D1 do
  if coeff(d1,r,0)=1 then
   D2:=D2 union {subs(r=r+1,d1)}:
   else
    D2:=D2 union {d1}:
 fi:

 od:


eq:=Z[C,C1,D,D1,1]-Z[C,C2,D,D2,1] +x:
 RETURN([eq,{[C,C2,D,D2,1]} ]):
fi:


fi:



end:








#fCDr(C,D,r,x,P): The algebraic equation satisfied by the ordinary generating function of the
#sequence enumerating Dyck paths with no upward runs in the range of C and no downward runs in the range of D. Try:
#fCDr({2*r+1},{2*r+1},r,x,P);
fCDr:=proc(C,D,r,x,P) local eq,gu,STILLTODO, ALREADYDONE,Z,a,STATE,var,var1,lu,gu1,i,i1:

ALREADYDONE:={}:

STILLTODO:={[C,C,D,D,0]}:

eq:={}:

while STILLTODO<>{} do
 STATE:=STILLTODO[1]:

ALREADYDONE:=ALREADYDONE union {STATE}:

 gu:=ChildCDr(STATE,r,Z,x):
 eq:=eq union {gu[1]}:

 STILLTODO:=(STILLTODO union gu[2]) minus ALREADYDONE:
od:


var:={seq(Z[op(a)], a in ALREADYDONE)}:

var1:=var minus {Z[C,C,D,D,0]}:

var:=[op(var1),Z[C,C,D,D,0]]:


gu:=subs(Z[C,C,D,D,0]=P,Groebner[Basis](eq,plex(op(var)))[1]);

gu:=factor(gu):

if not type(gu,`*`) then
 RETURN(gu):
fi:

lu:=SeqABCDr({},{},C,D,r,30):
lu:=add(lu[i]*x^(i-1),i=1..nops(lu)):

for i from 1 to nops(gu) do
gu1:=op(i,gu):
if {seq(coeff(taylor(subs(P=lu,gu1),x=0,31),x,i1),i1=0..30)}={0} then
gu1:=add(factor(coeff(gu1,P,i))*P^i,i=0..degree(gu1,P)):
 RETURN(gu1):
fi:
od:

FAIL:
end:




#CheckfCDr(C,D,r,K): Checks procedure fCDr(C,D,r,P,x) up to K terms. Try:
#CheckfCDr({2*r+2},{},50);
CheckfCDr:=proc(C,D,r,K) local P,x,eq,lu,mu,i:
eq:=fCDr(C,D,r,x,P):

lu:=SeqABCDr({},{},C,D,r,K):

mu:=add(lu[i]*x^(i-1),i=1..nops(lu)):

eq:=taylor(subs(P=mu,eq),x=0,K+1):

evalb({seq(coeff(eq,x,i),i=0..K)}={0}):

end: