######################################################################
## JONAS: Save this file as  JONAS. To use it,                       #
## stay in the                                                       #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read TEN: <Enter>                                 #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger                                       #
## Rutgers University ,                                              #
## zeilberg at math dot rutgers dot edu.                             #
######################################################################

with(combinat):


print(`This is JONAS, A Maple package, written by `):
print(`Doron Zeilberger`):

print(`accompanying his  article: `):

print(` "Another Proof that Euler Missed: Jonas Sjostrand's Amazingly Simple (and Lovely!)`):
print(` Proof of the No-Longer-Amazing Loehr-Warrington Lattice Paths Conjecture"`):
print(` by D. Zeilberger`):
print(` published in the Personal Journal of Ekhad and Zeilberger`):
print(` http://www.math.rutgers.edu/~zeilberg/pj.html `):
print(`available from the author's website.`):
print():
print(`First Written: Oct. 18, 2005. Tested for Maple 10. `):
print(`Version of Oct. 18, 2005.  `):
print():

print():
print(`The most current version is available on WWW at:`):
print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/JONAS .`):
print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`):
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 specific help for these also type "ezra(procedure_name);" `):
print():
 


ezra1:=proc()
if args=NULL then

print(`ALL the  procedures are: Aab, AttachLap, Bdok23, Bdok32, CarLaps, `):
print(`  CheckJ12J21, CheckJ21J12, CheckJ23J32,CheckJ32J23 `):
print(` Corpus, ExtractCycle,ExtractEdge, , ExtractLap, ExtractPath `):
print(` Fniab, JAab2, J12, J21, J23, J32, DownGreedyPath , LowestWeldPoint`):
print(` LWwords , ShiftDG, Targem   `):
print(`   `):
else
 ezra(args):
fi:
end:

ezra:=proc()
if args=NULL then

print(` JONAS: A Maple package for implementing the beautfiul bijections`):
print(` in Jonas Sjostrand's lovely proof of  the `):
print(` Loehr-Warrington Conjceture, as described in his article. `):
print(`Cylindrical Lattice Paths and the Loehr-Warrington TEN to teh Power N Conjecture`):
print(` arXiv.math.CO/0509521 `):
print(`and more succintly in my own article:`):
print(` Another Proof that Euler Missed: Jonas Sjostrand's Amazingly Simple (and Lovely!)`):
print(` Proof of the No-Longer-Amazing Loehr-Warrington Lattice Paths Conjecture`):
print(``):
print(`The MAIN procedures are:  Aab, Bdok23, Bdok32, CarLaps, `):
print(`  CheckJ12J21, CheckJ21J12, CheckJ23J32,CheckJ32J23 `):
print(` J12, J21, J23, J32 `):
print(`   `):

elif nargs=1 and args[1]=Aab then
print(` Aab(a,b,n): All the Loehr-Warrington words in {a-,-b} of length n(a+b)`):
print(` For example, try: Aab(3,2,3); `):

elif nargs=1 and args[1]=AttachLap then
print(`AttachLap(DG,L,a,b): Given an (a,b)- Sjostrand directed graph DG and a lap L`):
print(` (a word in Aab(a,b,1))`):
print(`attaches the lap L to the bottom of DG`):
print(`For example try: AttachLap(J12([-1,1],1,1),[-1,1],1,1);`):

elif nargs=1 and args[1]=Bdok23 then
print(`Bdok23(a,b,n): Checks that the range of Sjostrand's map 2->3 is what it should be`):
print(`For example, try: Bdok23(3,2,3);`):

elif nargs=1 and args[1]=Bdok32 then
print(`Bdok32(a,b,n): Checks that the range of Sjostrand's map 3->2 is what it should be`):
print(`For example, try: Bdok32(3,2,3);`):

elif nargs=1 and args[1]=CarLaps then
print(`CarLaps(a,b,n): Aab(a,b,1)^n`):

elif nargs=1 and args[1]=CheckJ12J21 then
print(`CheckJ12J21(a,b,n): checks that J12(J21(w))=w for all Sjostrand directed graphs`):
print(`with n(a+b) edges`):
print(`in {a,-b}. For example, try: CheckJ12J21(3,2,3)`):

elif nargs=1 and args[1]=CheckJ21J12 then
print(`CheckJ21J12(a,b,n): checks that J21(J12(w))=w for all Loehr-Warrington words of length n(a+b)`):
print(`in {a,-b}. For example, try: CheckJ21J12(3,2,3)`):

elif nargs=1 and args[1]=CheckJ23J32 then
print(`CheckJ23J32(a,b,n): checks that J23(J32(LapV))=LapV for all (a,b)-Sjostrand  Lap vectors`):
print(`with n components`):
print(`try: CheckJ23J32(3,2,3)`):

elif nargs=1 and args[1]=CheckJ32J23 then
print(`CheckJ32J23(a,b,n): checks that J32(J23(LapV))=LapV for all (a,b)-Sjostrand  Dirrected graphs`):
print(`with n(a+b) edges `):
print(`try: CheckJ32J23(3,2,3)`):

elif nargs=1 and args[1]=CheckJ12J21 then
print(`CheckJ12J21(a,b,n): checks that J12(J21(w))=w for all Sjostrand directed graphs`):
print(`with n(a+b) edges`):
print(`in {a,-b}. For example, try: CheckJ12J21(3,2,3)`):

elif nargs=1 and args[1]=Corpus then
print(`Corpus(a,b,Len,Sum1): The sequence of LW-words of `):
print(` length (a+b)*i+Sum1 mod (a+b)`):
print(`for i=0..Len. For example, try`):
print(`Corpus(3,2,3,0);`):

elif nargs=1 and args[1]=ExtractCycle then
print(`ExtractCycle(DG,Loc1,a,b): the `):
print(`left-greedy cycle from Loc1 to itself`):
print(`of tranversing the directed graph `):
print(`followed by the striped directed graph`):
print(`DG=[V,RightNeighs,LeftNeighs] `):

elif nargs=1 and args[1]=ExtractEdge then
print(`ExtractEdge(DG,Loc,a,b):extracts a right-greedy edge from directed graph DG`):
print(`situated at location Loc for the (a,b)-walks`):
print(`returns the edge, and remaing directed graph, and the new location`):
print(`For example, try: ExtractEdge(J12([-1,-1,1,1]),[0,0],1,1);`):


elif nargs=1 and args[1]=ExtractLap then
print(`ExtractLap(DG,a,b):Given a directed graph DG, extracts a left-greedy lap`):
print(`and also returns the reduced directed graph`):
print(`For example, try: ExtractLap(J12([-1,-1,1,1],1,1),1,1);`):

elif nargs=1 and args[1]=ExtractPath then
print(`ExtractPath(DG,Loc1,Loc2,a,b): the `):
print(`left-greedy path from Loc1 to Loc2`):
print(`of tranversing the directed graph `):
print(`followed by the striped directed graph`):
print(`DG=[V,RightNeighs,LeftNeighs] `):

elif nargs=1 and args[1]=Fniab then
print(`Fniab(n,i,a,b): the set of walks described in terms of edges that end with the  edge`):
print(`[i-a,i]:`):

elif nargs=1 and args[1]=JAab2 then
print(`JAab2(a,b,n): the set of forward-degree and backward-degree sequences`):

elif nargs=1 and args[1]=J12 then
print(`J12(w,a,b): `):
print(`converts the good walk w to the directed graph in the form`):
print(`[SetOfVertices, RightNeighs, LeftNeighs]`):
print(`where RightNeighs[v] (resp. LeftNeighs[v]) is a table indicating how many`):
print(`right edges to the right (resp. to the left) come out of vertex v`):
print(`For example, try:J12([-1,-1,1,1],1,1);`):

elif nargs=1 and args[1]=J21 then
print(`J21(DG,a,b): `):
print(`The Sjostrand Map 2->1 from directed Eulerian graphs to words`):
print(`For example, try:J21(J12([-1,-1,1,1],1,1),1,1);`):

elif nargs=1 and args[1]=J23 then
print(`J23(DG,a,b):The Sjostrand map`):
print(`For example, try: J23(J12([-1,-1,1,1]),1,1);`):

elif nargs=1 and args[1]=J32 then
print(` J32(LapV,a,b):Jonas Sjostrand's 3->2 map`):
print(`For example, try: J32([[-1,1],[-1,1]],1,1);`):

elif nargs=1 and args[1]=DownGreedyPath then
print(`DownGreedyPath(DG,Loc,a,b): the down-greedy way, starting at vertex Loc, of tranversion `):
print(`the directed graph DG=[V,UpNeighs,DownNeighs]`):
print(`where V is the set of vertices, UpNeighs is a table such that`):
print(`UpNeighs[v] is the number of up-going edges from v `):
print(`and  DownNeighs is a table such that`):
print(`DownNeighs[v] is the number of down-going edges from v. `):
print(`For example, try: DownGreedyPath(J12([-1,-1,1,1]),[0,0],1,1);`):

elif nargs=1 and args[1]=LowestWeldPoint then
print(`LowestWeldPoint(DG):Given a directed graph DG, finds the lowest weldpoint`):
print(`For example, try: LowestWeldPoint(J12([-1,-1,1,1],1,1));`):

elif nargs=1 and args[1]=LWwords then
print(`LWwords(a,b,Len1,Sum1): all the words in {a,-b} `):
print(`of length Len1 that sum to Sum1 and`):
print(`avoid a factor of the form a[-a](-b)`):
print(`(Here [A] means a word that adds up to A)`):
print(`For example, try: LWwords(3,2,10,0); LWwords(1,1,6,-2);`):

elif nargs=1 and args[1]=ShiftDG then
print(`ShiftDG(DG,s): shifts the directed graph DG s units to the right`):
print(`For example, try: ShiftDG(J12([1,1,-1,-1]),1);`):

elif nargs=1 and args[1]=Targem then
print(` Targem(w,a,b): translates the word w into a word in directed edges `):


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

fi:


end:


#############################
Aab:=proc(a,b,n) Walks(a,b,b*n,a*n):end:

with(combinat):

IsBad:=proc(w,a)
local i:
for i from 1 to nops(w) do
 if w[i]=a and convert([op(i..nops(w),w)],`+`)=0 then
  RETURN(true):
 fi:
od:
false:
end:

#Walks with steps with m steps a and n steps -b
Walks:=proc(a,b,m,n)
local gu,gu1,gu2,i1,i,w:
option remember:
if n=0 then RETURN({[a$m]}) fi:

if m=0 then RETURN({[(-b)$n]}) fi:

gu1:=Walks(a,b,m-1,n):
gu:={seq([op(gu1[i1]),a],i1=1..nops(gu1))}:

gu2:=Walks(a,b,m,n-1):

for i from 1 to nops(gu2) do
w:=gu2[i]:
if not IsBad(w,a) then
  gu:=gu union {[op(w),-b]}:
fi:
od:
gu:
end:

#LWwords(a,b,Len1,Sum1): 
#all the words in {a,-b} of length Len1 that sum to Sum1 and
#avoid a factor of the form a[-a](-b)
LWwords:=proc(a,b,Len1,Sum1)
local gu,gu1,gu2,i1,i,w:
option remember:
if Len1=0 then 
 if Sum1=0 then
  RETURN({[]}): 
 else
   RETURN({}): 
 fi:
fi:

gu1:=LWwords(a,b,Len1-1,Sum1-a):
gu:={seq([op(gu1[i1]),a],i1=1..nops(gu1))}:

gu2:=LWwords(a,b,Len1-1,Sum1+b):

for i from 1 to nops(gu2) do
w:=gu2[i]:
if not IsBad(w,a) then
  gu:=gu union {[op(w),-b]}:
fi:
od:
gu:
end:

Aab:=proc(a,b,n):LWwords(a,b,(a+b)*n,0);end:

#Corpus(a,b,Len,Sum1): 
#The sequence of LW-words of length (a+b)*i+Sum1 mod (a+b)
#for i=0..Len. For example, try
#Corpus(3,2,3,0);
Corpus:=proc(a,b,Len,Sum1) local Sum1a,i:
Sum1a:= Sum1 mod (a+b):
[seq(LWwords(a,b,(a+b)*i+Sum1a,Sum1),i=0..Len)]:
end:


###New Stuff
ez:=proc():print(`JAab(a,b,n)`):end:

#Targem(w,a,b): translates the word w into a word in directed edges
Targem:=proc(w,a,b) local i:
[seq(
[
[(i-1) mod (a+b), convert([op(1..i-1,w)],`+`)], 
[i mod (a+b) , convert([op(1..i,w)],`+`)]
], i=1..nops(w))]:
end:

#TargemG(w,a,b,Start1): translates the word w into a word in directed edges
#that starts at Start1
TargemG:=proc(w,a,b,Start1) local i:
[seq(
[
[(i-1) mod (a+b), Start1+convert([op(1..i-1,w)],`+`)], 
[i mod (a+b) , Start1+convert([op(1..i,w)],`+`)]
], i=1..nops(w))]:
end:

#Translates all the words in S
TargemS:=proc(S,a,b) local i: {seq(Targem(S[i],a,b),i=1..nops(S))}:end:

JAab:=proc(a,b,n): TargemS(Aab(a,b,n),a,b):end:


#J12(w,a,b): converts the good walk w to the directed graph in the form
#SetOfVertices, RightNeighs, LeftNeighs
#where RightNeighs[v] (resp. LeftNeighs[v]) is a table indicating how many
#right edges to the right (resp. to the left) come out of vertex v
#For example, try:J12([-1,-1,1,1],1,1);
J12:=proc(w,a,b) 
local w1,S,i, tsela,  RightN, LeftN:
w1:=Targem(w,a,b):
S:={seq(op(w1[i]),i=1..nops(w1))}:

for i from 1 to nops(S) do
RightN[S[i]]:=0:
LeftN[S[i]]:=0:
od:

for i from 1 to nops(w1) do
tsela:=w1[i]:

if tsela[1][2]<tsela[2][2] then
 RightN[tsela[1]]:= RightN[tsela[1]]+1:
else
 LeftN[tsela[1]]:= LeftN[tsela[1]]+1:
fi:

od:
[S,op(RightN),op(LeftN)]:
end:



#J12G(w,a,b,Start1): converts the good walk w to the directed graph 
#that starts at [0,Start], in the form
#SetOfVertices, RightNeighs, LeftNeighs
#where RightNeighs[v] (resp. LeftNeighs[v]) is a table indicating how many
#right edges to the right (resp. to the left) come out of vertex v
#For example, try:J12([-1,-1,1,1],1,1);
J12G:=proc(w,a,b,Start1) 
local w1,S,i, tsela,  RightN, LeftN:
w1:=TargemG(w,a,b,Start1):
S:={seq(op(w1[i]),i=1..nops(w1))}:

for i from 1 to nops(S) do
RightN[S[i]]:=0:
LeftN[S[i]]:=0:
od:

for i from 1 to nops(w1) do
tsela:=w1[i]:

if tsela[1][2]<tsela[2][2] then
 RightN[tsela[1]]:= RightN[tsela[1]]+1:
else
 LeftN[tsela[1]]:= LeftN[tsela[1]]+1:
fi:

od:
[S,op(RightN),op(LeftN)]:
end:



#LeftGreedyPathOld(DG,Loc,a,b): the left-greedy way of tranversing the directed graph 
#DG=[V,RightNeighs,LeftNeighs] starting at location Loc
LeftGreedyPathOld:=proc(DG,Loc,a,b) local V,RightN,LeftN, orekh, Loc1,i,DG1,w1,RightN1,LeftN1:
V:=DG[1]:
RightN:=DG[2]:
LeftN:=DG[3]:

for i from 1 to nops(V) do
RightN1[V[i]]:=RightN[V[i]]:
LeftN1[V[i]]:=LeftN[V[i]]:
od:


orekh:=add(RightN[V[i]],i=1..nops(V))+add(LeftN[V[i]],i=1..nops(V)):

if orekh=0 then
 RETURN([]):
fi:

if LeftN[Loc]>0 then
 LeftN1[Loc]:= LeftN1[Loc]-1:
 Loc1:=[(Loc[1]+1) mod (a+b),  Loc[2]-b]:
 DG1:=[V,op(RightN1),op(LeftN1)]:
 w1:=LeftGreedyPathOld(DG1,Loc1,a,b):
 RETURN([-b,op(w1)]):
fi:

 if RightN[Loc]=0 then
  print(`Got stuck`):
  RETURN(FAIL):
 fi:

 RightN1[Loc]:= RightN1[Loc]-1:
 Loc1:=[(Loc[1]+1) mod (a+b), Loc[2]+a]:
 DG1:=[V,op(RightN1),op(LeftN1)]:
 w1:=LeftGreedyPathOld(DG1,Loc1,a,b):
 RETURN([a,op(w1)]):
end:


#ExtractEdge(DG,Loc,a,b):extracts a left-greedy edge from directed graph DG
#situated at location Loc for the (a,b)-walks
#returns the edge, and remaing directed graph, and the new location
#For example, try: ExtractEdge(J12[-1,-1,1,1]),[0,0],1,1);
ExtractEdge:=proc(DG,Loc,a,b) local V,RightN,LeftN,RightN1,LeftN1,DG1,i,Loc1,orekh:

V:=DG[1]:
RightN:=DG[2]:
LeftN:=DG[3]:

for i from 1 to nops(V) do
 RightN1[V[i]]:=RightN[V[i]]:
 LeftN1[V[i]]:=LeftN[V[i]]:
od:


orekh:=add(RightN[V[i]],i=1..nops(V))+add(LeftN[V[i]],i=1..nops(V)):

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


if LeftN[Loc]>0 then
 LeftN1[Loc]:= LeftN1[Loc]-1:
 Loc1:=[(Loc[1]+1) mod (a+b) , Loc[2]-b]:
 DG1:=[V,op(RightN1),op(LeftN1)]:
 RETURN([-b,DG1,Loc1]):
fi:

if RightN[Loc]=0 then
 print(`Got stuck`):
 RETURN(FAIL):
fi:

RightN1[Loc]:= RightN1[Loc]-1:
 Loc1:=[(Loc[1]+1) mod (a+b), Loc[2]+a]:
 DG1:=[V,op(RightN1),op(LeftN1)]:
 RETURN([a,DG1,Loc1]):
end:

#DownGreedyPath(DG,Loc,a,b): the 
#left-greedy way of tranversing the directed graph 
#DG=[V,RightNeighs,LeftNeighs] starting at location Loc
DownGreedyPath:=proc(DG,Loc,a,b)  local gu,i,lu,orekh:

orekh:=add(DG[2][DG[1][i]],i=1..nops(DG[1]) )+
add(DG[3][DG[1][i]],i=1..nops(DG[1])):
if orekh=0 then
 RETURN([]):
fi:

gu:=ExtractEdge(DG,Loc,a,b):
lu:=DownGreedyPath(gu[2],gu[3],a,b):
[gu[1],op(lu)]:

end:


#ExtractPath(DG,Loc1,Loc2,a,b): the 
#left-greedy path from Loc1 to Loc2
#of tranversing the directed graph 
#followed by the striped directed graph
#DG=[V,RightNeighs,LeftNeighs] 
ExtractPath:=proc(DG,Loc1,Loc2,a,b)  local gu,lu:

if Loc1=Loc2 then
 RETURN([],DG):
fi:

gu:=ExtractEdge(DG,Loc1,a,b):
lu:=ExtractPath(gu[2],gu[3],Loc2,a,b):
[gu[1],op(lu[1])],lu[2]:

end:


#ExtractCycle(DG,Loc1,a,b): the 
#left-greedy cycle from Loc1 to Loc1
#of tranversing the directed graph 
#followed by the striped directed graph
#DG=[V,RightNeighs,LeftNeighs] 
ExtractCycle:=proc(DG,Loc1,a,b)  local gu,lu:

gu:=ExtractEdge(DG,Loc1,a,b):
lu:=ExtractPath(gu[2],gu[3],Loc1,a,b):
[gu[1],op(lu[1])],lu[2]:

end:


#ExtractLap(DG,a,b):Given a directed graph DG, extracts a left-greedy lap
#For example, try: ExtractLap(J12([-1,-1,1,1],1,1),1,1);
ExtractLap:=proc(DG,a,b) :
ExtractCycle(DG,LowestWeldPoint(DG),a,b):end:


#J23(DG,a,b):The Sjostrand map
#For example, try: J23(J12([-1,-1,1,1]),1,1);
J23:=proc(DG,a,b) local V,RightN,LeftN,orekh,DG1,Cyc,jon1,i,gu:

V:=DG[1]:
RightN:=DG[2]:
LeftN:=DG[3]:


orekh:=add(RightN[V[i]],i=1..nops(V))+add(LeftN[V[i]],i=1..nops(V)):
if orekh=0 then
 RETURN([]):
fi:

gu:=ExtractLap(DG,a,b):
Cyc:=gu[1]:
DG1:=gu[2]:
jon1:=J23(DG1,a,b):
[op(jon1),Cyc]:
end:

#CarLaps(a,b,n): Aab(a,b)^n
CarLaps:=proc(a,b,n)
local gu1,gu2,i,j:
if n=0 then
 RETURN({[]}):
fi:
gu1:=Aab(a,b,1):
gu2:=CarLaps(a,b,n-1):
{seq(seq([gu1[i],op(gu2[j])],i=1..nops(gu1)),j=1..nops(gu2))}:
end:


J21:=proc(DG,a,b): DownGreedyPath(DG,[0,0],a,b):end:


#CheckJ21J12(a,b,n): checks that J21(J12(w))=w for all Loehr-Warrington words of length n(a+b)
#in {a,-b}. For example, try: CheckJ21J12(3,2,3)
CheckJ21J12:=proc(a,b,n) local i,gu:
gu:=Aab(a,b,n):
evalb({seq(evalb(J21(J12(gu[i],a,b),a,b)=gu[i]),i=1..nops(gu))}={true}):
end:



#AreWeTheSame(DG1,DG2): are Directed graphs DG1 and DG2 the same
AreWeTheSame:=proc(DG1,DG2)
evalb({evalb(DG1[1]=DG2[1]),evalb(op(DG1[2])=op(DG2[2])),evalb(op(DG1[3])=op(DG2[3]))}
={true}):
end:


#CheckJ12J21(a,b,n): checks that J12(J21(DG))=DG for all Sjostrand directed graphs
#in {a,-b} with n(a+b) edges . For example, try: CheckJ12J21(3,2,3)
CheckJ12J21:=proc(a,b,n) local i,gu:
gu:=Aab(a,b,n):
gu:={seq(J12(gu[i],a,b),i=1..nops(gu))}:
evalb({seq(AreWeTheSame(J12(J21(gu[i],a,b),a,b),gu[i]),i=1..nops(gu))}={true}):
end:


#CheckJ23J32(a,b,n): checks that J23(J32(LapV))=LapV for all LapTuples of size n
#in {a,-b}. For example, try: CheckJ23J32(3,2,3)
CheckJ23J32:=proc(a,b,n) local i,gu:
gu:=CarLaps(a,b,n):
evalb({seq(evalb(J23(J32(gu[i],a,b),a,b)=gu[i]),i=1..nops(gu))}={true}):
end:

#CheckJ32J23(a,b,n): checks that J32(J32(DG))=DG for all Sjostrand directed graphs
#in {a,-b} with n(a+b) edges . For example, try: CheckJ32J23(3,2,3)
CheckJ32J23:=proc(a,b,n) local i,gu:
gu:=Aab(a,b,n):
gu:={seq(J12(gu[i],a,b),i=1..nops(gu))}:
evalb({seq(AreWeTheSame(J32(J23(gu[i],a,b),a,b),gu[i]),i=1..nops(gu))}={true}):
end:



#LowestWeldPoint(DG):Given a directed graph DG, finds the lowest weldpoint
#For example, try: LowestWeldPoint(J12([-1,-1,1,1],1,1));
LowestWeldPoint:=proc(DG) local V1,V,i,aluf,shi:

V1:={}:
V:=DG[1]:

for i from  1 to nops(V) do

if DG[3][V[i]]+DG[2][V[i]]>0 and V[i][1]=0 then
 V1:=V1 union {V[i]}:
fi:

od:

aluf:=V1[1]:
shi:=V1[1][2]:
for i from 2 to nops(V1) do
 if V1[i][2]<shi then
  aluf:=V1[i]:
  shi:=V1[i][2]:
 fi:
od:
aluf:
end:





Bdok23:=proc(a,b,n)
local gu,i:
gu:=Aab(a,b,n):
evalb({seq(J23(J12(gu[i],a,b),a,b),i=1..nops(gu))}=CarLaps(a,b,n)):
end:

Bdok32:=proc(a,b,n)
local gu,i:
gu:=CarLaps(a,b,n):
evalb({seq(J21(J32(gu[i],a,b),a,b),i=1..nops(gu))}=Aab(a,b,n)):
end:



#AttachLap(DG,L,a,b): Given an (a,b)- Sjostrand directed graph DG and a lap L (a word in Aab(a,b,1))
#attaches the lap L to the bottom of DG
#For example try: AttachLap(J12([-1,1],1,1),[-1,1],1,1);
AttachLap:=proc(DG,L,a,b) local orekh,V,RightN,LeftN,i,Lp,DG1,j:
V:=DG[1]: RightN:=DG[2]: LeftN:=DG[3]:
orekh:=add(RightN[V[i]],i=1..nops(V))+add(LeftN[V[i]],i=1..nops(V)):

if orekh=0 then
RETURN(J12(L,a,b)):
fi:

Lp:=LowestWeldPoint(DG)[2]:

for j from 0 while J12G(L,a,b,Lp-j*(a+b))[1] intersect V<>{} do od:
Lp:=Lp-(j-1)*(a+b):

DG1:=J12G(L,a,b,Lp):

Combine1(DG,DG1): 

end:


#Combine1(DG1,DG2): combines the directed graphs DG1 and DG2
Combine1:=proc(DG1,DG2) local V1,V2,RightN1,RightN2,LeftN1,LeftN2,RightN,LeftN,V,i:

V1:=DG1[1]: RightN1:=DG1[2]: LeftN1:=DG1[3]:
V2:=DG2[1]: RightN2:=DG2[2]: LeftN2:=DG2[3]:

V:=V1 union V2:

for i from 1 to nops(V) do

if member(V[i],V1) and member(V[i],V2) then
 LeftN[V[i]]:=  LeftN1[V[i]]+LeftN2[V[i]]:
 RightN[V[i]]:=  RightN1[V[i]]+RightN2[V[i]]:
elif not member(V[i],V1) and member(V[i],V2) then
 LeftN[V[i]]:=  LeftN2[V[i]]:
 RightN[V[i]]:= RightN2[V[i]]:
elif  member(V[i],V1) and not member(V[i],V2) then
 LeftN[V[i]]:=  LeftN1[V[i]]:
 RightN[V[i]]:= RightN1[V[i]]:
 else
  ERROR(`Something is wrong`):
fi:
od:

[V,op(RightN),op(LeftN)]:
end:




#J32(LapV,a,b):Jonas Sjostrand's 3->2 map
J32:=proc(LapV,a,b) local LapV1,L1,DG1:

if nops(LapV)=0 then
 ERROR(`LapV must have at least one element`):
fi:

L1:=LapV[nops(LapV)]:

if nops(LapV)=1 then
 RETURN(J12(L1,a,b)):
fi:

LapV1:=[op(1..nops(LapV)-1,LapV)]:

DG1:=J32(LapV1,a,b):
AttachLap(DG1,L1,a,b):

end: