######################################################################
## Dyck.txt Save this file as  Dyck.txt  to use it,                  #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read Dyck.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 March/April 2020: tested for Maple 18 `):
print(`Version  of April 2020  `):
print():
print(`This is Dyck.txt, A Maple package`):
print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's  article: `):
print(` Counting  Restricted Dyck Paths via (Numeric and Symbolic) Dynamic Programming `):
print(`Available from:`):
print(` http://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/kiss.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():


 

ezraBij:=proc()
if args=NULL then

print(` The bijections procedures are : AJ, invAJ, CheckAJ, dyck_to_motzkin, motzkin_to_dyck `):

else
ezra(args):

fi:

end:

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 a list of the procedures from FindRec.txt`):
print(` type "ezraFindRec();". For specific help type "ezraFindRec(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():
 




ezraStory:=proc()
if args=NULL then

print(` The Story procedures are :  Paper, PaperR, Theorem, TheoremR `):

else
ezra(args):
fi:

end:

ezraCheck:=proc()
if args=NULL then

print(` The Checking procedures are :  DyckPaths, DyckPathsABCD,  DyckPathsABCDr,`):
print(`IsGoodP, IsGoodN, IsGoodCP,IsGoodCN, SeqABCDbrute, SeqABCDrBrute, Targem1, Targem2 `):

else
ezra(args):
fi:

end:

ezra1:=proc()
if args=NULL then

print(` The supporting procedures are :  DrawDP, H, IsRange, IsRange1, NH, NV, RH, NRH, NRHr, NRV,NRVr,  OpeToRec, RV, V `):

else
ezra(args):

fi:


end:


ezra:=proc()
if args=NULL then

print(` Dyck.txt:  A Maple package for automatic enumeration of families of Dyck paths obeying certain kinds of restrictions  `):

print(`The MAIN procedures are`):
print(` Aeq,   InfoABCD, InfoABCDr, SeqABCD, SeqABCDr `):

elif nargs=1 and args[1]=DrawDP then
print(`DrawDP(L): plots the path L `):

elif nargs=1 and args[1]=Aeq then
print(`Aeq(L,x,P): The guessed algebraic equation in P(x) satisfied by the generating function of the sequence whose beginning`):
print(`starting at n=0 is L. Try:`):
print(`Aeq([seq(binomial(2*i,i)/(i+1),i=1..10)],x,P);`):


elif nargs=1 and args[1]=AJ then
print(`AJ(W): inputs a walk W of semi-length 2*n+1 that does not have peak-heights and valley-heights that are positive even integers`):
print(`outputs Alison J. Bu's mapping it to a word in {1,-1} with n 1's and n -1's. try:`):
print(`AJ([1,-1,1,-1]);`):


elif nargs=1 and args[1]=CheckAJ then
print(`CheckAJ(n): Checks Alison J. Bu's bijection for n. Try:`):
print(`CheckAJ(3);`):




elif nargs=1 and args[1]=DyckPaths then
print(`DyckPaths(n): The set of Dyck path of length 2n.: Try:`):
print(`DyckPathsP(5);`):

elif nargs=1 and args[1]=DrawDPs then
print(`DrawDPs(S): plots the set of paths in S . Try:`):
print(`DrawDPs(DyckPaths(3));`):


elif nargs=1 and args[1]=DyckPathsABCD then
print(`DyckPathsABCD(A,B,C,D,n): Given four finite sets of non-negative integers A, B, C,D, and a positive integer n, outputs`):
print(`the set of Dyck paths counted by SeqABCD(A,B,C,D,K). In other words:`):
print(`where no peak can be of height that belongs to A and`):
print(`no valley can be of height that belongs to  B:`):
print(`No length of upward run belongs to C and `):
print(`no length of downward run belongs to D. Try:`):
print(`DyckPathsABCD({1},{1},{2},{1},10);`):


elif nargs=1 and args[1]=DyckPathsABCDr then
print(`DyckPathsABCDr(A,B,C,D,r,n): Given four finite sets of arithmetical progressopms A, B, C,D, and a positive integer n, outputs`):
print(`the set of Dyck paths counted by SeqABCDr(A,B,C,D,r,K). In other words:`):
print(`where no peak can be of height that belongs to A and`):
print(`no valley can be of height that belongs to  B:`):
print(`No length of upward run belongs to C and `):
print(`no length of downward run belongs to D. Try:`):
print(`DyckPathsABCDr({},{},{2*r+1},{},6);`):


elif nargs=1 and args[1]=dyck_to_motzkin then
print(`dyck_to_motzkin(path): Robert Rougherty-Bliss's mapping from Dyck paths that avoid peak-heights of the form 2*r+3`):
print(`to Motzkin paths, giving a bijective proof of V. Retakh's result. Try:`):
print(`S:=DyckPathsABCDr({2*r+3},{},{},{},r,5): {seq(dyck_to_motzkin(s), s in S)};`):

elif nargs=1 and args[1]=H then
print(`H(m,n): The set of walks from the origin to [m,n] ending with a horizontal step. Try`):
print(` H(3,2); `):



elif nargs=1 and args[1]=InfoABCD  then
print(`InfoABCD(A,B,C,D,K,P,x,n,a,MaxC,M): inputs sets of positive integers A and B, and outputs the following list`):
print(`it regards the enumerating sequence of Dyck paths where`):
print(`where no H to  V corner [i,j] can be such that i-j belongs to A and`):
print(`no V to H corner [i,j] can be such that  i-j belongs to  B. `):
print(`No upward run belongs to C and no downward run belongs to D.`):
print(`(i)The algebraic equation F(P(x),x)=0 satisfied  by the ordinary generating function of the sequence`):
print(`(ii) The recurrence equation phrased in terms of a(n) and the initial conditions.`):
print(`(iii) Just for fun, the M-th term of the original sequence .`):
print(`(iv): For the OEIS the first K terms of the sequence, starting at n=0. :`):
print(`Note that if K is nos sufficiently large, only part of this will be outputted. Try:`):
print(`InfoABCD({1},{1},{1},{1},60,P,x,n,a,20,1000);`):


elif nargs=1 and args[1]=InfoABCDr  then
print(`InfoABCDr(A,B,C,D,r,K,P,x,n,a,MaxC,M): inputs sets of expressions in r, A and B, and outputs the following list`):
print(`it regards the enumerating sequence of Dyck paths where`):
print(`where no H to  V corner [i,j] can be such that i-j belongs to the range of A when r>=0 and`):
print(`no V to H corner [i,j] can be such that  i-j belongs to  the range of B when r>=0 and`):
print(`No upward run belongs to C and no downward run belongs to D. The output is`):
print(`(i)The algebraic equation F(P(x),x)=0 satisfied  by the ordinary generating function of the sequence`):
print(`(ii) The recurrence equation phrased in terms of a(n) and the initial conditions.`):
print(`(iii) Just for fun, the M-th term of the original sequence`):
print(`(iv) For the OEIS, the first K-a terms of the sequence, let's call this shifted sequence L. Try:`):
print(`InfoABCDr({1+2*r},{1+2*r},{},{},r, 60,P,x,n,a,20,1000);`):



elif nargs=1 and args[1]=invAJ then
print(`invAJ(W): inputs a walk W of semi-length n from [0,0] to [2*n,0] for some n`):
print(`outputs Alison J. Bu's mapping it to a Dyck path. try:`):
print(`invAJ([-1,1,-1,1]);`):

elif nargs=1 and args[1]=IsGoodCN then
print(`IsGoodCN(L,A): Inputs a Dyck path, and a set A of positive integers, decides whether`):
print(` none of the heights in valley locations belong to A. `):
print(` Try: `):
print(` L:=DyckPathsP(8)[13]; IsGoodCN(L,{1}); `):


elif nargs=1 and args[1]=IsGoodCP then
print(`IsGoodCP(L,A): Inputs a Dyck path, and a set A of positive integers, decides whether`):
print(` none of the heights in peak locations belong to A. `):
print(` Try: `):
print(` L:=DyckPaths(8)[13]; IsGoodCP(L,{1}); `):

elif nargs=1 and args[1]=IsGoodN then
print(`IsGoodN(L,A): Inputs a Dyck path, and a set A of positive integers, decides whether`):
print(` none of the downwards runs belong to A. `):
print(` Try: `):
print(` L:=DyckPathsP(8)[13]; IsGoodN(L,{1}); `):


elif nargs=1 and args[1]=IsGoodP then
print(`IsGoodP(L,A): Inputs a Dyck path, and a set A of positive integers, decides whether`):
print(` none of the upwards runs belong to A. `):
print(` Try: `):
print(` L:=DyckPaths(8)[13]; IsGoodP(L,{1}); `):



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

elif nargs=1 and args[1]=IsRange1 then
print(`IsRange1(f,k,n): inputs a linear expression f in k and a non-negative integer n`):
print(`decides whether f(k0)=n for some non-negative integer k0.`):
print(`Try:`):
print(`IsRange1(3*k+1,k,7);`):



elif nargs=1 and args[1]=motzkin_to_dyck then
print(`motzkin_to_dyck(path): Robert Rougherty-Bliss's mapping from Motzkin paths to Dyck paths that avoid peak-heights of the form 2*r+3`):
print(`giving the invese bijection of dyck_to_motzkin. Try:`):
print(`S:=DyckPathsABCDr({2*r+3},{},{},{},r,5): {seq(motzkin_to_dyck(dyck_to_motzkin(s)), s in S)};`):

elif nargs=1 and args[1]=NH then
print(`NH(m,n): The number of walks from the origin to [m,n] ending with a horizontal step. Try`):
print(` NH(3,2); `):


elif nargs=1 and args[1]=NRH then
print(`NRH(m,n,A,B): The number of walks from the origin to [m,n] ending with a horizontal step`):
print(`where no H to  V corner [i,j] can be such that i-j belongs to A and`):
print(`no V to H corner [i,j] can be such that  i-j belongs to B. Try:`):
print(`RNH(10,8,{3,5,7,9},{});`):

elif nargs=1 and args[1]=NRV then
print(`NRV(m,n,A,B): The number of walks from the origin to [m,n] ending with a vertical step`):
print(`where no H to  V corner [i,j] can be such that i-j belongs to A and`):
print(`no V to H corner [i,j] can be such that  i-j belongs to B. Try:`):
print(`NRV(10,10,{3,5,7,9},{});`):


elif nargs=1 and args[1]=NRVr then
print(`NRVr(m,n,A,B,r):  Inputs positive integers m>=n and sets A and B of linear expressions in the variable r`):
print(`Outputs`):
print(`The number of walks from the origin to [m,n] ending with a vertical step`):
print(`where no H to  V corner [i,j] can be such that i-j belongs to the range of A and`):
print(`no V to H corner [i,j] can be such that  i-j belongs to the range of B. Try:`):
print(`NRVr(10,10,{2*r+3},{},r);`):


elif nargs=1 and args[1]=NV then
print(`NV(m,n): The number of walks from the origin to [m,n] ending with a vertical step. Try`):
print(` NV(3,2); `):

elif nargs=1 and args[1]=OpeToRec then
print(`OpeToRec(ope,n,N,INI,a): inputs a recurrence operators ope in n and N and initial coditions, write it nicely`):
print(` in terms of a(n)=Combination of a(n-i), and the initial conditions. Try: `):
print(`OpeToRec(N^2-N-1,n,N,[1,1],a);`):

elif nargs=1 and args[1]=Paper then
print(`Paper(a0,K,P,x,n,a,MaxC,M): inputs a positve integer a0, and K,P,x,n,a,MaxC,M the same as in`):
print(`Theorem(A,B,C,D,K,P,x,n,a,MaxC, M) (q.v.), outputs an article with 2^(4*A) theorems outputted`):
print(`by Theorem(A,B,C,D,K,P,x,n,a,MaxC, M) where A,B,C,D, range over the subsets of {1, ..,a0}.`):
print(`Some of the theorems may be trivially equivalent due to redudancy of the conditions  and/or symmetry.Try:`):
print(`Paper(1,100,P,x,n,a,16,1000);`):


elif nargs=1 and args[1]=PaperR then
print(`PaperR(a0,b0,K,P,x,n,a,MaxC,M): inputs a positve integer A, and K,P,x,n,a,MaxC,M the same as in`):
print(`Theorem(A,B,C,D,K,P,x,n,a,MaxC, M) (look it up), outputs an article with 16 theorems outputted`):
print(`by TheoremR(A,B,C,D,r,K,P,x,n,a,MaxC, M) where A,B,C,D, range {{},{a0*r+b0}}`):
print(`Some of the theorems may be trivially equivalent due to redudancy of the conditions and/or symmetry.Try:`):
print(`PaperR(2,3,100,P,x,n,a,20,1000);`):

elif nargs=1 and args[1]=SeqABCD then
print(`SeqABCD(A,B,C,D,K): Given four finite sets of non-negative integers A, B,C,D and a positive integer K outputs`):
print(`the first K+1 terms (starting at n=0) of the sequence enumerating Dyck paths that `):
print(`where no peak can be of height that belongs to A and`):
print(`no valley can be of height that belongs to  B.`):
print(`No length of upward run belongs to C and `):
print(`no length of downward run belongs to D. Try: `):
print(`SeqABCD({1},{1},{1},{2},20);`):

elif nargs=1 and args[1]=SeqABCDbrute then
print(`SeqABCDbrute(A,B,C,D,K):  Same output as SeqABCD(A,B,C,D,K) but by complete brute force. FOR CHECKING PURPOSES ONLY`):
print(`Don't  make K too big! Try:`):
print(`SeqABCDbrute({1},{1},{1},{1},11);`):

elif nargs=1 and args[1]=SeqABCDr then
print(`SeqABCDr(A,B,C,D,r,K): Given four finite sets of linear expressions A, B,C,D in the variable r, and a positive integer K outputs`):
print(`the first K+1 terms (starting at n=0) of the sequence of Dyck paths `):
print(`where no peak can be of height that belongs to the range of A `):
print(`no valley can be of height that belongs to  the range of B.`):
print(`No length of upward run belongs to the range of C `):
print(`and no length of downward run belongs to the range of D. Try: `):
print(`SeqABCDr({2*r+3},{},{3*r+2},{3*r+1},r,20);`):


elif nargs=1 and args[1]=SeqABCDrBrute then
print(`SeqABCDrBure(A,B,C,D,r,K): Like SeqABCDrA,B,C,D,r,K) but by brute force. For checking purposed only. Try:`):
print(`SeqABCDrBrute({2*r+3},{},{3*r+2},{3*r+1},r,10);`):

elif nargs=1 and args[1]=Targem1 then
print(`Targem1(L): translates a Dyck path into a sequence of alternating positive and negative integers`):
print(`such that the partial sums are always non-negative. Try`):
print(`Targem1([1,-1,1,-1]);`):

elif nargs=1 and args[1]=Targem2 then
print(`Targem2(L): translates a Dyck path into a sequence of the partial sums when it changes sign. Try:`):
print(`Targem2([1,-1,1,-1]);`):

elif nargs=1 and args[1]=Theorem then
print(`Theorem(A,B,C,D,K,P,x,n,a,MaxC, M): Verbose version of InfoABCD(A,B,C,D,K,P,x,n,a,MaxC, M) (q.v.). Try:`):
print(`Theorem({1},{1},{1},{1},60,P,x,n,a,20,1000);`):


elif nargs=1 and args[1]=TheoremR then
print(`TheoremR(A,B,C,D,r,K,P,x,n,a,MaxC, M): Verbose version of InfoABCDr(A,B,C,D,r,K,P,x,n,a,MaxC, M) (q.v.). Try:`):
print(`TheoremR({2*r+3},{},{},{},r,60,P,x,n,a,20,1000);`):


elif nargs=1 and args[1]=V then
print(`V(m,n): The set of walks from the origin to [m,n] ending with a vertical step. Try`):
print(` V(3,2); `):

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

fi:


end:

#V(m,n): The set of walks from the origin to [m,n] ending with a vertical step
V:=proc(m,n) local gu,mu,mu1,k:
option remember:

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

if m=0 and n=0 then
 RETURN({[]}):
fi:

gu:={}:

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

gu:

end:



#H(m,n): The set of walks from the origin to [m,n] ending with a horizontal step
H:=proc(m,n) local gu,mu,mu1,k:
option remember:

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

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

gu:={}:

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

gu:

end:


#NV(m,n): The number of walks from the origin to [m,n] ending with a vertical step
NV:=proc(m,n) 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:

gu:=0:

for k from 1 to n do
 gu:=gu+NH(m,n-k):
od:

gu:

end:

#NH(m,n): The number of walks from the origin to [m,n] ending with a horizontal step
NH:=proc(m,n) 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:

gu:=0:

for k from 1 to m-n do
 gu:=gu+NV(m-k,n):
od:

gu:

end:


#RV(m,n,A,B): The set 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. Try:
#RV(10,10,{3,5,7,9},{});
RV:=proc(m,n,A,B) local gu,mu,mu1,k:
option remember:

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

if m=0 and n=0 then
 RETURN({[]}):
fi:

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

gu:={}:

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

gu:

end:



#RH(m,n): The set 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. Try:

RH:=proc(m,n,A,B) local gu,mu,mu1,k:
option remember:

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

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


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

gu:={}:

for k from 1 to m-n do
 mu:=RV(m-k,n,A,B):
 gu:= gu union {seq([op(mu1),k],mu1 in mu)}:
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:


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

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


###Start from SCHUTZENBERGER.txt
#empir(gu,degx,degP,x,P) 
#to "fit" an algebraic equation F(P(x),x)=0 of degree
#degP in P(x) and degx in n for P(x):=sum_i gu[i]*x^i
 
empir:=proc(gu,degx,degP,x,P)
local i1,i2,F,a,cand,lu,eq,var,mu,flo,pip,ka,mu1,Ftry:

if (1+degx)*(1+degP) > nops(gu)-15 then
 RETURN(`sequence too small`):
fi:
 
F:=0:
var:={}:
for i1 from 0 to degx do
 for i2 from 0 to degP do
  F:=F+a[i1,i2]*x^i1*P^i2:
  var:=var union {a[i1,i2]}:
 od:
od:
 
cand:=0:
 
for i1 from 0 to nops(gu)-1 do
 cand:=cand+op(i1+1,gu)*x^i1:
od:
 
lu:=subs(P=cand,F):
lu:=taylor(lu,x=0,nops(gu)-1):
 
eq:={}:
 
for i1 from 0 to nops(gu)-2 do
eq:=eq union {coeff(lu,x,i1)=0}
od:
 
mu:=solve(eq,var):

ka:=0:

for mu1 in mu do
 if op(1,mu1)=op(2,mu1) then
 ka:=ka+1:
 fi:
od:

if ka>1 then
 RETURN(FAIL):
fi:

F:=subs(mu,F):
 
if F=0 then
 RETURN(FAIL):
fi:
 
flo:=degree(F,P):
pip:=coeff(F,P,flo):
flo:=degree(pip,x):
pip:=coeff(pip,x,flo):
F:=normal(F/pip):

Ftry:=taylor(subs(P=add(gu[i1+1]*x^i1,i1=0..nops(gu)-1),F),x=0,nops(gu)+1):

if {seq(coeff(Ftry,x,i1),i1=0..nops(gu)-2)}<>{0} then
 RETURN(FAIL):
fi:

normal(F):
end:



Aeq:=proc(gu,x,P)
local degx,degP,L,lu:
 
for L from 1 to (nops(gu)-15)/3 do
for degP from 1 to L do
for degx from 0 to min(trunc((nops(gu)-15)/(1+degP))-1,L-degP) do
 
lu:=empir(gu,degx,degP,x,P):
 
if lu<>FAIL then
RETURN(lu):
fi:
od:
od:
od:
FAIL:
end:
 
###end from SCHUTZENBERGER.txt
 
   
 
##start from Findrec
ezraFindrec:=proc()
if args=NULL then

print(` FindRec.txt: A Maple package for empirically guessing partial recurrence`):
print(`equations satisfied by Discrete Functions of ONE Variable`):
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,MaxC): 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,2);`):

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):
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(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):
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)+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
ope:=ope[1]:

if {seq(add(subs(n=n1+i1,coeff(ope,N,i1))*f[n1+i1],i1=0..degree(ope,N)),n1=1..nops(f)-degree(ope,N))}<>{0} then
 RETURN(FAIL):
fi:

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


if SeqFromRec(ope,n,N,[op(1..degree(ope,N),f)],nops(f))<>f then
 RETURN(FAIL):
else
RETURN(ope):
fi:
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 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 and degree(ope,N)<>0 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,ka,i1:
ope1:=ope:

L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ka:={seq(subs(n=i1,lcoeff(ope1,N)),i1=1..K)}:

if member(0,ka) then
 RETURN(FAIL):
fi:

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

gu:=Ini:

for n1 from nops(Ini)+1 to K do
if member(0, { seq( subs(n=n1,denom(coeff(ope1,N,-j1)  ))  ,j1=1..L )}) then
 RETURN(FAIL):
else
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
fi:
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):
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)+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 from Findrec



#OpeToRec(ope,n,N,INI,a): inputs a recurrence operators ope in n and N and initial coditions, write it nicel
# in terms of a(n)=Combination of a(n-i), and the initial conditions. Try:
#OpeToRec(N^2-N-1,n,N,[1,1],a);
OpeToRec:=proc(ope,n,N,INI,a) local ORD,i,ope1:

ORD:=degree(ope,N):

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

ope1:=1-expand(ope1/coeff(ope1,N,0)):

[a(n)=add(factor(coeff(ope1,N,-i))*a(n-i),i=1..ORD), [seq(a(i)=INI[i],i=1..ORD)]]:

end:

#InfoABCD(A,B,C,D,K,P,x,n,a,MaxC,M): inputs sets of positive integers A and B, and outputs the following list
#it regards the enumerating sequence of Dyck paths where
#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 belongs to D
#(i): The first index n where it is non-zero, let's call it a
#(ii) The first K-a terms of the sequence, let's call this shifted sequence L
#(iii)The algebraic equation F(P(x),x)=0 satisfied  by the ordinary generating function of the sequence
#(iv) The pair [ope,INI], where ope is the linear recurrence equation op(n,N) annihilating the sequence L,  and INI are the needed
#initial conditions
#(v) Just for fun, the M-th term of the original sequence
#InfoABCD({1},{1},{1},{1},40,P,x,n,a,20,1000);
InfoABCD:=proc(A,B,C,D,K,P,x,n,a,MaxC, M) local L, L1, gu,eq,ope,INI,N:
L:=SeqABCD(A,B,C,D,K):
gu:=[]:

eq:=Aeq(L,x,P):

if eq=FAIL then
 RETURN([L]):
else
eq:=collect(eq,P):
eq:=subs(P=P(x),eq):
gu:=[op(gu),eq=0]:
fi:

L1:=[op(2..nops(L),L)]:

ope:=Findrec(L1,n,N,MaxC):

if ope=FAIL then
 RETURN([op(gu),L]):
else
INI:=[op(1..degree(ope,N),L1)]:
gu:=[op(gu),OpeToRec(ope,n,N,INI,a)]:
fi:

gu:=[op(gu),SeqFromRec(ope,n,N,INI,M)[M]]:

[op(gu),L]:
end:




#InfoABCDr(A,B,C,D,r,K,P,x,n,a,MaxC,M): inputs sets of expressions in r, A and B, and outputs the following list
#it regards the enumerating sequence of Dyck paths where
#where no H to  V corner [i,j] can be such that i-j belongs to the range of A when r>=0 and
#no V to H corner [i,j] can be such that  i-j belongs to  the range of B when r>=0 and
#no upward run belongs to C and no downward run belongs to D
#(i): The first index n where it is non-zero, let's call it a
#(ii) The first K-a terms of the sequence, let's call this shifted sequence L
#(iii)The algebraic equation F(P(x),x)=0 satisfied  by the ordinary generating function of the sequence
#(iv) The pair [ope,INI], where ope is the linear recurrence equation op(n,N) annihilating the sequence L,  and INI are the needed
#initial conditions
#(v) Just for fun, the M-th term of the original sequence
#InfoABCDr({1+2*r},{1+2*r},{},{},r, 40,P,x,n,N,1000);
InfoABCDr:=proc(A,B,C,D,r,K,P,x,n,a,MaxC, M) local N,L, L1, gu,eq,ope,INI:
L:=SeqABCDr(A,B,C,D,r,K):

gu:=[]:

eq:=Aeq(L,x,P):

if eq=FAIL then
 RETURN([L]):
else
eq:=subs(P=P(x),eq):
gu:=[op(gu),eq=0]:
fi:

L1:=[op(2..nops(L),L)]:

ope:=Findrec(L1,n,N,MaxC):

if ope=FAIL then
 RETURN([op(gu),L]):
else
INI:=[op(1..degree(ope,N),L1)]:
gu:=[op(gu),OpeToRec(ope,n,N,INI,a)]:
fi:

gu:=[op(gu),SeqFromRec(ope,n,N,INI,M)[M],L]:

gu:
end:





####Start Brute force


#DyckPaths(n): all the Dyck paths of length 2n. Try:
#DyckPaths(5);
DyckPaths:=proc(n) local k:
option remember:

if n<0  then {}:

elif n=0 then  {[]}:

else

{seq(op(DPk(n,k)),k=1..n)}:
 fi:
end:


#DPk(n,k): All the Dyck paths of length 2n of the form 1L1(-1) with L1 in DP(k-1). Try:
#DPk(5,2);
DPk:=proc(n,k) local gu1,gu2,mu1,mu2:
option remember:

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

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

 gu1:=DyckPaths(k-1):
 gu2:=DyckPaths(n-k):

{seq(seq([1,op(mu1),-1,op(mu2)],mu1 in gu1),mu2 in gu2)}:

end:

#IsGoodP(L,A): Inputs a Dyck path, and a set A of positive integers, decides whether
# none of the upwards runs belong to A.
#Try:
#L:=DP(8)[13]; IsGoodP(L,{1});
IsGoodP:=proc(L,A) local L1,i:
L1:=Targem1(L):

for i from 1 to nops(L1)/2 do
 if member(L1[2*i-1],A) then
   RETURN(false):
 fi:
od:

true:

end:


#IsGoodPr(L,A,r): Inputs a Dyck path L, and a set A of expressions in r, decides whether
# none of the upwards runs belong to the range of A.
#Try:
#L:=DP(8)[13]; IsGoodPr(L,{2*r+1},r);
IsGoodPr:=proc(L,A,r) local L1,i:
L1:=Targem1(L):

for i from 1 to nops(L1)/2 do

if IsRange(A,r,L1[2*i-1]) then
   RETURN(false):
 fi:

od:

true:

end:



#Targem1(L): translates a Dyck path into a sequence of alternating positive and negative integers
#such that the partial sums are always non-negative
Targem1:=proc(L) local i,j:

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

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

for i from 1 to nops(L) while L[i]=1 do od:

for j from i to nops(L) while L[j]=-1 do od:

[i-1,i-j,op(Targem1([op(j..nops(L),L)]))]:
end:


#Targem2(L): translates a Dyck path into a sequence of the partial sums when it changes sign
Targem2:=proc(L) local L1,i,L2,su:
L1:=Targem1(L):
L2:=[]:

su:=0:
for i from 1 to nops(L1) do
 su:=su+L1[i]:
 L2:=[op(L2),su]:
od:
L2:
end:

#IsGoodN(L,A): Inputs a Dyck path, and a set A of positive integers, decides whether
# none of the downwards runs belong to A.
#Try:
#L:=DP(8)[13]; IsGoodN(L,{1});
IsGoodN:=proc(L,A) local L1,i:
L1:=Targem1(L):

for i from 1 to nops(L1)/2 do
 if member(-L1[2*i],A) then
   RETURN(false):
 fi:
od:

true:

end:


#IsGoodNr(L,A,r): Inputs a Dyck path, and a set A of expressions in r, decides whether
# none of the downwards runs belong to A.
#Try:
#L:=DP(8)[13]; IsGoodNr(L,{2*r+1},r);
IsGoodNr:=proc(L,A,r) local L1,i:
L1:=Targem1(L):

for i from 1 to nops(L1)/2 do
if IsRange(A,r,-L1[2*i]) then
   RETURN(false):
 fi:
od:

true:

end:


#IsGoodCP(L,A): Inputs a Dyck path, and a set A of positive integers, decides whether
# none of the heights in a peak location belongs to A
#Try:
#L:=DP(8)[13]; IsGoodCP(L,{1});
IsGoodCP:=proc(L,A) local L1,i:
L1:=Targem2(L):

for i from 1 to nops(L1)/2 do
 if member(L1[2*i-1],A) then
   RETURN(false):
 fi:
od:

true:

end:


#IsGoodCPr(L,A,r): Inputs a Dyck path, and a set A of expressions in r, decides whether
# none of the heights in a peak location belongs to the range of A
#Try:
#L:=DP(8)[13]; IsGoodCPr(L,{2*r+1});
IsGoodCPr:=proc(L,A,r) local L1,i:
L1:=Targem2(L):

for i from 1 to nops(L1)/2 do
 if IsRange(A,r,L1[2*i-1]) then
   RETURN(false):
 fi:
od:

true:

end:


#IsGoodCN(L,A): Inputs a Dyck path, and a set A of positive integers, decides whether
# none of the heights in a valley location belongs to A
#Try:
#L:=DP(8)[13]; IsGoodCP(L,{1});
IsGoodCN:=proc(L,A) local L1,i:
L1:=Targem2(L):

for i from 1 to nops(L1)/2 do
 if member(L1[2*i],A) then
   RETURN(false):
 fi:
od:

true:

end:


#IsGoodCNr(L,A,r): Inputs a Dyck path, and a set A of expressions in r, decides whether
# none of the heights in a valley location belongs to the range of A
#Try:
#L:=DP(8)[13]; IsGoodCNr(L,{2*r+1},r);
IsGoodCNr:=proc(L,A,r) local L1,i:
L1:=Targem2(L):

for i from 1 to nops(L1)/2 do
 if IsRange(A,r,L1[2*i]) then
   RETURN(false):
 fi:
od:

true:

end:





#DyckPathsABCD(A,B,C,D,n): Given four finite sets of non-negative integers A, B, C,D, and a positive integer n, outputs
#the set of Dyck paths counted by SeqABCD(A,B,C,D,K). In other words:
#where no peak can be of height that belongs to A and`):
#no valley can be of height that belongs to  B:
#No length of upward run belongs to C and 
#no length of downward run belongs to D. Try
#DyckPathsABCD({1},{1},{2},{1},10);
DyckPathsABCD:=proc(A,B,C,D,n) local gu,mu,mu1:
mu:=DyckPaths(n):
gu:={}:
for mu1 in mu do

if (IsGoodP(mu1,C) and IsGoodN(mu1,D) and IsGoodCP(mu1,A) and IsGoodCN(mu1,B)) then
 gu:=gu union {mu1}:
fi:

od:

gu:

end:


#DyckPathsABCDr(A,B,C,D,r,n): Given four sets  A, B, C,D, of expressions in the variable r, and a positive integer n, outputs
#the set of Dyck paths counted by SeqABCDr(A,B,C,D,K). In other words:
#where no peak can be of height that belongs to the range of A and`):
#no valley can be of height that belongs to  the range of B:
#No length of upward run belongs to the range of C and 
#no length of downward run belongs to the range of D. Try
#DyckPathsABCDr({1},{1},{2},{1},10);
DyckPathsABCDr:=proc(A,B,C,D,r,n) local gu,mu,mu1:
mu:=DyckPaths(n):
gu:={}:
for mu1 in mu do

if (IsGoodPr(mu1,C,r) and IsGoodNr(mu1,D,r) and IsGoodCPr(mu1,A,r) and IsGoodCNr(mu1,B,r)) then
 gu:=gu union {mu1}:
fi:

od:

gu:

end:



#SeqABCDbrute(A,B,C,D,K):  Same output as SeqABCD(A,B,C,D,K) but by complete brute force. FOR CHECKING PURPOSES ONLY
#Don't  make K too big. Try:
#SeqABCDbrute({1},{1},{1},{1},11);
SeqABCDbrute:=proc(A,B,C,D,K) local i:

[seq(nops(DyckPathsABCD(A,B,C,D,i)),i=0..K)]:

end:


#SeqABCDrBrute(A,B,C,D,K):  Same output as SeqABCDr(A,B,C,D,r,K) but by complete brute force. FOR CHECKING PURPOSES ONLY
#Don't  make K too big. Try:
#SeqABCDrBrute({2*r+1},{},{},{},r,9);
SeqABCDrBrute:=proc(A,B,C,D,r,K) local i:

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

end:


#Targem3(L): Translates the Dyck path into the sequence of lattince points above the x-axis
Targem3:=proc(L) local i,L1,x:
L1:=[[0,0]]:
x:=0:

for i from 1 to nops(L) do
 x:=x+L[i]:
 L1:=[op(L1),[i,x]]:
od:

L1:

end:

#DrawDP(L): plots the path D
DrawDP:=proc(L): plot(Targem3(L)):end:



#DrawDPs(S): plots the set of paths in S . Try:
#DrawDPs(DyckPaths(3));
DrawDPs:=proc(S) local p,i,st,L,L1,j,sp: 

st:=0:

L:=Targem3(S[1]):

sp:=L[-1][1]:

p:=plot(L,axes=none,scaling=constrained):


for i from 2 to nops(S) do
L:=Targem3(S[i]):
st:=st+L[-1][1]+sp:
L1:=[seq([L[j][1]+st,L[j][2]],j=1..nops(L))]:
 p:=p,plot(L1):
od:

display(p):

end:

####end Brute force

##Start Stories


#Theorem(A,B,C,D,K,P,x,n,a,MaxC, M): Verbose version of InfoABCD(A,B,C,D,K,P,x,n,a,MaxC, M) (q.v.). Try:
#Theorem({1},{1},{1},{1},40,P,x,n,a,20,1000);
Theorem:=proc(A,B,C,D,K,P,x,n,a,MaxC, M) local gu,S:

gu:=InfoABCD(A,B,C,D,K,P,x,n,a,MaxC, M):


print(``):
print(`------------------------------------------------------------`):
print(``):
print(`On The Number of Dyck Paths Obeying Certain Restrictions`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
if A={} and B={} and C={} and D={} then
print(`Theorem: Let a(n) be the number of Dyck paths of semi-length n`):
else
print(`Theorem: Let a(n) be the number of Dyck paths of semi-length n obeying the following restrictions`):
if A<>{} then
 print(`The height of a peak can not belong to `, A):
fi:
if B<>{} then
 print(`The height of a valley can not belong to `, B):
fi:
if C<>{} then
 print(`No upward-run can belong to`, C):
fi:
if D<>{} then
 print(`No downward-run can belong to`, D):
fi:

fi:




S:=DyckPathsABCD(A,B,C,D,8):

if nops(S)<>{} then
print(`To make it crystal clear here is such a path of semi-length 8`):
 S:=S[rand(1..nops(S))()]:
print(``):
print(DrawDP(S));
fi:

print(``):

if nops(gu)=1 then
 print(`The first`, K, `terms of the sequence are `):
 print(gu[1]):
 RETURN():
fi:


if nops(gu)=2 then

print(`The ordinary generating function of the sequence a(n), let's call it `, P(x), ` satisfies the algebraic equation`):
print(``):
print(gu[1]):
print(``):
print(`For the sake of the OEIS, here are the first`, K, `terms `):
print(``):
print(gu[2]):
print(``):
RETURN():

fi:


print(``):
print(`The ordinary generating function of the sequence a(n), let's call it `, P(x), ` satisfies the algebraic equation`):
print(``):
print(gu[1]):
print(``):

print(``):
print(`The sequence a(n) satisfies the linear recurrence `):
print(``):
print(gu[2][1]):
print(``):
print(`subject to the initial conditions `):
print(``):
print(gu[2][2]):
print(``):
print(`Just for fun`, a(M), `equals `):
print(``):
print(gu[3]):
print(``):
print(`For the sake of the OEIS here are the first`, K, `terms. `):
print(``):
print(gu[4]):
print(``):
print(`------------------------------------------------------------`):

end:


#TheoremNumber(A,B,C,D,K,P,x,n,a,MaxC, M,Num): Verbose version of InfoABCD(A,B,C,D,K,P,x,n,a,MaxC, M) (q.v.). Try:
#TheoremNumber({1},{1},{1},{1},40,P,x,n,a,20,1000,3);
TheoremNumber:=proc(A,B,C,D,K,P,x,n,a,MaxC, M,Num) local gu,S:

gu:=InfoABCD(A,B,C,D,K,P,x,n,a,MaxC, M):


print(``):
print(`------------------------------------------------------------`):
print(``):
print(`Theorem Number`, Num ):
print(``):
if A={} and B={} and C={} and D={} then
print(`Let a(n) be the number of Dyck paths of semi-length n`):
else
print(`Theorem: Let a(n) be the number of Dyck paths of semi-length n obeying the following restrictions`):
if A<>{} then
 print(`The height of a peak can not belong to `, A):
fi:
if B<>{} then
 print(`The height of a valley can not belong to `, B):
fi:
if C<>{} then
 print(`No upward-run can belong to`, C):
fi:
if D<>{} then
 print(`No downward-run can belong to`, D):
fi:

fi:




S:=DyckPathsABCD(A,B,C,D,8):

if nops(S)<>{} then
print(`To make it crystal clear here is such a path of semi-length 8`):
 S:=S[rand(1..nops(S))()]:
print(``):
print(DrawDP(S));
fi:

print(``):

if nops(gu)=1 then
 print(`The first`, K, `terms of the sequence are `):
 print(gu[1]):
 RETURN():
fi:



if nops(gu)=2 then

print(`The ordinary generating function of the sequence a(n), let's call it `, P(x), ` satisfies the algebraic equation`):
print(``):
print(gu[1]):
print(``):
print(`For the sake of the OEIS, here are the first`, K, `terms `):
print(``):
print(gu[2]):
print(``):
RETURN():
fi:

print(``):
print(`The ordinary generating function of the sequence a(n), let's call it `, P(x), ` satisfies the algebraic equation`):
print(``):
print(gu[1]):
print(``):

print(``):
print(`The sequence a(n) satisfies the linear recurrence `):
print(``):
print(gu[2][1]):
print(``):
print(`subject to the initial conditions `):
print(``):
print(gu[2][2]):
print(``):
print(`Just for fun`, a(M), `equals `):
print(``):
print(gu[3]):
print(``):
print(`For the sake of the OEIS here are the first`, K, `terms. `):
print(``):
print(gu[4]):
print(``):
print(`------------------------------------------------------------`):

end:








#TheoremR(A,B,C,D,r,K,P,x,n,a,MaxC, M): Verbose version of InfoABCDr(A,B,C,D,r,K,P,x,n,a,MaxC, M) (q.v.). Try:
#TheoremR({1},{1},{1},{1},r,40,P,x,n,a,20,1000);
TheoremR:=proc(A,B,C,D,r,K,P,x,n,a,MaxC, M) local gu,S:

gu:=InfoABCDr(A,B,C,D,r,K,P,x,n,a,MaxC, M):

print(``):
print(`------------------------------------------------------------`):
print(``):
print(`On The Number of Dyck Paths Obeying Certain Restrictions`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
if A={} and B={} and C={} and D={} then
print(`Theorem: Let a(n) be the number of Dyck paths of semi-length n`):
else
print(`Theorem: Let a(n) be the number of Dyck paths of semi-length n obeying the following restrictions`):
if A<>{} then
 print(`The height of a peak can't be equal (for a non-negative integer r) to something of the form `, A):
fi:
if B<>{} then
 print(`The eight of a valley can't be  (for a non-negative integer r) equal to something of the form  `, B):
fi:
if C<>{} then
 print(`No upward-run can  be in (for a non-negative integer r)   equal to something of the form `, C):
fi:
if D<>{} then
 print(`No downward-run can be  in (for a non-negative integer r)   equal to something of the form `, D):
fi:

fi:




S:=DyckPathsABCD(A,B,C,D,8):

if nops(S)<>{} then
print(`To make it crystal clear here is such a path of semi-length 8`):
 S:=S[rand(1..nops(S))()]:
print(DrawDP(S));
fi:

print(``):



 if nops(gu)=1 then
  print(`The first `, K, ` terms of the sequence are `):
  print(gu[1]):
  RETURN():
 fi:



if nops(gu)=2 then
print(`The ordinary generating function of the sequence a(n), let's call it `, P(x), ` satisfies the algebraic equation`):
print(``):
print(gu[1]):
print(``):
print(`For the sake of the OEIS are the first `, K, ` terms `):
print(``):
print(gu[2]):
print(``):
RETURN():
fi:


print(``):
print(`The ordinary generating function of the sequence a(n), let's call it `, P(x), ` satisfies the algebraic equation`):
print(``):
print(gu[1]):
print(``):

print(`The sequence a(n) satisfies the linear recurrence `):
print(``):
print(gu[2][1]):
print(``):
print(`subject to the initial conditions `):
print(``):
print(gu[2][2]):
print(``):
print(`Just for fun`, a(M), `equals `):
print(``):
print(gu[3]):
print(``):
print(`For the sake of the OEIS here are the first`, K, `terms. `):
print(``):
print(gu[4]):
print(``):
print(`------------------------------------------------------------`):

end:




#TheoremRnumber(A,B,C,D,r,K,P,x,n,a,MaxC, M,Num): Verbose version of InfoABCDr(A,B,C,D,r,K,P,x,n,a,MaxC, M) (q.v.). Try:
#TheoremRnumber({1},{1},{1},{1},r,40,P,x,n,a,20,1000,3);
TheoremRnumber:=proc(A,B,C,D,r,K,P,x,n,a,MaxC, M,Num) local gu,S:

gu:=InfoABCDr(A,B,C,D,r,K,P,x,n,a,MaxC, M):

print(``):
print(`------------------------------------------------------------`):
print(``):
print(`Theorem Number`, Num):
print(``):
if A={} and B={} and C={} and D={} then
print(` Let a(n) be the number of Dyck paths of semi-length n`):
else
print(`Let a(n) be the number of Dyck paths of semi-length n obeying the following restrictions`):
if A<>{} then
 print(`The height of a peak can't be equal (for a non-negative integer r) to something of the form `, A):
fi:
if B<>{} then
 print(`The height of a valley can't be  (for a non-negative integer r) equal to something of the form  `, B):
fi:
if C<>{} then
 print(`No upward-run can  be in (for a non-negative integer r)   equal to something of the form `, C):
fi:
if D<>{} then
 print(`No downward-run can be  in (for a non-negative integer r)   equal to something of the form `, D):
fi:

fi:




S:=DyckPathsABCD(A,B,C,D,8):

if nops(S)<>{} then
print(`To make it crystal clear here is such a path of semi-length 8`):
 S:=S[rand(1..nops(S))()]:
print(DrawDP(S));
fi:

print(``):

 if nops(gu)=1 then
 print(`The first `, K, ` terms of the sequence are `):
 print(gu[1]):

  RETURN():
 fi:

if nops(gu)=2 then
print(`The ordinary generating function of the sequence a(n), let's call it`, P(x) , ` satisfies the algebraic equation`):
print(``):
print(gu[1]):
print(``):
print(`For the sake of the OEIS are the first `, K, ` terms `):
print(``):
print(gu[2]):
print(``):
RETURN():
fi:


print(``):
print(`The ordinary generating function of the sequence a(n), let's call it`, P(x) , ` satisfies the algebraic equation`):
print(``):
print(gu[1]):
print(``):

print(`The sequence a(n) satisfies the linear recurrence `):
print(``):
print(gu[2][1]):
print(``):
print(`subject to the initial conditions `):
print(``):
print(gu[2][2]):
print(``):
print(`Just for fun`, a(M), `equals `):
print(``):
print(gu[3]):
print(``):
print(`For the sake of the OEIS here are the first`, K, `terms. `):
print(``):
print(gu[4]):
print(``):
print(`------------------------------------------------------------`):

end:



#Paper(a0,K,P,x,n,a,MaxC,M): inputs a positve integer A, and K,P,x,n,a,MaxC,M the same as in
#Theorem(A,B,C,D,K,P,x,n,a,MaxC, M) (look it up), outputs an article with 2^(4*a0) theorems outputted
#by Theorem(A,B,C,D,K,P,x,n,a,MaxC, M) where A,B,C,D, range over the subsets of {1, ..,a0}.
#Some of the theorems may be trivially equivalent due to redudancy of the conditions.Try:
#Paper(1,100,P,x,n,a,16,1000);
Paper:=proc(a0,K,P,x,n,a,MaxC,M) local i,A,B,C,D,Num,S,t0:
S:=powerset({seq(i,i=1..a0)}):

t0:=time():
Num:=0:

print(``):
print(`---------------------------------------------------------------------------`):
print(``):
print(`Enumeration of The Number of Dyck Paths of Semi-Length n obeying various restrictions`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(` The number of Dyck paths with semi-length n is famously the Catalan numbers.`):
print(`In this article we will  explicitly enumerate Dyck paths with four kinds of restictions`):
print(` (i) the heights of the peaks are not allowed to take certain values`):
print(` (ii) the heights of the valleys are not allowed to take certain values`):
print(` (iii) the  upward runs can't have certain values`):
print(` (iv) the  downward runs can't have certain values`):
print(``):
print(`All the above restictions will be subsets of`, {seq(i,i=1..a0)}):

for A from 1 to nops(S) do
for B from 1 to nops(S) do
for C from 1 to nops(S) do
for D from 1 to nops(S) do
 Num:=Num+1:
 TheoremNumber(S[A],S[B],S[C],S[D],K,P,x,n,a,MaxC, M,Num):
od:
od:
od:
od:

print(``):
print(`This concludes this exciting paper with its`, Num, `theorems that took`, time()-t0, `to generate. `):
print(``):
print(`--------------------------------------------------`):


end:



#PaperR(a0,b0,K,P,x,n,a,MaxC,M): inputs a positve integer A, and K,P,x,n,a,MaxC,M the same as in
#Theorem(A,B,C,D,K,P,x,n,a,MaxC, M) (look it up), outputs an article with 16 theorems outputted
#by TheoremR(A,B,C,D,r,K,P,x,n,a,MaxC, M) where A,B,C,D, range {{},{a0*r+b0}}
#Some of the theorems may be trivially equivalent due to redudancy of the conditions.Try:
#PaperR(2,3,100,P,x,n,a,20,1000);
PaperR:=proc(a0,b0,K,P,x,n,a,MaxC,M) local r,A,B,C,D,Num,S,t0:
S:={{},{a0*r+b0}}:

t0:=time():
Num:=0:

print(``):
print(`---------------------------------------------------------------------------`):
print(``):
print(`Enumeration of The Number of Dyck Paths of Semi-Length n obeying various restrictions`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(` The number of Dyck paths with semi-length n is famously the Catalan numbers.`):
print(`In this article we will  explicitly enumerate Dyck paths with four kinds of restictions`):
print(`where   one of more of the conditions below are forbidden `):
print(` (i) the heights of the peaks are not allowed to to be of the from`, a0*r+b0):
print(` (ii) the heights of the valleys are not allowed to take certain values`, a0*r+b0):
print(` (iii) the  upward runs can't have certain values`, a0*r+b0):
print(` (iv) the  downward runs can't have certain values`, a0*r+b0):
print(``):

for A from 1 to nops(S) do
for B from 1 to nops(S) do
for C from 1 to nops(S) do
for D from 1 to nops(S) do
 Num:=Num+1:
 TheoremRnumber(S[A],S[B],S[C],S[D],r,K,P,x,n,a,MaxC, M,Num):
od:
od:
od:
od:

print(``):
print(`This concludes this exciting paper with its`, Num, `theorems that took`, time()-t0, `to generate. `):
print(``):
print(`--------------------------------------------------`):


end:



##End Stories




#Start Alison Bu's bijection
ez:=proc(): print(` AllW(n), AJ(W) , invAJ(W) , CheckAJ(n)  `): end:
with(combinat):


#invAJ(W): inputs a walk W of semi-length n from [0,0] to [2*n,0] for some n
#outputs Alison J. Bu's mapping it to a Dyck path. try:
#invAJ([-1,1,-1,1]);
invAJ:=proc(W) local n,DP,i,i1:
  n:=nops(W):
  DP:=[1]:
  if W[1]=-1 then DP:=[op(DP),-1,1]:
  else DP:=[op(DP), 1,1]: fi:
  for i from 2 to n do
    if add(W[i1],i1=1..i-1)=0 then
      if add(W[i1],i1=1..i)=-1 then DP:=[op(DP),-1,1]:
      else DP:=[op(DP),1,1]: fi:
    elif add(W[i1],i1=1..i-1)=-1 and add(W[i1],i1=1..i)=0 then DP:=[op(DP),-1,1]:
    elif abs(add(W[i1],i1=1..i-1))<abs(add(W[i1],i1=1..i)) then DP:=[op(DP),1,1]:
    else DP:=[op(DP),-1,-1]: fi: od:
   [op(DP),-1]:
 end:

#AJ(W): inputs a Dyck path with semi-length 2*n+2 with peak heights and valley-heights that are positive even integers
#outputs Alison J. Bu's mapping it to a word in U^n D^n. Try:
#AJ(invAJ([-1,1,-1,1]));

 AJ:=proc(DP) local n,W,i,j,i1:
  n:=(nops(DP)-2)/2:
  if DP[2]=1 then W:=[1]:
  else W:=[-1]: fi:
  for i from 2 to n do
    if DP[2*i]=1 then W:=[op(W),sign(add(W[j],j=1..i-1))]:
    else
      if DP[2*i+1]=-1 then W:=[op(W),-sign(add(W[j],j=1..i-1))]:
      else
         if add(W[i1],i1=1..i-1)=-1 then W:=[op(W),1]:
        elif add(W[i1],i1=1..i-1)=0 then W:=[op(W),-1]: fi:
   fi: fi: od:
  W:
end:


#AllW(n): all words in {1,-1} with n 1's and n -1's. Try:
#AllW(6);
AllW:=proc(n):
convert(permute([1$n, (-1)$n]),set):
end:


#CheckAJ(n): Checks Alison J. Bu's bijection for n
CheckAJ:=proc(n) local gu,gu1,r,guI,mu,muI, mu1:
gu:=DyckPathsABCDr({2*r+2},{2*r+2},{},{},r,2*n+1):
mu:=AllW(n):
guI:={seq(AJ(gu1),gu1 in gu)}:
muI:={seq(invAJ(mu1),mu1 in mu)}:

evalb(guI=mu) and evalb(muI=gu) and
evalb({seq(evalb(invAJ(AJ(gu1))=gu1),gu1 in gu)}={true}) and
evalb({seq(evalb(AJ(invAJ(mu1))=mu1),mu1 in mu)}={true}):

end:

#End Alison Bu's bijection

###Start Robert Dougherty-Bliss's bijective proof of V. Retakh's result

###
### To test this file, try check_motzkin_bijection_verbose(k) for small values
### of k.
###

# Bijectively turn a Dyck n-path with peaks in {1, 2, 4, 6, ...} into a Motzkin
# n-path.
dyck_to_motzkin := proc(path)
    local height, new_path, k, first, second:
    height := 0:
    new_path := []:

    for k from 1 to nops(path) by 2 do
        first := path[k]:
        second := path[k + 1]:

        height := height + first + second:

        if first = -1 and second = 1 then
            new_path := [op(new_path), 0]:
        elif first = 1 and second = 1 then
            new_path := [op(new_path), 1]:
        elif first = -1 and second = -1 then
            new_path := [op(new_path), -1]:
        else
            new_path := [op(new_path), 0]:
            if height <> 0 then
                print("Something has gone wrong!");
            fi:
        fi:
    od:

    new_path:
end:

# Bijectively turn a Motzkin n-path into a Dyck n-path with peaks in {1, 2, 4,
# 6, ...}.
motzkin_to_dyck := proc(path)
    local height, new_path, x:
    height := 0:
    new_path := []:

    for x in path do
        height := height + x:

        if x = 1 then
            new_path := [op(new_path), 1, 1]:
        elif x = -1 then
            new_path := [op(new_path), -1, -1]:
        elif height > 0 then
            new_path := [op(new_path), -1, 1]:
        else
            new_path := [op(new_path), 1, -1]:
        fi:
    od:

    new_path:
end:

# Return [A, B, C, R], where:
#   A = number of Dyck n-paths with peaks in {1, 2, 4, 6, 8, ...}
#   B = number of Motzkin n-paths.
#   C = Size of the image of A under `motzkin_bijection`
#   R = a list of pairs of elements and images who coincide under the map.
# In particular, the Motzkin map is a bijection when A = B = C.
check_motzkin_bijection := proc(n)
    local A, B, image, repeats, p, new_path:

    # Paths with peaks in {1, 2, 4, 6, ...}.
    A := DyckPathsABCDr({2 * r + 3}, {}, {}, {}, r, n):
    # Paths with no up-runs longer than 2.
    # (Same number as Motzkin paths, see the OEIS entry.)
    B := DyckPathsABCDr({}, {}, {r + 3}, {}, r, n):

    image := {}:
    repeats := {}:
    for p in A do
        new_path := dyck_to_motzkin(p):
        if new_path in image then
            repeats := {op(repeats), [p, new_path]}:
        fi:
        image := {op(image), new_path}:
    od:

    [nops(A), nops(B), nops(image), repeats]:
end:

# Return a list of paths where motzkin_to_dyck fails to invert dyck_to_motzkin.
# If this list is empty, then motzkin_to_dyck is injective for the given n.
check_motzkin_injective := proc(n)
    local A, failures, path,r:
    # Paths with peaks in {1, 2, 4, 6, ...}.
    A := DyckPathsABCDr({2 * r + 3}, {}, {}, {}, r, n):
    failures := []:

    for path in A do
        if motzkin_to_dyck(dyck_to_motzkin(path)) <> path then
            failures := [op(failures), path]:
        fi:
    od:

    return failures:
end:

check_motzkin_bijection_verbose := proc(max_n)
    local k:
    print(`[# of restricted Dyck paths, # of Motzkin paths, size of image]`);
    for k from 1 to max_n do
        print(check_motzkin_bijection(k));
    od;
end:

###end Robert Dougherty-Bliss's bijective proof of V. Retakh's result