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


print(`First Written: Oct. 2024: tested for Maple 2020 `):
print(`Version  : Oct. 2024  `):
print():
print(`This is GenPark.txt, A Maple package`):
print(`accompanying AJ Bu  and Doron Zeilberger's  article: `):
print(`" A Short Proof that the number of (a,b)-parking functions of length n is a*(a+bn)^(n-1)" `):

print():
print(`The most current version of this package is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/GenPark.txt   .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
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():
 





ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are:`):
print(` AnA, anr, bnr, cnkA, BnA, CnA, Findrec, PtorOpe, PtorOpeF `):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` GenPark.txt: A Maple package for experimenting with, and enumerating, generalized parking  functions  `):

print(`The MAIN procedures are`):
print(` anA, bnA, cnA, Fnab, GuessAP `):



elif nargs=1 and args[1]=AnA then
print(`AnA(n,A): Inputs a positive integer n, an increasing sequence A (with >=n terms) of pos. numbers, outputs the SET of increasing lists [a1,..., an] such that a[i]<=A[i]. Using the SECOND (clever) approach. For example try:`):
print(` AnA(5,[seq(i,i=1..30)]); `):


elif nargs=1 and args[1]=Anr then
print(`Anr(n,r): Same as AnA(n,[seq(r-1+i,i=1..n+3)]);. Try:`):
print(` Anr(5,4); `):

elif nargs=1 and args[1]=anA then
print(`anA(n,A): Inputs a positive integer n, an increasing sequence A (with >=n terms) of pos. numbers, outputs the number of increasing lists [a1,..., an] such that a[i]<=A[i]. Using the (clever) SECOND approach  For example try:`):
print(`[seq(anA(n,[seq(i,i=1..50)]),n=1..30)];`):


elif nargs=1 and args[1]=anr then
print(`anr(n,r): Same as anA(n,[seq(r-1+i,i=1..n+5)]). Try: `):
print(`[seq(anr(n,1),n=1..30)];`):


elif nargs=1 and args[1]=bnA then
print(`bnA(n,A): Inputs a positive integer n, an increasing sequence A (with >=n terms) of pos. numbers, outputs the number of  SCRAMBLINGS IF increasing lists [a1,..., an] such that a[i]<=A[i]. `):
print(`In other words the number of A-parking functions of length n.`):
print(`For example try:`):
print(`[seq(bnA(n,[seq(i,i=1..50)]),n=1..20)];`):


elif nargs=1 and args[1]=bnr then
print(`bnr(n,r): Same as bnA(n,[seq(r-i+1,i=1..2*n)]);`):
print(`For example try:`):
print(`[seq(bnr(n,2),n=1..20)];`):


elif nargs=1 and args[1]=BnA then
print(`BnA(n,A): Inputs a positive integer n, an increasing sequence A (with >=n terms) of pos. numbers, outputs the SET of  SCRAMBLINGS IF increasing lists [a1,..., an] such that a[i]<=A[i]. `):
print(`In other words the set of A-parking functions of length n.`):
print(`For example try:`):
print(`BnA(4,[seq(i,i=1..50)]);`):

elif nargs=1 and args[1]=cnkA then
print(`cnkA(n,k,A): Inputs a positive integer n, an increasing sequence A (with >=n terms) of pos. numbers, and k<=A[n], outputs the number of increasing lists [a1,..., an] such that a[i]<=A[i]  and a[n]=k. For example try:`):
print(`cnkA(10,3,[seq(i,i=1..50)]);`):


elif nargs=1 and args[1]=cnA then
print(`cnA(n,A): Inputs a positive integer n, an increasing sequence A (with >=n terms) of pos. numbers, outputs the number of increasing lists [a1,..., an] such that a[i]<=A[i]. Using the FIRST (naive) approach (via cnkA) For example try:`):
print(`[seq(cnA(n,[seq(i,i=1..50)]),n=1..30)];`):


elif nargs=1 and args[1]=CnA then
print(`CnA(n,A): Inputs a positive integer n, an increasing sequence A (with >=n terms) of pos. numbers, outputs the SET of increasing lists [a1,..., an] such that a[i]<=A[i]. Using the FIRST (naive) approach. For example try:`):
print(` CnA(5,[seq(i,i=1..30)]); `):

elif nargs=1 and args[1]=Findrec then
print(`Findrec(L,n,N): Given a list L starting at n=1 tries to find a linear recurrence equation with`):
print(`poly coffs.  of ofder 1`):
print(`e.g. try Findrec([seq(i!,i=1..30)],n,N);`):

elif nargs=1 and args[1]=Fnab then
print(`Fnab(n,a,b): the number of  genarlized (a,b)-parking functions of length n, i.e. all scramblings of sequences x[1],...,x[n] such that 1<=x[1]<=..<=x[n], and x[i]<=a+b*(i-1) for i=1..n. Computed by Dynamial programming.Try`):
print(`Fnab(4,a,b);`):


elif nargs=1 and args[1]=Gnab then
print(`Gnab(n,a,b): the conjectured, then proved, expression for Fnab(n,a,b), namely: (a*(a+b*n)^(n-1). Try`):
print(`Gnab(4,a,b);`):

elif nargs=1 and args[1]=GuessAP then
print(`GuessAP(a,b,n,R,K):  tries to guess  a formula for the number of increasing sequences of length n x[1]<=...<=x[n] of length n such that x[i]<=a*i-b (for 0<=b<a) using up to K data points. Try:`):
print(`GuessAP(1,0,n,R,30);`):

elif nargs=1 and args[1]=PtorOpe then
print(`PtorOpe(ope,n,N,R): solving the recurrence ope(n,N)x(n)=0, x(0)=1`):
print(`in terms of rising factorial R(n,k) (R is a symbol)`):
print(`For example try PtorOpe((n+1)*N-1,n,N,R);`):


elif nargs=1 and args[1]=PtorOpeF then
print(`PtorOpeF(ope,n,N): solving the recurrence ope(n,N)x(n)=0, x(0)=1`):
print(`in terms of rfactorials `):
print(`For example try PtorOpeF((n+1)*N-1,n,N);`):

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

fi:


end:



with(combinat):



##start from twoFone

#PtorOpe(ope,n,N,R): solving the recurrence ope(n,N)x(n)=0, x(0)=1
#in terms of rising factorial R(n,k) (R is a symbol)
#For example try PtorOpe((n+1)*N-1,n,N,R);
PtorOpe:=proc(ope,n,N,R)
local rat,c,t1,b1,t1c,b1c,i:

if degree(ope,N)<>1 then
 RETURN(FAIL):
fi:

rat:=normal(-coeff(ope,N,0)/coeff(ope,N,1)):
t1:=expand(numer(rat)):
b1:=expand(denom(rat)):

t1c:=coeff(t1,n,degree(t1,n)):
b1c:=coeff(b1,n,degree(b1,n)):

c:=t1c/b1c:

t1:=expand(t1/t1c):
b1:=expand(b1/b1c):


t1:=factor(t1):
b1:=factor(b1):

t1:=[solve(t1,n)]:
b1:= [solve(b1,n)]:

c^n* convert([seq(R(-t1[i],n),i=1..nops(t1))],`*`)/
convert([seq(R(-b1[i],n),i=1..nops(b1))],`*`):

end:

#PtorOpeF(ope,n,N): solving the recurrence ope(n,N)x(n)=0, x(0)=1
#in terms of factorials (R is a symbol)
#For example try PtorOpeF((n+1)*N-1,n,N,R);
PtorOpeF:=proc(ope,n,N)
local rat,c,t1,b1,t1c,b1c,i:

if degree(ope,N)<>1 then
 RETURN(FAIL):
fi:

rat:=normal(-coeff(ope,N,0)/coeff(ope,N,1)):
t1:=expand(numer(rat)):
b1:=expand(denom(rat)):

t1c:=coeff(t1,n,degree(t1,n)):
b1c:=coeff(b1,n,degree(b1,n)):

c:=t1c/b1c:

t1:=expand(t1/t1c):
b1:=expand(b1/b1c):


t1:=factor(t1):
b1:=factor(b1):

t1:=[solve(t1,n)]:
b1:= [solve(b1,n)]:

for i from 1 to nops(t1) do
 if type(t1[i],integer) and t1[i]>=0 then
   RETURN(FAIL):
 fi:
od:

for i from 1 to nops(b1) do
 if type(b1[i],integer) and b1[i]>=0 then
   RETURN(FAIL):
 fi:
od:

c^n* 
normal(convert([seq((-t1[i]+n-1)!/(-t1[i]-1)!,i=1..nops(t1))],`*`)/
convert([seq((-b1[i]+n-1)!/(-b1[i]-1)!,i=1..nops(b1))],`*`)):

end:



#end from twoFone

##start from Findrec 

##start from 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);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:


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

ope:=subs(var1,ope):

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

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

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

end:
 


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(L,n,N): Given a list L starting at n=1 tries to find a linear recurrence equation with
#poly coffs.  of ofder 1
#e.g. try Findrec([seq(i!,i=1..30)],n,N);
Findrec:=proc(f,n,N)
local DEGREE,ope:


for DEGREE from 0 to trunc(nops(f)/4)-4 do
 ope:=findrec([op(1..(2+DEGREE)*(2)+4,f)],DEGREE,1,n,N):
     if ope<>FAIL then
       RETURN(ope):
      fi:
 od:
FAIL:


end:

#end from Findrec

#cnrk: The number of lists [a1,..., an] such that a[i]<=r+i-1, a[1]<..<=a[n] and a[n]=k
cnrk:=proc(n,r,k) local i:
option remember:
if n=0 then
 RETURN(1):
fi:

if n=1 then
 if k<=r then
 RETURN(1):
else
 RETURN(0):
fi:
fi:


add(cnrk(n-1,r,i),i=1..min(r+n-2,k)):
end:

cnr:=proc(n,r) local k: option remember:add(cnrk(n,r,k),k=1..n+r-1):end:


#Cnrk: The set of lists [a1,..., an] such that a[i]<=r+i-1, a[1]<..<=a[n] and a[n]=k
Cnrk:=proc(n,r,k) local gu,i,lu,lu1:
option remember:
if n=0 then
 RETURN({[]}):
fi:

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

gu:={}:

for i from 1 to min(r+n-2,k) do
 lu:=Cnrk(n-1,r,i):
 gu:=gu union {seq([op(lu1),k],lu1 in lu)}:
od:
gu:
end:


Cnr:=proc(n,r) local k:
option remember:
{seq(op(Cnrk(n,r,k)),k=1..n+r-1)}:
end:


####

Anr:=proc(n,r) local k,gu,lu,lu1,i,gu1:
option remember:

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

if r>1 then
gu:=Anr(n,r-1):
gu:={seq(gu1+[1$n],gu1 in gu)}:
else
gu:={}:
fi:

for k from 1 to n do
 lu:=Anr(n-k,k+r-1):
  gu:=gu union {seq([1$k,op(lu1+[1$(n-k)]) ],lu1 in lu)}:
od:

gu:


end:

anr:=proc(n,r) local k:
option remember:
if n=0 then
 RETURN(1):
fi:
if n=1 then
 RETURN(r):
fi:

if r=0 then
 RETURN(0):
fi:

anr(n,r-1)+add(anr(n-k,k+r-1),k=1..n):

end:

#Bnr(n,r): all scamblings of the members of Anr(n,r)
Bnr:=proc(n,r) local gu,gu1:
gu:=Anr(n,r):

{seq(op(permute(gu1)),gu1 in gu)}:

end:



bnr:=proc(n,r) local k:
option remember:
if n=0 then
 RETURN(1):
fi:
if n=1 then
 RETURN(r):
fi:

if r=0 then
 RETURN(0):
fi:

bnr(n,r-1)+add(binomial(n,k)*bnr(n-k,k+r-1),k=1..n):

end:



######start General

#cnkA(n,k,A): The number of lists [a1,..., an] such that a[i]<=A[i] a[1]<..<=a[n] and a[n]=k
cnkA:=proc(n,k,A) local i:
option remember:

if nops(A)<=n then
 RETURN(FAIL):
fi:

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

if n=1 then
 if k<=A[1] then
 RETURN(1):
else
 RETURN(0):
fi:
fi:


add(cnkA(n-1,i,A),i=1..min(A[n-1],k)):
end:

cnA:=proc(n,A) local k: option remember:add(cnkA(n,k,A),k=1..A[n]):end:




anA:=proc(n,A) local k:
option remember:
if n=0 then
 RETURN(1):
fi:
if n=1 then
 RETURN(A[1]):
fi:

if min(op(A))<=0 then
 RETURN(0):
fi:

anA(n,A-[1$nops(A)])+add(anA(n-k,[op(k+1..nops(A),A)]-[1$(nops(A)-k)]),k=1..n):

end:


#I am here

AnA:=proc(n,A) local k,gu,lu,lu1,i,gu1:
option remember:

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


if min(op(A))>1 then
gu:=AnA(n,A-[1$nops(A)]):
gu:={seq(gu1+[1$n],gu1 in gu)}:
else
gu:={}:
fi:

for k from 1 to n do
 lu:=AnA(n-k,[op(k+1..nops(A),A)]-[1$(nops(A)-k)]):
  gu:=gu union {seq([1$k,op(lu1+[1$(n-k)]) ],lu1 in lu)}:
od:

gu:


end:

#CnAk: The set of lists [a1,..., an] such that a<=A and A[n]=k
CnAk:=proc(n,A,k) local gu,i,lu,lu1:
option remember:
if n=0 then
 RETURN({[]}):
fi:

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

gu:={}:

for i from 1 to min(A[n-1],k) do
 lu:=CnAk(n-1,A,i):
 gu:=gu union {seq([op(lu1),k],lu1 in lu)}:
od:
gu:
end:


CnA:=proc(n,A) local k:
option remember:
{seq(op(CnAk(n,A,k)),k=1..A[n])}:
end:



#BnA(n,A): all scamblings of the members of AnA(n,A)
BnA:=proc(n,A) local gu,gu1:
gu:=AnA(n,A):

{seq(op(permute(gu1)),gu1 in gu)}:

end:




bnA:=proc(n,A) local k:
option remember:
if n=0 then
 RETURN(1):
fi:
if n=1 then
 RETURN(A[1]):
fi:

if min(op(A))<=0  then
 RETURN(0):
fi:

bnA(n,A-[1$nops(A)]) +add(binomial(n,k)*bnA(n-k,[op(k+1..nops(A),A)]-[1$(nops(A)-k)] ),k=1..n):

end:


#GuessAP(a,b,n,R,K):  tries to guess  a formula for the number of increasing sequences of length n x[1]<=...<=x[n] of length n such that x[i]<=a*i-b (for 0<=b<a) using up to K data points, in terms of the raising factoral. R(a,n). Try:
#GuessAP(1,0,n,30,R);
GuessAP:=proc(a,b,n,R,K) local ope,L,K1,i,i1,N:

for K1 from 20  to K by 10 do
L:=[seq(anA(i,[seq(a*i1-b,i1=1..K1+4)]),i=1..K1)]:
ope:=Findrec(L,n,N):
if ope<>FAIL then
 RETURN(PtorOpe(ope,n,N,R)):
fi:
od:
FAIL:
end:


anARecConj:=proc(k,m,n,N) local i:
-(k+1)*((k+1)*n+k)/(k*n+k+1)*mul(((k+1)*n+k-i)/(k*n+k-i+1),i=1..m)*mul(((k+1)*n+(k-m)+i+1)/(k*n+k-m+i+1),i=m+1..k-1)+N:
end:

anARecTest:=proc(k,m,n) local L,Fr,n1,i,re,N:
L:=[seq(anA(n1,[seq(k*i-m,i=1..100)]),n1=1..50)]:
Fr:=Findrec(L,n,N,10):
re:=anARecConj(k,m,n,N): 
evalb(simplify(Fr-re)=0): 
end:


anAConj:=proc(k,m,N) local n:
[seq(binomial((k+1)*(n)+(k-1-m),n)*(k-m)/(k*n+k-m),n=1..N)]:
end:

anATest:=proc(k,m,N) local i,n:
evalb(anAConj(k,m,N)=[seq(anA(n,[seq(k*i-m,i=1..50)]),n=1..N)]):
end:

bnAConj:=proc(k,m,N) local n:
[seq((k-m)*(k*n+k-m)^(n-1),n=1..N)]:
end:

bnATest:=proc(k,m,N) local n,i:
evalb(bnAConj(k,m,N)=[seq(bnA(n,[seq(k*i-m,i=1..50)]),n=1..N)]):
end:











AJ:=proc(a,b,N) local n, i:[seq(bnA(n,[seq(a+b*(i-1),i=1..N)]),n=1..N)];end:
AJ1:=proc(a,b,N) local n: [seq(a*(a+b*n)^(n-1),n=1..N)]:end:


#Fnab(n,a,b): the number of  genarlized (a,b)-parking functions of length n, i.e. all scramblings of sequences x[1],...,x[n] such that 1<=x[1]<=..<=x[n], and x[i]<=a+b*(i-1) for i=1..n. Try
#Fnab(4,a,b);
Fnab:=proc(n,a,b) local r:
option remember:
if a=0 then
 RETURN(0):
fi:

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

Fnab(n,a-1,b)+add(binomial(n,r)*Fnab(n-r,a+b*r-1,b),r=1..n):
end:


#Gnab(n,a,b): (a*(a+b*n)^(n-1)
Gnab:=proc(n,a,b): a*(a+b*n)^(n-1):end:

Bdok:=proc(n,a,b) local r:a*(a+b*n)^(n-1)-add(binomial(n,r)*(a+b*r-1)*(a+b*n-1)^(n-r-1),r=0..n):end: