/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
AUTC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2 September 2015:
  Added comments and references
  Indented function code
  Removed some unused functions
  Replaced use of Chevalley-Weil formula by Eichler trace formula instead
  Changed output variables from a,b,c,... to x_0, ..., x_(g-1)
  Changed lots of internal variables to make code easier to read
  Set seed to 0 in RunExample functions and FindMatrixGenerators for reproducibility
  Changed output in cyclic trigonal case -- now it only outputs
    polynomials if genus eq 4
  
1 September 2015: Added functions DegreeThreeCharacters, IsPlaneQuintic

3 August 2015: Added function IsMultiplicityFree

28 April 2015: Added functions for calculations related to Streit,
"Homology, Belyi Functions, and Canonical Curves", Section 2

19 January 2014: Fix bug, SKG's need to be distinct elements of the group

20 June 2011: Collecting all my code snippets into one file.

I assume the user computes the group and its character table and classes at the beginning so that the order of this data is fixed once and for all.


*/




/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

EICHLER

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/

// Given a finite group G
// and a signature O
// list all #O-tuples that are surface kernel generators
// with signature O

AllSurfaceKernelGenerators:=function(G,O)
  if #O eq 3 then 
    L1:= { x: x in G | Order(x) eq O[1] };
    L2:= { y: y in G | Order(y) eq O[2]};
    triples:={@ [x,y,(x*y)^-1 ]: x in L1, y in L2 | Order(sub<G |x,y>) eq #G and Order((x*y)^-1) eq O[3] and #Set([x,y,(x*y)^-1 ]) eq 3@}; 
    return triples;
  end if;
  if #O eq 4 then 
    L1:= { x: x in G | Order(x) eq O[1] };
    L2:= { y: y in G | Order(y) eq O[2]};
    L3:= { z: z in G | Order(z) eq O[3]};
    quadruples:={@ [x,y,z,(x*y*z)^-1 ]: x in L1, y in L2, z in L3 | Order(sub<G |x,y,z>) eq #G and Order((x*y*z)^-1) eq O[4] @}; 
    return quadruples;
  end if;
  if #O eq 5 then 
    L1:= { x: x in G | Order(x) eq O[1]};
    L2:= { y: y in G | Order(y) eq O[2]};
    L3:= { z: z in G | Order(z) eq O[3]};
    L4:= { w: w in G | Order(w) eq O[4]};
    quintuples:={@ [x,y,z,w,(x*y*z*w)^-1 ]: x in L1, y in L2, z in L3, w in L4 | Order(sub<G |x,y,z,w>) eq #G and Order((x*y*z*w)^-1) eq O[5] @}; 
    return quintuples;
  end if;
  if #O eq 6 then 
    L1:= { x: x in G | Order(x) eq O[1]};
    L2:= { y: y in G | Order(y) eq O[2]};
    L3:= { z: z in G | Order(z) eq O[3]};
    L4:= { w: w in G | Order(w) eq O[4]};
    L5:= { v: v in G | Order(v) eq O[5]};
    sextuples:={@ [x,y,z,w,v,(x*y*z*w*v)^-1 ]: x in L1, y in L2, z in L3, w in L4, v in L5 | Order(sub<G |x,y,z,w,v>) eq #G and Order((x*y*z*w*v)^-1) eq O[6] @}; 
    return sextuples;
  end if;
end function;


// see Breuer Definition 11.1, page 38 and
// Lemma 11.5, page 39
// FixXuh is the size of the set of fixed points of an automorphism
// h of order $m$ with rotation constant $\zeta_{m}^{-u}$

FixXuh:= function(G,M,u,h)
  m:=Order(h);
  O:=[ Order(M[i]) : i in [1..#M] ];
  k:=Order(Centralizer(G,h));
  sum:=0;
  for i:=1 to #M do 
    if IsDivisibleBy(O[i],m) then 
      if IsConjugate(G,h,M[i]^(Round(O[i]*u/m))) then
        sum:=sum+1/O[i];
      end if;
    end if;
  end for;
  return k*sum;
end function;

//see Breuer Lemma 10.4 page 37

NumberOfFixedPoints:=function(G,M,h)
  N:=Normalizer(G,sub<G |h>);
  m:=Order(h);
  O:=[ Order(M[i]) : i in [1..#M] ];
  k:=Order(N);
  sum:=0;
  for i:=1 to #M do 
    if IsDivisibleBy(O[i],m) then 
      if IsConjugate(G,sub<G | h>,sub<G | M[i]^(Round(O[i]/m))>) then
        sum:=sum+1/O[i];
      end if;
    end if;
  end for;
  return k*sum;
end function;

// We test to see if there is a central involution
// with 2g+2 fixed points

IsHyperelliptic:=function(G,genus,M)
  Z:=Center(G);
  if Order(Z) eq 1 then 
     return false;
  end if;
  L:= {@ z : z in Z | Order(z) eq 2 @}; 
  if #L eq 0 then
    return false;
  end if;
  M:={@ z : z in L | NumberOfFixedPoints(G,M,z) eq 2*genus+2 @};
  if #M eq 0 then 
   return false;
  end if;
  if #M gt 0 then 
    return true,M;
  end if; 
end function;

// We test to see if there is an order 3 automorphism
// with g+2 fixed points

IsCyclicTrigonal:=function(G,genus,M)
  L:= {@ g :g in G | Order(g) eq 3 @}; 
  if #L eq 0 then
    return false;
  end if;
  M:={@ z : z in L | NumberOfFixedPoints(G,M,z) eq genus+2 @};
  if #M eq 0 then 
    return false;
  end if;
  if #M gt 0 then 
    return true,M;
  end if; 
end function;

// This function uses the Eichler trace formula to compute
// the value of the character of the G action on H^0(C,qK) for
// one conjugacy class represented by h
// M is a set of surface kernel generators 

EichlerEntry :=function(G,M,h,q) 
  m:=Order(h);
  z:=RootOfUnity(m);
  nu:=q mod m;
  if q eq 1 then
    sum:=1;
  end if;
  if q gt 1 then 
    sum:=0;
  end if;
  for u:=1 to m do 
    if GCD(u,m) eq 1 then
      sum:=sum+FixXuh(G,M,u,h)*(z^(u*nu))/(1-z^u);
    end if;
  end for;
  return sum;
end function;

// This function assembles the values of the previous function into a
// character
// See Farkas and Kra

Eichler:= function(G,genus,CCL,M,q)
  K<z>:=CyclotomicField(#G);
  if q eq 1 then
    ans:=[K!genus];
  end if;
  if q gt 1 then
    ans:=[K!(2*q-1)*(genus-1)];
  end if;
  for i:=2 to #CCL do   
    ans:=Append(ans,K!EichlerEntry(G,M,CCL[i][3],q));
  end for;
  return CharacterRing(G)!ans;
end function;  

// Given all the data fixed by the program
// this function creates matrix generators for a representation
// of G with the desired character

FindMatrixGenerators:=function(G,genus,T,CCL,M)
  SetSeed(0);
  chi:=Eichler(G,genus,CCL,M,1);
  K<z>:=CyclotomicField(#G);
  n:=#Generators(G);
  i:=1;
  while InnerProduct(T[i],chi) lt 1 do 
    i:=i+1;
  end while;
  ags:=ActionGenerators(GModule(T[i]));
  ags:=[ChangeRing(ags[k],K) : k in [1..n]];
  mats:=ags;
  for j:=1 to Floor(InnerProduct(T[i],chi)-1) do 
    mats:=[DiagonalJoin(mats[k],ags[k]): k in [1..n]];
  end for; 
  for i:=i+1 to #T do 
    if InnerProduct(T[i],chi) ge 1 then
      ags:=ActionGenerators(GModule(T[i]));
      ags:=[ChangeRing(ags[k],K) : k in [1..n]];
      for j:=1 to Floor(InnerProduct(T[i],chi)) do 
        mats:=[DiagonalJoin(mats[k],ags[k]): k in [1..n]];
      end for;
    end if; 
  end for;
  return mats,K,z;
end function;

// Given a list of numbers which are integers but which
// may lie in some other ring
// make it into a vector of integers

ToIntegerVector:=function(L,n)
  return [IntegerRing()!L[i]: i in [1..n]];
end function;

// We compare the decomposition of $I_3$ into irreducibles with
// that of $I_2 \otimes S_1$.  If the multiplicity of any irreducible
// in I_3 exceeds its multiplicity in $I_2 \otimes S_1$, then it is
// clear that the canonical ideal is not generated by quadrics.

IsClearlyNotGeneratedByQuadrics:=function(G,genus,T,CCL,cm,M)
  chiS1:=Eichler(G,genus,CCL,M,1);
  decompS1:=ToIntegerVector(Decomposition(T,chiS1),#CCL);
  decompH0CK:=decompS1;
  decompI1:=[decompS1[i]-decompH0CK[i]: i in [1..#CCL]]; 
  decompS2:=ToIntegerVector(Decomposition(T,Symmetrization(chiS1,[2])),#CCL);
  decompH0C2K:=ToIntegerVector(Decomposition(T,Eichler(G,genus,CCL,M,2)),#CCL);
  decompI2:=[decompS2[i]-decompH0C2K[i]: i in [1..#CCL]];
  decompS3:=ToIntegerVector(Decomposition(T,Symmetrization(chiS1,[3])),#CCL);
  decompH0C3K:=ToIntegerVector(Decomposition(T,Eichler(G,genus,CCL,M,3)),#CCL);
  decompI3:=[decompS3[i]-decompH0C3K[i]: i in [1..#CCL]];
  chiI2:=Symmetrization(chiS1,[2]) - Eichler(G,genus,CCL,M,2);
  decompI2timesS1:=ToIntegerVector(Decomposition(T,chiI2*chiS1),#CCL);
  printf "I_1      =";Sprint(decompI1); 
  printf "S_1      =";Sprint(decompS1);
  printf "H^0(C,K) =";Sprint(decompH0CK);
  printf "I_2      =";Sprint(decompI2);
  printf "S_2      =";Sprint(decompS2);
  printf "H^0(C,2K)=";Sprint(decompH0C2K);
  printf "I_3      =";Sprint(decompI3);
  printf "S_3      =";Sprint(decompS3);
  printf "H^0(C,3K)=";Sprint(decompH0C3K);
  printf "I2timesS1=";Sprint(decompI2timesS1);
  printf "Is clearly not generated by quadrics? ";
  for i:=1 to #CCL do 
    if decompI2timesS1[i] lt decompI3[i] then
      return true,chiI2;
    end if;
  end for;
  return false,chiI2;
end function;



/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PLANE QUINTICS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

*/

// Find all the degree three characters of a finite group
// Input a character table T
// List degree 3 characters, then degree 1+1+1, then degree 1+2

DegreeThreeCharacters:=function(T)
  D1:=[];
  D2:=[];
  D3:=[];
  for i:=1 to #T do
    if T[i][1] eq 1 then 
      D1:=Append(D1,T[i]);
    end if;
  end for;
  for i:=1 to #T do
    if T[i][1] eq 2 then 
      D2:=Append(D2,T[i]);
    end if;
  end for;
  for i:=1 to #T do
    if T[i][1] eq 3 then 
      D3:=Append(D3,T[i]);
    end if;
  end for;
  for i:=1 to #D1 do
    for j:=i to #D1 do
      for k:=j to #D1 do
        D3:=Append(D3,D1[i]+D1[j]+D1[k]);
      end for;
    end for;
  end for;
  for i:=1 to #D1 do
    for j:=1 to #D2 do
      D3:=Append(D3,D1[i]+D2[j]);
    end for;
  end for;
  return D3;
end function;

// Test if the character chi_I_2 of a genus 6 Riemann surface
// is Sym^2(psi) for some degree 3 character psi
// Input a character table T and the character chi

PlaneQuinticObstructed:=function(T,chi)
  L:=DegreeThreeCharacters(T);
  for i:=1 to #L do
    if Symmetrization(L[i],[2]) eq chi then
      return false,L[i];
    end if;
  end for;
  return true,[];
end function;


// Given the character chi_I_2,
// find all degree 3 characters psi such that
// Sym^2(psi) = chi_I_2
PlaneCharacters:=function(T,chiI2)
  L:=DegreeThreeCharacters(T);
  M:=[];
  for i:=1 to #L do 
    if Symmetrization(L[i],[2]) eq chiI2 then 
      M:=Append(M,L[i]);
    end if;
  end for;
  return M;  
end function;

// Given a character psi such that Sym^2 psi = chi_I_2,
// and the character alpha of something in I_3 not in I_2 otimes S_1, 
// find a character theta in Sym^5 psi such that 
// psi*theta = alpha

QuinticCharacters:=function(T,psi)
  sym5psi:=Symmetrization(psi,[5]);
  A:=[];
  for i:=1 to #T do 
    if IntegerRing()!(T[i][1]) eq 1 then 
      D:=Decomposition(T,sym5psi);
      if IntegerRing()!(D[i]) ge 1 then 
          A:=Append(A,i);
      end if;
    end if;
  end for;
  return A;  
end function;

PlaneQuinticMatrixGenerators:=function(G,T,psi)
  SetSeed(0);
  K<z>:=CyclotomicField(#G);
  n:=#Generators(G);
  i:=1;
  while InnerProduct(T[i],psi) lt 1 do 
    i:=i+1;
  end while;
  ags:=ActionGenerators(GModule(T[i]));
  ags:=[ChangeRing(ags[k],K) : k in [1..n]];
  mats:=ags;
  for j:=1 to Floor(InnerProduct(T[i],psi)-1) do 
    mats:=[DiagonalJoin(mats[k],ags[k]): k in [1..n]];
  end for; 
  for i:=i+1 to #T do 
    if InnerProduct(T[i],psi) ge 1 then
      ags:=ActionGenerators(GModule(T[i]));
      ags:=[ChangeRing(ags[k],K) : k in [1..n]];
      for j:=1 to Floor(InnerProduct(T[i],psi)) do 
        mats:=[DiagonalJoin(mats[k],ags[k]): k in [1..n]];
      end for;
    end if; 
  end for;
  return mats,K,z;
end function;

// In the following function, we attempt to analyze the plane quintic
// The character psi such that Sym^2 psi is the action on 
// I_2 was found in "IsPlaneQuintic"
// Note: the equation defining a plane quintic need not be invariant
// It could be a covariant
// (and indeed, this occurs)




/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

DECOMPOSE G ACTION

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/

// Compute the dot product of two lists M and N
Dot := function(M,N)
  n:=#M;
  ans:=0;  
  for i:=1 to n do
     ans:=ans+M[i]*N[i];
  end for;
return ans;
end function;  


// Turn a sequence of length n into a list
ToList:=function(Y,n)
  L:=[];
  for i:=1 to n do
    L:=Append(L,Y[i]);
  end for;
return L;
end function;


// This function writes the basis of the image of the
// matrix M (in the IsotypicalSubspace function below) as polynomials
MyElements:=function(Z,Bd)
  Y:=Basis(Z);
  degZ:=Degree(Z);
  n:=#Y;
  elts:=[];
  for i:=1 to n do
     elts:=Append(elts,Dot(ToList(Y[i],degZ),Bd));
  end for;
return elts;
end function;

// Given a matrix group, tranpose everything
TransposeGroup:=function(G)
  n:=NumberOfRows(G.1);
  K:=CoefficientRing(G.1);
  GLnK:=GeneralLinearGroup(n,K);
  return sub<GLnK | [Transpose(G.i) : i in [1..#Generators(G)]]>;
end function;


MyActionGenerators:=function(Gmod)
  L:=ActionGenerators(Gmod);
  return [ Transpose(L[i]) : i in [1..#L]];
end function;


SquareCoefficient:=function(i,j)
  if i eq j then 
    return 1/2;
  end if;
  return 1;
end function;


CubeCoefficient:=function(i,j,k)
  if #Set([i,j,k]) eq 1 then 
    return 1/6;
  end if;
  if #Set([i,j,k]) eq 2 then 
    return 1/2;
  end if;
  if #Set([i,j,k]) eq 3 then 
    return 1;
  end if;
end function;


FifthPowerCoefficient:=function(i,j,k,l,m)
  S:={*i,j,k,l,m *};
  a:=1;
  for t:=Minimum(S) to Maximum(S) do 
    a:=a/(Factorial(Multiplicity(S,t)));
  end for;
  return a;
end function;

MySymmetrization:=function(M,d)
  n:=NumberOfRows(M);
  if d eq 2 then
    C:=[];
    for i:=1 to n do
      for j:=i to n do
        Cij:=[]; 
        for p:=1 to n do 
          for q:=p to n do 
            Cij:=Append(Cij,SquareCoefficient(p,q)*(M[p][i]*M[q][j]+M[q][i]*M[p][j]));
          end for;
        end for;
        C:=Append(C,Cij); 
      end for;
    end for;
    return Transpose(Matrix(C));
  end if;
  if d eq 3 then
    C:=[];
    for i:=1 to n do
      for j:=i to n do
        for k:=j to n do 
          Cij:=[]; 
          for p:=1 to n do 
            for q:=p to n do 
              for r:=q to n do 
                Cij:=Append(Cij,CubeCoefficient(p,q,r)*(M[p][i]*M[q][j]*M[r][k]+M[p][i]*M[r][j]*M[q][k]+M[q][i]*M[p][j]*M[r][k]+M[q][i]*M[r][j]*M[p][k]+M[r][i]*M[p][j]*M[q][k]+M[r][i]*M[q][j]*M[p][k]));
              end for;
            end for;
          end for;
          C:=Append(C,Cij); 
        end for;
      end for;
    end for;
    return Transpose(Matrix(C));
  end if;
  if d eq 5 then
    C:=[];
    for i:=1 to n do
      for j:=i to n do
        for k:=j to n do 
          for l:=k to n do 
            for m:=l to n do 
              Cij:=[]; 
              for p:=1 to n do 
                for q:=p to n do 
                  for r:=q to n do 
                    for s:=r to n do 
                      for t:=s to n do
                        Cij:=Append(Cij,FifthPowerCoefficient(p,q,r,s,t)*(M[p][i]*M[q][j]*M[r][k]*M[s][l]*M[t][m] + M[p][i]*M[q][j]*M[r][k]*M[t][l]*M[s][m] + M[p][i]*M[q][j]*M[s][k]*M[r][l]*M[t][m] + M[p][i]*M[q][j]*M[s][k]*M[t][l]*M[r][m] + M[p][i]*M[q][j]*M[t][k]*M[r][l]*M[s][m] + M[p][i]*M[q][j]*M[t][k]*M[s][l]*M[r][m] + M[p][i]*M[r][j]*M[q][k]*M[s][l]*M[t][m] + M[p][i]*M[r][j]*M[q][k]*M[t][l]*M[s][m] + M[p][i]*M[r][j]*M[s][k]*M[q][l]*M[t][m] + M[p][i]*M[r][j]*M[s][k]*M[t][l]*M[q][m] + M[p][i]*M[r][j]*M[t][k]*M[q][l]*M[s][m] + M[p][i]*M[r][j]*M[t][k]*M[s][l]*M[q][m] + M[p][i]*M[s][j]*M[q][k]*M[r][l]*M[t][m] + M[p][i]*M[s][j]*M[q][k]*M[t][l]*M[r][m] + M[p][i]*M[s][j]*M[r][k]*M[q][l]*M[t][m] + M[p][i]*M[s][j]*M[r][k]*M[t][l]*M[q][m] + M[p][i]*M[s][j]*M[t][k]*M[q][l]*M[r][m] + M[p][i]*M[s][j]*M[t][k]*M[r][l]*M[q][m] + M[p][i]*M[t][j]*M[q][k]*M[r][l]*M[s][m] + M[p][i]*M[t][j]*M[q][k]*M[s][l]*M[r][m] + M[p][i]*M[t][j]*M[r][k]*M[q][l]*M[s][m] + M[p][i]*M[t][j]*M[r][k]*M[s][l]*M[q][m] + M[p][i]*M[t][j]*M[s][k]*M[q][l]*M[r][m] + M[p][i]*M[t][j]*M[s][k]*M[r][l]*M[q][m] + M[q][i]*M[p][j]*M[r][k]*M[s][l]*M[t][m] + M[q][i]*M[p][j]*M[r][k]*M[t][l]*M[s][m] + M[q][i]*M[p][j]*M[s][k]*M[r][l]*M[t][m] + M[q][i]*M[p][j]*M[s][k]*M[t][l]*M[r][m] + M[q][i]*M[p][j]*M[t][k]*M[r][l]*M[s][m] + M[q][i]*M[p][j]*M[t][k]*M[s][l]*M[r][m] + M[q][i]*M[r][j]*M[p][k]*M[s][l]*M[t][m] + M[q][i]*M[r][j]*M[p][k]*M[t][l]*M[s][m] + M[q][i]*M[r][j]*M[s][k]*M[p][l]*M[t][m] + M[q][i]*M[r][j]*M[s][k]*M[t][l]*M[p][m] + M[q][i]*M[r][j]*M[t][k]*M[p][l]*M[s][m] + M[q][i]*M[r][j]*M[t][k]*M[s][l]*M[p][m] + M[q][i]*M[s][j]*M[p][k]*M[r][l]*M[t][m] + M[q][i]*M[s][j]*M[p][k]*M[t][l]*M[r][m] + M[q][i]*M[s][j]*M[r][k]*M[p][l]*M[t][m] + M[q][i]*M[s][j]*M[r][k]*M[t][l]*M[p][m] + M[q][i]*M[s][j]*M[t][k]*M[p][l]*M[r][m] + M[q][i]*M[s][j]*M[t][k]*M[r][l]*M[p][m] + M[q][i]*M[t][j]*M[p][k]*M[r][l]*M[s][m] + M[q][i]*M[t][j]*M[p][k]*M[s][l]*M[r][m] + M[q][i]*M[t][j]*M[r][k]*M[p][l]*M[s][m] + M[q][i]*M[t][j]*M[r][k]*M[s][l]*M[p][m] + M[q][i]*M[t][j]*M[s][k]*M[p][l]*M[r][m] + M[q][i]*M[t][j]*M[s][k]*M[r][l]*M[p][m] + M[r][i]*M[p][j]*M[q][k]*M[s][l]*M[t][m] + M[r][i]*M[p][j]*M[q][k]*M[t][l]*M[s][m] + M[r][i]*M[p][j]*M[s][k]*M[q][l]*M[t][m] + M[r][i]*M[p][j]*M[s][k]*M[t][l]*M[q][m] + M[r][i]*M[p][j]*M[t][k]*M[q][l]*M[s][m] + M[r][i]*M[p][j]*M[t][k]*M[s][l]*M[q][m] + M[r][i]*M[q][j]*M[p][k]*M[s][l]*M[t][m] + M[r][i]*M[q][j]*M[p][k]*M[t][l]*M[s][m] + M[r][i]*M[q][j]*M[s][k]*M[p][l]*M[t][m] + M[r][i]*M[q][j]*M[s][k]*M[t][l]*M[p][m] + M[r][i]*M[q][j]*M[t][k]*M[p][l]*M[s][m] + M[r][i]*M[q][j]*M[t][k]*M[s][l]*M[p][m] + M[r][i]*M[s][j]*M[p][k]*M[q][l]*M[t][m] + M[r][i]*M[s][j]*M[p][k]*M[t][l]*M[q][m] + M[r][i]*M[s][j]*M[q][k]*M[p][l]*M[t][m] + M[r][i]*M[s][j]*M[q][k]*M[t][l]*M[p][m] + M[r][i]*M[s][j]*M[t][k]*M[p][l]*M[q][m] + M[r][i]*M[s][j]*M[t][k]*M[q][l]*M[p][m] + M[r][i]*M[t][j]*M[p][k]*M[q][l]*M[s][m] + M[r][i]*M[t][j]*M[p][k]*M[s][l]*M[q][m] + M[r][i]*M[t][j]*M[q][k]*M[p][l]*M[s][m] + M[r][i]*M[t][j]*M[q][k]*M[s][l]*M[p][m] + M[r][i]*M[t][j]*M[s][k]*M[p][l]*M[q][m] + M[r][i]*M[t][j]*M[s][k]*M[q][l]*M[p][m] + M[s][i]*M[p][j]*M[q][k]*M[r][l]*M[t][m] + M[s][i]*M[p][j]*M[q][k]*M[t][l]*M[r][m] + M[s][i]*M[p][j]*M[r][k]*M[q][l]*M[t][m] + M[s][i]*M[p][j]*M[r][k]*M[t][l]*M[q][m] + M[s][i]*M[p][j]*M[t][k]*M[q][l]*M[r][m] + M[s][i]*M[p][j]*M[t][k]*M[r][l]*M[q][m] + M[s][i]*M[q][j]*M[p][k]*M[r][l]*M[t][m] + M[s][i]*M[q][j]*M[p][k]*M[t][l]*M[r][m] + M[s][i]*M[q][j]*M[r][k]*M[p][l]*M[t][m] + M[s][i]*M[q][j]*M[r][k]*M[t][l]*M[p][m] + M[s][i]*M[q][j]*M[t][k]*M[p][l]*M[r][m] + M[s][i]*M[q][j]*M[t][k]*M[r][l]*M[p][m] + M[s][i]*M[r][j]*M[p][k]*M[q][l]*M[t][m] + M[s][i]*M[r][j]*M[p][k]*M[t][l]*M[q][m] + M[s][i]*M[r][j]*M[q][k]*M[p][l]*M[t][m] + M[s][i]*M[r][j]*M[q][k]*M[t][l]*M[p][m] + M[s][i]*M[r][j]*M[t][k]*M[p][l]*M[q][m] + M[s][i]*M[r][j]*M[t][k]*M[q][l]*M[p][m] + M[s][i]*M[t][j]*M[p][k]*M[q][l]*M[r][m] + M[s][i]*M[t][j]*M[p][k]*M[r][l]*M[q][m] + M[s][i]*M[t][j]*M[q][k]*M[p][l]*M[r][m] + M[s][i]*M[t][j]*M[q][k]*M[r][l]*M[p][m] + M[s][i]*M[t][j]*M[r][k]*M[p][l]*M[q][m] + M[s][i]*M[t][j]*M[r][k]*M[q][l]*M[p][m] + M[t][i]*M[p][j]*M[q][k]*M[r][l]*M[s][m] + M[t][i]*M[p][j]*M[q][k]*M[s][l]*M[r][m] + M[t][i]*M[p][j]*M[r][k]*M[q][l]*M[s][m] + M[t][i]*M[p][j]*M[r][k]*M[s][l]*M[q][m] + M[t][i]*M[p][j]*M[s][k]*M[q][l]*M[r][m] + M[t][i]*M[p][j]*M[s][k]*M[r][l]*M[q][m] + M[t][i]*M[q][j]*M[p][k]*M[r][l]*M[s][m] + M[t][i]*M[q][j]*M[p][k]*M[s][l]*M[r][m] + M[t][i]*M[q][j]*M[r][k]*M[p][l]*M[s][m] + M[t][i]*M[q][j]*M[r][k]*M[s][l]*M[p][m] + M[t][i]*M[q][j]*M[s][k]*M[p][l]*M[r][m] + M[t][i]*M[q][j]*M[s][k]*M[r][l]*M[p][m] + M[t][i]*M[r][j]*M[p][k]*M[q][l]*M[s][m] + M[t][i]*M[r][j]*M[p][k]*M[s][l]*M[q][m] + M[t][i]*M[r][j]*M[q][k]*M[p][l]*M[s][m] + M[t][i]*M[r][j]*M[q][k]*M[s][l]*M[p][m] + M[t][i]*M[r][j]*M[s][k]*M[p][l]*M[q][m] + M[t][i]*M[r][j]*M[s][k]*M[q][l]*M[p][m] + M[t][i]*M[s][j]*M[p][k]*M[q][l]*M[r][m] + M[t][i]*M[s][j]*M[p][k]*M[r][l]*M[q][m] + M[t][i]*M[s][j]*M[q][k]*M[p][l]*M[r][m] + M[t][i]*M[s][j]*M[q][k]*M[r][l]*M[p][m] + M[t][i]*M[s][j]*M[r][k]*M[p][l]*M[q][m] + M[t][i]*M[s][j]*M[r][k]*M[q][l]*M[p][m]
 ));
                       end for;
                      end for;
                    end for;
                  end for;
                end for;
              C:=Append(C,Cij); 
            end for;
          end for; 
        end for;
      end for;
    end for;
    return Transpose(Matrix(C));
  end if;
end function;

// Given a finite group G, a polynomial ring S, a representation rho
// of G into GL(S_1), a degree d, a character table T, and an index j
// for a row of the character table, compute a basis of the isotypical
// subspace of S_d with character chi_j

IsotypicalSubspace:=function(G,rho,T,S,d,j)
  n:=NumberOfRows(rho(G.1));
  n:=Binomial(n+d-1,d);
  K:=CoefficientRing(S);
  Bd:=[];
  if d eq 2 then
    for p:=1 to Rank(S) do 
      for q:=p to Rank(S) do 
        Bd:=Append(Bd,S.p*S.q);
      end for;
    end for;
  end if;
  if d eq 3 then 
    for p:=1 to Rank(S) do
      for q:=p to Rank(S) do
        for r:=q to Rank(S) do
          Bd:=Append(Bd,S.p*S.q*S.r);
        end for;
      end for;
    end for;
  end if;
  if d eq 5 then 
    for p:=1 to Rank(S) do
      for q:=p to Rank(S) do
        for r:=q to Rank(S) do 
           for s:=r to Rank(S) do
             for t:=s to Rank(S) do
               Bd:=Append(Bd,S.p*S.q*S.r*S.t*S.s);
             end for;
          end for;
        end for;
      end for;
    end for;
  end if;
  nm:=Inverse(NumberingMap(G));
  N:=Order(G);
  M:=ZeroMatrix(K,n,n);
  for i:=1 to N do
    M:=M + ComplexConjugate(T[j](nm(i)))*Transpose(MySymmetrization(rho(nm(i)),d));
  end for;
  printf "CharacterRow ";print j;
  printf "Dimension ";print Dimension(Image(M));
  printf "Multiplicity "; print Dimension(Image(M))/T[j][1];
  return MyElements(Image(M),Bd),M;
end function;


AnalyzePlaneQuintic:=function(G,T,psi)
  answer:=[];
  MatrixGens:=PlaneQuinticMatrixGenerators(G,T,psi);
  K<z>:=CoefficientRing(MatrixGens[1]);
  S<y_0,y_1,y_2>:=PolynomialRing(K,3);
  GL3K:=GeneralLinearGroup(3,K);
  MatrixGens:=[GL3K!MatrixGens[i] : i in [1..#MatrixGens]];
  rho:=hom<G->GL3K | MatrixGens>;
  thetas:=QuinticCharacters(T,psi);
  for j:=1 to #thetas do 
    answer:=Append(answer,IsotypicalSubspace(G,rho,T,S,5,thetas[j]));
  end for;
  return answer;
end function;

// Given a matrix group G acting on a polynomial ring S,
// and two bases of polynomials L1 and L2 that generate
// isomorphic irreducible representations, change the basis
// of L2 so that G acts by the same matrix on L1 and the new
// basis of L2
FindParallelBases:=function(G,L1,L2)
  S:=Parent(L1[1]);
  d:=Degree(L1[1]);
  n:=#L1;
  Degd,F:=GModule(TransposeGroup(G),S,d);
  Sub1:=sub<Degd | [F(L1[i]) : i in [1..n]]>;
  Sub2:=sub<Degd | [F(L2[i]) : i in [1..n]]>;
  K:=CoefficientRing(S); 
  GLnK:=GeneralLinearGroup(n,K);
  W1:=GModule(sub<GLnK | MyActionGenerators(Sub1)>);
  W2:=GModule(sub<GLnK | MyActionGenerators(Sub2)>);
  P:=AHom(W2,W1).1;
  P:=ChangeRing(P,S);
  L2:=Matrix(S,n,1,[L2[i]: i in [1..n]]);
  return Transpose(P)*L2;
end function;


/* Given a group G, genus, and signature E,
compute the classes and character table of G
to fix them once and for all;
find a set of surface kernel generators;
test whether the corresponding curve is hyperelliptic,
cyclic trigonal;
compute the character of G on differentials and the canonical
ideal;
determine if the curve is a plane quintic;
and if the curve is not special in these ways, compute
matrix surface kernel generators and then 
candidate quadrics in the canonical ideal.
*/

RunExample:=function(G,genus,E)
  print "Set seed to 0.";
  SetSeed(0);
  T:=CharacterTable(G);
  print T;
  CCL:=Classes(G);
  print CCL;
  cm:=ClassMap(G);
  AllSKG:=AllSurfaceKernelGenerators(G,E);
  SKG:=AllSKG[1];
  print "Surface kernel generators: ",SKG;
  hyperellipticBoolean:=IsHyperelliptic(G,genus,SKG);
  print "Is hyperelliptic? ",hyperellipticBoolean;
  if hyperellipticBoolean then
    return "Curve is hyperelliptic";
  end if;
  cyclicTrigonalBoolean:=IsCyclicTrigonal(G,genus,SKG);
  print "Is cyclic trigonal? ",cyclicTrigonalBoolean;
  if cyclicTrigonalBoolean then
    if genus gt 7 then 
      return "Curve is cyclic trigonal";
    end if;
  end if;
  print "Multiplicities of irreducibles in relevant G-modules:";
  clearlyNotGeneratedByQuadricsBoolean,chiI2:=IsClearlyNotGeneratedByQuadrics(G,genus,T,CCL,cm,SKG);
  print clearlyNotGeneratedByQuadricsBoolean;
  if genus eq 6 then
    planeQuinticBoolean,psi:=PlaneQuinticObstructed(T,chiI2);
    print "Plane quintic obstruction? ",planeQuinticBoolean;
  end if;
  MatrixGens,K,z:=FindMatrixGenerators(G,genus,T,CCL,SKG);
  print "Matrix generators for action on H^0(C,K):";
  printf "Field K "; print K;
  print MatrixGens;
  GLgK:=GeneralLinearGroup(genus,K);
  rho:=hom<G->GLgK | [<G.i,GLgK!MatrixGens[i]> : i in [1..#Generators(G)]]>;
  print "Matrix Surface Kernel Generators:";
  SetPrintLevel("Magma");
  print [rho(SKG[i]): i in [1..#SKG]];
  SetPrintLevel("Default");
  print "Finding quadrics:";
  if genus eq 4 then
    S<x_0,x_1,x_2,x_3>:=PolynomialRing(K,genus);
  end if;
  if genus eq 5 then
    S<x_0,x_1,x_2,x_3,x_4>:=PolynomialRing(K,genus);
  end if;
  if genus eq 6 then
    S<x_0,x_1,x_2,x_3,x_4,x_5>:=PolynomialRing(K,genus);
  end if;
  if genus eq 7 then
    S<x_0,x_1,x_2,x_3,x_4,x_5,x_6>:=PolynomialRing(K,genus);
  end if;
  if genus eq 8 then
    S<x_0,x_1,x_2,x_3,x_4,x_5,x_6,x_7>:=PolynomialRing(K,genus);
  end if;
  if genus eq 9 then
    S<x_0,x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8>:=PolynomialRing(K,genus);
  end if;
  if genus eq 10 then
    S<x_0,x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9>:=PolynomialRing(K,genus);
  end if;
  chiS1:=Eichler(G,genus,CCL,SKG,1);
  I2:=Decomposition(T,Symmetrization(chiS1,[2])-Eichler(G,genus,CCL,SKG,2));
  I2:=[IntegerRing()!I2[i]: i in [1..#CCL]];
  Quadrics:=[];
  for j:=1 to #CCL do
    if I2[j] ge 1 then
      printf("I2 contains a " cat Sprint(I2[j]*T[j][1]) cat "-dimensional subspace of ");
      IT,projmat:=IsotypicalSubspace(G,rho,T,S,2,j);
      print IT;
      Quadrics:=Append(Quadrics,IT);
    end if;
  end for;
  Cubics:=[];
  if clearlyNotGeneratedByQuadricsBoolean then
    print "Finding cubics:";
    I3:=Decomposition(T,Symmetrization(chiS1,[3])-Eichler(G,genus,CCL,SKG,3));
    I3:=[IntegerRing()!I3[i]: i in [1..#CCL]];
    for j:=1 to #CCL do
      if I3[j] ge 1 then
        printf("I3 contains a " cat Sprint(I3[j]*T[j][1]) cat "-dimensional subspace of ");
        IT:= IsotypicalSubspace(G,rho,T,S,3,j);
        print IT;
        Cubics:=Append(Cubics,IT);
      end if;
    end for;
  return MatrixGens,[rho(SKG[i]): i in [1..#SKG]],Quadrics,Cubics;  
  end if;
  return MatrixGens,[rho(SKG[i]): i in [1..#SKG]],Quadrics,[];
end function;

RunGivenSKG:=function(G,genus,SKG)
  E:=[ Order(SKG[i]): i in [1..#SKG]];
  print "Set seed to 0.";
  SetSeed(0);
  T:=CharacterTable(G);
  print T;
  CCL:=Classes(G);
  print CCL;
  cm:=ClassMap(G);
  print "Surface kernel generators: ",SKG;
  hyperellipticBoolean:=IsHyperelliptic(G,genus,SKG);
  print "Is hyperelliptic? ",hyperellipticBoolean;
  if hyperellipticBoolean then
    return "Curve is hyperelliptic";
  end if;
  cyclicTrigonalBoolean:=IsCyclicTrigonal(G,genus,SKG);
  print "Is cyclic trigonal? ",cyclicTrigonalBoolean;
  if cyclicTrigonalBoolean then
    if genus gt 7 then 
      return "Curve is cyclic trigonal";
    end if;
  end if;
  print "Multiplicities of irreducibles in relevant G-modules:";
  clearlyNotGeneratedByQuadricsBoolean,chiI2:=IsClearlyNotGeneratedByQuadrics(G,genus,T,CCL,cm,SKG);
  print clearlyNotGeneratedByQuadricsBoolean;
  if genus eq 6 then
    planeQuinticBoolean,psi:=PlaneQuinticObstructed(T,chiI2);
    print "Plane quintic obstruction? ",planeQuinticBoolean;
/*    if planeQuinticBoolean eq false then
      return AnalyzePlaneQuintic(G,T,psi);
   end if;*/
  end if;
  MatrixGens,K,z:=FindMatrixGenerators(G,genus,T,CCL,SKG);
  print "Matrix generators for action on H^0(C,K):";
  printf "Field K "; print K;
  print MatrixGens;
  GLgK:=GeneralLinearGroup(genus,K);
  rho:=hom<G->GLgK | [<G.i,GLgK!MatrixGens[i]> : i in [1..#Generators(G)]]>;
  print "Matrix Surface Kernel Generators:";
  SetPrintLevel("Magma");
  print [rho(SKG[i]): i in [1..#SKG]];
  SetPrintLevel("Default");
  print "Finding quadrics:";
  if genus eq 4 then
    S<x_0,x_1,x_2,x_3>:=PolynomialRing(K,genus);
  end if;
  if genus eq 5 then
    S<x_0,x_1,x_2,x_3,x_4>:=PolynomialRing(K,genus);
  end if;
  if genus eq 6 then
    S<x_0,x_1,x_2,x_3,x_4,x_5>:=PolynomialRing(K,genus);
  end if;
  if genus eq 7 then
    S<x_0,x_1,x_2,x_3,x_4,x_5,x_6>:=PolynomialRing(K,genus);
  end if;
  if genus eq 8 then
    S<x_0,x_1,x_2,x_3,x_4,x_5,x_6,x_7>:=PolynomialRing(K,genus);
  end if;
  if genus eq 9 then
    S<x_0,x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8>:=PolynomialRing(K,genus);
  end if;
  if genus eq 10 then
    S<x_0,x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9>:=PolynomialRing(K,genus);
  end if;
  chiS1:=Eichler(G,genus,CCL,SKG,1);
  I2:=Decomposition(T,Symmetrization(chiS1,[2])-Eichler(G,genus,CCL,SKG,2));
  I2:=[IntegerRing()!I2[i]: i in [1..#CCL]];
  Quadrics:=[];
  for j:=1 to #CCL do
    if I2[j] ge 1 then
      printf("I2 contains a " cat Sprint(I2[j]*T[j][1]) cat "-dimensional subspace of ");
      IT,projmat:=IsotypicalSubspace(G,rho,T,S,2,j);
      print IT;
      Quadrics:=Append(Quadrics,IT);
    end if;
  end for;
  Cubics:=[];
  if clearlyNotGeneratedByQuadricsBoolean then
    print "Finding cubics:";
    I3:=Decomposition(T,Symmetrization(chiS1,[3])-Eichler(G,genus,CCL,SKG,3));
    I3:=[IntegerRing()!I3[i]: i in [1..#CCL]];
    for j:=1 to #CCL do
      if I3[j] ge 1 then
        printf("I3 contains a " cat Sprint(I3[j]*T[j][1]) cat "-dimensional subspace of ");
        IT:= IsotypicalSubspace(G,rho,T,S,3,j);
        print IT;
        Cubics:=Append(Cubics,IT);
      end if;
    end for;
  return MatrixGens,[rho(SKG[i]): i in [1..#SKG]],Quadrics,Cubics;  
  end if;
  return MatrixGens,[rho(SKG[i]): i in [1..#SKG]],Quadrics,[];
end function;



// Assume the matrices have entries coercible into a cyclotomic field
// of order Order(G)


RunGivenSKMatrixGenerators:=function(OrdG,genus,MatrixSKG)
  print "Set seed to 0.";
  SetSeed(0);
  K<z>:=CyclotomicField(OrdG);
  GLgK:=GeneralLinearGroup(genus,K);
  G:=sub<GLgK | [GLgK!MatrixSKG[i] : i in [1..#MatrixSKG]]>;
  MatrixSKG:=[G!MatrixSKG[i] : i in [1..#MatrixSKG]];
  E:=[Order(MatrixSKG[i]) : i in [1..#MatrixSKG]];
  T:=CharacterTable(G);
  print T; 
  CCL:=Classes(G);
  print CCL;
  cm:=ClassMap(G);
  rho:=hom<G->GLgK | [<G.i,GLgK!MatrixSKG[i]> : i in [1..#Generators(G)]]>;
  hyperellipticBoolean:=IsHyperelliptic(G,genus,MatrixSKG);
  print "Is hyperelliptic? ",hyperellipticBoolean;
  if hyperellipticBoolean then
    return "Curve is hyperelliptic";
  end if;
  cyclicTrigonalBoolean:=IsCyclicTrigonal(G,genus,MatrixSKG);
  print "Is cyclic trigonal? ",cyclicTrigonalBoolean;
  if cyclicTrigonalBoolean then
    if genus gt 7 then 
      return "Curve is cyclic trigonal";
    end if;
  end if;
  print "Multiplicities of irreducibles in relevant G-modules:";
  clearlyNotGeneratedByQuadricsBoolean,chiI2:=IsClearlyNotGeneratedByQuadrics(G,genus,T,CCL,cm,MatrixSKG);
  print clearlyNotGeneratedByQuadricsBoolean;
  if genus eq 6 then
    planeQuinticBoolean,psi:=PlaneQuinticObstructed(T,chiI2);
    print "Plane quintic obstruction? ",planeQuinticBoolean;
/*    if planeQuinticBoolean eq false then
      return AnalyzePlaneQuintic(G,T,psi);
   end if; */
  end if;
  print "Matrix Surface Kernel Generators:";
  printf "Field K "; print K;
  SetPrintLevel("Magma");
  print [MatrixSKG[i]: i in [1..#MatrixSKG]];
  SetPrintLevel("Default");
  print "Finding quadrics:";
  if genus eq 4 then
    S<x_0,x_1,x_2,x_3>:=PolynomialRing(K,genus);
  end if;
  if genus eq 5 then
    S<x_0,x_1,x_2,x_3,x_4>:=PolynomialRing(K,genus);
  end if;
  if genus eq 6 then
    S<x_0,x_1,x_2,x_3,x_4,x_5>:=PolynomialRing(K,genus);
  end if;
  if genus eq 7 then
    S<x_0,x_1,x_2,x_3,x_4,x_5,x_6>:=PolynomialRing(K,genus);
  end if;
  chiS1:=Eichler(G,genus,CCL,MatrixSKG,1);
  I2:=Decomposition(T,Symmetrization(chiS1,[2])-Eichler(G,genus,CCL,MatrixSKG,2));
  I2:=[IntegerRing()!I2[i]: i in [1..#CCL]];
  Quadrics:=[];
  for j:=1 to #CCL do
    if I2[j] ge 1 then
      printf("I2 contains a " cat Sprint(I2[j]*T[j][1]) cat "-dimensional subspace of ");
      IT,projmat:=IsotypicalSubspace(G,rho,T,S,2,j);
      print IT;
      Quadrics:=Append(Quadrics,IT);
    end if;
  end for;
  Cubics:=[];
  if clearlyNotGeneratedByQuadricsBoolean then
    print "Finding cubics:";
    I3:=Decomposition(T,Symmetrization(chiS1,[3])-Eichler(G,genus,CCL,MatrixSKG,3));
    I3:=[IntegerRing()!I3[i]: i in [1..#CCL]];
    for j:=1 to #CCL do
      if I3[j] ge 1 then
        printf("I3 contains a " cat Sprint(I3[j]*T[j][1]) cat "-dimensional subspace of ");
        IT:= IsotypicalSubspace(G,rho,T,S,3,j);
        print IT;
        Cubics:=Append(Cubics,IT);
      end if;
    end for;
  return MatrixSKG,MatrixSKG,Quadrics,Cubics;  
  end if;
  return MatrixSKG,MatrixSKG,Quadrics,[];
end function;


IsMultiplicityFree:=function(G,genus,E)
  T:=CharacterTable(G);
  CCL:=Classes(G);
  cm:=ClassMap(G);
  AllSKG:=AllSurfaceKernelGenerators(G,E);
  SKG:=AllSKG[1];
  print "Is hyperelliptic? ",IsHyperelliptic(G,genus,SKG);
  if IsHyperelliptic(G,genus,SKG) eq true then
    return "Curve is hyperelliptic";
  end if;
  print "Is cyclic trigonal? ",IsCyclicTrigonal(G,genus,SKG);
  chiS1:=Eichler(G,genus,CCL,SKG,1);
  decompS1:=Decomposition(T,chiS1);
  m:=IntegerRing()!Maximum(decompS1);
  print decompS1;
  print "Is multiplicity free?";
  return m le 1;
end function;



/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

STREIT

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

*/


// Given a finite group G and an element h,
// let H = <h>.  Compute the left cosets Hg.

LeftCosetsByElement:=function(G,h)
  n:=Order(G);
  k:=Order(h);
  f:=NumberingMap(G);
  finv:=Inverse(f);
  S:={@ i : i in [1..n] @};
  for i:=1 to Floor(n/k) do
    g:=finv(S[1]);
    S:=IndexedSetToSet(S);
    Tk:={ G | };
    for t:=1 to k do
      g:=g*h;
      Tk:=Include(Tk,g);
      S:=Exclude(S,f(g));
    end for;
    Tk:=SetToIndexedSet(Tk);
    if i eq 1 then
      T:={@Tk@};
    end if;
    if i ge 2 then 
      T:=Include(T,Tk);
    end if;
    S:=SetToIndexedSet(S);
  end for;
  return SetToIndexedSet(T);
end function;

// Given a finite group G, 
// a numbering map f: G -> \mathbb{N}, and
// an element h generating H=<h>, 
// compute Streit's $W_{G,H}$

WGH:=function(G,f,h)
  n:=Order(G);
  k:=Order(h);
  H:=[h^t: t in [1..k]];
  L:=LeftCosetsByElement(G,h);
  M:=[];
  for i:=1 to #L do 
    M[i]:=[0 : j in [1..Order(G)]];
    rhogWH:=[ f(L[i][1]*H[j]): j in [1..#H]];
    for j:=1 to #rhogWH do
      M[i][rhogWH[j]] :=1;
    end for;
  end for;
  return M;
end function;

// Given a finite group G, 
// a numbering map f: G -> \mathbb{N}, and
// a set M of Belyi surface kernel generators
// compute Streit's $W$

StreitW:=function(G,f,M)
  V:=VectorSpace(RationalField(),Order(G));
  S:=sub<V | WGH(G,f,M[1])>;
  for i:=2 to #M do 
    S:=S+sub<V | WGH(G,f,M[i])>;
  end for;
  return S;
end function;