##########################################################################

## GAP functions
Read("C:/texts/papers/dessins/GAP/braid.txt");
S := FreeGroup("s1", "s2", "s3");;  # the braid group B_4
s1 := S.1;; s2 := S.2;; s3 := S.3;;
\1 := One(S);;
tau := s2^-1*s3*s1^-1*s2;;  # The twist of two imaginary branches
G := FreeGroup("a1", "a2", "a3", "a4");;
a1 := G.1;; a2 := G.2;; a3 := G.3;; a4 := G.4;;
rel_slope := a1*a2*a3*a4/(a3*a4*a1*a2);
## Input data; to be initialized before calling "init_group"
beta  := [];;      # The braids $\beta_i$
gamma := [];;      # The braids $\gamma_i$
braid := [];;      # The monodromies $\mu_i$
N := 0;;           # The reference fiber $x_N$
K := 0;;           # The improper fiber $x_K$
kappa := One(G);;  # The slope
## Output (for further manual analysis)
Paths := [];;  # Will contain the paths from $x_N$ to $x_i$ (in paths[i+1])
Gens  := [];;  # Will contain the generators in $x_i$ in terms of $x_N$
Rels  := [];;  # Will contain the relations at $x_i$ in terms of $x_N$
Slope := One(G);;  # The slope at $x_N$
Infty := [];;      # The relation at infinity
dPlane := \1;      # The braid about $\{\Im x>0\}$
rPlane := [];;     # The relations about $\{\Im x>0\}$
g := [];;          # The group $\pi_1$ to be computed
h := [];;          # The commutant $[\pi_1,\pi_1]$
gens := [];;       # Generators of g
rels := [];;       # Relations of g

## A home made cyclic word reduction
reduce := function(g)
 local r, r2, try_factor, a;
 try_factor := function()
  for a in [a1, a2, a3, a4, a1^-1, a2^-1, a3^-1, a4^-1] do
   r2 := a*r*a^-1;
   if Length(r2) < Length(r) then r := r2; return true; fi;
#   if r2 < r then r := r2; return true; fi;
  od;
  return false;
 end;
 if IsList(g) then return List(g, e -> reduce(e)); fi;
 r := g;
 while try_factor() do ; od;
 return r;
end;

## Create g and h using the data computed
create_group := function()
 g := G/Union(Union(Rels), Infty, rPlane);;
 if ValueOption("Simplified") = true then g := SimplifiedFpGroup(g); fi;
 gens := GeneratorsOfGroup(g);
 rels := RelatorsOfFpGroup(g);
 Print(gens, ", ", Length(rels), " relations\n");
 Print(Sum(rels, g -> Length(g)), ": ", List(rels, g -> Length(g)), "\n");
 h := DerivedSubgroup(g);
 return g;
end;

## Translate a list from $x_i$ to $x_N$
translate := function(i, list)
 return Braid(Paths[i+1]^-1, G, list);
end;

## This function computes $\pi_1$ (in g) and its commutant (in h)
#   - "n" is the number of the reference fiber
#   - option "NoPlane": do not use the relation about $\{\Im x>0\}$
#   - option "Reduced": change to the cyclically reduced relations
#   - option "Simplified": call "SimplifiedFpGroup" before returning
init_group := function(n)
 local i, path_n, red, fl;

 red := function(g)
  if ValueOption("Reduced") = true then return reduce(g); fi;
  return g;
 end;
 fl := l -> Filtered(l, g -> g <> One(G));

 N := n;
# Initialize paths from $x_0$ to $x_i$
 Paths := [1..Length(beta)+1];
 Paths[1] := \1;
 for i in [1..Length(beta)] do
  Paths[i+1] := Paths[i]*gamma[i]*beta[i];
 od;
# Recompute the Paths from $x_N$ to $x_i$
 path_n := Paths[N+1];
 Paths := List(Paths, p -> path_n^-1*p);
 Gens  := List([1..Length(beta)], i -> translate(i, GeneratorsOfGroup(G)));
# Braid relations at $x_i$
 Rels := List(braid, b -> [fl(BraidRelations(b, G))[2]]);
# Modified braid relations at the improper fiber $x_K$
 Slope := One(G);
 if K > 0 then
  Slope := translate(K, kappa);
  Rels[K] := fl(List(GeneratorsOfGroup(G),
                      g -> Braid(braid[K], G, g)*kappa/(kappa*g)));
 fi;
# All relations recomputed at $x_N$
 Rels := List([1..Length(Rels)], i -> red(translate(i, Rels[i])));
# The relation at infinity recomputed at $x_N$
 Infty := [red((a1*a2*a3*a4)^2/Slope)];
# The relation $\pth_{r+1}=\id$ about $\partial\{\Im x>0\}$
 dPlane := Paths[Length(Paths)]*path_n*(s1*s2*s3)^4; rPlane := [];
 if ValueOption("NoPlane") <> true then
  rPlane := fl(red(BraidRelations(dPlane, G)));
 fi;
 return create_group();
end;

## Shift singular fibers by dn
shift := function(dn)
 local shift_index, shift_list;
 shift_index := i -> (i - 1 + dn) mod Length(beta) + 1;
 shift_list := l -> List([1..Length(beta)], i -> l[shift_index(i)]);
 beta  := shift_list(beta);
 gamma := shift_list(gamma);
 braid := shift_list(braid);
 if K <> 0 then K := (K - 1 - dn) mod Length(beta) + 1; fi;
end;

## These are from the book
qgroup := function(arg) return FactorGroupFpGroupByRels(g, arg); end;
qsize := function(arg) return Size(CallFuncList(qgroup, arg)); end;
P := [];;
simplify := function(arg)
 P := PresentationNormalClosure(g, Subgroup(g, arg));
 SimplifyPresentation(P);
end;
