########################################################################
#
#  Software for "Writing elements of PSL(2,q) as commutators"
#
#     by Darryl McCullough and Marcus Wanderley
#
#  Version of August 8, 2011
#
#  Contact:  dmccullough@math.ou.edu
#
#  Written for GAP 4.4.10
#
########################################################################

########################################################################
#
#  This is just a brute-force calculation of the traces of generating
#  pairs for SL(2,F_q), feasible only for small q (up to the lower 
#  teens).
#
#  q must be at least 3. 
#
#  The function FindTraces( q ) gives a list of the traces.
#
########################################################################

########################################################################
#
#  initialization
#
#  The routine  SetQ(q)  reinitializes the order of the field.
#
########################################################################

# initialize global constants

# the order of the field is  q = p^deg

q := 13;
p := Characteristic(Z(q));
deg := DegreeFFE(Z(q));

u := Z(q);
v := Z(q^2)^(q-1);

zero := 0*Z(q);
one := Z(q)^0;
two := one + one;
three := two + one;
four := three + one;
five := four + one;

#  Define some basic functions for finite field elements ( FFE's ).
#
#    LogForm -- Converts finite field elements from 
#       exponential to logarithmic form.
#
#    ExpForm -- Converts finite field elements from 
#       logarithmic to exponential form.
#
#    FindSquareRoot -- Finds a square root of a FFE (the result
#       may lie in a quadratic extension of the field).


LogForm := function( expFormFFE ) 
   # returns 0 if FFE = zero, and n if FFE = u^n with 0 < n < q
 
   if expFormFFE = 0*Z(q) then
      return 0;
   elif LogFFE( expFormFFE, Z(q) ) = 0 then
      return q-1;
   else 
      return LogFFE( expFormFFE, Z(q) ); 
   fi; 
   end;


ExpForm := function( logFormFFE ) 
   # returns 0*Z(q) if logFormFFE = 0, and Z(q)^n if logFormFFE = n
   # with 0 < n < q
 
   if logFormFFE = 0 then
      return 0*Z(q);
   else  
      return Z(q)^logFormFFE; 
   fi; 
   end;

FindSquareRoot := function( FFE )

   if FFE = zero then
      return zero;
   fi;
 
   if p = 2 then 
      return FFE^( q/2 );
   else 
      return Z(q^2)^( LogForm( FFE ) * (q+1)/2 );
   fi;
   end;


# Initialize lists of field elements and group elements.
#
#    primeField  is a list of FFE's of the prime field in exponential
#       form, ordered as  [ zero, one, two, ... , p-1 ].
#
#    FFEList  is a list of all FFE's in exponential form.  Certain
#       executions are much faster using this list than using a call 
#       to GAP's built-in  GF(q).
#
#    slElements  is a list of the elements of SL(2,q).
#


FindPrimeField := function()  
   # findPrimeField  returns a list of the exponential form of 
   # the prime field elements [ zero, one, ... , p-1 ].
 
   local 
      FFE, primeFieldFFEVector, loopCounter;
 
   FFE := zero;
   primeFieldFFEVector := [];
   for loopCounter in [0..p-1] do
      Add( primeFieldFFEVector, FFE );
      FFE := FFE + one;
   od;

   return primeFieldFFEVector;
   end;

primeField := FindPrimeField();


FindFFEList := function()
   local
      logFFE, FFEList;
 
   FFEList := [];
   for logFFE in [0..q-1] do 
      Add( FFEList, ExpForm( logFFE ) ); 
   od;
 
   return FFEList;
   end;

FFEList := FindFFEList();


FindSLElements := function()
   local 
      listOfSLElements,
      a, b, c, d, adMinusOne;
 
   listOfSLElements := [];
   for a in FFEList do
      for d in FFEList do
         adMinusOne := a*d - one;
         for b in FFEList do
            for c in FFEList do
               if b*c = adMinusOne then
                  Add( listOfSLElements, [ [ a, b ], [ c, d ] ] );
               fi;
            od;
         od;
      od;
   od;
 
   return listOfSLElements;
   end;

slElements := FindSLElements();

SLInverse := function( matrix )
   return [ [ matrix[2][2], -matrix[1][2] ], [ -matrix[2][1], matrix[1][1] ] ];
   end;

# Initialize tn, the set of traces of elements of order n in PSL(2,q),
# for 2 \leq n \leq 5.  # First we must find the roots mu1 and -mu2 of
# the polynomial x^2 - x - 1 #


FindMu1 := function()       
   # mu1 is one root of x^2 - x - 1 = 0
   # -mu2 will be the other
 
   if p = 2 then
      return Z(4); 
   else 
      return ( one + FindSquareRoot( five )) * Inverse( two ) ; 
   fi;
   end;

FindMu2 := function()       
   # mu2 is one root of x^2 + x - 1 = 0, 
   # -mu1 is the other
 
   return FindMu1() - one;
   end;

mu1 := FindMu1();
mu2 := FindMu2();

FindT2 := function()
   return Set( [ zero ] );
   end;

FindT3 := function()
   return Set( [ one, -one ] );
   end;

FindT4 := function()
   if p = 2 then
      return Set( [] );
   else 
      return Set( [ FindSquareRoot( two ), -FindSquareRoot ( two ) ] );
   fi;
   end;

FindT5 := function()
   return Set( [ mu1, -mu1, mu2, -mu2 ] );
   end;

t2 := FindT2();
t3 := FindT3();
t4 := FindT4();
t5 := FindT5();


SetQ := function( n );  
   # reinitialize all constants for a new value of q

   q            := n;  
   p            := Characteristic(Z(q));
   deg          := DegreeFFE(Z(q));
   zero         := 0*Z(q);
   one          := Z(q)^0;
   two          := one + one;
   three        := two + one;
   four         := three + one;
   five         := four + one;
   u            := Z(q);
   v            := Z(q^2)^(q-1);
   mu1          := FindMu1();
   mu2          := FindMu2();
   primeField   := FindPrimeField();
   FFEList      := FindFFEList();
   slElements   := FindSLElements();
   slElements   := FindSLElements();
   t2           := FindT2();
   t3           := FindT3();
   t4           := FindT4();
   t5           := FindT5();
   end;

########################################################################
#
# input-output
#
#   routines for displaying elements, vectors, and pairs of matrices
#
########################################################################


PrintFFE := function( expFormFFE )   
   # print a field element as an integer, if in the prime field, or
   # as u^n, if not
   local 
      integerForm;

   #  if it's in the prime field, look it up in the primeField list
   if DegreeFFE( expFormFFE ) = 1 then
      integerForm := Position( primeField, expFormFFE ) - 1;
      Print(integerForm);

   #  if it's not in the prime field, write it as u^n
   else
      integerForm := LogForm( expFormFFE );
      if integerForm = 1 then
         Print("u");
      else 
         Print("u^", integerForm);
      fi;
   fi;
   end;


PrintVector := function( listOfFFElements )
   local 
      comma, FFE;
   Print("[ ");
   comma := "";
   for FFE in listOfFFElements do
      Print( comma );
      PrintFFE( FFE );
      comma := ", ";
   od;
   Print(" ]");
   end;


########################################################################
#
#  Next we define the function  IsGeneratingPair( A, B ).
#
#  It works with the F_q-triple associated to the matirx pair ( A, B ):
#  (alpha, beta, gamma) = ( Trace(A), Trace(B), Trace(AB) ) .
#
#  (A crude check of whether Order( < A, B, > ) = Order( SL( 2, q ) )
#  is far too inefficient.)
#
########################################################################


Fricke := function( a, b, c ) 
   # the Q-value of an F-triple
 
   return a^2 + b^2 + c^2 - a * b * c - two;
   end;


IsSingularTriple := function( alpha, beta, gamma )

   if Fricke( alpha, beta, gamma ) = two then
      return true;
   else 
      return false;
   fi;
   end;


IsDihedralTriple := function( alpha, beta, gamma )

   if alpha in t2 and beta in t2 then
      return true;
   fi;
   if alpha in t2 and gamma in t2 then 
      return true; 
   fi;
   if beta in t2 and gamma in t2 then 
      return true; 
   fi;
 
   return false;
   end;


IsA4Triple := function( alpha, beta, gamma )

   if alpha in Union( t2, t3 ) 
      and beta in Union( t2, t3 )
      and gamma in Union( t2, t3 )
      and Fricke( alpha, beta, gamma ) = zero then 
      return true; 
   else
      return false;
   fi;
   end;


IsS4Triple := function( alpha, beta, gamma )

   if alpha in Union( t2, t3, t4) 
      and beta in Union( t2, t3, t4)
      and gamma in Union( t2, t3, t4)
      and Fricke( alpha, beta, gamma ) = one then 
      return true; 
   else
      return false;
   fi;
   end;

# the next several functions are auxilary routines for
# the  IsA5Triple  routine which follows them

InterchangeMu1AndMinusMu2 := function(FFE)
   #  A routine needed for one of the steps of the 
   #  normalization process of the  IsA5Triple  routine.
   #  This one interchanges  mu1 <--> -mu2  and  mu2 <--> -mu1 .

   if FFE = mu1 then 
      return -mu2;
   elif FFE = -mu1 then 
      return mu2;
   elif FFE = mu2 then 
      return -mu1;
   elif FFE = -mu2 then 
      return mu1; 
   fi;
 
   return FFE;
   end;

MakeAlphaIntoMu1 := function( alpha, beta, gamma )
   # This is the first technical procedure needed in the 
   # normalization process of the  IsA5Triple  routine.
   # It uses any mu1 or -mu1 in the triple to make alpha = mu1
   # by Nielsen moves.
 
   local
      oldAlpha;
 
   if alpha = -mu1 then 
      alpha := -alpha; 
      beta := -beta;
   fi;
   if beta  = mu1 then 
      oldAlpha := alpha; 
      alpha := beta; 
      beta := oldAlpha;
   fi;
   if beta  = -mu1 then 
      oldAlpha := alpha; 
      alpha := -beta; 
      beta := -oldAlpha;
   fi;
   if gamma = mu1 then 
      oldAlpha := alpha; 
      alpha := gamma; 
      gamma := oldAlpha;
   fi;
   if gamma = -mu1 then 
      oldAlpha := alpha; 
      alpha := -gamma; 
      gamma := -oldAlpha;
   fi;
   return [ alpha, beta, gamma ];
   end;

MakeBetaIntoMuOne := function( alpha, beta, gamma )
   # This is the second technical procedure needed in the 
   # normalization process of the  IsA5Triple  routine.
   # It uses any mu1 or -mu1 in the second or third coordinate
   # to make beta = mu1.

   local
      oldBeta;  
 
   if beta = -mu1 then 
      beta := -beta; 
      gamma := -gamma;
   fi;
   if gamma =  mu1 then 
      oldBeta := beta;  
      beta :=  gamma; 
      gamma := oldBeta;
   fi;
   if gamma = -mu1 then 
      oldBeta := beta; 
      beta := -gamma; 
      gamma := -oldBeta;
   fi;
 
   return [ alpha, beta, gamma ];
   end;

MakeBetaIntoMuTwo := function( alpha, beta, gamma )
   # This is the third technical procedure needed in the 
   # normalization process of the  IsA5Triple  routine.
   # It uses any mu2 or -mu2 in the second or third coordinate
   # to make beta = mu2.
 
   local
      oldBeta;  
 
   if beta = -mu2 then 
      beta := -beta; 
      gamma := -gamma;
   fi;
   if gamma = mu2 then 
      oldBeta := beta; 
      beta := gamma;  
      gamma := oldBeta;
   fi;
   if gamma = -mu2 then 
      oldBeta := beta; 
      beta := -gamma; 
      gamma := -oldBeta;
   fi;
 
   return [ alpha, beta, gamma ];
   end;

MakeBetaIntoOne := function( alpha, beta, gamma )
   # This is the fourth technical procedure needed in the 
   # normalization process of the  IsA5Triple  routine.
   # It uses any 1 or -1 in the second or third coordinate
   # to make beta = 1.

   local
      oldBeta;  
 
   if beta = -one then 
      beta := -beta; 
      gamma := -gamma;
   fi;
   if gamma = one then 
      oldBeta := beta; 
      beta := gamma; 
      gamma := oldBeta;
   fi;
   if gamma = -one then 
      oldBeta := beta; 
      beta := -gamma; 
      gamma := -oldBeta;
   fi;
 
   return [ alpha, beta, gamma ];
   end;

IsA5Triple := function( alpha, beta, gamma )
   local 
      partiallyNormalizedTriple, normalizedTriple;

   # First, put the triple into normalized form for testing.

   # The first step is to use any  mu1, -mu1, mu2, or -mu2  to
   # make  alpha = mu1.
 
   partiallyNormalizedTriple := MakeAlphaIntoMu1( alpha, beta, gamma );
   alpha := partiallyNormalizedTriple[ 1 ];
   beta := partiallyNormalizedTriple[ 2 ];
   gamma := partiallyNormalizedTriple[ 3 ];
 
   # If no mu1 or -mu1 was found, interchange  mu1  and  -mu2  and
   # go through the first step of normalization again.
 
   if not alpha = mu1 then
      alpha := InterchangeMu1AndMinusMu2(alpha);
      beta := InterchangeMu1AndMinusMu2(beta) ;
      gamma := InterchangeMu1AndMinusMu2(gamma);

      partiallyNormalizedTriple := MakeAlphaIntoMu1( alpha, beta, gamma );
      alpha := partiallyNormalizedTriple[ 1 ];
      beta := partiallyNormalizedTriple[ 2 ];
      gamma := partiallyNormalizedTriple[ 3 ];
   fi;
 
   # If there was originally any  mu1, -mu1, mu2, or -mu2, 
   # then  alpha = mu1  now.
 
   # The second step of normalization is to use any mu1 (or -mu1)
   # in the second or third coordinate to make beta = mu1.
 
   partiallyNormalizedTriple := MakeBetaIntoMuOne( alpha, beta, gamma );
   alpha := partiallyNormalizedTriple[ 1 ];
   beta := partiallyNormalizedTriple[ 2 ];
   gamma := partiallyNormalizedTriple[ 3 ];
 
   # If beta has not been changed to a mu1, look for an mu2 or -mu2
   # in the second or third coordinate.
 
   if not beta = mu1 then
      partiallyNormalizedTriple := MakeBetaIntoMuTwo( alpha, beta, gamma );
      alpha := partiallyNormalizedTriple[ 1 ];
      beta := partiallyNormalizedTriple[ 2 ];
      gamma := partiallyNormalizedTriple[ 3 ];
   fi;
 
   # If beta has not been changed to a mu1 or a mu2, look for 
   # a 1 or -1 in the second or third coordinate and try to make
   # beta = 1.
 
   if (not beta = mu1) and (not beta = mu2) then
      partiallyNormalizedTriple := MakeBetaIntoOne( alpha, beta, gamma );
      alpha := partiallyNormalizedTriple[ 1 ];
      beta := partiallyNormalizedTriple[ 2 ];
      gamma := partiallyNormalizedTriple[ 3 ];
   fi;
 
   normalizedTriple := [ alpha, beta, gamma ];
   
   # Now test the normalized triple to see if it is one of
   # the ten A5 cases given in the paper 
   # Exceptional subgroups of SL(2,F).
 
   if normalizedTriple = [ mu1, mu1, one ] then 
       return true;
   fi;
   if normalizedTriple = [ mu1, mu1, mu1 ] then 
       return true;
   fi;
   if normalizedTriple = [ mu1, mu2, zero ] then 
       return true;
   fi;
   if normalizedTriple = [ mu1, mu2, one ] then 
       return true; 
   fi;
   if normalizedTriple = [ mu1, one, zero ] then 
       return true; 
   fi;
   if normalizedTriple = [ mu1, one, one ] then 
       return true; 
   fi;
   if p = 5 and normalizedTriple = [ mu1, mu1, -one ] then 
       return true; 
   fi;
   if p = 5 and normalizedTriple = [ mu1, mu1, zero ] then 
       return true; 
   fi;
   if p = 5 and normalizedTriple = [ mu1, mu2, -one ] then 
       return true; 
   fi;
   if p = 5 and normalizedTriple = [ mu1, mu2, mu2 ] then 
       return true; 
   fi;
   return false;
   end;


IsProperSubfieldTriple := function( alpha, beta, gamma )
   if DegreeFFE( [ alpha, beta, gamma ] ) < deg then
      return true;
   else 
      return false;
   fi;
   end;


IsLinearSubgroupTriple := function(alpha, beta, gamma)
   local 
      pi;
 
   if IsSingularTriple( alpha, beta, gamma )
      or IsDihedralTriple( alpha, beta, gamma ) 
      or IsA4Triple( alpha, beta, gamma )
      or IsS4Triple( alpha, beta, gamma )
      or IsA5Triple( alpha, beta, gamma ) then
      return false;
   fi;

   # first test for an ordinary linear subgroup, or a PGL of 
   # less-than-minimal index

   if IsProperSubfieldTriple( alpha, beta, gamma ) 
      then return true;
   fi;

   # Now test for the PGL-type. It the degree of the PGL-subgroup
   # is not half of the degree of the field, then the entries of the
   # triple lie in a proper subfield, which has already been detected.
   # In the remaining case,  deg  must be even and each of  alpha^2,  
   # beta^2,  gamma^2, and  alpha * beta * gamma  must lie in the 
   # field of order p^(deg/2).
   if ( not p = 2 ) and IsEvenInt( deg ) then
      pi := u^( (p^(deg/2) + 1) / 2 );
      if (IsInt( ( deg/2 )/DegreeFFE( alpha ) ) 
            or IsInt( ( deg/2 )/DegreeFFE( pi * alpha ) ) ) 
         and
         (IsInt( ( deg/2 )/DegreeFFE( beta )) 
            or IsInt( ( deg/2 )/DegreeFFE( pi *  beta ) ) ) 
         and
         (IsInt( ( deg/2 )/DegreeFFE( gamma )) 
            or IsInt(( deg/2 )/DegreeFFE( pi * gamma ) ) ) 
         and
         (IsInt( ( deg/2 )/DegreeFFE( alpha * beta * gamma ) ) ) then
         return true;
      fi;
   fi;

   return false;
   end;
                  

IsProperSubgroupTriple := function( alpha, beta, gamma )

   if q = 3 then
      if IsA4Triple( alpha, beta, gamma ) then
         return false;
      else
         return true;
      fi;
   elif q = 4 or q = 5 then
      if IsA5Triple( alpha, beta, gamma ) then
         return false;
      else
         return true;
      fi;
   else
      if IsSingularTriple( alpha, beta, gamma )
         or IsDihedralTriple( alpha, beta, gamma )
         or IsA4Triple( alpha, beta, gamma ) 
         or IsS4Triple( alpha, beta, gamma ) 
         or IsA5Triple( alpha, beta, gamma ) 
         or IsLinearSubgroupTriple( alpha, beta, gamma ) then
         return true;
      fi;
   fi;
   return false;
   end;

IsGeneratingPair := function( A, B )
   return not IsProperSubgroupTriple( Trace(A), Trace(B), Trace( A * B ) );
   end;

##########################################################################
#
#  And now the main goal...
#
##########################################################################

FindTraces := function( fieldOrder )
   local
      i, j, a, b, tr, 
      length, traceList, genPairs, pair,
      element, outputList;

   if not IsPrimePowerInt( fieldOrder ) then
      Print( "\n", fieldOrder, " is not a prime power.\n" );
      return;
   fi;

   SetQ( fieldOrder );

   traceList := [];

   length := Length( slElements );

   for i in [1..length-1] do
      for j in [i+1..length] do
         a := slElements[i];
         b := slElements[j];
         if IsGeneratingPair( a, b ) then
             tr := Trace( a * b * SLInverse( a ) * SLInverse( b ) );
             if ( not ( tr in traceList ) ) then
                Add( traceList, tr );             
             fi;
         fi;
      od;
   od;

   # pretty print the resulting traces, use primefield order when q = p
   # and multiplicative order otherwise

   outputList := [ ];
   if deg = 1 then
      for element in primeField do
         if element in traceList then
            Add( outputList, element);
         fi;
      od;
   else
      for element in FFEList do
         if element in traceList then
            Add( outputList, element);
         fi;
      od;
   fi;

   PrintVector( outputList );
   Print("\n");
   end;