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


print(`First Written: March 20, 2024: tested for Maple 2020 `):
print(`Version  : March 20, 2024 `):
print():
print(`This is AbelBijection.txt, A Maple package`):
print(`accompanying Gil Kalai and Doron Zeilberger's  article: `):
print(`  Bijective and Automated Approaches to Abel Sums`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/mamarim/mamarimhtml/abelKZ.html `):
print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/AbelBijection.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(`For a list of the supporting functions type: ezra1();`):
print(`For a list of the checking functions type: ezraC();`):
print():
print(`-----------------------------------------------------`):






ezraS:=proc()
if args=NULL then

print(`The Story procedures are:  `):
print(``):

else
ezra(args):

fi:


end:




ezraC:=proc()
if args=NULL then

print(`The Checking procedures are:  CheckCO,  CheckCOrand, CheckGK, CheckGKrand, CheckGK1, Verify1 `):
print(``):

else
ezra(args):

fi:


end:

ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are:  AF, Cyc, Orb, TriplesAB, TriplesCD, RandAB, RandCD, RF `):
print(``):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` Abel.txt: A Maple package for  investigating Abel-type identities`):

print(`The MAIN procedures are:  B1,B2,GK, iGK `):


print(``):


elif nargs=1 and args[1]=AF then
print(`AF(m,n): All the functions from [m] to [n]. Try:`):
print(`AF(4,3);`):


elif nargs=1 and args[1]=B1 then
print(`B1(n,p): The left side of Eq. (1) in the Kalai-Zeilberger paper. Try:`):
print(`B1(10,4);`):


elif nargs=1 and args[1]=B2 then
print(`B2(n,p): The right side of Eq. (1) in the Kalai-Zeilberger paper. Try`):
print(`B2(10,4); Also try:`):
print(`{seq(seq(B1(n,p)-B2(n,p),n=0..10),p=0..10)};`):



elif nargs=1 and args[1]=CheckCO then
print(`CheckCO(n,p): checks the crucial observation for all functions from [n+p] to [n]. Try:`):
print(`CheckCO(4,0);`):

elif nargs=1 and args[1]=CheckCOrand then
print(`CheckCOrand(n,p,K): checks the crucial observation for K random functions from [n+p] to [n]. Try:`):
print(`CheckCOrand(40,0,10);`):


elif nargs=1 and args[1]=CheckGK then
print(`CheckGK(n,p): Check mapping GK for all functions from [n+p] to [n]. Try:`):
print(`CheckGK(4,0);`):

elif nargs=1 and args[1]=CheckGK1 then
print(`CheckGK1(f,p): Check mapping GK for the mapping f from [n+p] to [n]. Try:`):
print(`f:=RF(11,10): CheckGK1(f,1);`):

elif nargs=1 and args[1]=Cyc then
print(`Cyc(f,i): The cycle that ends the orbit starting at i. try:`):
print(`Cyc([4,1,3,2,5],3);`):


elif nargs=1 and args[1]=CheckGKrand then
print(`CheckGKrand(n,p,K): Check mapping GK for K random  functions from [n+p] to [n]. Try:`):
print(`CheckGKrand(10,2,100);`):

elif nargs=1 and args[1]=GK then
print(`GK(TL,p): Given a left-triple [A,B,f] where f maps [n+p] to [n], outputs the correponding right-triple [C,D,f] using the Gil Kalai mapping. Try:`):
print(`TL:=RandAB(10,0); GK(TL,0);`):

elif nargs=1 and args[1]=iGK then
print(`iGK(TR,p): Given a right-triple TR=[C,D,f] where f maps [n+p] to [n], outputs the correponding left-triple [A,B,f] using the Gil Kalai inverse mapping. Try:`):
print(`TR:=RandCD(10,0); iGK(TR,0);`):

elif nargs=1 and args[1]=Orb then
print(`Orb(f,i): inputs a function f from [n] to [n] and outputs the orbit that starts at i, try:`):
print(`Orb([1,1,1,1],2);`):

elif nargs=1 and args[1]=TriplesAB then
print(`TriplesAB(f,p): inputs a function f from [n+p] to [n] (where n:=nops(f)-p) and outputs the set of triples`):
print(`[A,B,f] such that A and B are disjoint, A union B = [n] and such that f(A) is a subset of A and f(B union {i,i=n+1..n+p}) is a subset of B. Try:`):
print(`TriplesAB([1,2,3],0);`):


elif nargs=1 and args[1]=TriplesCD then
print(`TriplesCD(f,p): inputs a function f from [n+p] to [n] (where n:=nops(f)-p) and outputs the set of triples`):
print(`[C,D,f] such that C union D =[n], C and D are disjoint, and such that f(D union {n+1, ...,n+p})=D. Try:`):
print(`TriplesCD([1,2,3],0);`):

elif nargs=1 and args[1]=RandAB then
print(`RandAB(n,p): a random triplet [f,A,B] such that AUB=[n], A intersect B={} and f(A) is a subset of A and f(B union {n+1,...,n+p}) is a subset of B.`):
print(`Try:`):
print(`RandAB(5,4);`):


elif nargs=1 and args[1]=RandCD then
print(`RandCD(n,p): a random triplet [f,C,D] such that CUD=[n], C intersect D={} f(D union {n+1,...n+p})=D`):
print(`Try:`):
print(`RandCD(5,0);`):

elif nargs=1 and args[1]=RF then
print(`RF(m,n): a random function from [m] to [n]. Try`):
print(`RF(7,4);`):

elif nargs=1 and args[1]=Verify1 then
print(`Verify1(n1,p1): verifies Eq. (1) of the paper for all n<=n1, p<=p1. Try:`):
print(`Verify1(5,3);`):

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

fi:


end:


with(combinat):

#Verify1(n1,p1): verifies Eq. (1) of the paper for all n<=n1, p<=p1. Try:
#Verify1(5,3);
Verify1:=proc(n1,p1) local k,n,p:
evalb({seq(seq(add(binomial(n,k)*k^k*(n-k)^(n-k+p),k=0..n)- add(binomial(n,k)*n^k*(n-k)!*stirling2(p+n-k,n-k),k=0..n),n=0..n1),p=0..p1)}={0}):
end:

#RandAB(n,p): a random triplet [f,A,B] such that f maps [n+p] to [n], AUB=[n], A intersect B={} and f(A) is a subset of A and f(B union {i,i=n+1..n+p)) is a subset of B.
#Try:
#RandAB(5,2);
RandAB:=proc(n,p) local ra,i,A,B,ra1,ra2,f:
ra:=rand(0..1):

A:={}:
B:={}:
for i from 1 to n do
if ra()=0 then
 A:=A union {i}:
else
 B:=B union {i}:
fi:
od:

ra1:=rand(1..nops(A)):
ra2:=rand(1..nops(B)):

f:=[]:

for i from 1 to n do
if member(i,A) then
 f:=[op(f),A[ra1()]]:
else
f:=[op(f),B[ra2()]]:
fi:
od:

for i from n+1 to n+p do
f:=[op(f),B[ra2()]]:
od:

[A,B,f]:

end:



#RandCD(n,p): a random triplet [f,C,D] such that CUD=[n], C intersect D={} and f(D)=D
#Try:
#RandCD(5,2);
RandCD:=proc(n,p) local ra,i,C,D,ra1,ra2,f,c:
ra:=rand(0..1):

C:={}:
D:={}:
for i from 1 to n do
if ra()=0 then
 C:=C union {i}:
else
 D:=D union {i}:
fi:
od:

ra1:=rand(1..n):
ra2:=randperm(D):
c:=0:
f:=[]:

for i from 1 to n do
if member(i,C) then
 f:=[op(f),ra1()]:
else
c:=c+1:
f:=[op(f),ra2[c]]:
fi:
od:

for i from n+1 to n+p do
f:=[op(f),ra2[rand(1..nops(ra2))()]]:
od:

[C,D,f]:

end:








#TriplesAB(f,p): inputs a function f from [n+p] to [n] (where n:=nops(f)-p) and outputs the set of triples
#[A,B,f] 
#TriplesAB([1,2,3],0);
TriplesAB:=proc(f,p) local n,S,A,T,i,gu,a,b,B:

n:=nops(f)-p:
T:={seq(i,i=1..n)}:


S:=powerset(n):


gu:={}:
for A in S  do
 B:=T minus A:
 if {seq(f[a], a in A)} minus A={} and
      {seq(f[b], b in B union {seq(i,i=n+1..n+p)})} minus B={}  then
   gu:=gu union {[A,B,f]}:
 fi:
od:
gu:


end:


#TriplesCD(f,p): inputs a function f from [n+p] to [n] (where n:=nops(f)=p) and outputs the set of triples [C,D,f] such
#that C and D are disjoint, C union D is [n]
#and such that f(D union {n+1..n+p})=D. Try
#TriplesCD([1,2,3],0);
TriplesCD:=proc(f,p) local n,S,A,T,i,gu,b,B:

n:=nops(f)-p:
T:={seq(i,i=1..n)}:

S:=powerset(n):


gu:={}:
for A in S  do
 B:=T minus A:

 if {seq(f[b], b in B union {seq(i,i=n+1..n+p)} )}=B  then
   gu:=gu union {[A,B,f]}:
 fi:
od:
gu:


end:

#AF(m,n): All the functions from [m] to [n]. Try
#AF(4,3);
AF:=proc(m,n) local gu,gu1,i:
option remember:
if m=0 then
RETURN({[]}):
fi:
gu:=AF(m-1,n):
{seq(seq([op(gu1),i],gu1 in gu),i=1..n)}:

end:

#CheckCO(n,p): checks the crucial observation for all functions from [n+p] to [n]. Try:
#CheckCO(4,0);
CheckCO:=proc(n,p) local gu,f:
gu:=AF(n+p,n):
evalb({ seq(nops(TriplesAB(f,p))-nops(TriplesCD(f,p)), f in gu)}={0}):
end:

#RF(m,n): a random function from [m] to [n]. Try
#RF(7,4);
RF:=proc(m,n) local ra,i:
ra:=rand(1..n):
[seq(ra(),i=1..m)]:
end:



#CheckCOrand(n,p,K): checks the crucial observation for K random functions from [n+p] to [n]. Try:
#CheckCOrand(40,0,10);
CheckCOrand:=proc(n,p,K) local i,gu,f:
gu:={seq(RF(n+p,n),i=1..K)}:
evalb({ seq(nops(TriplesAB(f,p))-nops(TriplesCD(f,p)), f in gu)}={0}):
end:



#Orb(f,i): inputs a function f from [n] to [n] and outputs the orbit that starts at i, try:
#Orb([1,1,1,1],2);
Orb:=proc(f,i) local gu:
gu:=[i]:

while not member(f[gu[nops(gu)]],{op(gu)}) do
 gu:=[op(gu),f[gu[nops(gu)]]]:
od:
[op(gu),f[gu[-1]]]:

end:

#Cyc(f,i): The cycle that ends the orbit starting at i. try:
#Cyc([4,1,3,2,5],3);
Cyc:=proc(f,i) local gu,j,akha:
gu:=Orb(f,i):

akha:=gu[nops(gu)]:

for j from 1 to nops(gu) while gu[nops(gu)-j]<>akha do od:
[op(nops(gu)-j..nops(gu),gu)]:
end:



#GK(TL,p): Given a left-triple [A,B,f] where f is a map from [n+p] to [n], outputs the correponding right-triple [C,D,f] using the Gil Kalai mapping. Try:
#TL:=RandAB(10,0); GK(TL);
GK:=proc(TL,p) local B,f,C,D,b,n,i:
B:=TL[2]: f:=TL[3]:
n:=nops(f)-p:

D:={seq(op(Cyc(f,b)), b in B)}:

D:=D union {seq(op(Orb(f,i)),i=n+1..n+p)} minus {seq(i,i=n+1..n+p)}:

C:={seq(i, i=1..n)} minus D:
[C,D,f]:
end:



#iGK(TR,p): Given a right-triple TR=[C,D,f] where f mapls [n+p] to [n], outputs the correponding left-triple [A,B,f] using the Gil Kalai inverse mapping. Try:
#TR:=RandCD(10,0); iGK(TR,0);
iGK:=proc(TR,p) local A,B,f,D,n,i:
D:=TR[2]: f:=TR[3]:
n:=nops(f)-p:
B:={seq(f[i],i=n+1..n+p)}:
for i from 1 to n do
 if {op(Cyc(f,i))} minus D={} then
    B:=B union {i}:
 fi:
od:
A:={seq(i,i=1..n)} minus B:
[A,B,f]:
end:


#CheckGK1(f,p): Check mapping GK for the mapping f from [n+p] to [n] where n:=nops(f)-p. Try:
#f:=RF(11,10): CheckGK1(f,1);
CheckGK1:=proc(f,p) local Gu1,Gu2, gu1, gu2:
Gu1:=TriplesAB(f,p):
Gu2:=TriplesCD(f,p):



if not evalb({seq(iGK(gu2,p),gu2 in Gu2)}=Gu1) then
 RETURN(false):
fi:




if not evalb({seq(GK(gu1,p),gu1 in Gu1)}=Gu2) then
 RETURN(false):
fi:



if {seq(evalb(iGK(GK(gu1,p),p)=gu1), gu1 in Gu1)}<>{true} then
 RETURN(false):
fi:



if {seq(evalb(GK(iGK(gu2,p),p)=gu2), gu2 in Gu2)}<>{true} then
 RETURN(false):
fi:

true:

end:

#CheckGK(n,p): checks GK for all functions f from  [n+p] to [n]. try:
#CheckGK(5,0);
CheckGK:=proc(n,p) local gu,f:
gu:=AF(n+p,n):
for f in gu do
 if not CheckGK1(f,p) then
   lprint(f):
   RETURN(false):
 fi:
od:
true:
end:



#CheckGKrand(n,p,K): checks GK forK random functions functions f from  [n] to [n]. try:
#CheckGKrand(10,2,10);
CheckGKrand:=proc(n,p,K) local gu,f,i:
gu:={seq(RF(n+p,n),i=1..K)}:

for f in gu do
 if not CheckGK1(f,p) then
   RETURN(false):
 fi:
od:
true:
end:



B1:=proc(n,p) local k: add(binomial(n,k)*k^k*(n-k)^(n-k+p),k=0..n):end:
B2:=proc(n,p) local k: add(binomial(n,k)*n^k*(n-k)!*combinat[stirling2](p+n-k,n-k),k=0..n):end: