ListOfTypes:=function(Ksquare: MaximalCoefficient:=3*(Ksquare+2))
  list:=[* *];
  for m1 in [ 1..MaximalCoefficient] do 
  for m2 in [m1..MaximalCoefficient] do
  for m3 in [m2..MaximalCoefficient] do
  for m4 in [m3..MaximalCoefficient] do
  for m5 in [m4..MaximalCoefficient] do
  for m6 in [m5..MaximalCoefficient] do
  for m7 in [m6..MaximalCoefficient] do
    Theta:=5-(m1^-1)-(m2^-1)-(m3^-1)-(m4^-1)-(m5^-1)-(m6^-1)-(m7^-1);
    if Theta ne 0 then Alpha:=Ksquare*4^-1*Theta^-1;
    if Alpha in IntegerRing() and Alpha ge 1 then
       seqm:=[Integers() | m1,m2,m3,m4,m5,m6,m7 ];
    if forall{m : m in seqm | 2*Alpha/m in IntegerRing()} then bads:=0;
       for m in seqm do
         if Alpha/m notin IntegerRing() then bads +:=1;
       end if; end for;
    if bads le 4-Ksquare/2 then Append(~list,[seqm,[Alpha]]);
    end if; end if; end if; end if;
  end for; end for; end for; end for; end for; end for; end for;
return list;
end function;

ElsOfOrd:=function(group,order)
  Els:={ };
  for g in group do if Order(g) eq order then Include(~Els, g);
  end if; end for;
  return Els;
end function;

ExistSphericalGenerators:=function(group,seq)
  test:=false;
  for x1 in ElsOfOrd(group,seq[1]) do
  for x2 in ElsOfOrd(group,seq[2]) do
  for x3 in ElsOfOrd(group,seq[3]) do
  for x4 in ElsOfOrd(group,seq[4]) do
  for x5 in ElsOfOrd(group,seq[5]) do
  for x6 in ElsOfOrd(group,seq[6]) do
    if Order(x1*x2*x3*x4*x5*x6) eq seq[7]
       and #sub<group|x1,x2,x3,x4,x5,x6> eq #group
       then test:=true; break x1;
    end if;
  end for; end for; end for; end for; end for; end for;
return test;
end function;

ListGroups:=function(Ksquare)
  list:=[* *]; L:=ListOfTypes(Ksquare);
    for t1 in [1..#L] do
    for t2 in [t1..#L] do
      T1:=L[t1]; T2:=L[t2]; ord:=8*T1[2][1]*T2[2][1]/Ksquare;
    if ord in IntegerRing() and ord le 2000 and ord notin {1024, 1152} then
    for G in SmallGroups(IntegerRing()!ord: Warning := false) do
    if ExistSphericalGenerators(G,T1[1]) then
    if ExistSphericalGenerators(G,T2[1]) then
      Append(~list,[* T1[1],T2[1],G *]);
    end if; end if; end for; end if; end for; end for;
return list;
end function;

FSGUpToConjugation:=function(gr,seq)
  Heaven:={@ @}; Hell:={@ @};
    for x1 in ElsOfOrd(gr,seq[1]) do
    for x2 in ElsOfOrd(gr,seq[2]) do
    for x3 in ElsOfOrd(gr,seq[3]) do
    for x4 in ElsOfOrd(gr,seq[4]) do
    for x5 in ElsOfOrd(gr,seq[5]) do
    for x6 in ElsOfOrd(gr,seq[6]) do
      if Order(x1*x2*x3*x4*x5*x6) eq seq[7] then
      if #sub<gr|x1,x2,x3,x4,x5,x6> eq #gr then
      if [x1,x2,x3,x4,x5,x6,(x1*x2*x3*x4*x5*x6)^-1] notin Hell then
         Include(~Heaven,[x1,x2,x3,x4,x5,x6,(x1*x2*x3*x4*x5*x6)^-1]);
        for g in gr do
         Include(~Hell, [x1^g,x2^g,x3^g,x4^g,x5^g,x6^g,((x1*x2*x3*x4*x5*x6)^-1)^g]);
        end for;
      end if; end if; end if;
    end for; end for; end for; end for; end for; end for;
return Heaven;
end function;

CheckSings:=function(gr,seq1,seq2,Ksquare)
  SetGens1:=FSGUpToConjugation(gr,seq1);
  SetGens2:=FSGUpToConjugation(gr,seq2);
  test:=false;
    for gens1 in SetGens1 do
    for gens2 in SetGens2 do
      Bool:=true; Nodes:=0;
      for gen1 in gens1 do
      for gen2 in gens2 do
      for d1 in [1..Order(gen1)-1] do
      for d2 in [1..Order(gen2)-1] do
        if IsConjugate(gr,gen1^d1,gen2^d2) then
           if Order(gen1^d1) ge 3 then Bool:=false; break gen1;
           elif Order(gen1^d1) eq 2
             then Nodes +:= Order(gr)/(2*d1*d2*#Conjugates(gr,gen1^d1));
             if Nodes gt 8-Ksquare then break gen1; end if;
        end if; end if;
      end for; end for; end for; end for;
      if Bool then
        if 8-Ksquare eq Nodes then test:=true; break gens1;
      end if; end if; 
    end for; end for;
return test;
end function;

ExistingNodalSurfaces:=function(Ksquare)
  M:=[* *];
  for triple in ListGroups(Ksquare) do
     G:=triple[3]; T1:=triple[1]; T2:=triple[2]; 
  if CheckSings(G,T1,T2,Ksquare) then Append(~M, [* G,T1,T2 *]);
  end if; end for;
return M;
end function;

AutGr:=function(gr)
  Aut:=AutomorphismGroup(gr); A:={ Aut!1 };
  repeat
    for g1 in Generators(Aut) do
    for g2 in A do
      Include (~A,g1*g2);
    end for; end for;
  until  #A eq #Aut;
return A;
end function;

HurwitzMove:=function(seq,idx)
  return Insert(Remove(seq,idx),idx+1,seq[idx]^seq[idx+1]);
end function;

HurwitzOrbit:=function(seq)
orb:={ }; Purgatory:={ seq };
  repeat
    ExtractRep(~Purgatory,~gens); Include(~orb, gens);
      for k in [1..#seq-1] do hurgens:=HurwitzMove(gens,k);
        if hurgens notin orb then Include(~Purgatory, hurgens);
      end if; end for;
  until IsEmpty(Purgatory);
  orbcut:={  };
  for gens in orb do test:=true;
    for k in [1..#seq-1] do
      if Order(gens[k]) gt Order(gens[k+1]) then test:=false; break k;
      end if;
    end for;
    if test then Include(~orbcut, gens);
    end if;
  end for;
return orbcut;
end function;

CheckSingsEl:=function(gr,seqgrelt1,seqgrelt2,Ksquare)
  Answer:=true; Nodes:=0;
    for g1 in seqgrelt1 do
    for g2 in seqgrelt2 do
    for d1 in [1..Order(g1)-1] do
    for d2 in [1..Order(g2)-1] do
    if IsConjugate(gr,g1^d1,g2^d2) then
    if Order(g1^d1) ge 3 then Answer:=false; break g1;
    elif Order(g1^d1) eq 2 then
      Nodes +:=Order(gr)/(2*d1*d2*#Conjugates(gr,g1^d1)); 
    if Nodes gt 8-Ksquare then Answer:=false; break g1;
    end if; end if; end if; end for; end for; end for; end for;
return Answer and Nodes eq 8-Ksquare;
end function;

SphGens:=function(gr,seq)
  Gens:={ };
  for x1 in ElsOfOrd(gr,seq[1]) do
  for x2 in ElsOfOrd(gr,seq[2]) do
  for x3 in ElsOfOrd(gr,seq[3]) do
  for x4 in ElsOfOrd(gr,seq[4]) do
  if Order(x1*x2*x3*x4) eq seq[5] then
  if sub<gr|x1,x2,x3,x4> eq gr then
    Include(~Gens, [x1,x2,x3,x4,(x1*x2*x3*x4)^-1]);
  end if; end if; end for; end for; end for; end for;
return Gens;
end function;

FindAllComponents:=function(gr,seq1,seq2,Ksquare)
  Comps:={@ @}; Heaven:={ }; Hell:={ }; Aut:=AutGr(gr);
  NumberOfCands:=#SphGens(gr,seq1)*#SphGens(gr,seq2);
  for gen1 in SphGens(gr,seq1) do
  for gen2 in SphGens(gr,seq2) do
    if gen1 cat gen2 notin Hell then
      Include(~Heaven, [gen1,gen2]);
      orb1:=HurwitzOrbit(gen1); orb2:=HurwitzOrbit(gen2);
      for g1 in orb1 do for g2 in orb2 do for phi in Aut do
         Include(~Hell, phi(g1 cat g2));
         if #Hell eq NumberOfCands then break gen1;
         end if;
      end for; end for; end for;
    end if;
  end for; end for;
  for gens in Heaven do
  if CheckSingsEl(gr,gens[1],gens[2],Ksquare) then
    Include(~Comps, gens);
  end if; end for;
return Comps;
end function;

PolyGroup:=func<seq | Group<a1,a2,a3,a4,a5 |
           a1^seq[1], a2^seq[2], a3^seq[3], a4^seq[4], a5^seq[5], a1*a2*a3*a4*a5>>;

Pi1s:=function(gr,seq1,seq2,Ksquare)
  fundgrps:={@ @};
  T1:=PolyGroup(seq1); T2:=PolyGroup(seq2); T1xT2:=DirectProduct(T1,T2);
  in2:=hom< T2->T1xT2 | [T1xT2.6, T1xT2.7, T1xT2.8, T1xT2.9, T1xT2.10]>;
  grxgr,inG:=DirectProduct(gr,gr); m:=NumberOfGenerators(gr);
  L:=[ ]; for i in [1..m] do Append(~L,grxgr.i*grxgr.(i+m)); end for;
  Diag:=hom<gr->grxgr|L>(gr);
  for gen in FindAllComponents(gr,seq1,seq2,Ksquare) do
    f2:=hom< T2->gr |gen[2][1],gen[2][2],gen[2][3],gen[2][4],gen[2][5]>;
    f:=hom<T1xT2->grxgr|
                  inG[1](gen[1][1]),inG[1](gen[1][2]),inG[1](gen[1][3]),
                                    inG[1](gen[1][4]),inG[1](gen[1][5]),
                  inG[2](gen[2][1]),inG[2](gen[2][2]),inG[2](gen[2][3]),
                                    inG[2](gen[2][4]),inG[2](gen[2][5])>;
    H:=Rewrite(T1xT2,Diag@@f); rels:=[ ];
    for i in [1..5] do
    if IsEven(Order(gen[1][i])) then a:=IntegerRing()!(Order(gen[1][i])/2);
    for j in [1..5] do
    if IsEven(Order(gen[2][j])) then b:=IntegerRing()!(Order(gen[2][j])/2);
      test,h:= IsConjugate(gr,gen[1][i]^a,gen[2][j]^b);
      if test then for c in Centralizer(gr,gen[1][i]^a) do
         Append(~rels, T1xT2.i^a * ((T1xT2.(j+5)^b)^(in2((h^-1*c) @@ f2))));
      end for; end if;
    end if; end for; end if; end for;
    Include(~fundgrps, Simplify(quo<H | rels>));
  end for;
return fundgrps;
end function;