// Intrinsics used by code from paper 12:
// _Computer aided discovery of a fast algorithm for testing conjugacy
//  in braid groups_, by Volker Gebhardt

// This is actually the file "braid.m" referred to in the paper,
// renamed to be consistent with the other papers.  There may be
// some cosmetic differences between the code here and the code
// shown in the paper.  Also, note that some extra intrinsics not
// mentioned in the paper are also included here.


// Section 2: Background: braid groups and testing conjugacy
// Subsection 2.2: Permutation braids and normal form

intrinsic MyNormalForm(x::GrpBrdElt) -> RngIntElt, SeqEnum
{ Return the infimum of x and a sequence of permutations defining the
  simple factors in the normal form of x. (x must be positive.) }

  require Id(Parent(x)) le x : "The argument must be a positive element.";
  delta := FundamentalElement(Parent(x));
  k := 0;
  seq := [];
  while not IsId(x) do
    d := GCD(x, delta);
    if d eq delta then
       k +:= 1;
     else
       Append(~seq, InducedPermutation(d)^-1);
    end if;
    x := d^-1 * x;
  end while;
  return k, seq;
end intrinsic;

// Subsection 2.3: Super summit sets and conjugacy testing

intrinsic MySuperSummitRepresentative(x::GrpBrdElt) -> GrpBrdElt, GrpBrdElt
{ Return a representative of the super summit set of x and an element
  conjugating x to this representative }

  n := NumberOfStrings(Parent(x));
  conj := Id(Parent(x));
  // maximise infimum
  count := n*(n - 1)/2 - 1;
  inf := Infimum(x);
  while count gt 0 and CanonicalLength(x) gt 0 do
    x, c := Cycle(x);  count -:= 1;
    conj := conj * c;
    if Infimum(x) gt inf then
      count := n*(n - 1)/2 - 1;
      inf := Infimum(x);
    end if;
  end while;
  // minimise supremum
  count := n*(n - 1)/2 - 1;
  sup := Supremum(x);
  while count gt 0 and CanonicalLength(x) gt 0 do
    x, c := Decycle(x);  count -:= 1;
    conj := conj * c;
    if Supremum(x) lt sup then
       count := n*(n - 1)/2 - 1;
       sup := Supremum(x);
    end if;
  end while;
  return x, conj;
end intrinsic;

intrinsic MyIsConjugate(x::GrpBrdElt, y::GrpBrdElt) -> BoolElt, RngIntElt, GrpBrdElt
{ Return whether x and y are conjugate, the number of super summit
  elements computed for the test, and an element conjugating x to y. }

  x, c_x := MySuperSummitRepresentative(x);
  S := {@ x @};
  conj := {@ c_x @};
  y, c_y := MySuperSummitRepresentative(y);
  if y eq x then
    return true, 1, c_x * c_y^-1;
  end if;
  // close S with respect to conjugation by simple elements
  pos := 1;
  while pos le #S do
    for s in { MinimalElementConjugatingToSuperSummit(S[pos], a)
               : a in Generators(Parent(x)) }
    do
      ns := LeftNormalForm(S[pos]^s);
      if ns notin S then
        Include(~S, ns);
        Include(~conj, LeftNormalForm(conj[pos]*s));
        if y eq ns then
          return true, #S, conj[pos] * s * c_y^-1;
        end if;
      end if;
    end for;
    pos +:= 1;
  end while;
  // S is closed with respect to conjugation by simple elements, that
  // is, S is the super summit set of x. Since S does not contain y,
  // the super summit sets of x and y are distinct.
  return false, #S, _;
end intrinsic;

// This next intrinsic (MyIsConjugate_Old) is not used in the paper
intrinsic MyIsConjugate_Old(x::GrpBrdElt, y::GrpBrdElt) -> BoolElt, RngIntElt
{ Return whether x and y are conjugate and the number of super summit
  elements computed for the test. }

  x := MySuperSummitRepresentative(x);
  S := {@ x @};
  y := MySuperSummitRepresentative(y);
  if y eq x then
    return true, 1;
  end if;
  // close S with respect to conjugation by simple elements
  pos := 1;
  while pos le #S do
    for s in { MinimalElementConjugatingToSuperSummit(S[pos], a)
               : a in Generators(Parent(x)) }
    do
      ns := LeftNormalForm(S[pos]^s);
      if ns notin S then
        Include(~S,ns);
        if y eq ns then
          return true, #S;
        end if;
      end if;
    end for;
    pos +:= 1;
  end while;
  // S is closed with respect to conjugation by simple elements, that
  // is, S is the super summit set of x. Since S does not contain y,
  // the super summit sets of x and y are distinct.
  return false, #S;
end intrinsic;


// Section 3: Coming across another class invariant
// Subsection 3.2: Identifying the smaller invariant

intrinsic Circuits(S::SetIndx) -> SetIndx
{ Given the super summit set S of an element x, return a set containing
  the vertex sets of the circuits of the graph induced on S by cycling. }

  C := {@ @};
  seen := [ false : s in S ];
  while exists(pos){ i : i in [1..#seen] | not seen[i] } do
    // follow trajectory until we arrive at previuosly seen element
    T := [ S | ];
    while not seen[pos] do
      seen[pos] := true;
      Append(~T, S[pos]);
      pos := Index(S, Cycle(S[pos]));
    end while;
    // if periodic part of trajectory is new, add it to C
    if S[pos] in T then
      Include(~C, {@ T[i] : i in [Index(T, S[pos])..#T] @});
    end if;
  end while;
  return C;
end intrinsic;

// Subsection 3.3: Looking for a way of computing ultra summit sets

intrinsic SatisfiesConvexity(S::SetIndx) -> BoolElt
{ Return whether a set S of braid group elements satisfies the "convexity"
  property.  S must be closed under conjugation with the fundamental element }

  B := Universe(S);
  D := [ B | s : s in Sym(NumberOfStrings(B)) ];  // simple elts.
  seen := [ false : i in [1..#S] ];
  // mark elements reachable from S[1] by a chain of elements
  //   of S linked by conjugation with simple elements
  new := { 1 };
  seen[1] := true;
  while #new gt 0 do
    ExtractRep(~new, ~i);
    u := S[i];
    for d in D do
      idx := Index(S, u^d);
      // u^d is in S if and only if idx is positive
      if idx gt 0 and not seen[idx] then
        Include(~new, idx);
        seen[idx] := true;
      end if;
    end for;
  end while;
  return forall{ i : i in [1..#S] | seen[i] };
end intrinsic;

intrinsic SatisfiesGCD(S::SetIndx) -> BoolElt
{ Return whether a set S of braid group elements satisfies the "gcd" property }

  B := Universe(S);
  D := [ B | s : s in Sym(NumberOfStrings(B)) ];  // simple elts.
  for s in S do
    for i := 1 to #D-1 do
      if s^D[i] in S then
        for j := i+1 to #D do
          if s^D[j] in S then
            if s^LeftGCD(D[i], D[j]) notin S then
              return false;
            end if;
          end if;
        end for;
      end if;
    end for;
  end for;
  return true;
end intrinsic;


// Section 4: On the way to a proof
// Subsection 4.1: Linking cycling and conjugation

intrinsic SimpleElementAdjacency(C::SetIndx) -> Mtrx
{ Given an indexed set of sets of elements of a braid group, return the
  the adjacency matrix for conjugation with simple elements }

  if #C eq 0 then
    return Matrix(Integers(), 0, 0, []);
  end if;
  B := Universe(C[1]);
  D := [ B | s : s in Sym(NumberOfStrings(B)) ];  // simple elts.
  A := Matrix(Integers(), #C, #C, [ 0 : i in [1..(#C)^2] ]);
  Ux := &join(C);
  for i := 1 to #C do
    for b in D do
      for p in C[i] do
        pb := p^b;
        if pb in Ux and exists(j){ k : k in [1..#C] | pb in C[k] } then
          A[i,j] +:= 1;
        end if;
      end for;
    end for;
  end for;
  return A;
end intrinsic;

intrinsic CheckTransport(x::GrpBrdElt, s::GrpBrdElt) -> RngIntElt
{ Return the number of simple elements t satisfying Cycle(x^s) eq Cycle(x)^t }

  B := Parent(x);
  D := [ B | s : s in Sym(NumberOfStrings(B)) ];  // simple elts.
  return #{ t : t in D | Cycle(x)^t eq Cycle(x^s) };
end intrinsic;

// Subsection 4.2: Defining the transport map

intrinsic MyTransport(x::GrpBrdElt, s::GrpBrdElt) -> GrpBrdElt
{ Return the transport of s along x -> Cycle(x).  x and x^s must be super
  summit elements }

  k := -Infimum(x);
  delta := FundamentalElement(Parent(x));
  lfx  := GCD(delta, delta^k*x);
  lfxs := GCD(delta, delta^k*(x^s));
  return (lfx^d)^-1 * s * lfxs^d where d := delta^k;
end intrinsic;

// Subsection 4.3: Proving properties of the transport map

intrinsic MovingFactors(x::GrpBrdElt, s::GrpBrdElt) -> SeqEnum
{ Given an element x = \delta^k A_1..A_r and a simple element s such that x
  and g^s are super summit elements, return a sequence [s_1,...,s_(r+1)]
  such that the i-th simple factor in the normal form of x^s is
  s_i^-1 A_i s_(i+1) }

  require IsSuperSummitRepresentative(x)
      and IsSuperSummitRepresentative(x^s) :
      "x and x^s should be super summit elements";
  B := Parent(x);
  delta := FundamentalElement(B);
  mov := [ s^d where d := delta^Infimum(x) ];
  xs := delta^(-Infimum(x)) * x^s;
  x := delta^(-Infimum(x)) * x;
  r := CanonicalLength(x);
  for i := 1 to r do
    Ai := GCD(delta, x);
    Asi := GCD(delta, xs);
    Append(~mov, NormalForm(Ai^-1 * mov[#mov] * Asi));
    x := Ai^-1 * x;
    xs := Asi^-1 * xs;
  end for;
  return mov;
end intrinsic;

intrinsic TestMovingFactors(x::GrpBrdElt) -> BoolElt
{ Return whether moving factors for conjugates of x satisfy implications
  for partial ordering and gcd }

  require IsSuperSummitRepresentative(x) :
          "x should be a super summit element";
  B := Parent(x);
  r := CanonicalLength(x);
  Dx := [ B | s : s in Sym(NumberOfStrings(B))
                | IsSuperSummitRepresentative(x^(B!s)) ];
  for s in Dx do
    for t in Dx do
      u := GCD(s, t);
      movs := MovingFactors(x, s);
      movt := MovingFactors(x, t);
      movu := MovingFactors(x, u);
      if exists{ i : i in [1..r+1] | movu[i] ne GCD(movs[i], movt[i]) } then
        return false;
      end if;
    end for;
  end for;
  return true;
end intrinsic;


// Section 5: Computing minimal simple elements

intrinsic Pullback(x::GrpBrdElt, s::GrpBrdElt) -> GrpBrdElt
{ Given a super summit element x and a simple element s, return the minimal
  simple element u such that x^u is a super summit element and s is a divisor
  of the transport of u from x to c(x). }

  B := Parent(x);
  require s in B : "Elements do not belong to a common braid group";
  LeftNormalForm(~x);
  require IsSuperSummitRepresentative(x) :
          "The first argument must be a super summit element";
  require IsSimple(s) : "The second argument must be simple";
  delta := FundamentalElement(B);
  tau := delta^Infimum(x);
  tau1 := tau*delta;
  cfp := [ B | cf : cf in CFP(x)[3] ];  // the simple factors
  A1s := cfp[1]^-1*delta;
  m := s^tau; // this should be moving into first simple factor
  // ...ensure that the product of the first simple factor of
  //    conjugate with m is simple
  u1 := (A1s^-1 * LCM(A1s, m))^t where t := tau1^-1;
  // ...and that m is a divisor of the second simple factor
  u2 := m;
  for i := 2 to #cfp do
    u2 := cfp[i]^-1 * LCM(cfp[i], u2);
  end for;
  u := MinimalElementConjugatingToSuperSummit(x, LeftLCM(u1, u2));
  return LeftNormalForm(u);
end intrinsic;


// The remaining intrinsics in this file are not part of the paper,
// but may be of interest to those wishing to experiment with Braid
// groups and ultra summit sets

intrinsic Orbit(g::GrpBrdElt, op::Intrinsic) -> SeqEnum, SeqEnum, RngIntElt
{ The orbit of g under the operation op, the conjugating elements along the
  orbit and the index of the first element in the periodic part. }

  conj := Id(Parent(g));
  Orb := [ g ];
  Conj := [ conj ];
  repeat
    g, c := op(g);
    conj *:= c;
    index := Index(Orb, g);
    Append(~Orb, g);
    Append(~Conj, conj);
  until index gt 0;
  return Orb, Conj, index;
end intrinsic;

intrinsic MyMinimalElementConjugatingToUltraSummit(x::GrpBrdElt, a::GrpBrdElt
     : Shortcut := false) -> GrpBrdElt, BoolElt
{ Given an ultra summit element x and a simple element a, return the minimal
  element c such that a le c and x^c is an ultra summit element.  If the flag
  Shortcut is set to true, the computation will be aborted if c is seen to
  be non-minimal; in this case, the function returns id, true. }

  B := Parent(x);
  require a in B : "Elements do not belong to a common braid group";
  require IsUltraSummitRepresentative(x) :
          "The first argument must be an ultra summit element";
  orb := Orbit(x, Cycle);
  require IsSimple(a) : "The second argument must be simple";
  a := MinimalElementConjugatingToSuperSummit(x, a);
  // compute periodic points under iterated transport along the
  // trajectory of x under cycling
  tran := {@ B | @};
  b := a;
  idx := 0;
  repeat
    Include(~tran, b);
    for i := 1 to #orb - 1 do
      b := Transport(orb[i], b);
      if IsId(b) then
        idx := -1;
        break;
      end if;
    end for;
    if idx eq 0 then
      idx := Index(tran, b);
    end if;
  until idx ne 0;
  if idx gt 0 then
    // case 1: if a stable point with divisor a exists, this is the
    //         desired element c
    if exists(f){ tran[i] : i in [idx..#tran] | a le tran[i] } then
      return f, false;
    end if;
    // case 2: if there is no stable point with divisor a and the iterated
    //         transport has non-trivial stable points, the desired element
    //         will not be minimal
    if Shortcut then
      return B!1, true;
    end if;
  end if;
  // case 3: pull a back along the trajectory of x until pullback becomes
  //         stable; a suitable stable point will bring us into case 1
  tran := {@ B | @};
  b := a;
  repeat
    Include(~tran, b);
    for i := #orb - 1 to 1 by -1 do
      b := Pullback(orb[i], b);
    end for;
    idx := Index(tran, b);
  until idx gt 0;
  l := #tran+1 - idx; // we want a multiple of the period length l...
  idx +:= idx mod l;  // ...this is it
  // compute stable points under iterated transport along the period
  b := tran[idx];
  tran := {@ B | @};
  repeat
    Include(~tran, b);
    for i := 1 to #orb - 1 do
      b := Transport(orb[i], b);
    end for;
    idx := Index(tran, b);
  until idx gt 0;
  if not exists(f){ tran[i] : i in [idx..#tran] | a le tran[i] } then
    assert false;  // this is a problem...
  end if;
  return f, false;
end intrinsic;

intrinsic MyUltraSummitRepresentative(x::GrpBrdElt) -> GrpBrdElt
{ Given an element x of a braid group, return an element of the ultra summit
  set of x. }

  x, c := SuperSummitRepresentative(x);
  orb, conj , i := Orbit(x, Cycle);
  return orb[i], c*conj[i];
end intrinsic;

intrinsic MyIsConjugate_Ultra(x::GrpBrdElt, y::GrpBrdElt) -> BoolElt, RngIntElt, GrpBrdElt
{ Return whether x and y are conjugate, the number of super summit
  elements computed for the test and an element conjugating x to y. }

  x, c_x := MyUltraSummitRepresentative(x);
  U := {@ x @};
  conj := {@ c_x @};
  y, c_y := MyUltraSummitRepresentative(y);
  if y eq x then
    return true, 1, c_x * c_y^-1;
  end if;
  // close U with respect to conjugation by simple elements
  pos := 1;
  while pos le #U do
    for s in { MyMinimalElementConjugatingToUltraSummit(U[pos], a)
               : a in Generators(Parent(x)) }
    do
      nu := LeftNormalForm(U[pos]^s);
      if nu notin U then
        Include(~U, nu);
        Include(~conj, LeftNormalForm(conj[pos]*s));
        if y eq nu then
          return true, #U, conj[pos] * s * c_y^-1;
        end if;
      end if;
    end for;
    pos +:= 1;
  end while;
  // U is closed with respect to conjugation by simple elements, that
  // is, U is the ultra summit set of x. Since U does not contain y,
  // the ultra summit sets of x and y are distinct.
  return false, #U, _;
end intrinsic;