#%Format: latex
#% This file is marked for TeX processing
#[caption]"Proof of \autoref{th.Alexander}"
#================================================================
#: Auxiliary functions and objects
#----------------------------------------------------------------
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^\perp$ if "transposed" option is specified
 if ValueOption("transposed") = true then return [v[2], -v[1]]; fi;
 return v;
end;
vectors := [   #: Canonical generators, see \eqref{eq.aT}
  function(N, p, t) return v\*([t^0, 0*t]); end,   #: type $\vI$
  function(N, p, t) return v\*([t^0, t+t^0]); end, #: type $\vII$
  function(N, p, t)                                #: type $\vIII_+$ or $\III_0$
   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]); fi;
   return fail;
  end,
  function(N, p, t)                                #: type $\vIII_-$
   if (p <> 3) and (N mod 3 = 0) then return v\*([-t, t^(e(N, p)/3)]); 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 $\bu_T$ as a list
 return List([1..Length(vectors)], i -> vectors[i](N, p, t));
end;
Region := function (N, p, t)   #"Region" #: Powers of $\xy$ (a region of width $N$)
 t := BurauPath((x*y)^-1, t);  return List([0..N-1], i -> t^i);
end;
Coeff := [5, 5, 4, 4, 3];      #"Coeff"  #: Coefficients for $k\suptext{all}$ in Euler's formula
#================================================================
#: "CreateSkeleton(N, p, psi, vector[: full])" creates the skeleton%
#:    $\skeleton$ of the universal subgroup $H\subset\BG3$ corresponding%
#:    to the triple $(N, p, \psi)$, starting from the generator "vector".%
#:    Also stored are "N", "p", "psi",%
#:    and the automorphism "U" such that $\skeleton/U$ is the%
#:    skeleton of the universal subgroup of $\Bu3$.%
#:    If "full" is specified, the extended skeleton $\Skt H$ is computed
#----------------------------------------------------------------
CreateSkeleton := function(N, p, psi, vector)   #"CreateSkeleton"
 local pp, i, G, H, U, 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 U in G then return Image(hom, U); fi;
  for i in G do if (vec^U)^i = vec then
   rt := RightTransversal(G, H);
   return PermList(List(rt, g -> PositionCanonical(rt, i*g)));
  fi; od;
  return Image(hom, U^3);
 end;
 pp := p;  if p = 0 then p := 7; fi;     #: A safe value
 pd := p^DegreeIndeterminate(psi, t);
 root := RootsOfUPol(GF(pd), PolynomialModP(psi, p))[1];
 U := 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, U^3); fi;
 hom := FactorCosetAction(G, H);
 return Skeleton(Index(G, H), Image(hom, m1), Image(hom, m2),
          "N", N, "p", pp, "psi", psi, "U", proj(), "vector", vector,
          "orbit", List(Vectors(N, p, root: transposed), IsFound),
          "expected", pd);
end;
#: Data to be printed
DataFields := ["N", "p", "psi", "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_.psi := [get_.N[1], "psi"];
get_.id := [get_.N[1], "id"];  #%
get_.vector := [get_.N[1], "vector"];
get_.orbit := [get_.N[1], "orbit"];
#================================================================
#: \hyperlink{step.1}{Step~$1$} of the algorithm (see \autopageref{ss.step.1}):%
#:   call \kwd{Pass1}.
#----------------------------------------------------------------
#: "tryPathList(N, p, list)" does the main job:%
#:   it tries a list of elements $g_1,\ldots,g_r\in\MG$ given in "list",%
#:   returning a list of pairs $(p,\mp_\xi)$ that fail the test.
#: The input values "N" and "p=0,2,3" are used for computing $M$ and detecting%
#:   allowed types $T\ne0$.
#----------------------------------------------------------------
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)  #: Add to the list all common factors as integral polynomials
  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()    #: Fails if $\CR=0$; otherwise, calls "addp" on all prime factors of $\CR$
  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)); #: $g_i\in\MG$ as Burau matrices
 vec := Filtered(Vectors(N, p, t), v -> v <> fail); #: admissible types as vectors
 reg := Region(N, p, t: normalized);                #: a trivial region
 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  #: the main loop
  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;   #: if one of the resultants is zero, return "fail"
    fi;
   od;
  od; od;
 od; od;
 return found;
end;
#----------------------------------------------------------------
#: "tryPathLists(N, ...)" calls \kwd{tryPathList} iteratedly on all primes $p=0,1,2$%
#:   and all sets $\{g_1,\ldots,g_r\}$ passed as arguments, computing intersection%
#:   over all sets and union over all primes.%
#:   The result is united with the master list passed as option "result=..." (if any).
#----------------------------------------------------------------
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          #: loop over $p=0,2,3$
  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  #: loop over all arguments
   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;             #: after one successful step, compute intersections
   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;
#[skiplines]
\[do\] := [];  #: just initialize the global
#[endskip]%
#----------------------------------------------------------------
#: "Pass1()" initializes all variables and calls \kwd{tryPathLists} on the hand tuned%
#:   lists contained in \GAPref{Burau\_start}.%
#:   For further processing, the pairs $(p,\mp_\xi)$ with $N\le5$ (computed manually)%
#:   from \GAPref{Burau\_fixed} are also added to the result, which is saved%
#:   to \GAPref{log/pass1}.%
#:   Output is logged to "log/pass1.log".
#----------------------------------------------------------------
Pass1 := function()   #"Pass1"
 local Res, rc;
 Res := Set([]);
 OutputLogTo("log/pass1.log");                      #: log output, just in case
 PushOptions(rec(result := Res));
 \[do\] := tryPathLists; #%
 Read("Burau_start.txt");  #: \hyperlink{step.1}{Step~$1$} of the algorithm
 \[do\] := function(N, p, psi, vec) AddSet(Res, [p, N, psi, [vec]]); end;
 Read("Burau_fixed.txt");                           #: add pairs with $N\le5$
 PopOptions();  OutputLogTo();
 Print(">>Total found: ", Length(Res), "\n");
 OpenOut("log/pass1.txt");                          #: save the result
 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;
#================================================================
#: \hyperlink{step.2}{Step~$2$} of the algorithm (see \autopageref{ss.step.2}):%
#:   call \kwd{Pass2}.
#----------------------------------------------------------------
#: "tryEuler(N, p, psi)" checks if $N<6+k\suptext{all}$, see \eqref{eq.Euler},%
#:   with an improved estimate on $k\suptext{all}$
#----------------------------------------------------------------
tryEuler := function(N, p, psi)   #"tryEuler"
 local vec, reg, pd, root, v, u, r, test;
 pd := p^DegreeIndeterminate(psi, t);
 root := RootsOfUPol(GF(pd), PolynomialModP(psi, p))[1];
 vec := Vectors(N, p, root); #%
 reg := Region(N, p, root);
 test := function(u) #: return "Coeff[u]" if "u" and "v" can share a region, or $0$ otherwise
  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(N, p, psi)" checks a candidate $(p,\mp_\xi)$ for genus zero
#----------------------------------------------------------------
tryGenus := function(N, p, psi)   #"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, psi], ", size=", S.Count, ", genus=", g,
        ", orbit=", get(S, "orbit"), "\n");
  S.info := SkeletonInfo(S);
  if (g = 0) and IsOutputTextStream(out) then            #: save if $g=0$
#: Output format:%
#:   $[N,p,\mp_\xi,"SkeletonInfo",\skeleton\ne\skeleton/U,\text{type $T$}]$ %
#:   \# $[\text{types}]$
   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, psi) then            #: \kwd{tryEuler} first
  Print(prefix, [N, p, psi], ", Euler <= 0\n");
  return;
 fi;
 list := ValueOption("vectors"); #%
 if list = fail then list := [1..5]; fi;
 for i in list do                                        #: try all types $T\ne0$,
  if found[i] then continue; fi;                         #: but avoid repetition
  SS := CreateSkeleton(N, p, psi, i);                    #: compute $\skeleton$
  found := UnionBlist(found, SS.orbit);
  SS.id := SS.U <> One(SS.U);  print(SS);  prefix := " ";        #: try $\skeleton$
  if SS.id and (g > 0) then print(FactorSkeleton(SS, SS.U)); fi; #: try $\skeleton/U$
 od;
end;
#----------------------------------------------------------------
#: "Pass2()" reads candidates $(p,\mp_\xi)$ from \GAPref{log/pass1}, tests them for%
#:   genus zero, and saves those passed to \GAPref{log/genus0}.%
#:   Output is logged to "log/pass2.log".
#----------------------------------------------------------------
Pass2 := function()  #"Pass2"
 \[do\] := tryGenus;
 OpenOut("log/genus0.txt"); #%
 OutputLogTo("log/pass2.log");
 Read("log/pass1.txt");
 OutputLogTo(); #%
 CloseStream(out);
end;
#================================================================
#[endinput]
