##########################################################################
## BijTools: Save this file as  BijTools to use it,                      #
## stay in the same directory, get into Maple                            #
## (by typing: maple <Enter> )                                           #
## and then type:  read BijTools: <Enter>                                #
## Then follow the instructions given there                              #
##                                                                       #
## Written by Philip Matchett Wood, Rutgers University,                  #
##  matchett@math.rutgers.edu.                                           # 
##########################################################################

print(`First Written: Jan. 09, 2007: tested for Maple 10 `):
print(`Version of Jan. 09, 2007:  `):
print():
print(`This is BijTools, A Maple package`):
print(`accompanying the article: `):
print(`A translation method for finding combinatorial bijections`):
print(` by  Philip Matchett Wood  and Doron Zeilberger`):
print(`Available from:`):
print(`http://www.math.rutgers.edu/~matchett/Publications/`):

print( ):
print(`The most current version of this program is  available on WWW at:`):
print(`http://www.math.rutgers.edu/~matchett/Publications/`):
print(`Please report all bugs to: matchett at math dot rutgers dot edu .`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(`type "Help();". For specific help type "Help(procedure_name);" `):
print():


## NOTE: the Help function is defined at the end of this file.

# the function Hlp lists the main functions---less detailed than Help. 
Hlp:=proc():
print(`IndicesSeq(B), EntriesSeq(B), MultSet(S1,S2),`):
print(`IdBij(S), BuildBij(L1,L2), InvBij(B)`):
print(`ComposeBij(B1,B2), MultBij(B1,B2), AddBij(B1,B2), SumBij(L), SubtBij(B1,B2),`):
print(`BijOut(B,a,b,AnOrder), BijEqual(B1,B2), IsBij(B,Domain,Range),`):
print(`SplitProdBij(B1xB2), SwitchProdBij(B1xB2)`): 
end:

with(ListTools):

# HELPER FUNCTIONS

# IndicesSeq(B): returns the indices in a plain sequence,  instead of as a
# sequence of single-element lists.
IndicesSeq:=proc(B)
seq( op([indices(B)][i]) ,i=1..nops([indices(B)])):
end:

# EntriesSeq(B): returns the entries in a plain sequence,  instead of as a
# sequence of single-element lists.
EntriesSeq:=proc(B)
seq( op([entries(B)][i]) ,i=1..nops([entries(B)])):
end:

# MultSet(S1,S2): given two set S1 and S2, returns the Cartesian product 
#  $ S1 \\cross  S2 $
MultSet:=proc(S1,S2) local x,y,S1xS2:
S1xS2:={}:
for x in S1 do
	for y in S2 do
		S1xS2:=S1xS2 union {[x,y]}:
	od:
od:

return S1xS2:
end:

# BASIC OPERATIONS for BUILDING BIJECTIONS

# IdBij(S): creates the identity bijection from  S  to  S  (where S is a set).
IdBij:=proc(S) local x,B:
B:=table([seq(S[i] = S[i],i=1..nops(S))]):
B:
end:

# BuildBij(L1,L2):  creates the bijection from L1 to L2 (as sets) so that
# B(L1[i]) = L2[i].
BuildBij:=proc(L1,L2) local B:
# check that _sets_ have same size
if nops({op(L1)}) <> nops({op(L2)}) then
	print(`BuildBij: sets have different sizes.`):
	return FAIL:
fi:

B:=table([seq(L1[i] = L2[i],i=1..nops(L1))]):
B:
end:

# InvBij(B): takes a bijection B and inverts it.
InvBij:=proc(B) local Domain,Range,Bout:
Domain:=[IndicesSeq(B)]:
Range:=[EntriesSeq(B)]:

# NOTE: this depends on entries and indices giving the lists in the same
# order, which they are supposed to do.
Bout:=table([seq(Range[i] = Domain[i], i=1..nops(Range))]):

Bout:
end:



# OPERATIONS for COMBINING BIJECTIONS

# ComposeBij(B1,B2): returns the bijection B1 applied to B2.
ComposeBij:=proc(B1,B2) local B1ofB2, x, Domain:
Domain:={IndicesSeq(B2)}:

# check if the range of B2 is the domain of B1
if {IndicesSeq(B1)} <> {EntriesSeq(B2)} then 
	print(`ComposeBij: range of B2 not equal to domain of B1`):
	return FAIL:
fi:

B1ofB2:=table():
for x in Domain do
	B1ofB2[x]:=B1[B2[x]]:
od:

B1ofB2:
end:


# MultBij(B1,B2): creates the cartesian product bijection given two bijections.
## NOTE: this may need fixing to multiply 3 or more sets (avoid [[1,2],3] and
## such stuff).
# Further Note: don't multiply by zero, as it may screw things up.
MultBij:=proc(B1,B2) local i,j,B1xB2,I1,I2:

if (nops(B1) =0 or nops(B2) = 0) then
	print(`MultBij: cannot multiply by the empty set (zero elements)`):
	return FAIL:
fi:

I1:=[IndicesSeq(B1)]:
I2:=[IndicesSeq(B2)]:
#I1:=[seq( op([indices(B1)][i]) ,i=1..nops([indices(B1)]))]:
#I2:=[seq( op([indices(B2)][i]) ,i=1..nops([indices(B2)]))]:

B1xB2:=table(
	[seq(
		seq( [I1[i],I2[j]]=[B1[I1[i]],B2[I2[j]]] ,i=1..nops(I1))
	, j=1..nops(I2))]
):

B1xB2:
end:

# SplitProdBij(B1xB2): separates the two component bijections B1 and B2 of a
# bijection that is the product of B1 and B2.
# NOTE: this function checks to make sure that B1xB2 is indeed a product
# bijection to start with.
SplitProdBij:=proc(B1xB2) local p,x,y,i,j,B1,B2,S1,S2:

# Test if domain of B1xB2 is valid for a product bijection
for p in {IndicesSeq(B1xB2)} do
	if nops(p) <> 2 then
		print(`SplitProdBij: B1xB2 not a product bijection, bad domain!!`):
		return FAIL:
	fi:
od:
# Test if range of B1xB2 is valid for a product bijection
for p in {IndicesSeq(B1xB2)} do
	if nops(B1xB2[p]) <> 2 then
		print(`SplitProdBij: B1xB2 not a product bijection, bad range!!`):
		return FAIL:
	fi:
od:

# defining the domains for the component bijections
S1:={seq( op(indices(B1xB2)[i])[1],i=1..nops([indices(B1xB2)]) )}:
S2:={seq( op(indices(B1xB2)[j])[2],j=1..nops([indices(B1xB2)]) )}:

# Checking that the bijection B1xB2 is the same on [a,x] for every
# x and the same on [x,b] for every x as well. 
for x in S1 do 
	if nops( {seq(B1xB2[[x,y]][1],y in S2)} ) <> 1 then
		print(`SplitProdBij: B1xB2 not a product bijection---inconsistent!!`):
		return FAIL:
	fi:
od:
for y in S2 do 
	if nops( {seq(B1xB2[[x,y]][2],x in S1)} ) <> 1 then
		print(`SplitProdBij: B1xB2 not a product bijection---inconsistent!!`):
		return FAIL:
	fi:
od:

# Defining B1 and B2.
B1:=table(
	[seq(S1[i]=B1xB2[[S1[i],S2[1]]][1],i=1..nops(S1)) ]
):
B2:=table(
	[seq(S2[j]=B1xB2[[S1[1],S2[j]]][2],j=1..nops(S2)) ]
):

B1,B2:

end:


# SwitchProdBij(B1xB2): given a product bijection, say B1:A --> B and
# B2:C-->D, returns the corresponding bijection from A x C to D x B 
# (NOTE: the switch---the range of the returned bijection is _not_ B x D.)
# Also, note that, unlike SplitProdBij, we do _not_ need to check for
# "inconsistency".
SwitchProdBij:=proc(B1xB2) local p,B, B2xB1:

# Test if domain of B1xB2 is valid for a product bijection
for p in {IndicesSeq(B1xB2)} do
	if nops(p) <> 2 then
		print(`SplitProdBij: B1xB2 not a product bijection, bad domain!!`):
		return FAIL:
	fi:
od:
# Test if range of B1xB2 is valid for a product bijection
for p in {IndicesSeq(B1xB2)} do
	if nops(B1xB2[p]) <> 2 then
		print(`SplitProdBij: B1xB2 not a product bijection, bad range!!`):
		return FAIL:
	fi:
od:

B2xB1:= table([
	seq( p = [ B1xB2[p][2], B1xB2[p][1] ], p in [IndicesSeq(B1xB2)] )
]):

return B2xB1:
end:



# AddBij(B1,B2):  returns the bijection on the disjoint union.  If the sets
# (index sets) are not disjoint, returns FAIL.
AddBij:=proc(B1,B2) local i,j,B1plusB2:

if ({indices(B1)} intersect {indices(B2)} <> {}) or 
   ({entries(B1)} intersect {entries(B2)} <> {})  then 
	print(`AddBij: sets not disjiont!`):
	return FAIL:
fi:

B1plusB2:=table(
	[op(op(op(B1))),op(op(op(B2)))] 
):
B1plusB2:
end:


# SumBij(L):  returns the bijection on the disjoint union of a List of
# bijections.  If the sets
# (index sets) are not disjoint, returns FAIL.
SumBij:=proc(L) local i,j,SumBi,B:
SumBi:=table():
for B in L do
	SumBi:=AddBij(SumBi,B):
od:
return SumBi:
end:


# SubtBij(B1,B2)---use combinatorial inversion to find the bijection B1-B2...
# here $B1: A \to B$ and $B2: C \to D$ where $C \subset A$ and $D \subset B$.
SubtBij:=proc(B1,B2) local As,Bs,Cs,Ds,NewBij,P,a,b,x,B2Inv:
# check that $C \subset A$ and $D \subset B$
As:={IndicesSeq(B1)}:
Bs:={EntriesSeq(B1)}:
Cs:={IndicesSeq(B2)}: 
Ds:={EntriesSeq(B2)}:
if not( (Cs subset As) and (Ds subset Bs) ) then
	print(`Cannot subtract bijections---subset condition fails.`):
	return FAIL:
fi:

B2Inv:=InvBij(B2):
NewBij:=table():
for x in (As minus Cs) do 
	P:=[x,B1[x]]:
	while (P[nops(P)]) in Ds do
		a:= InvBij(B2)[P[nops(P)]]:
		b:= B1[a]:
		P:=[op(P),a,b]:
	od:
	NewBij[x]:=P[nops(P)]:
od:
NewBij:
end:


# OUTPUT function


# BijOut(Bijection,a,b,AnOrder)
# A general output function for bijecions that depend on a single parameter.
# Outputs the bijection for n=a to n=b.
# AnOrder is a total order relation defined on the domain of the bijection
#
BijOut:=proc(Bij,a,b,AnOrder) local x,B, i:

for i from a to b do
    print(`CASE n=`,i,`:`):
    B:=Bij(i):
    if (i <> 1 and i <>-1) then
        for x in sort([IndicesSeq(B)],AnOrder) do
            print(x,`  |--->  `,B[x])
        od:
    else
        for x in sort([IndicesSeq(B)]) do
            print(x,`  |--->  `,B[x])
        od:
    fi:
    print(`================================================`):
od:
end:

# EQUALITY CHECKING function

# BijEqual(B1,B2):  takes two bijections (i.e., tables in Maple) for arguments
#   returns  true   if bijections are equal
#   returns  false  otherwise.
BijEqual:=proc(B1,B2) local D1,D2,R1,R2, x:

D1:= {IndicesSeq(B1)}:
D2:= {IndicesSeq(B2)}:
R1:= {EntriesSeq(B1)}:
R2:= {EntriesSeq(B2)}:

# checking domain and range have the same size
if (nops(D1) <> nops(R1)) or (nops(D2)<>nops(R2)) then
	return false: # i.e., one of B1 or B2 is not a bijection
elif (D1<>D2) or (R1<>R2) then
	return false: # i.e., B1 and B2 do not have same domain and range
fi:

for x in D1 do # note D1=D2 at this point in the function
	if B1[x] <> B2[x] then 
		return false:
	fi:
od:

return true:
end:

# BIJECTION CHECKING function

# IsBij(B,Domain,Range):  takes one bijection (i.e., a table in Maple) and
# two sets (the domain and the range) for arguments
#   returns  true   B is a bijection from Domain to Range
#   returns  false  otherwise.
IsBij:=proc(B,Domain, Range) local D1,R1:

if nops(Domain)<> nops(Range) then
	print(`ERROR in IsBij: given domain and range have different`):
	print(`cardinalities, so a bijection is impossible.`):
	return ERROR:
fi:

D1:= {IndicesSeq(B)}:
R1:= {EntriesSeq(B)}:


# checking domain and range of B match the given sets Domain and Range
if (D1 <> Domain) or (R1 <> Range) then
	return false: # i.e., B does not map Domain to Range
fi:

return true:
end:


###############################################################################
########## Help function defined below # ######################################

Help:=proc()
if args=NULL then

print(`BijTools: A Maple package implementing the`):
print(`basic tools for using the translation method described in `):
print(`A translation method for finding combinatorial bijections`):
print(`by Philip Matchett Wood and Doron Zeilberger.`):
print(`The MAIN procedures are`):
print(`IndicesSeq, EntriesSeq`):
print(`IdBij, BuildBij, InvBij`):
print(`ComposeBij, MultBij, AddBij, SumBij, SubtBij,`):
print(`BijOut, BijEqual, IsBij,`):
print(`SplitProdBij, SwitchProdBij`): 
print():
print(`For help with a specific procedure, type Help(Procedure);`):
print(`For example, for help with BuildBij, type: Help(BuildBij);`):

elif nargs=1 and args[1]=IndicesSeq then
print(`IndicesSeq(B); outputs the domain of a bijection B as a sequence`):
print(`Note that a bijection is simply a table in Maple.`):

elif nargs=1 and args[1]=EntriesSeq then
print(`EntriesSeq(B); outputs the range of a bijection B as a sequence`):
print(`Note that a bijection is simply a table in Maple.`):

elif nargs=1 and args[1]=MultSet then
print(`MultSet(S1,S2): given two set S1 and S2, returns the Cartesian product`):
print(`$ S1 \\cross  S2 $.  For example, try typing S1:={1,2}; S2:={a,b};  `):
print(`followed by  S1xS2:=MultSet(S1,S2);`):


elif nargs=1 and args[1]=IdBij then
print(`IdBij(S): given a set S, returns the bijection sending x |-> x`):
print(` for every x in S`):
print(`For example, try  B:= IdBij({1,2,3}):`):
print(`Note that a bijection is just a table, so to see how B is defined`):
print(`try typing  op(B);  or   B[2];`):

elif nargs=1 and args[1]=BuildBij then
print(`BuildBij(L1,L2): given two lists L1 and L2, returns a bijection`):
print(`sending the k-th element of L1 to the k-th element of L2.`):
print(`For example, try  B:=BuildBij([a,b,c], [1,2,3]):`):
print(`Note that a bijection is just a table, so to see how B is defined`):
print(`try typing  op(B);  or   B[c];`):

elif nargs=1 and args[1]=InvBij then
print(`InvBij(B): given a bijection B returns the inverse bijection.`):
print(`For example, try  B:=BuildBij([a,b,c], [1,2,3]):  followed by`):
print(`Binv:=InvBij(B):`):
print(`Note that a bijection is just a table, so to see how Binv is defined`):
print(`try typing  op(Binv);  or   Binv[2];`):


elif nargs=1 and args[1]=ComposeBij then
print(`ComposeBij(B1,B2): given two bijections B1 and B2, this function`):
print(`checks that the range of B2 equals the domain of B1 and then`):
print(`returns the bijection B1 applied to B2.`):
print(`For example, try  B1:=BuildBij([1,2,3], [a,b,c]):  followed by`):
print(`B2:=BuildBij([a,b,c],[101,102,103]):  followed by`):
print(`B2ofB1:= ComposeBij(B2,B1):`):
print(`Note that a bijection is just a table, so to see how B2ofB1 is defined`):
print(`try typing  op(B2ofB1);  or   B2ofB1[2];`):

elif nargs=1 and args[1]=MultBij then
print(`MultBij(B1,B2): given two bijections B1: S1->T1 and B2: S2->T2,`):
print(`this function returns the natural bijection F: S1 x S2 -> T1 x T2`):
print(`where  x  denotes Cartesian product.`):
print(`For example, try  B1:=BuildBij([1,2,3], [a,b,c]):  followed by`):
print(`B2:=BuildBij([w,y,z],[23,25,26]):  followed by`):
print(`B1xB2:= MultBij(B1,B2):`):
print(`Note that a bijection is just a table, so to see how B1xB2 is defined`):
print(`try typing  op(B1xB2);  or   B1xB2[[2,z]];`):

elif nargs=1 and args[1]=AddBij then
print(`AddBij(B1,B2): given two bijections  B1: S1->T1 and B2: S2->T2,`):
print(`this function returns the natural bijection F: S1 union S2 -> T1 union T2`):
print(`where  union  denotes disjoint union.  If S1 intersects S2 or`):
print(`T1 intersects T2 is non-empty, then the function returns FAIL.`):
print(`For example, try  B1:=BuildBij([1,2,3], [a,b,c]):  followed by`):
print(`B2:=BuildBij([w,y,z],[23,25,26]):  followed by`):
print(`B1uB2:= AddBij(B1,B2):`):
print(`Note that a bijection is just a table, so to see how B1uB2 is defined`):
print(`try typing  op(B1uB2);  or   B1uB2[2];   or B1uB2[z]; `):

elif nargs=1 and args[1]=SumBij then
print(`SumBij(L): Given a list of bijections L, successively applies the`):
print(`function AddBij.  For example, if B1, B2, and B3 are bijections,`):
print(`then SumBij([B1,B2,B3]); returns AddBij(AddBij(B3,B2),B1).`):

elif nargs=1 and args[1]=SubtBij then
print(`SubtBij(B1,B2): given two bijections B1 and B2, this function`):
print(`uses combinatorial inversion to find the bijection B1-B2.`):
print(`here $B1: A --> B$ and $B2: C --> D$ where`):
print(`$C \\subset A$ and $D \\subset B$, and so SubtBij returns the`):
print(`bijection sending $A \\setminus C --> B \\setminus D$ given by `):
print(`the "alternating pairs algorithm".`):
print(`For example, try B1:=BuildBij([a1,a2,a3],[b1,b2,b3]):`):
print(`followed by  B2:=BuildBij([a1],[b3]):`):
print(`followed by  B1minusB2:=SubtBij(B1,B2):`):
print(`Note that a bijection is just a table, so to see how B1minusB2 is defined`):
print(`try typing  op(B1minusB2);  or   B1minusB2[a3];   or B1minusB2[a2]; `):


elif nargs=1 and args[1]=BijOut then
print(`BijOut(B,a,b,AnOrder): A general output function for bijecions that`):
print(`depend on a single parameter.  Outputs the bijection for n=a to n=b.`):
print(`AnOrder is a total order relation defined on the domain of the bijection.`):
print(`Note that B *must* be a Maple function of an integer.`):
print(`For example, type  Help(phi5rec)  and follow the instuructions`):
print(`after typing   read TransMethodZeck: `):

elif nargs=1 and args[1]=BijEqual then
print(`BijEqual(B1,B2): given two bijections B1 and B2, returns `):
print(`true if the bijections are the same and false otherwise.`):

elif nargs=1 and args[1]=IsBij then
print(`IsBij(B,Domain,Range): given a bijection B and two sets Domain and Range`):
print(`this function returns true if B is a bijection from Domain to Range`):
print(`and returns false otherwise`):



elif nargs=1 and args[1]=SplitProdBij then
print(`SplitProdBij(B1xB2): separates the two component bijections B1 and B2`): 
print(`of a bijection that is the product of B1 and B2.`):
print(`NOTE: this function checks to make sure that B1xB2 is indeed`):
print(`a product bijection to start with.`):

elif nargs=1 and args[1]=SwitchProdBij then
print(`SwitchProdBij(B1xB2): given a product bijection,`): 
print(`say B1:A --> B and  B2:C-->D, returns the corresponding bijection`): 
print(`from A x C to D x B`):
print(`NOTICE the switch: the range of the returned bijection is _not_ B x D.`):

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

fi:
end:


#EOF