#%Format: latex
#% This file is marked for TeX processing
#[caption]"Computation for the proof of \autoref{th.Alexander}"
#================================================================
#: "CreateSkeleton(N, p, kappa, vector[: full])" creates the skeleton%
#:    $\skeleton$ of the universal subgroup of $\BG3$ corresponding%
#:    to the triple $(N, p, \kappa)$, starting from the generator "vector".%
#:    Also stored are "N", "p", "kappa",%
#:    and the automorphism "T" such that $\skeleton/T$ is the%
#:    skeleton of the universal subgroup of $\Bu3$
#----------------------------------------------------------------
Reread("SkPi1.txt");   #: (see \file{SkPi1})
e := function(N, p)    #: The $\e_p(N)$ function
 if (p = 2) or (N mod 4 = 0) then return N; fi;
 if N mod 2 = 0 then return N/2; fi;
 return N*2;
end;
v\* := function(v)     #: Convert $\bv$ to $\bv^*$ if "transposed" option is specified
 if ValueOption("transposed") = true then return [v[2], -v[1]]; fi;
 return v;
end;
vectors := [                                        #: Canonical generators:
  function(N, p, t) return v\*([0*t, t^0]); end,    #: type $\vI$
  function(N, p, t) return v\*([t + t^0, t]); end,  #: type $\vII$
  function(N, p, t)                                 #: type $\vIII_+$
   if (p = 3) and (N mod 3 <> 0) then return v\*([t^0, -t]); fi;
   if (p <> 3) and (N mod 3 = 0) then return v\*([t^(e(N, p)/3), t^0]); fi;
   return fail;
  end,
  function(N, p, t)                                 #: type $\vIII_-$
   if (p <> 3) and (N mod 3 = 0) then return v\*([-t^(e(N, p)/3), t^0]); fi;
   return fail;
  end,
  function(N, p, t)                                 #: type $\vIV$
   if e(N, p) mod 2 <> 0 then return v\*([t^((e(N, p) - 1)/2), t^0]); fi;
   return fail;
  end];
Vectors := function(N, p, t)   #"Vectors"  #: All generators
 return List([1..Length(vectors)], i -> vectors[i](N, p, t));
end;
Region := function (N, p, t)   #"Region" #: Powers of $\xy$
 t := BurauPath((x*y)^-1, t);  return List([0..N-1], i -> t^i);
end;
Coeff := [5, 5, 4, 4, 3];      #"Coeff"  #: Coefficients for Euler's formula
CreateSkeleton := function(N, p, kappa, vector)   #"CreateSkeleton"
 local pp, i, G, H, T, pd, root, vec, m, m1, m2, os, rt, hom, IsFound, proj;
 IsFound := function(v)
  if v = fail then return false; fi;
  for i in [0..pd-1] do if v^(m^i) in os.orbit then return true; fi; od;
  return false;
 end;
 proj := function()
  if T in G then return Image(hom, T); fi;
  for i in G do if (vec^T)^i = vec then
   rt := RightTransversal(G, H);
   return PermList(List(rt, g -> PositionCanonical(rt, i*g)));
  fi; od;
  return Image(hom, T^3);
 end;
 pp := p;  if p = 0 then p := 7; fi;     #: A safe value
 pd := p^DegreeIndeterminate(kappa, t);
 root := RootsOfUPol(GF(pd), PolynomialModP(kappa, p))[1];
 T := root*IdentityMat(2); #%
 m := Z(pd)*IdentityMat(2);
 m1 := BurauPath(x, root: transposed); #%
 m2 := BurauPath(y, root: transposed);
 G := Group(m1, m2); #%
 vec := vectors[vector](N, p, root: transposed);
 os := OrbitStabilizer(G, vec); #%
 H := os.stabilizer;
 if ValueOption("full") <> true then H := ClosureGroup(H, T^3); fi;
 hom := FactorCosetAction(G, H);
 return Skeleton(Index(G, H), Image(hom, m1), Image(hom, m2),
          "N", N, "p", pp, "kappa", kappa, "T", proj(), "vector", vector,
          "orbit", List(Vectors(N, p, root: transposed), IsFound),
          "expected", pd);
end;
#: Data to be printed
DataFields := ["N", "p", "kappa", "info", "id", "vector"];
UserData := function(S) return List(DataFields, i -> get(S, i)); end;
get_.N := [function(S, arg) return get(S.parent, arg[2]); end, "N"];
get_.p := [get_.N[1], "p"];  #%
get_.kappa := [get_.N[1], "kappa"];
get_.id := [get_.N[1], "id"];  #%
get_.vector := [get_.N[1], "vector"];
get_.orbit := [get_.N[1], "orbit"];
#================================================================
#: Pass 1 functions
#----------------------------------------------------------------
tryPathList := function(N, p, list)  #"tryPathList"
 local mat, vec, reg, cp, m1, m2, v1, v2, i, zero, addp, found, d, f, r, s;
 addp := function(pp)
  if N mod pp = 0 then return; fi;
  for f in Set(Factors(Gcd(PolynomialModP(cp, pp), PolynomialModP(d, pp)))) do
   AddSet(found, [pp, N, UnivariatePolynomial(Integers,
       List(CoefficientsOfUnivariatePolynomial(f), IntFFE), t)]);
  od;
 end;
 zero := function()
  d := DeterminantMatDestructive([(vec[v1]^mat[m1])^reg[i], vec[v2]^mat[m2]]);
  r := Resultant(cp, d, t);
  if p <> 0 then if r mod p = Zero(r) then addp(p); fi; return false; fi;
  if r = Zero(r) then return true; fi;
  for s in Set(Filtered(Factors(AbsInt(r)), s -> not(s in [1, 2, 3]))) do
   addp(s);
  od;
  return false;
 end;
 mat := List(list, m -> BurauPath(m, t: normalized));
 vec := Filtered(Vectors(N, p, t), v -> v <> fail);
 reg := Region(N, p, t: normalized);
 found := Set([]);
 cp := Value(CyclotomicPolynomial(Integers, e(N, p)),
             [Indeterminate(Integers)], [t]);
 for m1 in [1..Length(mat)] do for m2 in [1..Length(mat)] do
  for v1 in [1..Length(vec)] do for v2 in [1..Length(vec)] do
   for i in [1..Length(reg)] do
    if (i <> 1) or (m1 <> m2) or (v1 <> v2) then
     if zero() then return fail; fi;
    fi;
   od;
  od; od;
 od; od;
 return found;
end;
tryPathLists := function(arg)   #"tryPathLists"
 local N, i, p, isfound, found, res, Res;
 N := arg[1];  Res := Set([]);
 Print(">>N=", N, ":\n");
 for p in [0, 2, 3] do
  if (p <> 0) and (N mod p = 0) then continue; fi;
  found := Set([]); isfound := false;
  Print(" p=", p, ":\n");
  for i in [2..Length(arg)] do
   if arg[i] = [] then continue; fi;  Print("  ", arg[i], ": ");
   res := tryPathList(N, p, arg[i]);
   if res = fail then Print("fail\007\n"); continue; fi;
   if isfound then IntersectSet(found, res); else found := res; fi;
   isfound := true;
   Print("found=", Length(res), "/", Length(found), "\n");
   if found = [] then break; fi;
  od;
  if not isfound then Print("Failed\n"); return fail; fi;
  UniteSet(Res, found);
 od;
 Print("Total=", Length(Res), ": ", Set(List(Res, i -> i[1])), "\n");
 res := ValueOption("result");
 if IsSet(res) then UniteSet(res, Res); fi;
end;
\[do\] := [];
Pass1 := function()   #"Pass1"
 local Res, rc;
 Res := Set([]);
 OutputLogTo("log/pass1.log");
 PushOptions(rec(result := Res));
 \[do\] := tryPathLists; #%
 Read("Burau_start.txt");
 \[do\] := function(N, p, kappa, vec)
  AddSet(Res, [p, N, kappa, [vec]]);
 end;
 Read("Burau_fixed.txt");
 PopOptions();
 OutputLogTo();
 Print(">>Total found: ", Length(Res), "\n");
 OpenOut("log/pass1.txt");
 PrintTo(out, "\# Total found: ", Length(Res), "\n");
 for rc in Res do
  PrintTo(out, "\\[do\\](", rc[2], ", ", rc[1], ", ", rc[3]);
  if Length(rc) > 3 then PrintTo(out, ": vectors := ", rc[4]); fi;
  PrintTo(out, ");\n");
 od;
 CloseStream(out);
end;
#================================================================
#: Pass 2 functions
#----------------------------------------------------------------
tryEuler := function(N, p, kappa)   #"tryEuler"
 local vec, reg, pd, root, v, u, r, test;
 pd := p^DegreeIndeterminate(kappa, t);
 root := RootsOfUPol(GF(pd), PolynomialModP(kappa, p))[1];
 vec := Vectors(N, p, root); #%
 reg := Region(N, p, root);
 test := function(u)
  if vec[u] = fail then return 0; fi;
  for r in reg do
   if Determinant([v^r, vec[u]]) = Zero(root) then return Coeff[u]; fi;
  od;
  return 0;
 end;
 for v in vec do
  if v = fail then continue; fi;
  if N < 6 + Sum(List([1..5], test)) then return true; fi;
 od;
 return false;
end;
tryGenus := function(N, p, kappa)   #"tryGenus"
 local print, list, SS, g, found, i, prefix;
 print := function(S)
  g := Genus(S);
  if not IsBound(S.id) then S.id := fail; fi;
  Print(prefix, [N, p, kappa], ", size=", S.Count, ", genus=", g,
        ", orbit=", get(S, "orbit"), "\n");
  S.info := SkeletonInfo(S);
  if (g = 0) and IsOutputTextStream(out) then
   PrintTo(out, "\\[do\\](", UserData(S), "); \# ", S.orbit, "\n");
  fi;
 end;
 found := List(Vectors(N, p, t), i -> i = fail);  prefix := "*";
 if (N > 10) and not tryEuler(N, p, kappa) then
  Print(prefix, [N, p, kappa], ", Euler <= 0\n");
  return;
 fi;
 list := ValueOption("vectors"); #%
 if list = fail then list := [1..5]; fi;
 for i in list do
  if found[i] then continue; fi;
  SS := CreateSkeleton(N, p, kappa, i);
  found := UnionBlist(found, SS.orbit);
  SS.id := SS.T <> One(SS.T);  print(SS);  prefix := " ";
  if SS.id and (g > 0) then print(FactorSkeleton(SS, SS.T)); fi;
 od;
end;
Pass2 := function()  #"Pass2"
 \[do\] := tryGenus;
 OpenOut("log/genus0.txt"); #%
 OutputLogTo("log/pass2.log");
 Read("log/pass1.txt");
 OutputLogTo(); #%
 CloseStream(out);
end;
#================================================================
#: Further processing
#: "ReadSkeleton(...)" is intended to read a skeleton from "log/genus0.txt"
#----------------------------------------------------------------
ReadSkeleton := function(data)
 local S;
 S := CreateSkeleton(data[1], data[2], data[3], data[6]);
 S.info := data[4]; #%
 S.id := data[5]; #%
 S.orbit := data[6]; #%
 S.data := data;
 Print([data[1], data[2], data[3]], ", Count=", S.Count, "\n");
 return S;
end;
ListUnique := function(list)   #"ListUnique"
 local i, j, current, res;
 i := 1; #%
 res := [];
 while i < Length(list) do
  current := [list[i]];#%
  j := i + 1;
  while j <= Length(list) do
   if IsIsomorphic(list[i], list[j]) then
    Add(current, list[j]);#%
    Remove(list, j);#%
    continue;
   fi;
   j := j + 1;
  od;
  if Length(current) > 1 then Add(res, current); fi;
  i := i + 1;
 od;
 return res;
end;
AddUnique := function(list, S)
 local SS;
 for SS in list do if IsIsomorphic(S, SS) then return; fi; od;
 Add(list, S, 1);
end;
SplitProduct := function(S1, S2, r)
 local S, T, list, filter, p, pd, pp, s1, s2;
 filter := function(comp)
  s1 := ExtractComponent(S, comp);#%
  s2 := s1;
  if comp[1]^T in comp then
   p := MappingPermListList([1..Length(comp)], comp);
   s2 := FactorSkeleton(s1, p*RestrictedPerm(T, comp)*p^-1);
   s1.info := SkeletonInfo(s2);
  fi;
  if Genus(s2) > 0 then return fail; fi;
  s1.expected := pd;#%
  s1.p := pp;
  Unbind(s1.parent);  #: to save memory
  return s1;
 end;
 pd := S1.expected*S2.expected; #%
 pp := S1.p*S2.p;
 S := Prod(S1, S2); #%
 T := prod(S1.T, S1.Count, S2.T, S2.Count)^r;
 list := Filtered(List(Components(S), filter), s -> s <> fail);
 ListUnique(list);    #: To provide identification
 for p in [1..Length(list)] do list[p].data := [S1.data, S2.data, p, r]; od;
 return list;
end;
ReadSkeleton2 := function(data1, data2, iso1, iso2)
 return SplitProduct(ReadSkeleton(data1), ReadSkeleton(data2), 3: full);
end;
Isomorphisms := function()   #"Isomorphisms"
 local list, S, full;
 list := []; #%
 full := ValueOption("full") = true;
 \[do\] := function(data)
  S := ReadSkeleton(data);
  if full or (Genus(S) = 0) then Add(list, S); fi;
  if not full and (S.T <> One(S.T)) then Add(list, FactorSkeleton(S, S.T)); fi;
 end;
 Read("log/genus0.txt");
 return List(ListUnique(list), l -> List(l, UserData));
end;
Intersections := function()   #"Intersections"
 local list, S, i, j, k, res1, res2, iso1, iso2, process, temp;
 list := [];
 \[do\] := function(data)
  if [data[1], data[2]] = [1, 3] then return; fi;
  Add(list, ReadSkeleton(data: full));
 end;
 process := function(l)
  temp := List(l, u -> PositionProperty(list, v -> IsIsomorphic(u, v)) <> fail);
  for k in [1..Length(l)] do l[k] := l[k].info; od;
  return temp;
 end;
 Read("log/genus0.txt");
 OpenOut("log/intersection.txt");
 for i in [1..Length(list)] do
  Print(UserData(list[i]), "\n");
  for j in [i+1..Length(list)] do
   res1 := SplitProduct(list[i], list[j], 1);#%
   iso1 := process(res1);
   if res1 = [] then continue; fi;
   Print(">>", UserData(list[j]), "\n");
   res2 := SplitProduct(list[i], list[j], 3);#%
   iso2 := process(res2);
   Print("  ", res1, "\n", "  ", res2, "\n");
   PrintTo(out, "\\[do\\](", UserData(list[i]), ", ", UserData(list[j]), ", ",
           iso1, ", ", iso2, ");\n\# ", res1, "\n\# ", res2, "\n");
  od;
 od;
 CloseStream(out);
end;
Inclusions := function()
 local list, S, S1, S2, i, Exceptions, less, remove, print;
 list := []; #%
 Exceptions := [[7, 43], [7, 29], [9, 19]];
 less := function(S1, S2)
  i := Length(S1.sub) - Length(S2.sub);
  if i = 0 then i := S1.Count - S2.Count; fi;
  if i = 0 then i := S2.p - S1.p; fi;
  return i > 0;
 end;
 remove := function(S, list, l)
  if l = 0 then list := S.sup; fi;
  for S1 in list do
   if l > 0 then
    i := Position(S1.sub, S);
    if i <> fail then Remove(S1.sub, i); fi;
   fi;
   remove (S, S1.sup, l+1);
  od;
 end;
 print := function(S, l)
  if (l = 0) and IsBound(S.printed) then return; fi;
  S.printed := 1;
  i := Concatenation(String("", l), String(Position(list, S), -2),
                     String(SkeletonInfo(S), l-42));
  Print(i, "\n");#%
  PrintTo(out, i, " <- ", S.data, "\n");#%
  for S in S.sub do print(S, l+1); od;
 end;
 \[do\] := function(data1, data2, iso1, iso2)
  if [data2[1], data2[2]] in Exceptions then return; fi;
  S1 := ReadSkeleton(data1: full); S2 := ReadSkeleton(data2: full);
  AddUnique(list, S1); AddUnique(list, S2);
  for S in SplitProduct(S1, S2, 3) do AddUnique(list, S); od;
  for S in SplitProduct(S1, S2, 1) do AddUnique(list, S); od;
 end;
 Read("log/intersection.txt");
 Print("Read ", Length(list), " skeletons. Sorting...\n");
 for S in list do S.sub := []; S.sup := []; od;
 for S1 in list do for S2 in list do
  if S1.Count <> S2.Count then
   if Morphisms(S1, S2) <> [] then Add(S1.sub, S2); Add(S2.sup, S1); fi;
  fi;
 od; od;
 Sort(list, less); #%
 for S in list do Sort(S.sub, less); od;
 for S in list do remove(S, fail, 0); od;
 OpenOut("log/inclusion.txt");
 for S in list do print(S, 0); od;
 CloseStream(out);
end;
Intersections2 := function()   #"Intersections"
 local list1, list2, i, j, k;
 \[do\] := function(data1, data2, iso1, iso2)
  AddSet(list1, data1); AddSet(list1, data2);
  Add(list2, [data1, data2]); Add(list2, [data2, data1]);
 end;
 list1 := Set([]);  list2 := [];
 Read("log/intersection.txt");
 for i in [1..Length(list1)] do
  for j in [i+1..Length(list1)] do
   for k in [j+1..Length(list1)] do
    if ([list1[i], list1[j]] in list2) and
       ([list1[i], list1[k]] in list2) and
       ([list1[j], list1[k]] in list2) then
     Print([list1[i], list1[j], list1[k]], "\n");
    fi;
   od;
  od;
 od;
end;
#================================================================
#: Computation of $\pi_1$
#----------------------------------------------------------------
Skipped := function(S)
 local except;
 Unbind(S.info); #%
 Print("-> ", SkeletonInfo(S));
 except := ValueOption("except");
 if except <> fail then
  if SkeletonInfo(S) in except then Print(": skipped\n"); return true; fi;
 fi;
 Print("\n"); #%
 return false;
end;
analyzePi1 := function(S)
 local g, ab, size, expected, isinf;
 isinf := function(gr)
  g := gr;#%
  ab := AbelianInvariants(g);#%
  Print("  ", ab, ",");
  if 0 in ab then Print("  skipped\n"); return true; fi;
  return false;
 end;
 if Skipped(S) then return; fi;
 g := Pi1(S: simplified); #%
 expected := S.expected;
 if isinf(g) then return; fi;
 if isinf(DerivedSubgroup(g)) then return; fi;
 size := Size(g);#%
 Print("  ", size, "/", expected, ";");#%
 expected := size/expected;
 if expected <> 1 then
  Print("  >>");#%
  if isinf(DerivedSubgroup(g)) then return; fi;
  Print("  ", Size(g));#%
  if Size(g) <> expected then Print("  !!!"); fi;
 fi;
 Print("\n");
end;
Exceptions := Set([ #: Artificial exceptions:
 [ 200, 0, 0, "5^40", 20 ], # [10, 3], #: coset table calculation fails
 [ 840, 0, 0, "1^28 7^116", 28 ],  # [7, 29], #: too slow
 [ 1848, 0, 0, "1^42 7^258", 42 ], # [7, 43], #: too slow
 [ 300, 0, 0, "5^60", 300 ], #: coset table calculation fails
 [ 640, 0, 0, "8^80", 64 ],  #: coset table calculation fails
[]]);
Pi1s := function()   #"Pi1s"
 OutputLogTo("log/pi1.txt");
 \[do\] := function(data)
  if data[5] = fail then return; fi;
  analyzePi1(ReadSkeleton(data): full, except := Exceptions);
 end;
 Read("log/genus0.txt");
 \[do\] := function(data1, data2, iso1, iso2)
  local S, list, i;
  if not(false in iso2) then return; fi;
  list := ReadSkeleton2(data1, data2, iso1, iso2);
  for i in [1..Length(list)] do if iso2[i] then list[i] := fail; fi; od;
  list := Filtered(list, S -> S <> fail);
  Print("*Processing ", Length(list), " group(s)\n");
  for S in list do
   Print("*", SkeletonInfo(S), " ");#%
   analyzePi1(S: except := Exceptions);
  od;
 end;
 Read("log/intersection.txt");
 OutputLogTo();
end;

analyzePi1 := function(S)
 local res, g, gg, ab, absize, size, expected, isinf;
 isinf := function(gr)
  g := gr;#%
  ab := AbelianInvariants(g);#%
  Add(res, ab);
  if 0 in ab then Add(res, "skipped"); return true; fi;
  return false;
 end;
 gg := Pi1(S: simplified); #%
 expected := S.expected; #%
 res := [];
 if isinf(gg) then return res; fi;
 absize := Product(ab);
 if isinf(DerivedSubgroup(g)) then return res; fi;
 size := Size(gg)/absize;#%
 Add(res, size);#%
 Add(res, expected);#%
 expected := size/expected;
 if expected <> 1 then
  if isinf(DerivedSubgroup(g)) then return res; fi;
  Add(res, Size(g));#%
  if Size(g) <> expected then Add(res, "!!!"); fi;
 fi;
 if not IsString(res[Length(res)]) then Add(res, "OK"); fi;
 return res;
end;
Pi1s := function()   #"Pi1s"
 local Res, res, fname, i, j, list, shorten, process, done, save;
 Res := []; #%
 fname := "log/BurauPi1.txt";
 if IsReadableFile(fname) then Res := ReadAsFunction(fname)(); fi;
 shorten := function(data) return [data[1], data[2], data[3], data[6]]; end;
 done := function(ID)
  for j in Res do if j[1] = ID then return true; fi; od;
  return false;
 end;
 save := function()
  OpenOut(fname);
  PrintTo(out, "return [\n", Res[1]);
  for j in [2..Length(Res)] do PrintTo(out, ",\n", Res[j]); od;
  PrintTo(out, "\n];\n");
  CloseStream(out);
 end;
 process := function(ID, S)
  for i in Res do if i[1] = ID then return; fi; od;
  res := analyzePi1(S);
  Print("*", SkeletonInfo(S), ": ", res, "\n");
  Add(Res, [ID, res]); #%
  save();
 end;
 \[do\] := function(data1, data2, iso1, iso2)
  if not(false in iso2) then return; fi;
  if done([shorten(data1), shorten(data2)]) then return; fi;
  list := ReadSkeleton2(data1, data2, iso1, iso2);
  for i in [1..Length(list)] do if iso2[i] then list[i] := fail; fi; od;
  list := Filtered(list, S -> S <> fail);
  Print("*Processing ", Length(list), " group(s)\n");
  for i in [1..Length(list)] do
   process([shorten(data1), shorten(data2), i], list[i]);
  od;
  Add(Res, [[shorten(data1), shorten(data2)], "done"]);
  save();
 end;
 Read("log/intersection.txt");
 \[do\] := function(data)
  if data[5] = fail then return; fi;
  if done(shorten(data)) then return; fi;
  process(shorten(data), ReadSkeleton(data: full));
 end;
 Read("log/genus0.txt");
end;

#================================================================
#[endinput]
