/********************************************************************************************************************************
	This code was used to do the calculations for the main theorems of the paper "On rigid varieties isogenous to a product of 
  curves". The focus is on  The code was written in a way s.t. it can be easily adapted for other usecases. 
********************************************************************************************************************************/


// *******************************************************************************************************************************
// ************************************************** GENERAL USEFULL FUNCTIONS **************************************************
// *******************************************************************************************************************************


// Constants for directory/file naming:
_FOLDERNAME := "Computation_Bihol_Classes";
_CHI_SUBFOLDERNAME := "Absol_Faithfull_Chi";
_KERNELS_SUBFOLDERNAME := "With_Kernels";
_TRIANGLEGRPFILENAME := "TriangleGrpsInfo.txt";

/* 
  input:
    - str: (MonStgElt) The string to split.
    - delimiter: (MonStgElt) Single-character delimiter.
  output:
    (SeqEnum) Sequence of non-empty substrings obtained by splitting str by delimiter.
  description:
    Splits a string into parts separated by the given delimiter, ignoring empty parts.
*/ 
function SplitString(str, delimiter)
  parts := [];
  current := "";
  for i in [1..#str] do
    if str[i] eq delimiter then
      if #current gt 0 then
        Append(~parts, current);
      end if;
      current := ""; 
    else
      current cat:= str[i]; 
    end if;
  end for;
  if #current gt 0 then
    Append(~parts, current);
  end if;
  return parts;
end function;

/*
  input:
    - maxGrpOrder: (RngIntElt) Maximum group order checked.
    - triangleGroups: (SetEnum) Set of tuples <grpOrder, grpIdentifier> for groups with hyperbolic generating triples.
  description:
    Saves computed triangle groups to a file to reduce computational costs for future runs.
*/
procedure WriteTriangleGroups(maxGrpOrder, triangleGroups)
  filename := _FOLDERNAME cat "/" cat _TRIANGLEGRPFILENAME ;
  PrintFile(filename, IntegerToString(maxGrpOrder) cat "\n" cat &cat[ IntegerToString(tup[1]) cat "," cat IntegerToString(tup[2]) 
      cat ";" : tup in triangleGroups] : Overwrite := true);
// Besser mit Join arbeiten: Join([...], ", "); 
end procedure;


/* 
  output:
    (RngIntElt), (SetEnum) Maximum group order checked, and the set of tuples < grpOrder, grpIdentifier > identifying 
        the groups admitting hyperbolic generating triples
  description:
    Reads the saved list of triangle groups from file.
*/ 
ReadTriangleGroups := function()
  triangleGroups := {};
  try 
    file := Open(_FOLDERNAME cat "/" cat _TRIANGLEGRPFILENAME, "r");
    maxGrpOrderString := Gets(file);
    triangleGroupsString := Gets(file);
    maxGrpOrder := StringToInteger(maxGrpOrderString);
    triangleGroups := { < StringToInteger(num) : num in SplitString(tup, ",") > : tup in SplitString(triangleGroupsString, ";")};
  catch e // If there is a problem with the directory/files, inform the user in the console
    print("-- Probems reading the TriangleGroups File --");
    return 0, { };
  end try;
  return maxGrpOrder, triangleGroups;
end function;

/*
  input:
    - tup: (Tup/SeqEnum) Sequence to be converted.
  output:
    (SetEnum) The input sequence converted into a set.
  description:
    Converts a tuple/sequence into a set.
*/
ToSet := function(tup) 
  return { el : el in tup }; 
end function;

/* 
  input: 
    - tup : (Tup/SeqEnum/SetEnum) The collection.
    - elem : Element to count.
  output:
    (RngIntElt) Number of occurrences of elem in tup.
  description:
    Counts the number of appearances of a specific element.
*/
CntElem := function(tup, elem) 
  cnt := 0; 
  for i in tup do 
    if i eq elem then 
      cnt := cnt + 1; 
    end if; 
  end for; 
  return cnt; 
end function;

/*
  input:
    - T: (SeqEnum) Type of a generating triple.
  output:
    Rational number 1 - 1/m1 - 1/m2 - 1/m3.
  description:
    Theta value used to determine whether a generating triple yields a curve of genus > 2.
*/
Theta := function(T)
  t := 1 - 1/T[1] - 1/T[2] - 1/T[3];
  return t; 
end function;

/* 
  input: 
    - mset : (Tup/(Multi-)Set) A multiset, set, or sequence of integers.
  output: 
    (SeqEnum) Ordered sequence of the elements.
  description:
    Converts an unordered multiset or set into an ascending sorted sequence.
*/ 
OrderedSeq := function(mset)
  seq := [ ];
	while #mset ne 0 do 
    min := Minimum(mset);
		Append(~seq, min); 
    Exclude(~mset, min); 
	end while;
  return seq; 
end function;

/* 
  input: 
    - gOrd: (RngIntElt) Order of group G.
    - genus: (RngIntElt) Genus of quotient C/G.
    - d: (RngIntElt) Divisor parameter for type constraints.
  output: 
     (SetEnum) Set of possible generating triple types.
  description:
    Computes all triple types [m1, m2, m3] that satisfy the Riemann-Hurwitz formula for the given genus.
*/ 
ListOfTypes := function(gOrd,genus,d) 
  types := {}; 
  orders := [k : k in Divisors(Gcd(gOrd,d)) | 2 le k and k le 4*genus +2]; 
  for i1 in [1..#orders] do
    for i2 in [i1..#orders] do
      for i3 in [i2..#orders] do 
	      typ := [orders[i1],orders[i2], orders[i3]];
	      if 2*genus-2 eq gOrd*Theta(typ) then 
	        Include(~types,typ);
	      end if;
	    end for;
    end for;
  end for;
  return types;
end function;

/* 
  input: 
    - gOrd: (RngIntElt) Order of a group G.
    - genus: (RngIntElt) Genus of the quotient curve.
  output: 
    (BoolElt) True if the group satisfies Conder's genus bounds; otherwise false.
  description:
    For genus up to 301, uses Conder's empirical upper bounds for maximal group orders.
    For larger genus, applies the general Hurwitz bound (84(g - 1)).
*/ 
SatisfiesConderBound := function(gOrd, genus)
  CondersList := [48,168,120,192,150,504,336,320,432,240,120,360,1092,504,720,1344,
  168,720,228,480,1008,192,216,720,750,624,1296,672,264,720,372,1536,1320,544,672,
  1728,444,912,936,960,410,1512,516,1320,2160,384,408,1920,1176,1200,2448,832,456,
  1296,1320,1344,1368,928,504,1440,732,1488,1512,3072,1014,576,804,2448,660,672,
  710,2160,876,1776,1800,912,1176,1872,948,1536,3888,1312,696,4032,1200,2064,712,
  2640,744,3600,2184,768,2232,768,3420,3840,1164,2352,1320,2400,1010,4896,1236,
  2496,2016,1696,888,5184,1308,2640,2664,2688,936,3420,936,1856,9828,960,1734,2880,
  5808,2928,1080,1488,1500,3024,1524,10752,3096,2080,1310,3960,1596,3216,6480,2176,
  1128,1152,1668,2184,1144,1152,1352,6912,12180,3504,3528,2368,1224,3600,1812,3648,
  2448,1320,1550,3744,1884,3792,1288,3840,1320,3888,1956,2624,6600,1344,1368,8064,
  4056,2720,6840,2064,1416,1856,3000,2640,1432,2848,1464,8640,2172,4368,4392,1600,
  1512,4464,1512,1536,4536,2640,1910,5760,2316,4656,3420,4704,1608,2640,2388,4800,
  4824,3232,1656,7344,2050,4944,1672,4368,2904,5040,2532,3392,1720,1728,1800,6480,
  2604,5232,5256,5280,1800,5328,2676,5376,5400,3616,1848,6840,2748,6072,3276,3712,
  1896,11232,1896,1920,5688,2176,1944,9600,2892,5808,11664,3904,4116,2624,2964,5952,
  2008,12000,2510,6048,12144,6096,4896,12288,2088,6192,3108,4160,3132,2112,2136,7920,
  2200,6384,2152,3216,2184,5040,3252,4896,6552,4384,3750,2304,3324,6672,6696,2688,
  2810,2304,3396,2400,10260,4368,2328,13824,13872,4640,6984,4672,2376,7056,2376,7104,
  3888,4768,4056,9000]; 
  if genus le 301 then 
    return gOrd le CondersList[genus-1];
  end if; 
  return gOrd le 84*(genus-1);
end function;

/* 
  output: 
    (SetEnum) Set of tuples <group order, group identifier>.
  description:
    Returns a predefined set of groups of order less than 1024 that are known to admit Beauville structures
    of dimension 2 or 3, based on the classification by Carta and Fairbairn.
*/ 
GetAllGroupsDimTwoAndThree := function() 
  return { <600, 145>, <960, 5710>, <168, 42>, <512, 1641>, <360, 118>, <729, 35>, <432,735>, <512, 1632>, <480, 218>, 
  <672, 1043>, <972, 141>, <660, 13>, <300, 22>,<256, 306>, <729, 87>, <729, 44>, <256, 295>, <640, 21455>, <512, 335>, 
  <972,757>, <720, 769>, <960, 5699>, <960, 5716>, <720, 770>, <243, 13>, <864, 2666>, <640, 21454>, <960, 5693>, 
  <960, 5712>, <672, 1046>, <125, 3>, <960, 5709>, <625, 7>, <640, 21536>, <672, 1254>, <256, 298>, <512, 325>, 
  <640, 21456>, <343,3>, <336, 114>, <864, 4445>, <128, 36>, <841, 2>, <432, 763>, <960, 5721>, <512,351>, <864, 4680>, 
  <729, 9>, <512, 1637>, <960, 5694>, <900, 88>, <960, 639>, <625, 4>, <729, 40>, <121, 2>, <1008, 884>, <729, 38>, 
  <864, 4449>, <600, 54>, <960, 5704>, <480, 951>, <289, 2>, <729, 41>, <512, 1572>, <864, 4187>, <25, 2>, <512, 1649>, 
  <432, 749>, <729, 65>, <672, 1255>, <504, 156>, <864, 4446>, <720,415>, <720, 764>, <720, 413>, <729, 34>, 
  <972, 179>, <729, 45>, <960, 5696>, <960, 5719>, <960, 5707>, <275, 3>, <864, 4708>, <720, 409>, <720, 412>, <243,4>, 
  <120, 34>, <400, 213>, <864, 4186>, <961, 2>, <240, 90>, <672, 1044>, <800,1065>, <240, 91>, <361, 2>, <1008, 881>, 
  <729, 37>, <480, 952>, <512, 1643>, <864, 4665>, <972, 877>, <529, 2>, <729, 48>, <512, 1634>, <360, 119>, 
  <512, 1644>, <972, 135>, <336, 208>, <480, 950>, <720, 411>, <625, 2>, <1008, 517>, <972, 138>, <512, 1642>, 
  <432, 623>, <480, 219>, <640, 787>, <336, 209>, <960, 5723>, <480, 948>, <729, 68>, <720, 766>, <512, 1574>, 
  <576, 8652>, <960, 5725>, <392, 39>, <972, 183>, <243, 3>, <960, 5711>, <1008, 882>, <720, 763>, <1008,883>, 
  <720, 767>, <504, 157>, <775, 3>, <672, 1048>, <320, 1635>, <960, 5724>, <360, 120>, <49, 2>, <169, 2>, <324, 160>, 
  <240, 189> };
end function;


// *******************************************************************************************************************************
// *********************************************** SPECIAL TREATMENT OF THREEFOLDS ***********************************************
// *******************************************************************************************************************************


/* 
  input: 
    - chi: (RngIntElt) Holomorphic Euler characteristic.
  output: 
    (SeqEnum) List of sequences [ [GroupOrder], [T1], [T2], [T3] ] satisfying the conditions.
  description:
    Computes all admissible triple type data that could yield a rigid threefold with Euler number chi.
*/
AdNDunmixed := function(chi)
  boundGOrd := Floor(168*Sqrt(-21*chi));                  
  setAdUnmixed := {};
  ordersDimTwoAndThree := {grpId[1] : grpId in GetAllGroupsDimTwoAndThree()};
  for GOrd in [1..boundGOrd] do
    generas := {};
    divi := Divisors(-chi*GOrd);  
    for n1 in [1..#divi] do 
      for n2 in [n1..#divi] do
        d1 := divi[n1]; 
        d2 := divi[n2]; 
        d3 := -chi*GOrd/(d1*d2);
        if d3 in IntegerRing() and d3 ge d2 then 
          d3 := IntegerRing()!d3;
          if SatisfiesConderBound(GOrd,d1+1) and SatisfiesConderBound(GOrd,d2+1) and SatisfiesConderBound(GOrd,d3+1) then 
            Include(~generas,[d1+1,d2+1,d3+1]); // g-1 must divide -chi*GOrd
          end if; 
        end if; 
      end for; 
	  end for;
    for gens in generas do
      g1 := gens[1]; 
      g2 := gens[2]; 
      g3 := gens[3];						           
      for T1 in ListOfTypes(GOrd,g1,(g2-1)*(g3-1)) do   	     
        for T2 in ListOfTypes(GOrd,g2,(g1-1)*(g3-1)) do	         
          for T3 in ListOfTypes(GOrd,g3,(g1-1)*(g2-1)) do
            D := [[GOrd]] cat OrderedSeq([T1,T2,T3]);
            if GOrd in ordersDimTwoAndThree or GOrd ge 1023 then 	               	
              Include(~setAdUnmixed,D);
            end if; 
          end for; 
        end for; 
      end for;
    end for; 
  end for;
  return OrderedSeq(setAdUnmixed); 
end function;

/* 
  input: 
    - G : (Grp) A finite group.
    - order : (RngIntElt) Order of the searched elements 
  output: 
    (SetEnum) Set of elememts of G of given order 
  description:
    Computes all elements of G of given order
*/ 
ElemsOfOrder := function(G, order)
  return { g : g in G | Order(g) eq order};
end function;

/* 
  input: 
    - G : (Grp) A finite group.
    - T : (SeqEnum) The type of the wanted generating triple
  output: 
    (BoolElt) True if G admits a generating triple of type T.
  description:
    Checks whether the group G has a generating triple of the given type.
*/ 
TriangleCurvesEx := function(G, T)
  genTripleExists := false;
  for a in ElemsOfOrder(G,T[1]) do 
    for b in ElemsOfOrder(G,T[2]) do 
      if Order((a*b)^-1) eq T[3] and #sub<G|a,b> eq #G then 
        genTripleExists := true; 
        break a;
      end if; 
    end for;
  end for;
  return genTripleExists; 
end function;


/* 
  input: 
    - chi: (RngIntElt) Holomorphic Euler characteristic.
  output: 
    (SeqEnum) List of [ [GroupOrder, GroupId], [T1], [T2], [T3] ] where each Ti is the type of a generating triple.
  description:
    Filters admissible groups and triple types that actually exist in the database.
*/ 
AdGroups := function(chi)
  Results:=[];
  for D in AdNDunmixed(chi) do
	  GOrd := D[1][1]; 
    T1 := D[2]; 
    T2 := D[3]; 
    T3 := D[4];
    if GOrd ge 1024 then 
      printf "!! Warning skipping GOrd >= 1024 !! \nHere GOrd = %o possible, but will not be checked \n", GOrd;
      continue D;
    end if; 
    for H in GetAllGroupsDimTwoAndThree() do
      if H[1] eq GOrd then 
        G := SmallGroup(GOrd,H[2]);
        if TriangleCurvesEx(G,T1) and TriangleCurvesEx(G,T2) and TriangleCurvesEx(G,T3) then 
          Append(~Results, [[GOrd,H[2]],T1,T2,T3]);
        end if;
      end if;
   	end for;
  end for;
  return Results;
end function;

// *******************************************************************************************************************************
// ***************************************************** BRAIDMOVE FUNCTIONS *****************************************************
// *******************************************************************************************************************************

BraidMove := function(V,i)
   c := V[i]*V[i+1]*V[i]^-1;
   return Insert(Remove(V,i+1),i,c); 
end function;

BraidMove2 := function(S,i)
	c := S[i]*S[i+1]*S[i]^-1;
  return Insert(Remove(S,i+1),i,c); 
end function;

BraidOrbit2 := function(S)
  orb := { }; 
  Trash := {S}; 
  repeat
    ExtractRep(~Trash, ~gens); 
    Include(~orb, gens);
		for i in [1,2] do 
      newgens := BraidMove2(gens,i);
	    if newgens notin orb then 
			  Include(~Trash, newgens);
		  end if; 
	  end for;
  until IsEmpty(Trash);
  return orb;
end function;

/* 
  input: 
    - G : (Grp) A finite group.
    - T : (SeqEnum) The type of the wanted generating triple
  output: 
    (SetEnum) Set of generating tiples (one from each orbit)
  description:
    Computes all orbits of gen triples of type T under the action of Braidmoves 
*/ 
TriangleModB3 := function(G,T);
  gens := {};
  for a in ElemsOfOrder(G,T[1]) do 
	for b in ElemsOfOrder(G,T[2]) do 
	  c := (a*b)^-1;
	  if Order(c) eq T[3] and #sub<G|a,b> eq #G then 
		Include(~gens,[a,b,c]);
	  end if; 
	end for;
  end for;
  Repres := {};
  while not IsEmpty(gens) do 
    V := Rep(gens); 
    Include(~Repres,V); 
    orb := BraidOrbit2(V);
    for W in orb do 
      Exclude(~gens, W);
    end for; 
  end while; 
  return Repres;
end function;

/* 
  input: 
    - G : (Grp) A finite group.
    - V : (SeqEnum) Generating triple of G of the form (h1, ..., hr)
  output: 
    (SetEnum) Stabilizing-set of V
  description:
    Determines the stabilizer set (union of conjugates of the subgroups generated by th elements) of V  
*/ 
StabSet := function(G,V) 
  stabS := {Id(G)};
  for i in [1..#V] do 
    g := V[i]; 
	  for n in [1..Order(g)-1] do
      try 
	      stabS := stabS join Conjugates(G,g^n);
      catch e // sometimes problems with assignment of g to G 
        H := Parent(V[i]);
        isIsom, isom := IsIsomorphic(H, G);
        g := isom(V[i]);
        stabS := stabS join Conjugates(G,g^n);
        // printf "Problem im StabSet besteht mit G = %o  H = %o \n",IdentifyGroup(G), IdentifyGroup(H);
      end try;
	  end for; 
  end for;
  return stabS; 
end function;

/* 
  input: 
    - G : (Grp) A finite group.
    - normSubGr :  The kernel K such that imageSet is subset of G/K
    - pi : (Map) The map from G to G/K
    - imageSet : (SetEnum) Subset of a quotient of G of which the preimage will be computed
  output: 
    (SetEnum) the preimage of imageSet
  description:
    Determines the preimage of a subset of a quotient of the group G
*/ 
GetPreimageSet := function(G, normSubGr, pi, imageSet)
  preImageReps := imageSet @@ pi;
  return {rep*ns : ns in normSubGr, rep in preImageReps};
end function;



/*
  input: 
    - chi: (RngIntElt) Holomorphic Euler characteristic.
  output: 
    (SeqEnum) List of [ [GroupOrder, GroupId], [T1], [T2], [T3] ] where each Ti is the type of a generating triple.
  description:
    Computes all possible sequences, which admit gen triples with trivial intersection of stabilizers  
*/
ExistsThreefold := function(chi)
  dataGenTriples := AdGroups(chi);
  Results := [];
  for dataSeq in dataGenTriples do 
    GOrd := dataSeq[1][1];
    n := dataSeq[1][2]; 
    T1 := dataSeq[2]; 
    T2 := dataSeq[3]; 
    T3 := dataSeq[4];
    G := SmallGroup(GOrd,n); 
    for V1 in TriangleModB3(G,T1) do 
      for V2 in TriangleModB3(G,T2) do 
        for V3 in TriangleModB3(G,T3) do 
          if #(StabSet(G,V1)  meet StabSet(G,V2) meet StabSet(G,V3)) eq 1 then 
            Append(~Results, dataSeq); 
            break V1;		
          end if;
        end for;
      end for;
    end for;
  end for;
  return Results;
end function;

/*  
  input: 
    - G : (Grp) A finite group.
  output: 
    (BoolElt) True if the group G is a triangle group (Theta > 0)
  description:
    Checks if G is a triangle group
*/
IsTriangleGroup := function(G)
  genTripleExists := false;
  for g in G do 
    for h in G do 
      type := [Order(g), Order(h), Order(g*h)];
      if type[1] le type[2] and type[2] le type[3] and Theta(type) gt 0 and #sub<G|g,h> eq #G then 
        genTripleExists := true; 
        break g;
      end if; 
    end for;
  end for;
  return genTripleExists; 
end function;

/*  
  input: 
    - maxGrpOrder : (RngIntElt) Grouporder up to which the groups are checked
  output: 
    (SetEnum) Group-identifying tuples <GroupOrder, GroupId>
  description:
    Reads the already computed triangle groups and computes further if needed
*/
GetGoodTriangleGroups := function(maxGrpOrder)
  grpOrder, grpIdenitifier := ReadTriangleGroups();
  if maxGrpOrder gt grpOrder then // We need to compute further
    for n in [(grpOrder+1)..maxGrpOrder] do  
      for idx in [1..NumberOfSmallGroups(n)] do
        if IsTriangleGroup(SmallGroup(n, idx)) then
          Include(~grpIdenitifier, <n,idx>);
        end if;
      end for;
    end for;
    WriteTriangleGroups(maxGrpOrder, grpIdenitifier);
  end if;
  return grpIdenitifier;
end function;


/*  
  input: 
    - G : (Grp) A finite group.
    - T : (SeqEnum) The type of the wanted generating triple
  output: 
    (SetEnum) Generating tiples of the (Aut x B_3)-action (one from each orbit)
  description:
    Computes all orbits of gen triples of type T under the action of braidmoves and automorphisms 
*/
// Alternatively, use the results of "Database of topological types of actions on curves"
TriangleModB3Aut := function(G,T);
  AutG := AutomorphismGroup(G); 
  f,P := PermutationRepresentation(AutG);
  gens := {};
  for a in ElemsOfOrder(G,T[1]) do 
	  for b in ElemsOfOrder(G,T[2]) do 
		  c:=(a*b)^-1;
			if Order(c) eq T[3] and #sub<G|a,b> eq #G then 
				Include(~gens,[a,b,c]);
			end if; 
	  end for;
  end for;
  Repres := {};
  while not IsEmpty(gens) do 
    V := Rep(gens); 
    Include(~Repres,V); 
    orb := BraidOrbit2(V);
    for W in orb do 
      for p in P do
        phi := p @@ f; 	
        Exclude(~gens, [phi(W[1]),phi(W[2]),phi(W[3])]);
      end for; 
    end for;
  end while; 
  return Repres;
end function;


// *******************************************************************************************************************************
// ******************************************************** HODGE NUMBERS ********************************************************
// *******************************************************************************************************************************


/*  
  input: 
    - X : character of a group G 
    - g : element of a group G
  output: 
    - value rho(g), where "rho" is a representation affording the character X
  description:
    coefficients of the characteristic polynomial are determined from the power sums of the eigenvalues X(g^i)
*/ 
CharPol := function(X,g)
  d := Degree(X); 
  L := [];
	for i in [1..d] do
		Append(~L,X(g^i));
	end for;
  return Polynomial(PowerSumToCoefficients(L)); 
end function;

/*  
  input: 
    - pol : characteristic polynomial 
    - n : dimension of the presentation
  output: 
    - (SeqEnum) of elements [a, mult] where exp(a*2*pi*i/n) is a root of pol(X,g) with multiplicity mult
  description:
    Compute the diagonal of the representation
*/ 
ListNa := function(pol,n)
	L := []; 
  F := CyclotomicField(n); 
  z := F.1;            // z is equal to exp(2*pi*i/n)
	for r in Roots(Polynomial(F,pol),F) do
		for a in [0..n-1] do 
			if r[1] eq z^a then            
				Append(~L,[a,r[2]]); 
        break a; 
			end if;
		end for; 
  end for;
  return L; 
end function;

/*  
  input: 
    - V : the generating vector of G of the form (h1,...,hr,a1,b1,...,ah,bh)
    - T : the type of V of the form [ genus(C/G), ord(h1), ..., ord(hr) ]
    - X : irreducible character of a group G 
    - G : the group
  output: 
    - the multiplicity of X in the character of the representation obtained by V 
  description:
    Implementation of the Chevalley-Weil formula
*/ 
ChevWeil := function(V,T,X,G)
  h := T[1]; 
  ram := [V[i] : i in [1..#T-1]];  // ram is the ramification part of V i.e. (h1,...,hr)
  mult := (h-1)*X(Id(G));
  for g in ram do 
    n := Order(g); 
    for t in ListNa(CharPol(X,g),n) do
      mult := t[1]*t[2]/n + mult;
    end for; 
  end for;
  return mult; 
end function;

/*  
  input: 
    - G : the group
    - V : the generating vector of G o fthe form (h1,...,hr,a1,b1,...,ah,bh)
    - T : the type of V of the form [ genus(C/G), ord(h1), ..., ord(hr) ]
  output: 
    - character of the representation "phi" from a corresponding generating vector 
  description:
    Determines the character of the representation (using Chevalley-Weil formula)
*/
Chi_phi := function(G,V,T)
  CT := CharacterTable(G); 
  h := T[1]; 
  char := h*CT[1]; // CT[1] is the trivial character of G
	for i in [2..#CT] do 
	  X := CT[i]; 
    m := ChevWeil(V,T,X,G); 
    char := char + m*X;
	end for;
  return char; 
end function;

/*  
  input: 
    - G0 : the maximal diagonal acting subgroup of G  
    - X1 : the character of the action on the first factor
    - X2 : the character of the action on the second factor
    - X3 : the character of the action on the third factor
  output: 
    - (SeqEnum) the hodge numbers [h30, h20, h10, h11, h21]
  description:
    Computes the Hodgenumbers [h30, h20, h10, h11, h21]
*/
HodgeDiamondInvarG0 := function(G0,X1,X2,X3) 
  p := PrincipalCharacter(G0); 
  X1bar := ComplexConjugate(X1);
  X2bar := ComplexConjugate(X2); 
  X3bar := ComplexConjugate(X3); 
  chi10 := X1 + X2 + X3;
  chi20 := X1*X2 + X1*X3 + X2*X3;
  chi30 := X1*X2*X3;
  chi11 := X1bar*X2 +  X1*X2bar +  X1bar*X3 +  X1*X3bar + X2bar*X3 + X2*X3bar + 3*p;
  chi21 := X1bar*X2*X3 +  X1*X2bar*X3  + X1*X2*X3bar + 2*(X1 + X2 + X3); 
  h10 := InnerProduct(chi10,p); 
  h20 := InnerProduct(chi20,p); 
  h30 := InnerProduct(chi30,p);
  h11 := InnerProduct(chi11,p); 
  h21 := InnerProduct(chi21,p);
  return [h30, h20, h10, h11, h21]; 
end function;

/*  
  input: 
    - G : the maximal diagonal actimg subgroup of G  
    - V1 : the first generating triple of G
    - V2 : the second generating triple of G
    - V3 : the third generating triple of G
    - T1 : the type of V1 of the form [ genus(C/G), ord(h1), ..., ord(hr) ]
    - T2 : the type of V2 of the form [ genus(C/G), ord(h1), ..., ord(hr) ]
    - T3 : the type of V3 of the form [ genus(C/G), ord(h1), ..., ord(hr) ]
  output: 
    - (SeqEnum) the hodge numbers [h30, h20, h10, h11, h21]
  description:
    Computes the Hodgenumbers [h30, h20, h10, h11, h21] in the case of an unmixed, absol. faithfull action and n = 3
*/
HodgeDiamondUnmixed := function(G,V1,V2,V3,T1,T2,T3) 
  X1 := Chi_phi(G,V1,T1); 
  X2 := Chi_phi(G,V2,T2); 
  X3 := Chi_phi(G,V3,T3);     
  return HodgeDiamondInvarG0(G,X1,X2,X3);
end function;

/*
  input: 
    - G: the group  
    - kernelComps: (SeqEnum) of KernelComponent-records
    - elem: the element of generating triples in the tuple form as returned by GetBiholomReprs
  output: 
    - (SeqEnum) the hodge numbers [h30, h20, h10, h11, h21]
  description:
    Computes the Hodgenumbers [h30, h20, h10, h11, h21] in the case of an unmixed action and n = 3
*/
ComputeHodgeNumbers := function(G, kernelComps, elem)
  charSeq := [];
  for kerIndx in [1..#elem] do
    quotSubGr, piQuotSubGr := G / kernelComps[kerIndx]`kernel`subgroup;
    for V in elem[kerIndx] do
      parentTriple := Parent(V[1]);
      isIsom, isom := IsIsomorphic(parentTriple, quotSubGr);
      triple := [isom(V[1]),isom(V[2]),isom(V[3])];
      Append(~charSeq, LiftCharacter(Chi_phi(quotSubGr, triple, [0,Order(triple[1]),Order(triple[2]),Order(triple[3])]), piQuotSubGr, G));
    end for;
  end for;
  return HodgeDiamondInvarG0(G,charSeq[1],charSeq[2],charSeq[3]);
end function;

// *******************************************************************************************************************************
// ******************************************* ORBIT REPRESENTATIVES OF CHAPTER 4 ACTION *****************************************
// *******************************************************************************************************************************


// rec<KernelComponent | kernel := (Kernelrecord), typeRecs := [ TypeRec1, ... ];
KernelComponent := recformat<kernel : Rec, typeRecs : SeqEnum, qAutsMap : Map, quotSubGr : GrpPC, piQuotSubGr : HomGrp>;
// rec<TypeRec | bAutReprs := [ BAutRepr1, ... ], multType := (how often this type occurs), type := [n1, n2, n3]>;
TypeRec := recformat<bAutReprs : SeqEnum, multType : RngIntElt, type : SeqEnum>;
// rec<BAutRepr | triple := [V1, V2, V3], permReprOrbits := [PermReprOrbit1, ... ], autKAct := (The kern-stabilizing action on the orbit reprs)>;
BAutRepr := recformat<triple : SeqEnum, permReprOrbits : SeqEnum, autKAct : Map>;
// rec<PermReprOrbit | rep := (Fixed repr of orbit), orbit := (orbit of the Aut(G/Ni) / BAut(G/Ni, Si) action) >;
PermReprOrbit := recformat<rep : GrpPermElt, orbit : SetEnum>;

/* 
  input: 
    - G: (Grp) A finite group.
    - kernelComps: (SeqEnum) of KernelComponent-records
  description:
    Outputs relevant information of the kernel components (for debugging and quick information on console)
*/
procedure OutputKernelComps(G, kernelComps)
  printf "-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ \n";
  for kernelComp in kernelComps do
    quotSubGr, piQuotSubGr := G / kernelComp`kernel`subgroup;
    printf "- Kernel: %o   QuotSubGr: %o   #Types: %o \n", IdentifyGroup(kernelComp`kernel`subgroup), IdentifyGroup(quotSubGr), #kernelComp`typeRecs;
    for typeRec in kernelComp`typeRecs do
      printf "-- Type: %o with multiplicity %o   #TriplesMod3Aut: %o \n", typeRec`type, typeRec`multType, #typeRec`bAutReprs;
      cnt := 0;
      for bAutRepr in typeRec`bAutReprs do
        cnt +:= 1;
        printf "--- %o.Triple has %o Aut()/BAut(S%o) orbits \n", cnt, #bAutRepr`permReprOrbits ,cnt;
      end for;
    end for;
  end for;
  printf "-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ \n";
end procedure;

/* 
  input: 
    - qAuts: the automorphisms Aut(G/Ni) (as permutation represenation)
    - braidAuts: the automorphisms of braid type BAut(G/Ni, Si) (as permutation represenation)
  output: 
    - (SeqEnum) of PermReprOrbit-records
  description:
    Computes the orbits of Aut(G/Ni) / BAut(G/Ni, Si) and fixes an orbit representative
*/
BAutOrbits := function(qAuts, braidAuts)
  permReprOrbits := []; 
  qAutsSet := Set(qAuts);
  while not IsEmpty(qAutsSet) do
    elem := Rep(qAutsSet); 
    orbitElem := { bAut * elem : bAut in braidAuts}; // orbit of elem
    Append(~permReprOrbits, rec<PermReprOrbit | rep := elem, orbit := orbitElem>);
    qAutsSet diff:= orbitElem; 
  end while;
  return permReprOrbits;
end function;

/* 
  input: 
    - tripleReps: the triple representatives
    - qAutsMap: the map to get the automorphisms form their permutation represenation
    - qAuts: the automorphisms Aut(G/Ni) (as permutation represenation)
    - mapRestrKerStabAutQ: the relevant automorphisms of the quotient-group
  output: 
    - (SeqEnum) of PermReprOrbit-records
  description:
    Computes the orbits of Aut(G/Ni) / BAut(G/Ni, Si) and fixes an orbit repr 
*/
GetBAutReprs := function(tripleReps, qAutsMap, qAuts, mapRestrKerStabAutQ)
  bAutReprs := [];
  for triple in tripleReps do 
    braidOrbs := BraidOrbit2(triple); 
    braidAuts := sub<qAuts | {p : p in qAuts | [(p @@ qAutsMap)(elem) : elem in triple] in braidOrbs}>;
    permReprOrbits := BAutOrbits(qAuts, braidAuts);
    permReprs := {pRO`rep : pRO in permReprOrbits};
    cartProd := [x : x in CartesianProduct(permReprs, Domain(mapRestrKerStabAutQ))];
    // The action of the kernel stabilizing automorphisms on the orbit representatives 
    mapForStabilizingAuts := map< cartProd -> permReprs | [<x, y> : x in cartProd | exists(y){ pRO`rep : pRO in permReprOrbits | 
                x[1]*mapRestrKerStabAutQ(x[2]) in pRO`orbit }] >; 
    Append(~bAutReprs, rec<BAutRepr | triple := triple, permReprOrbits := permReprOrbits, autKAct := mapForStabilizingAuts>);
  end for;
  return bAutReprs;
end function;

/* 
  input: 
    - G: (Grp) A finite group.
    - autK: the kernel stabilizing automorphisms (as permutation represenation)
    - kernel: (Rec) the kernel-record
    - types: (SeqEnum/multiset) of according types of the form [n1,n2,n3] 
  output: 
    - (SeqEnum) of KernelComponent-records 
  description:
    Computes all the relevant information/oorbits of the KernelComponents we need for the GetBiholomReprs-function
*/
GetKernelComp := function(G, autK, kernel, types)
  gAutsMap, gAuts := PermutationRepresentation(AutomorphismGroup(G)); 
  quotSubGr, piQuotSubGr := G / kernel`subgroup;
  autQ := AutomorphismGroup(quotSubGr);
  qAutsMap, qAuts := PermutationRepresentation(autQ);
  // Get the automorphisms on the quotient G/K of the kernel K stabilizing automorphisms of G
  kerStabAutQ := { qAutsMap(autQ !hom<quotSubGr -> quotSubGr | 
                    [<el,piQuotSubGr(Inverse(gAutsMap)(g)(Inverse(piQuotSubGr)(el)))> : el in Generators(quotSubGr)]>) : g in autK };
  // Get the map from the kernel K stabilizing automorphisms of G to the automorphisms on the quotient G/K
  mapRestrKerStabAutQ := map< autK -> kerStabAutQ | [ <g, qAutsMap(autQ !hom<quotSubGr -> quotSubGr | 
                    [<el,piQuotSubGr((Inverse(gAutsMap)(g))(Inverse(piQuotSubGr)(el)))> : el in Generators(quotSubGr)]>)> : g in autK] >;
  kernelComp := rec<KernelComponent | kernel := kernel, typeRecs := [], qAutsMap := qAutsMap>;
  // Compute for all Type-records the Aut(G/Ni) / BAut(G/Ni, Si) representatives
  branchingTypes := Multiset(types);
  for type in Set(branchingTypes) do
    // Get the triples from database/computation
    modB3AutReps := TriangleModB3Aut(quotSubGr, type);
    bAutReprs := GetBAutReprs(modB3AutReps, qAutsMap, qAuts, mapRestrKerStabAutQ);
    Append(~(kernelComp`typeRecs), rec<TypeRec | bAutReprs := bAutReprs, multType := Multiplicity(branchingTypes, type), type := type>);
  end for; 
  return kernelComp;
end function;

/* 
  input: 
    - triple: the triples of permuation representations
    - permRep: the permutation representation of the automorphism that will be applied on the triple
    - qAutsMap: the map to get the automorphism from the permutation representation
  output: 
    - (SeqEnum) The generating triple after application of the automorphism
  description:
    Computes image of the generating triple
*/
ApplyPermRepToTriple := function(triple, permRep, qAutsMap)
  autQ := permRep @@ qAutsMap;
  return [autQ(triple[1]), autQ(triple[2]), autQ(triple[3])];
end function;


/* 
  input: 
    - kernelComps: (SeqEnum) of KernelComponent-records
    - autK: the automorphisms stabilizing the kernels of kernelComps
  output: 
    - (SetEnum) of tuples of generating triples for each KernelComponent and each Type
  description:
    Computes a generating triple for each biholomorphism class. The output is seperated in tuples to distinguish the KernelComponents and Types    
*/
GetBiholomReprs := function(kernelComps, autK) 
  genTripleTups := {};
  unionTypeRecs := [ [tR : tR in kC`typeRecs] : kC in kernelComps]; 
  // Sadly we have to do the computations with indices because Records are not hashable
  for bAutReprIndxTup in CartesianProduct([ CartesianProduct([ Multisets({1..#tR`bAutReprs}, tR`multType) : tR in kC`typeRecs]) : kC in kernelComps]) do
    qBAutReprs := {};
    // The information in the tuple
    bAutReprs := <<< unionTypeRecs[i][j]`bAutReprs[indx] : indx in ToSet(bAutReprIndxTup[i][j])> : j in [1..#kernelComps[i]`typeRecs]> : i in [1..#kernelComps]>;

    cP := CartesianProduct(  [ CartesianProduct ( [ CartesianProduct( [  Multisets({ pR`rep : pR in unionTypeRecs[i][j]`bAutReprs[indx]`permReprOrbits}, 
            CntElem(bAutReprIndxTup[i][j], indx)) : indx in ToSet(bAutReprIndxTup[i][j]) ]) : j in [1..#kernelComps[i]`typeRecs] ]) : i in [1..#kernelComps] ] );
    qBAut := { x : x in cP };
    // Get for each orbit of the kernel stabilizing automorphisms a representative
    while(not IsEmpty(qBAut)) do
      elem := Rep(qBAut);
      Include(~qBAutReprs, elem);
      autKOrbitElem := { <<< {* bAutReprs[i][j][mSIndx]`autKAct(mSRepr,aK) : mSRepr in elem[i][j][mSIndx] *} : mSIndx in [1..#elem[i][j]] > : 
                               j in [1..#kernelComps[i]`typeRecs] > : i in [1..#kernelComps] > : aK in autK }; 
      qBAut diff:= autKOrbitElem;
    end while;
    genTripleTups join:= { << ApplyPermRepToTriple(bAutReprs[i][j][mSIndx]`triple, mSRepr, kernelComps[i]`qAutsMap) : mSRepr in qBAutRepr[i][j][mSIndx], 
          mSIndx in [1..#qBAutRepr[i][j]], j in [1..#kernelComps[i]`typeRecs] > : i in [1..#kernelComps] > : qBAutRepr in qBAutReprs};
  end for;
  return genTripleTups;
end function;


// *******************************************************************************************************************************
// *********************************************** COMPUTATIONS FOR CHAPTER 6 TABLES *********************************************
// *******************************************************************************************************************************

/* 
  input:
    - [optional] logOutput : (BoolElt) indicating whether the result are logged in a file
    - [optional] outputTriples: (BoolElt) indicating whether the generating triples are logged/output
  description:
    Computes spherical systems of generators for each biholomorphism class with hodge numbers for the group (Z_5)^2 with kernels 
    K1 = <Id(G)>, K2 = <Id(G)>, K3 = <G.2> (unmixed action, n=3)
*/
procedure CheckSmallestGroupWithKernel(: logOutput := false, outputTriples := false)
  // For file logging, to stop automatic linebreaks
  SetAutoColumns(false);
  SetColumns(0);
  G := SmallGroup(25,2); 
  type1 := [5,5,5]; 
  type2 := [5,5,5]; 
  type3 := [5,5,5];
  printf "\n----------------------- %o with kernels ----------------------- \n", IdentifyGroup(G);
  logText := "Type: " cat Sprint(<type1, type2, type3>) cat " \n";
  startTime := Realtime(); // For compution-time logging
  gAutsMap, gAuts := PermutationRepresentation(AutomorphismGroup(G)); 
  // the kernels form the proof of section 6
  normSubGrps :=  [ nSG : nSG in NormalSubgroups(G) | nSG`order eq 1 or nSG`subgroup eq sub<G|G.2>]; 
  kerSubGrps := { kerComp`subgroup : kerComp in normSubGrps }; // the subgroups to the kernels
  kerStabilizingAuts := sub<gAuts | {p : p in gAuts | forall{K : K in kerSubGrps | sub<G | {(p @@ gAutsMap)(h) : h in K}> in kerSubGrps} }>;
  Q1, pi1 := G / normSubGrps[1]`subgroup;
  Q2, pi2 := G / normSubGrps[2]`subgroup;
  branchingTypes := {* type1,type2,type3 *}; 
  // Two trivial kernels and one non trivial kernels
  kernelComponents := [GetKernelComp(G, kerStabilizingAuts, normSubGrps[1], [type1,type2]), GetKernelComp(G, kerStabilizingAuts, normSubGrps[2], [type3])]; 
  logText cat:= "Kernels: " cat Sprint(<IdentifyGroup(kerComp`kernel`subgroup) : mult in [1..typeRec`multType] ,typeRec in kerComp`typeRecs, 
      kerComp in kernelComponents>) cat " \n";
  OutputKernelComps(G, kernelComponents);
  IsoClasses := {* *};
  tripleCombs := GetBiholomReprs(kernelComponents, kerStabilizingAuts);
  for elem in tripleCombs do 
    if #(GetPreimageSet(G, normSubGrps[1]`subgroup, pi1, StabSet(Q1, elem[1][1])) meet
         GetPreimageSet(G, normSubGrps[1]`subgroup, pi1, StabSet(Q1, elem[1][2])) meet
         GetPreimageSet(G, normSubGrps[2]`subgroup, pi2, StabSet(Q2, elem[2][1]))) eq 1 then 
      hodgeNumbers := ComputeHodgeNumbers(G, kernelComponents,  elem);
      Include(~IsoClasses, hodgeNumbers);
      if outputTriples then 
        outputText := "Hodgenumber: " cat Sprint(hodgeNumbers) cat "\n" cat "Triple: " cat Sprint(elem) cat "\n \n";
        logText cat:= outputText;
        printf outputText;
      end if;
    end if;
  end for; 
  printf "Hodgenumbers: %o \n", IsoClasses;
  logText cat:= "Hodgenumbers: " cat Sprint(IsoClasses) cat "\n \n";
  endTime := Realtime(); 
  durationTime := Round(endTime - startTime);
  logText cat:= "Computationtime:   " cat IntegerToString(Round(durationTime / 60)) cat " min " cat IntegerToString(durationTime mod 60) cat " sec \n";
  if logOutput then
    folderstructure := _FOLDERNAME cat "/" cat _KERNELS_SUBFOLDERNAME;
    if System("test -d " cat folderstructure) ne 0 then
      System("mkdir " cat folderstructure);
    end if;
    PrintFile(folderstructure cat "/" cat "Grouporder" cat IntegerToString(#G) cat ".txt", logText : Overwrite := true);
  end if;
end procedure;

/* 
  input:
    - chiMax : (RngIntElt) the upper bound of checked values of chi
    - [optional] logOutput : (BoolElt) indicating whether the result are logged in a file
    - [optional] chiMin: (RngIntElt) the lower bound of checked values of chi
    - [optional] outputTriples: (BoolElt) indicating whether the generating triples are logged/output
  description:
    Computes spherical systems of generators for each biholomorphism class with hodge numbers for the given values of chi 
    (absolutely faithfull, unmixed action, n=3)
*/
procedure CheckThreefoldsAbsFaithfull(chiMax: logOutput := false, chiMin := 1, outputTriples := false)
  // For file logging, to stop automatic linebreaks
  SetAutoColumns(false);
  SetColumns(0);
  // Avoid problems with sign 
  chiMin := Abs(chiMin);
  chiMax := Abs(chiMax);
  for i in [chiMin..chiMax] do
    printf "\n---------------------------  chi = -%o  --------------------------- \n", i;
    logText := "";
    startTime := Realtime();
    threefoldData := ExistsThreefold(-i);
    for data in threefoldData do 
      logText cat:= "------------------ Group: " cat Sprint(<data[1][1],data[1][2]>) cat " ------------------ \n";
      G := SmallGroup(data[1][1],data[1][2]); 
      type1 := data[2]; 
      type2 := data[3]; 
      type3 := data[4];
      logText cat:= "Type: " cat Sprint(<type1, type2, type3>) cat " \n";
      gAutsMap, gAuts := PermutationRepresentation(AutomorphismGroup(G)); 
      // the trivial subgroup, because the action is absolutely faithfull
      normSubGrp :=  [ nSG : nSG in NormalSubgroups(G) | nSG`order eq 1][1]; 
      kerStabilizingAuts := gAuts; // All stabilize the trivial subgroup
      qSubGrp, piQ := G / normSubGrp`subgroup;
      branchingTypes := {* type1,type2,type3 *}; 
      kernelComponents := [GetKernelComp(G, kerStabilizingAuts, normSubGrp, branchingTypes)]; // Only one kernel
      OutputKernelComps(G, kernelComponents);
      IsoClasses := {* *};
      tripleCombs := GetBiholomReprs(kernelComponents, kerStabilizingAuts);
      for elem in tripleCombs do 
        // Problems with right group assignment in GetPreimageSet(G, normSubGrp`subgroup, piQ, StabSet(qSubGrp, trip)) 
        V := elem[1];
        try 
          if not V[1][1] in qSubGrp then error "Problem"; end if;
        catch e
          H := Parent(V[1][1]); 
          isIsom, isom := IsIsomorphic(H, qSubGrp);
          V := [[isom(V[1][1]),isom(V[1][2]),isom(V[1][3])], [isom(V[2][1]),isom(V[2][2]),isom(V[2][3])], [isom(V[3][1]),isom(V[3][2]),isom(V[3][3])]];
        end try;
        if #&meet([ GetPreimageSet(G, normSubGrp`subgroup, piQ, StabSet(qSubGrp, trip)) : trip in V]) eq 1 then 
          hodgeNumbers := HodgeDiamondUnmixed(qSubGrp ,V[1],V[2],V[3],[0] cat type1,[0] cat type2,[0] cat type3);
          Include(~IsoClasses, hodgeNumbers);
          if outputTriples then 
            outputText := "Hodgenumber: " cat Sprint(hodgeNumbers) cat "\n" cat "Triple: " cat Sprint(V) cat "\n";
            logText cat:= outputText;
            print(outputText);
          end if;
        end if;
      end for; 
      printf "Hodgenumbers: %o \n", IsoClasses;
      logText cat:= "Hodgenumbers: " cat Sprint(IsoClasses) cat "\n \n";
    end for;
    endTime := Realtime(); 
    durationTime := Round(endTime - startTime);
    logText cat:= "Computationtime:   " cat IntegerToString(Round(durationTime / 60)) cat " min " cat IntegerToString(durationTime mod 60) cat " sec \n";
    if logOutput then
      folderstructure := _FOLDERNAME cat "/" cat _CHI_SUBFOLDERNAME ;
      if System("test -d " cat folderstructure) ne 0 then
        System("mkdir " cat folderstructure);
      end if;
      PrintFile(folderstructure cat "/" cat "chi" cat IntegerToString(i) cat ".txt", logText : Overwrite := true);
    end if;
  end for;
end procedure;

// Directory-check
if System("test -d " cat _FOLDERNAME) ne 0 then
  System("mkdir " cat _FOLDERNAME);
end if;
// Calling our two cases 
CheckSmallestGroupWithKernel(: logOutput := true);
CheckThreefoldsAbsFaithfull(-5: logOutput := true);