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

with(combinat):

Digits:=100:

print(`First Written: July 2022: tested for Maple 20  `):
print():
print(`This is AreaBounce.txt, A Maple package`):
print(`accompanying AJ Bu and Doron Zeilberger's  article: `):
print(` Experimenting with the  Area  Statistics of Generalized Dyck Paths`):
print(`-------------------------------`):

print(`-------------------------------`):
print():
print(`The most current version is available on WWW at:`):
print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/AreaBounce.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`-------------------------------`):
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):

print(`-------------------------------`):
print(`For a list of the supporting functions type: ezra1();`):
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print():
 print(`-------------------------------`):


ezra1:=proc()
if args=NULL then

print(` The supporting procedures are: `):
print(`area, bounce, DPn, DPnk, qbin `):

else
 ezra(args):
fi:
end:

ezra:=proc()

if args=NULL then

print(` AreaBounce.txt: A Maple package for experimenting with the Area and Bounce Statistics on Dyck Paths`):
print(`The main procedures are : CheckSym, CheckSymF, Fn, Fnk, Gn, Gnk, Thm34`):



elif nargs=1 and args[1]=area then
print(`area(P): The area of the Dyck Path P, try: `):
print(` area([1,1,1,-1,-1,-1]);`):


elif nargs=1 and args[1]=bounce then
print(` bounce(P): The bounce statistic of the Dyck Path P, try: `):
print(` bounce([1,1,1,-1,-1,-1]);`):


elif nargs=1 and args[1]=CheckSym then
print(`CheckSym(n): Checks that Fn(n,q,t)=Fn(n,t,q). Try:`):
print(`CheckSym(5);`):



elif nargs=1 and args[1]=CheckSymF then
print(`CheckSymF(n): Checks that Gn(n,q,t)=Gn(n,t,q). Same output as CheckSym(n), but faster. Try:`):
print(`CheckSymF(5);`):


elif nargs=1 and args[1]=DPnk then
print(`DPnk(n,k): The set of Dyck-paths of semi-length n that start with k 1's. Try:`):
print(` DPnk(5);`):

elif nargs=1 and args[1]=DPn then
print(`DPn(n): The set of Dyck-paths of semi-length n. Try:`):
print(` DPn(5);`):


elif nargs=1 and args[1]=Fn then
print(`Fn(n,q,t): The weight-enumertor of all Dyck-paths of semi-length n according to the weight q^area(P)*t^bounce(P), In other words the weight-enumerator of DPn(n). It is done directly. Try:`):
print(`Fn(5,q,t);`):


elif nargs=1 and args[1]=Fnk then
print(`Fnk(n,k,q,t): The weight-enumertor of all Dyck-paths of semi-length n that start with k 1's, according to the weight q^area(P)*t^bounce(P), In other words the weight-enumerator of DPnk(n,k). It is done directly. Try:`):
print(`Fn(5,q,t);`):


elif nargs=1 and args[1]=Gn then
print(`Gn(n,q,t): The weight-enumertor of all Dyck-paths of semi-length n according to the weight q^area(P)*t^bounce(P), In other words the weight-enumerator of DPn(n). Done via the recurrence of Haglund's theorem. `):
print(` Same output as Fn(n,q,t), but faster. Try:`):
print(`Gn(5,q,t);`):


elif nargs=1 and args[1]=Fnk then
print(`Fnk(n,k,q,t): The weight-enumertor of all Dyck-paths of semi-length n that start with k 1's, according to the weight q^area(P)*t^bounce(P), In other words the weight-enumerator of DPnk(n,k). It is done directly. Try:`):
print(`Fn(5,q,t);`):


elif nargs=1 and args[1]=Gnk then
print(`Gnk(n,k,q,t): The weight-enumertor of all Dyck-paths of semi-length n that start with k 1's, according to the weight q^area(P)*t^bounce(P), In other words the weight-enumerator of DPnk(n,k). Done via the recurrence of Haglund's theorem. `):
print(` Same output as Fnk(n,k,q,t), but faster. Try:`):
print(`Gnk(5,3,q,t);`):


elif nargs=1 and args[1]=qbin then
print(`qbin(q,a,b): The   q-binomial coefficient  (1-q^a)*...(1-q^(a-b+1))/((1-q)*...*(1-q^b)). Try: `):
print(`qbin(q,5,3);`):


elif nargs=1 and args[1]=Thm34 then
print(`Thm34(n,k):  verifies empirically Haglund's Theorem 3.4. Try: `):
print(`Thm34(5,4);`):


print(``):


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

fi:



end:

Help:=proc():
print(`BouncePath(P), NorthPath(P), bounceInfo(P), Agf(n,t),qbin(q,a,b)`):
print(`Thm34(n,k)`):
end:

BouncePath:=proc(P) local n,height,L,pi1,heightB,co,i,c:
	n:=nops(P)/2:
	L:=[]:
	height:=0:
	heightB:=0:
	pi1:=subs({-1=0},P):
	
	while heightB<n do
		co:=0:
		if height=heightB then 
			for i from 2*heightB+1 to 2*n while pi1[i]=1 do 
				co:=co+1:
				
			od:

		else
			c:=height-heightB: 
			for i from 2*heightB+1 to 2*n while (c<>0 or pi1[i]=1) do
				co:=co+1:
				if pi1[i]=0 then
					c:=c-1:
				fi:
			
			od:

		fi:
		
	L:=[op(L),1$co,(-1)$co]: 
	height:=height+add(pi1[i],i=2*heightB+1..min(2*heightB+2*co,2*n)): 
		     heightB:=heightB+co:
		  

	od:
	L: 
end:


#NorthPath(P): inputs a path P, and removes up to the end of the first North run and also outputs 
#the length of the first North run
NorthPath:=proc(P) local n,pi1,co,i:
	n:=nops(P):
	pi1:=P:
	
	while pi1[1]=-1 do
		pi1:=[op(2..nops(pi1),pi1)]: 
	od:
	co:=0:
	
	for i from 1 while pi1[i]=1 do
		co:=co+1:
	od:
	
	[op(co+1..nops(pi1),pi1)],co:
end:

#bounceInfo(P): Inputs a Dyck path P
#outputs the bounce, the number of bounces, the length of each bounce, and the peaks
bounceInfo:=proc(P) local n,height,pi1,NP,lengths,peaks,P1,i:
	P1:=BouncePath(P):
	n:=nops(P1)/2:
	pi1:=P1:
	height:=[0]:
	lengths:=[]:
	peaks:=[]:
	while height[-1]<>n do
		NP:=NorthPath(pi1):
		pi1:=NP[1]: 
		height:=[op(height),height[-1]+NP[2]]: 
		lengths:=[op(lengths),NP[2]]:
		peaks:=[op(peaks), [height[-2],height[-1]]]:
	od: 
	add(n-height[i],i=2..nops(height)-1), nops(height)-1,lengths,peaks:
end:

#bounce(P): Inputs a Dyck path P
#outputs the bounce
bounce:=proc(P) local n,height,pi1,NP,P1,i:
	P1:=BouncePath(P):
	n:=nops(P1)/2:
	pi1:=P1:
	height:=[0]:
	while height[-1]<>n do
		NP:=NorthPath(pi1):
		pi1:=NP[1]: 
		height:=[op(height),height[-1]+NP[2]]: 
	od: 
	add(n-height[i],i=2..nops(height)-1):
end:


#area(P): inputs  a Dyck path P
#outputs the area
area:=proc(P) local n,width,pi1,height,area,NP,i:
	n:=nops(P)/2:
	width:=0:
	pi1:=P:
	height:=0:
	area:=0:
	
	while height<n do
		NP:=NorthPath(pi1):
		pi1:=NP[1]: 
		for i from 1 to nops(pi1) while pi1[i]=-1 do 
			width:=width+1: 
		od: 
		height:=height+NP[2]: 
		area:=area+height*width: 
		width:=0: 
	od:
	area-add(i,i=1..n):
end:


DPn:=proc(n) local i, P:
	[seq([seq( op([ 1$P[2*i-1], (-1)$abs(P[2*i])]),i=1..nops(P)/2  )] , P in V(n,n))]:
end:


Finfo:=proc(n,q,t) local paths, P:
	paths:=DPn(n):
	for P in paths do
		print(P, area(P), bounce(P)):
	od:

	add(q^area(P)*t^bounce(P), P in paths):
end:


#Fn(n,q,t): The weight-enumertor of all Dyck-paths of semi-length n according to the weight q^area(P)*t^bounce(P), done directly. Try:
#Fn(5,q,t);
Fn:=proc(n,q,t) local paths, P:
	paths:=DPn(n):

	add(q^area(P)*t^bounce(P), P in paths):
end:

Agf:=proc(n,t) local k, su:
option remember:
 if n=0 then 
    RETURN(1):
 fi:

 su:=0:

 for k from 0 to n-1 do
   su:=expand(su+t^k*Agf(k,t)*Agf(n-1-k,t)):
od:

su:
end:

#qbin(q,a,b): The q-binomial coefficients 
qbin:=proc(q,a,b) local i:
	normal(mul(1-q^i,i=a-b+1..a)/mul(1-q^i,i=1..b)):
end:

#DPnk(n,k): the set of paths from [0,0] to [n,n] that start with exactly k "1" steps (1 is North -1 is East)
DPnk:=proc(n,k) local P,j,i:
	if n=0 and k=0 then 
		RETURN([[]]):
	fi:

	[seq( seq( [1$k,(-1)$i,seq( op([1$P[2*j-1],(-1)$abs(P[2*j])]), j=1..nops(P)/2)],P in V1(n-k,n-i,k-i)),i=1..k)]:
end:

#Fnk(n,k,q,t):
Fnk:=proc(n,k,q,t) local P,L:
  L:=DPnk(n,k): 
  add(q^area(P)*t^bounce(P),P in L):
end:

#Thm34(n,k,q,t): checks theorem 3.4 in The q, t-Catalan Numbers and the Space of Diagonal Harmonics
#by Haglund
Thm34:=proc(n,k) local r,q,t:
	evalb(simplify(Fnk(n,k,q,t)=add(qbin(q,r+k-1,r)*t^(n-k)*q^binomial(k,2)*Fnk(n-k,r,q,t),r=0..n-k))):
end:

##### From Dyck.txt  #######

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

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

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

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

gu:={}:

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

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

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

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

gu:={}:

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

gu:

end:


#Gnk(n,k,q,t): checks theorem 3.4 in The q, t-Catalan Numbers and the Space of Diagonal Harmonics
Gnk:=proc(n,k,q,t) local r:
option remember:
if n=0 then
 if k=0 then
  RETURN(1):
else
 RETURN(0):
fi:
fi:

if k<=0 then
 RETURN(0):
fi:

	expand(normal(add(qbin(q,r+k-1,r)*t^(n-k)*q^binomial(k,2)*Gnk(n-k,r,q,t),r=0..n-k))):

end:

Gn:=proc(n,q,t) local k: add(Gnk(n,k,q,t),k=1..n):end:


#CheckSym(n): Checks that Fn(n,q,t)=Fn(n,t,q)
CheckSym:=proc(n) local q,t,gu:
gu:=Fn(n,q,t):
evalb(subs({q=t,t=q},gu)-gu=0):
end:



#CheckSymF(n): Checks that Gn(n,q,t)=Gn(n,t,q)
CheckSymF:=proc(n) local q,t,gu:
gu:=Gn(n,q,t):
evalb(subs({q=t,t=q},gu)-gu=0):
end: