# F/out(A,B)[3] gives a finite presentation of the vertical
# amalgam decomposition (F9,F9;F81) of the index 4 subgroup
# Gamma0 of a locally quasi-primitive (10,10)-group Gamma
# described by the matrices A and B.

amE := function(A,B)
local i,j,cA,cB,S;
cA := DimensionsMat(A)[1];
cB := DimensionsMat(A)[2];
S := [];
for i in [1..cA] do
  for j in [1..cB] do
    Add(S,[i,A[i][j]+cA,B[i][j],j]);
  od;
od;
return S;
end;

spt := function(VT,ET,E,t)
local i;
if Size(VT) = t then 
  return ET;
else
  i := 1;
  while ((E[i][1] in VT and E[i][2] in VT) or
         (not E[i][1] in VT and not E[i][2] in VT)) do
    i := i+1;
  od;
  if E[i][1] in VT then
     return spt(Union([VT,[E[i][2]]]), Union([ET,[E[i]]]), E, t);
  else
    return spt(Union([VT,[E[i][1]]]), Union([ET,[E[i]]]), E, t);
  fi;
fi;
end;

dij := function(A,B)
local cA,cB,cvp,E,st,w,i,j,k,Q,l,mcvp,m;
cA := DimensionsMat(A)[1];
cB := DimensionsMat(A)[2];
cvp := [];
cvp[1] := [0,1,0,0];
E := amE(A,B);
st := spt([1],[],E,2*cA);
w := [];
for j in [1..Size(E)] do 
  if E[j] in st then 
    w[j] := 1;
  else
    w[j] := 1000;
  fi;
od;
for i in [2..2*cA] do
  cvp[i] := [1000000,i,0,0];
od;
Q := [1..2*cA];
while Size(Q) > 0 do
  mcvp := [];
  for l in Q do
    Add(mcvp,cvp[l]);
  od;
  m := Minimum(mcvp);
  RemoveSet(Q,m[2]);
  for k in [1..Size(E)] do
    if E[k][1] = m[2] then
      if cvp[E[k][2]][1] > m[1]+w[k] then
         cvp[E[k][2]][1] :=  m[1]+w[k];
         cvp[E[k][2]][3] :=  m[2];
         cvp[E[k][2]][4] := k;
      fi;
    fi;
    if E[k][2] = m[2] then
      if cvp[E[k][1]][1] > m[1]+w[k] then
         cvp[E[k][1]][1] :=  m[1]+w[k];
         cvp[E[k][1]][3] :=  m[2];
         cvp[E[k][1]][4] := k;
      fi;
    fi; 
  od;
od;
return cvp;
end;

ws := function(v,L,d,E,c)
local l,r;
if v = 1 then
  return L;
else
  if v = E[d[v][4]][1] then
    l := E[d[v][4]][3];
    r := E[d[v][4]][4];
  else
    l := c+1-E[d[v][4]][3];
    r := c+1-E[d[v][4]][4];
  fi;
  Add(L[1],l);
  Add(L[2],r);
  return ws(d[v][3],L,d,E,c);
fi;
end;

inv := function(L,c)
local i;
for i in [1..Size(L)] do
  L[i] := c+1-L[i];
od;
return L;
end;

gens := function(A,B)
local E,w,cA,cB,st,i,d,G;
cA := DimensionsMat(A)[1];
cB := DimensionsMat(A)[2];
w := [[],[]];
E := amE(A,B);
st := spt([1],[],E,2*cA);
G := Difference(E,st);
d := dij(A,B);
for i in [1..Size(G)] do
  w[1][i] := Concatenation(inv(Reversed(ws(G[i][1],[[],[]],d,E,cB)[1]), cB),
                           [G[i][3]],
                           ws(G[i][2],[[],[]],d,E,cB)[1]);
  w[2][i] := Concatenation(inv(Reversed(ws(G[i][1],[[],[]],d,E,cB)[2]), cB),
                           [G[i][4]],
                           ws(G[i][2],[[],[]],d,E,cB)[2]);
od;
return w;
end;

F := FreeGroup("s1","s2","s3","s4","s5","s6","s7","s8","s9",
               "t1","t2","t3","t4","t5","t6","t7","t8","t9");
s1 := F.1;;
s2 := F.2;;
s3 := F.3;;
s4 := F.4;;
s5 := F.5;;
s6 := F.6;;
s7 := F.7;;
s8 := F.8;;
s9 := F.9;;
t1 := F.10;;
t2 := F.11;;
t3 := F.12;;
t4 := F.13;;
t5 := F.14;;
t6 := F.15;;
t7 := F.16;;
t8 := F.17;;
t9 := F.18;;

U := FreeGroup("u1","u2","u3","u4","u5","u6","u7","u8","u9");
u1 := U.1;;
u2 := U.2;;
u3 := U.3;;
u4 := U.4;;
u5 := U.5;;
u6 := U.6;;
u7 := U.7;;
u8 := U.8;;
u9 := U.9;;

z1 := function(x,W)
local y;
if x = 1 then y := W[1]*W[1]^-1;
elif x = 2 then y := W[9]^-1;
elif x = 3 then y := W[8]^-1;
elif x = 4 then y := W[7]^-1;
elif x = 5 then y := W[6]^-1;
elif x = 6 then y := W[5]^-1;
elif x = 7 then y := W[4]^-1;
elif x = 8 then y := W[3]^-1;
elif x = 9 then y := W[2]^-1;
elif x = 10 then y := W[1]^-1;
fi;
return y;
end;

z2 := function(x,W)
return (z1(11-x,W))^-1;
end;

ns := function(L,W)
return z1(L[1],W)*z2(L[2],W);
end;

wd := function(L,W)
local w,i;
w := W[1]*W[1]^-1;
if Size(L) = 2 then
  return ns(L,W);
else
  for i in [1,3..Size(L)-1] do
    w := w*ns([L[i],L[i+1]],W);
  od;
  return w;
fi;
end;

out := function(A,B)
local i, g, h, S, T, ST;
g := gens(A,B)[1];
h := gens(A,B)[2];
S := [];
T := [];
ST := [];
for i in [1..Size(g)] do 
  S[i] := wd(g[i],[u1,u2,u3,u4,u5,u6,u7,u8,u9]);
  T[i] := wd(h[i],[u1,u2,u3,u4,u5,u6,u7,u8,u9]);
  ST[i] := wd(g[i],[s1,s2,s3,s4,s5,s6,s7,s8,s9])*
          (wd(h[i],[t1,t2,t3,t4,t5,t6,t7,t8,t9]))^-1;
od;
return [S,T,ST];
end;


# Example:
# A := [ 
# [ 4, 4, 3, 8, 1, 1, 9, 2, 6, 6 ], 
# [ 3, 5, 1, 6, 3, 4, 3, 4, 8,10 ], 
# [ 6,10, 4, 2, 6, 2, 4, 1, 9, 2 ], 
# [ 7, 7, 2, 3, 2, 5, 5, 3, 1, 1 ], 
# [10, 6, 8, 4, 4,10, 8, 6, 2, 9 ],
# [ 1, 1, 5, 7, 7, 3, 2, 9, 5, 3 ], 
# [ 8, 9,10, 9, 8, 6, 6,10, 4, 4 ], 
# [ 9, 2, 9, 5, 9, 7, 1, 5,10, 7 ], 
# [ 5, 3, 6, 1,10, 8, 7, 8, 7, 8 ], 
# [ 2, 8, 7,10, 5, 9,10, 7, 3, 5 ] ];

# B := [ 
# [ 8, 3, 2, 9, 7, 4, 1, 6,10, 5 ], 
# [10, 3, 5, 1, 6, 7, 8, 9, 2, 4 ], 
# [ 4, 7, 6, 3, 8, 5,10, 9, 2, 1 ], 
# [ 9,10, 2, 1, 4, 7, 6, 5, 8, 3 ], 
# [ 2, 3, 6, 5, 4, 1,10, 9, 8, 7 ], 
# [ 6, 1, 2, 5, 4, 3,10, 9, 8, 7 ], 
# [ 4, 3,10, 5, 8, 7, 6, 9, 1, 2 ], 
# [10, 9, 4, 1, 6, 3, 2, 5, 8, 7 ], 
# [ 4, 9, 2,10, 3, 5, 6, 7, 8, 1 ], 
# [ 7, 3, 2, 6,10, 8, 5, 1, 4, 9 ] ];