########################################################################
#
#  Software for "Nielsen equivalence of generating pairs of SL(2,q)"
#
#     by Darryl McCullough and Marcus Wanderley
#
#  Version of August 8, 2011
#
#  Contact:  dmccullough@math.ou.edu
#
#  Written for GAP 4.4.5
#
########################################################################

########################################################################
#
#  The program's output appears on the console and is also written to
#     a file named  outputfile . To change the name of the output file, 
#     alter the String constant in this command:

outputfile := String("outputfile");

#
#
#
########################################################################

########################################################################
#
#  Some routines to try are:
#
#     FindMarkoffClasses( 13 );  finds the Markoff classes when q = 13. 
#        Conjecture C is the assertion that there is one Markoff class
#        for each field element other than 2.
#
#     TestConjectureA( 13 );  will test Conjecture A when q = 13.
#        The algorithm, described in section 9 of the paper, assumes 
#        that the field satisfies Conjecture C. For each possible value
#        of Q, it constructs a graph and then makes random walks in that
#        graph, looking for a certain kind of circuit.  Finding such a 
#        circuit verifies that there is one Nielsen class for that value 
#        of Q. The program may appear to hang when q = 1 mod 4 and the 
#        Q-value is -2, as (consistent with Conjecture A) no such circuit 
#        exists in that case, and the program takes a while before giving 
#        up.  The variables  triesLimit and  maxCircuitLength  in the 
#        routine  FindCircuit  control the duration of the search.
#        triesLimit is the number of random walks that will be taken.
#        maxCircuitLength  initially tells  FindCircuit  how many steps 
#        to take in each walk before giving up in a search for a good 
#        circuit. Once a good circuit has been found, the routine gives 
#        up a walk when the length of the current shortest good circuit 
#        has been exceeded.
#
#        On our desktop machines,  FindMarkoffClasses   and TestConjectureA
#        become impractical around q = 100.
#
#     FindAllNielsenClasses( 5 );  will find the Nielsen classes when
#        q = 5. It directly calculates the sizes of the Aut(F_2)-orbits
#        on generating pairs. On our machines, FindAllNielsenClasses is 
#        already slow when q = 7, and almost impractical for q = 9.
#
########################################################################

########################################################################
#
#  initialization
#
#  Use the routine  SetQ(q)  to reinitialize the order of the field.
#
########################################################################

# initialize global constants

# the order of the field is  q = p^deg
# here we set it to 13, the smallest field in
# the nonexceptional range, but any prime power
# can be used

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;

matrixI      := [ [  one, zero  ], [ zero,  one ] ];
matrixMinusI := [ [ -one, zero  ], [ zero, -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).
#
#    squares  is a list of the FFE's that are squares.
#
#    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();


FindSquares := function()
   local 
      listOfSquares, loopCounter;
 
   listOfSquares := [ zero ];
   for loopCounter in [1..q-1] do
      Add( listOfSquares, u^( 2 * loopCounter ) ); 
   od;
 
   listOfSquares := DuplicateFreeList( listOfSquares );
   return listOfSquares;
   end;

squares := FindSquares();


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();


# 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();
   matrixI      := [ [  one, zero  ], [ zero,  one ] ];
   matrixMinusI := [ [ -one, zero  ], [ zero, -one ] ];
   primeField   := FindPrimeField();
   FFEList      := FindFFEList();
   squares      := FindSquares();
   slElements   := FindSLElements();
   t2           := FindT2();
   t3           := FindT3();
   t4           := FindT4();
   t5           := FindT5();
   end;


# define the Fricke function on triples and the
# trace invariant on pairs

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

TraceInvariant := function( A, B )
   # tr( [ A, B ] )
 
   return Fricke( Trace(A), Trace(B), Trace(A * B) );
   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
   # also, appends the same output to  outputfile
   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;
      if integerForm < 10 then
         AppendTo(outputfile, "  ", integerForm, "  ");
         Print("  ", integerForm, "  ");
      else 
         AppendTo( outputfile, " ", integerForm, "  ");
         Print(" ", integerForm, "  ");
      fi;

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


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


PrintPair := function( pair )
   # for a pair (A, B) of elements of SL(2,q), print A, B, and AB, 
   # also print their traces and tr( [A, B] )

   local 
      A, B, AB;

   A  := pair[1];
   B  := pair[2];
   AB := A * B;

   AppendTo( outputfile, "\n"); 
   Print("\n"); 

   #  print the first row
   PrintVector( A[ 1 ] ); 
   AppendTo( outputfile, "    "); 
   Print("    "); 
   PrintVector( B[ 1 ] );
   AppendTo( outputfile, "    "); 
   Print("    "); 
   PrintVector( AB[ 1 ] );
   AppendTo( outputfile, "   Traces: ");
   Print("   Traces: ");
   PrintVector( [ Trace( A ), 
      Trace( B  ), 
      Trace( AB ) ] );
   AppendTo( outputfile, "\n");
   Print("\n");
   
   #  print the second row
   PrintVector( A[ 2 ] ); 
   AppendTo( outputfile, " ,  "); 
   Print(" ,  "); 
   PrintVector( B[ 2 ] );
   AppendTo( outputfile, " ,  "); 
   Print(" ,  "); 
   PrintVector( AB[ 2 ] );
   AppendTo( outputfile, "   tr( [A, B] ) = ");
   Print("   tr( [A, B] ) = ");
   PrintFFE( TraceInvariant( A, B ));
   AppendTo( outputfile, "\n");
   Print("\n");
   end;


########################################################################
#
# slices
#
#   routine to print out the Q-values on a slice of the space of triples
#     it was used to produce the examples in section 12 of the paper
#
########################################################################


PrintSlice := function( logFormAlpha )
   local 
      qValues, logFormBeta, logFormGamma;

   for logFormGamma in [q-1,q-2..0] do
      qValues := [];
      for logFormBeta in [0..q-1] do
         Add( qValues, Fricke( ExpForm( logFormAlpha ), 
         ExpForm( logFormBeta ), ExpForm( logFormGamma )));
      od;
      PrintVector( qValues );
      Print("\n");
   od;
   end;


########################################################################
#
#  essential_characters
#
#     A sequence of routines to test whether an F_q-triple is 
#        inessential, i. e., is the Fricke trace of a proper subgroup.
#
#     IsProperSubgroupTriple( alpha, beta, gamma )  tells whether
#        ( alpha, beta, gamma ) is inessential.
#
########################################################################


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;


########################################################################
#
# trace_orbits
#
#   routines to calculate the Markoff equivalence classes of traces
#
########################################################################


RR := function( triple )
   return [ triple[1], triple[3], triple[2] ];
   end;

SS := function( triple )
   return [ triple[2], triple[1], triple[3] ];
   end;

TT := function( triple )
   return [ triple[1], triple[2], triple[1] * triple[2] - triple[3] ];
   end;


ListEssentialTriples := function( logFormFFE )    
   # make a list of all essential triples

   local 
      i, j, k, 
      triplesList;  

   triplesList := [ ]; 
   for i in [0..q-1] do
      for j in [0..q-1] do
         for k in [0..q-1] do
            if Fricke( ExpForm(i), ExpForm(j), ExpForm(k) ) = ExpForm( logFormFFE ) and
            (not IsProperSubgroupTriple( ExpForm(i), ExpForm(j), ExpForm(k) ) ) then
               Add( triplesList, [ ExpForm(i), ExpForm(j), ExpForm(k) ] );
            fi;
         od;
      od;
   od;

   return triplesList;
   end;


FindListContainingElement := function( listOfLists, elementOfOneList ) 
   # For a list of lists, and an element of one of the member lists,
   # find the position of the member list that contains the element.

   local 
      memberList;
                       
   for memberList in listOfLists do    
      if elementOfOneList in memberList then 
         return Position( listOfLists, memberList );
      fi; 
   od;
   end;


CoalesceTripleLists := function( triplesList )
   local 
      resultList, triple,
      xMember, imageMember;
   
   # Start with a list of triples closed under RR, SS, and TT,
   # in our case the list of all essential triples with a
   # certain Q-value. 
   # We want to produce a list of the orbits under the action
   # of the group generated by RR, SS, and TT.
   
   # Make a list of lists with each triple made into a one-element
   # member list.
   resultList := [];                          
   for triple in triplesList do                    
      Add( resultList, [ triple ] );   
   od;                                        

   # Go through the input list of triples.
   # If x and RR(x) are not in the same member list, concatenate 
   # the member lists that contain them.                                  

   for triple in triplesList do                        
      xMember := FindListContainingElement( resultList, triple ); 
      imageMember := FindListContainingElement( resultList, RR( triple ) );
      if not( xMember = imageMember ) then   
         Append( resultList[ xMember ], resultList[ imageMember ] ); 
         Unbind( resultList[ imageMember ] );  
      fi;                                    

      # Do the same for SS(x) and TT(x).

      imageMember := FindListContainingElement( resultList, SS( triple ) );
      if not( xMember = imageMember ) then   
         Append( resultList[ xMember], resultList[ imageMember ] );
         Unbind( resultList[ imageMember ] );
      fi;

      imageMember := FindListContainingElement( resultList, TT( triple ) );
      if not( xMember = imageMember ) then
         Append( resultList[ xMember ], resultList[ imageMember ] );
         Unbind( resultList[ imageMember ] );
      fi; 
   od;

   return Compacted(resultList); 
   end;


FindMarkoffClasses := function( fieldOrder )
   local 
      logFormFFE, frickeValue, 
      listOfOrbits, orbit;

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

   if not ( fieldOrder = q ) then
      SetQ( fieldOrder );
   fi;
   AppendTo( outputfile, "\nThe value of q is ", fieldOrder, "\n");
   Print( "\nThe value of q is ", fieldOrder, "\n");

   for logFormFFE in [1..fieldOrder] do

      # if the field has prime order, use zero, one, etc in order
      if IsPrimeInt( fieldOrder ) then
         frickeValue := LogForm( primeField[ logFormFFE ]); 
      else 
         frickeValue := logFormFFE - 1;
      fi;

      listOfOrbits := CoalesceTripleLists( ListEssentialTriples( frickeValue ) );

      AppendTo( outputfile, "\nFor Fricke( alpha, beta, gamma ) =" );
      Print( "\nFor Fricke( alpha, beta, gamma ) =" );
      PrintFFE( ExpForm( frickeValue ) );
      AppendTo( outputfile, ", the number of orbits is ", 
         Size( listOfOrbits ), ".\n" );
      Print( ", the number of orbits is ", Size( listOfOrbits ), ".\n" );

      if not Size( listOfOrbits ) = 0 then
         AppendTo( outputfile, "  The orbits are of length " );
         Print( "  The orbits are of length " );
         for orbit in listOfOrbits do 
            AppendTo( outputfile, " ", Size( orbit ) );
            Print( " ", Size( orbit ) );
         od;
      fi;
   od;
   AppendTo( outputfile, "\n\n" );
   Print( "\n\n" );

   end;

########################################################################
#
# pair_orbits
#
#   Routines to calculate the Nielsen equivalence classes of pairs.
#   
#   The vector of generating pairs with a given Q-value grows rapidly,
#      and these routines only work for very small q  ( q \leq 9 on our
#      desktop machines ).
#
########################################################################


RRPair := function( L )
   return [ Inverse(L[1]), L[1] * L[2] ];
   end;

SSPair := function( L )
   return [ L[2] , L[1] ];
   end;

TTPair := function( L )
   return [ Inverse(L[1]), L[2] ];
   end;

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

GeneratingPairs := function( n )   
   # Return a list of all generating pairs with trace invariant equal
   # to expForm( n )

   local 
      pairsList,
      A, B;

   pairsList := [ ]; 
   for A in slElements do
     for B in slElements do
        if TraceInvariant( A, B ) = ExpForm( n ) and
           IsGeneratingPair( A, B ) then
           Add( pairsList, [ A, B ] );
        fi;
     od;
   od;

   return pairsList;
   end;

CoalescePairLists := function( pairsList )
   # Same routine as CoalesceTripleLists, except for generating pairs.

   local 
      resultList, pair, 
      xMember, imageMember;

   resultList := [];         
   for pair in pairsList do 
       Add( resultList, [ pair ] ); 
   od;                         

   for pair in pairsList do               
      xMember := FindListContainingElement( resultList, pair );    
      imageMember := FindListContainingElement( resultList, RRPair( pair ) );
      if not( xMember = imageMember ) then     
         Append( resultList[ xMember ], resultList[ imageMember ] ); 
         Unbind( resultList[ imageMember ] );       
      fi;                      
      imageMember := FindListContainingElement( resultList, SSPair( pair ) );
      if not( xMember = imageMember ) then     
         Append( resultList[ xMember ], resultList[ imageMember ] );
         Unbind( resultList[ imageMember ] );
      fi;

      imageMember := FindListContainingElement( resultList, TTPair( pair ) );
      if not( xMember = imageMember ) then
         Append( resultList[ xMember ], resultList[ imageMember ] );
         Unbind( resultList[ imageMember ] );
      fi; 
   od;

   return Compacted(resultList); 
   end;

FindNielsenClasses := function( logFormFFE )
   # find the Nielsen classes of generating pairs with
   # Q-value equal to logFormFFE

   local 
      orbitsList, orbit;

   orbitsList := CoalescePairLists( GeneratingPairs( logFormFFE ) );

   Print( "\nFor Fricke( alpha, beta, gamma ) =" );
   PrintFFE( ExpForm( logFormFFE ) );
   Print( ", the number of orbits is ", Size( orbitsList ), ".\n" );

   if not Size( orbitsList ) = 0 then
      Print( "  The orbits are of length " );
      for orbit in orbitsList do 
         Print( " ", Size( orbit ) );
      od;
   fi;

   Print( "\n" );
   end;

FindAllNielsenClasses := function( fieldOrder )
   # find the Nielsen classes of generating pairs when the
   # field has order  fieldOrder
   local
      logFormFFE, qValue;
   if not IsPrimePowerInt( fieldOrder ) then
      Print( "\n", fieldOrder, " is not a prime power.\n" );
      return;
   fi;

   if not ( fieldOrder = q ) then
      SetQ( fieldOrder );
   fi;

   for logFormFFE in [1..fieldOrder] do

      # if the field has prime order, use zero, one, etc in order
      if IsPrimeInt( fieldOrder ) then
         qValue := LogForm( primeField[ logFormFFE ]); 
      else 
         qValue := logFormFFE - 1;
      fi;
      FindNielsenClasses( qValue );
   od;
   end;


########################################################################
#
# circuits
#
#   routines to test Nielsen equivalence by finding 
#      circuits in the graph of pairs
#
########################################################################


FindGeneratingPair := function( x, beta, gamma )
   # find a generating pair ( A, B ) so that
   # A = [ [ x, 0 ], [ 0, x^{-1} ] ] and
   # Fricke trace of (A, B) is ( x + x^{-1}, beta, gamma )

   # B will be of the form [ [ a, 1 ], [ c, d ] ]

   local 
      a, c, d;

   if x = one or x = -one then
      Print( "x (the first parameter) cannot be 1 or -1.\n" );    
      return;
   fi;

   if IsProperSubgroupTriple( x + x^-1, beta, gamma ) then
      Print( "\n( alpha, beta, gamma ) is not an essential triple.\n" );
      return [ matrixI, matrixI ];
   else
      a := ( gamma - beta * (x^-1) ) / ( x - x^-1 );
      c := ( gamma * beta * (x + x^-1) - gamma^2 - beta^2 ) / ( x - x^-1)^2 - one;
      d := ( beta * x - gamma ) / ( x - x^-1 );
      return [ [ [ x, zero ], [ zero, x^-1 ] ], [ [ a, one ], [ c, d ] ] ];
   fi;
   end;

FindConjugater := function( M, x ) 
   # Find  C  so that CMC^{-1} = [ x     0 ]
   #                             [ 0   x^-1]
   # assuming that the trace of  M  is x + x^{-1}
   # with  x  not  0,  1,  or  -1.
 
   # M is [ [ a, b ], [ c, d ] ].
   # Apart the from exceptional cases when a = x or a = x^{-1},
   # C will be [ [ X, Y ], [ Z, one ] ].

   local 
      a, b, c, d,
      X, Y, Z;
 
   if not Trace( M ) = x + x^-1 then 
      Print("\nTrace(M) is not equal to x + x^-1.\n");
      return matrixI;
   fi;
 
   if x = zero or x = one or x = -one then
      Print( "\nx cannot be  0, 1, or -1.\n" );
      return;
   fi;
 
   a := M[1][1];
   b := M[1][2];
   c := M[2][1];
   d := M[2][2];
 
   if a = x and b = zero and c = zero and d = x^-1 then
       return matrixI;
   elif a = x^-1 and b = zero and c = zero and d = x then
       return [ [ zero, one ], [ -one, zero ] ];
   elif a = x and c = zero and not b = zero then
       return [ [b^-1 * (x - x^-1), one], [zero, b * ( x - x^-1)^-1] ];
   elif a = x and b = zero and not c = zero then
       return [ [-c * (x - x^-1)^-1, zero], [one, -c^-1 * ( x - x^-1)] ];
   elif a = x^-1 and not b = zero and c = zero then
       return [ [zero, -b * (x^-1 - x)^-1], [ b^-1 * (x^-1 - x), one] ];
   elif a = x^-1 and b = zero and not c = zero then
       return [ [one, c^-1 * (x - x^-1)], [-c * ( x - x^-1)^-1, zero] ];
   else
      # a is not x and is not x^{-1}, so we can take
      Z := c * (d - x)^-1;
      Y := c^-1 * ( a - x ) * (d - x) * (x - x^-1)^-1;
      X := -c * (a - x)^-1 * Y;
      return [ [ X, Y ], [ Z, one ] ];
   fi;
   end;


IsConjugatePair := function( X, P )
   # Test whether the pair  X  is conjugate to the pair  P  assuming that
   # P[1] is diagonal and has the same trace as  X[1].
 
   # The algorithm is to conjugate X to a pair with  X[1]  diagonal, then 
   # return true or false according to whether or not the product of the 
   # (1,2)-entries of  P[2]  and the new  X[2]  is a square.
 
   local 
      C, conjugatedX2;
 
   C := FindConjugater( X[1], P[1][1][1] );
   conjugatedX2 := C * X[2] * Inverse(C);
 
   return conjugatedX2[1][2] * ( P[2][1][2] ) in squares;
   end;
 
 
FindCircuit := function( startingPair )  
   # Test whether the starting pair is Nielsen equivalent to a
   # nonconjugate pair with the same Fricke trace,

   # The algorithm assumes that the first matrix of 
   # the initial pair is diagonal.

   # The algorithm is to take random walks in the graph of generating 
   # pairs having the same Q-values, where the edges correspond to the 
   # Nielsen moves T, U, and V, and see if we return to a nonconjugate 
   # pair with the same trace.

   local  
      currentRecord, triesMade, 
      triesLimit, maxCircuitLength,
      stepNames, walkHistory, stepCount, 
      currentPair, lastStep, nextStep,
      currentRecordPair, currentRecordWalk, 
      step;
   currentRecordPair := [];
   currentRecordWalk := [];
   maxCircuitLength := 10000;

   # The variable  currentRecord  tells the current shortest loop 
   # that returns to a nonconjugate matrix with the same trace.
   # we initialize it to maxCircuitLength, so that walks will give
   # up in finite time 
   currentRecord := maxCircuitLength;

   # The variable  triesMade  tells how many attempts have already been 
   # made to break the current record.
   # The variable  triesLimit  tells how many tries to make to beat the 
   # current record before giving up.
   triesMade := 0;
   triesLimit := 100;

   # check to make sure the first matrix is not parabolic, since
   # findConjugater assumes it is not
   if Trace( startingPair[1] ) = two or Trace( startingPair[1] ) = -two then
      Print("The first matrix of the starting pair should not be parabolic.\n");
      return;
   fi;

   AppendTo( outputfile, "\n\nFor ( q, Q-value ) = ( ", q, ",");
   Print( "\n\nFor ( q, Q-value ) = ( ", q, ",");
   PrintFFE( TraceInvariant ( startingPair[1], startingPair[2] ));
   AppendTo( outputfile, "), the starting pair is: \n" );
   Print( "), the starting pair is: \n" );
   PrintPair( startingPair ); 
   AppendTo( outputfile, "\n" );
   Print( "\n" );

   # the walk is recorded as a list of characters called walkHistory, whose 
   # ith entry is a Nielsen move  R, S, or T
   stepNames := [ "R", "S", "T" ];

   # each trip through the next repeat loop is a single walk in the graph
   repeat
      stepCount := 0;
      walkHistory := [];
      currentPair := startingPair;
      lastStep := Random( [1, 2, 3] );

      # each trip through the next repeat loop is a single step of
      # the walk
      repeat
         # first choose a random next step different from the last step
         repeat
            nextStep := Random( [ 1, 2, 3 ] );
         until not nextStep = lastStep;

         # take the next step, and update
         if nextStep = 1 then
            currentPair := RRPair( currentPair );
         elif nextStep = 2 then   
            currentPair := SSPair( currentPair );
         else
            currentPair := TTPair( currentPair );
         fi;
         Add( walkHistory, stepNames[ nextStep ] );
         lastStep := nextStep;
         stepCount := stepCount + 1;

         # stop the walk when either it fails--- i. e. the number of steps
         # taken exceeds the current record, or when it succeeds--- it 
         # reaches a pair that has the same Fricke trace as the starting 
         # pair, but is not conjugate to the starting pair:

      until ( stepCount >= currentRecord ) or
         ( (Trace( currentPair[1] ) = Trace( startingPair[1] ) and
            Trace( currentPair[2] ) = Trace( startingPair[2] ) and
            Trace( currentPair[1] * currentPair[2] ) =
            Trace( startingPair[1] * startingPair[2] ) ) and
            ( not IsConjugatePair( currentPair, startingPair ) ) );
      triesMade := triesMade + 1;
      
      # see if current record has been broken, and if so,
      # update the current record information: the new record, 
      # the walk that achieved it, and the ending pair

      if ( stepCount < currentRecord ) then
         currentRecord := stepCount;
         currentRecordPair := currentPair;
         currentRecordWalk := walkHistory;
         triesMade := 0;        
      fi;

   until ( triesMade >= triesLimit );

   # test whether a good circuit was ever found, and report
   if currentRecord = maxCircuitLength then
      if ( (q mod 4) = 1 ) and ( TraceInvariant( startingPair[1], 
            startingPair[2] ) = -two ) then
         AppendTo( outputfile, "No circuit found. This is consistent with Conjecture A");
         AppendTo( outputfile, " since q = 1 mod 4 and Q-value = -2.");
         Print("No circuit found. This is consistent with Conjecture A");
         Print(" since q = 1 mod 4 and Q-value = -2.");
      else    
         AppendTo( outputfile, "No circuit found. Conjecture A not verified.");
         Print( "No circuit found. Conjecture A not verified.");
      fi;
   else
      if currentRecord = 1 then
         AppendTo( outputfile, "The shortest walk found had ", 
         currentRecord, " step: ");
         Print( "The shortest walk found had ", currentRecord, " step: ");
      else
         AppendTo( outputfile, "The shortest walk found had ", 
         currentRecord, " steps: ");
         Print( "The shortest walk found had ", currentRecord, " steps: ");
      fi;
      for step in currentRecordWalk do
         AppendTo( outputfile, step," ");
         Print( step," ");
      od;
      AppendTo( outputfile, "\n\nThe ending pair is:\n");
      Print( "\n\nThe ending pair is:\n");
      PrintPair(currentRecordPair);
   fi;
   end;


########################################################################
#
# triple_circuits
#
#   routines to test Conjecture P by finding paths in
#      the graph of pairs, basically the same routines as
#      for testing Nielsen equivalence
#
########################################################################


FindTripleCircuit := function(aIntFFE, bIntFFE, cIntFFE)
   # same routine as findCircuit, except using the space of triples
   # and seeking a circuit that gets to a weakly Markoff equivalent 
   # but not Markoff equivalent triple

   local 
      alpha, beta, gamma, 
      currentRecord, triesMade, 
      triesLimit, maxCircuitLength,
      stepNames, walkHistory, 
      startingTriple, currentTriple,
      lastStep, nextStep, stepCount, 
      step;

   maxCircuitLength := 10000;
   currentRecord := maxCircuitLength;
   triesMade := 0;
   triesLimit := 2000;
   stepNames := [ "R", "S", "T" ];

   alpha := ExpForm( aIntFFE );
   beta  := ExpForm( bIntFFE );
   gamma := ExpForm( cIntFFE );
   startingTriple := [ alpha, beta, gamma ];
   Print( "\nThe starting triple is: " );
   PrintVector( startingTriple ); 
   Print( "\n" );

   repeat
      stepCount := 0;
      walkHistory := [];
      currentTriple := startingTriple;
      lastStep := Random( [1, 2, 3] );

      # first look for a circuit that returns to (-a, b, -c)
      repeat
         repeat     
            nextStep := Random( [1, 2, 3] );
         until not nextStep = lastStep;
         if nextStep = 1 then
            currentTriple := RR(currentTriple);
         elif nextStep = 2 then   
            currentTriple := SS(currentTriple);
         else
            currentTriple := TT(currentTriple);
         fi;
         Add( walkHistory, stepNames[nextStep] );
         lastStep := nextStep;
         stepCount := stepCount + 1;

      until ( stepCount >= currentRecord ) or
         ( currentTriple[1] = -startingTriple[1] and
         currentTriple[2] = startingTriple[2] and
         currentTriple[3] =  -startingTriple[3] ) ;

      # if not already above currentRecord, look for a circuit 
      # that returns to (a, -b, -c)
      if stepCount < currentRecord then
         repeat
            repeat     
               nextStep := Random( [1, 2, 3] );
            until not nextStep = lastStep;
            if nextStep = 1 then
               currentTriple := RR(currentTriple);
            elif nextStep = 2 then   
               currentTriple := SS(currentTriple);
            else
               currentTriple := TT(currentTriple);
            fi;
            Add(walkHistory, stepNames[nextStep]);
            lastStep := nextStep;
            stepCount := stepCount + 1;
   
         until ( stepCount >= currentRecord ) or
            ( currentTriple[1] = startingTriple[1] and
            currentTriple[2] = -startingTriple[2] and
            currentTriple[3] =  -startingTriple[3] ) ;
      fi;
      triesMade := triesMade + 1;

      # check to see if it's a new record
      if stepCount < currentRecord then 
         currentRecord := stepCount;   
         triesMade := 0;           
         for step in walkHistory do    
            Print(step," ");
         od;
         Print("\n\n");
         PrintVector(currentTriple);
         Print("\n\n",stepCount," steps.\n\n");
      fi;

   until ( triesMade >= triesLimit );

   if currentRecord = maxCircuitLength then
      AppendTo( outputfile, "No circuit found. Conjecture P not verified.");
      Print( "No circuit found. Conjecture P not verified.");
   fi;
   end;


########################################################################
#
#  conjectures
#
#     routines for testing conjectures
#
#     The first routine lists the values of Q for which Conjecture A
#       is trivial.
#
#     The last two routines test Conjectures A and C for a range of
#       field orders.
#
########################################################################


ListTrivialCases := function( fieldOrder )
   # list FFE's for which 2 - x is a square
   # Conjecture A is trivial for these cases

   # this is not used in any other routines,
   # but was used in checking the circuits
   # method

   local
      logFFE, FFE,
      elementList,
      qValue;

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

   if not ( q = fieldOrder ) then
      SetQ( fieldOrder );
   fi;

   elementList := [];
   for logFFE in [0..q-1] do
      FFE := ExpForm( logFFE );
      if ( not( FFE = two) ) and not ( two - FFE in squares ) then
         Add( elementList, logFFE );
      fi;
   od;

   for logFFE in [0..q-1] do
      if IsPrimeInt( q ) then      
         qValue := LogForm( primeField[ logFFE + 1 ] ); 
      else 
         qValue := logFFE;
      fi;
      if qValue in elementList then
         PrintFFE ( ExpForm( qValue ) );
         Print("\n");
      fi;
   od;
   end;


FindEssentialTripleInQLevel := function( qValue )
   # Given an FFE qValue other than 2, find an essential triple with 
   # trace invariant equal to qValue.

   # Use a brute-force method: look in the slice with alpha = u + u^-1
   # and just test triples until one with the correct Fricke trace is 
   # found.

   local
      alpha, logBeta, logGamma;

   if qValue = two then
      Print( "\nNo essential triple can have Q-value = 2.\n" );
      return [];
   fi;

   if ( q < 5 ) then
      Print("\nq must be at least 5.\n");
      return [];
   fi;

   alpha := u + u^-1;
   for logBeta in [0..q-1] do
      for logGamma in [0..q-1] do
         if ( Fricke( alpha, ExpForm( logBeta) , ExpForm( logGamma ) ) =
            qValue ) and
            ( IsProperSubgroupTriple( alpha, ExpForm( logBeta ), ExpForm( logGamma ) ) = 
            false ) then
            return [ alpha, ExpForm(logBeta) , ExpForm(logGamma) ];
         fi;                        
      od;
   od;
   AppendTo( outputfile, "\n\nNo essential triples found with ",
      "(q, Q-value) = ( ", q);
   Print( "\n\nNo essential triples found with (q, Q-value) = ( ", q);
   AppendTo( outputfile, "," );
   Print( "," );
   PrintFFE( qValue );
   AppendTo( outputfile, ")\n" );
   Print(")\n");
   return [];
   end;


TestConjectureA := function( fieldOrder )
   # tests Conjecture A for field of order  fieldOrder  , assuming
   # Conjecture C

   local
      qLogValue, qValue,
      basePair, baseTriple;

   if ( fieldOrder < 5 ) then
      Print("\n", "q must be at least 5.\n");
      return;
   fi;

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

   if ( IsEvenInt( fieldOrder ) ) then
      Print( "\n", fieldOrder, " is not odd.\n" );
      return;
   fi;

   if not ( fieldOrder = q ) then
      SetQ( fieldOrder );
   fi;

   AppendTo( outputfile, "\n\nTesting Conjecture A for q = ", fieldOrder,"\n" );

   for qLogValue in [0..q-1] do

      # if the field has prime order, use zero, one, etc in order
      if IsPrimeInt( fieldOrder ) then
         qValue := primeField[ qLogValue + 1 ]; 
      else 
         qValue := ExpForm( qLogValue );
      fi;

      if ( not ( qValue = two ) ) then
         baseTriple := FindEssentialTripleInQLevel( qValue ); 
         if not baseTriple = [] then
            basePair := FindGeneratingPair(u, baseTriple[ 2 ], baseTriple[ 3 ] );
            FindCircuit( basePair );
         fi;
      fi;
   od;
   end;

TestConjectureAForRange := function( m, n )
   # test Conjecture A for odd prime powers q with 3 \leq m \leq q \leq n, 
   # assuming that Conjecture C is already verified for q
   # the results are written to the screen and to file outputfile

   local
      fieldOrder;

   if ( m > n ) then
      Print("\nStart of range must be less than or equal to end of range.\n");
      return;
   fi;

   for fieldOrder in [m..n] do
      TestConjectureA( fieldOrder );
   od;
   end;

TestConjectureCForRange := function( m, n )
   # test Conjecture C for prime powers q with 3 \leq m \leq q \leq n
   # the results are written to the screen and to file outputfile

   local
      fieldOrder;

   if ( m > n ) then
      Print("\nStart of range must be less than or equal to end of range.\n");
      return;
   fi;

   for fieldOrder in [m..n] do
      FindMarkoffClasses( fieldOrder );
   od;
   end;