action := function(LFT, F,q)	
    

    Fq := {@ f : f in F @};
    perm := [];
    for f in Fq do 
    	denom := LFT[2][1]*f + LFT[2][2];
    	if denom eq 0 then
            Append(~perm, q + 1); 
	else
	    Append(~perm, Position(Fq, (LFT[1][1]*f + LFT[1][2])*denom^-1 ));
        end if;
     end for;
     if LFT[2][1] eq 0 then
	Append(~perm, q + 1);
     else
        Append(~perm, Position(Fq,LFT[1][1]*LFT[2][1]^-1));
     end if;
     return perm;
end function;



FindGstar := function(ABCRList, Q)
    Gstar := [];
    Fq := FiniteField(Q);
    for L in ABCRList do
        if not (L eq []) then
	    A := action(L[1], Fq, Q);
	    B := action(L[2], Fq, Q);
	    R := action(L[4], Fq, Q);
	    Append(~A, Q+2);
	    Append(~A, Q+3);
	    Append(~B, Q+2);
	    Append(~B, Q+3);
	    if IsSquare(Determinant(L[4])) then
		Append(~R, Q+3);
		Append(~R, Q+2);
	    else
		Append(~R, Q+2);
		Append(~R, Q+3);
	    end if;
	    S := SymmetricGroup(Q+3);
	    a := S!A;
	    b := S!B;
	    theta := S!R;
            p := theta*a^-1;
 	    q := theta;
	    r := theta*b;
	    Gs := PermutationGroup< Q+3 | p, q, r>;
       	    Gstar := Append(Gstar,<Gs,p,q,r>);
	end if;
   end for;
   return Gstar;
end function;


Findq := function(n)
// outputs array of first n q's chronologically by p where
// where p is  7 (q=p), p is prime and \pm 1 mod 7 (q=p)
// or where p = \pm 2 or \pm 3 mod 7 (q=p^3).

  total := n;
  counter := 1;
  cubed := false;
  Nextp := 1;
  qList := [];
   
  while (counter le total) do
     Nextp := NextPrime(Nextp);
     if ((Nextp mod 7 eq 1 mod 7) or (Nextp mod 7 eq 6 mod 7) or (Nextp eq 7)) then
        qList := Append(qList, Nextp);
     else 
        qList := Append(qList, (Nextp)^3);
     end if;
     counter := counter + 1;
  end while;

  return(qList);

end function;

NFindq := function(n, N)
// Finds q less than or equal to a specified number N up to n elements

  total := n;
  counter := 1;
  cubed := false;
  Nextp := 1;
  qList := [];
   
  while ((counter le total) and (Nextp lt N)) do
     Nextp := NextPrime(Nextp);
     if (((Nextp mod 7 eq 1 mod 7) or (Nextp mod 7 eq 6 mod 7) or (Nextp eq 7)) and (Nextp le N)) then
        qList := Append(qList, Nextp);
        counter := counter + 1;
     else 
        if (Nextp^3 le N) then  
           qList := Append(qList, (Nextp)^3);
           counter := counter + 1;
        end if;
     end if;
  end while;

  return(qList);

end function;

SortFindq := function(n)
//sorts found q
   qList := Findq(n);
   return(Sort(qList));
end function;

NSortFindq := function(n, N)
//sorts found q less than N, up to n elements.
   qList := NFindq(n, N);
   return(Sort(qList));
end function;

Findqle := function(N)
//finds all q less than or equal to N
  
  counter := 1;
  cubed := false;
  Nextp := 1;
  qList := [];
   
  while (Nextp lt N) do
     Nextp := NextPrime(Nextp);
     if (((Nextp mod 7 eq 1 mod 7) or (Nextp mod 7 eq 6 mod 7) or (Nextp eq 7)) and (Nextp le N)) then
        qList := Append(qList, Nextp);
        counter := counter + 1;
     else 
        if (Nextp^3 le N) then  
           qList := Append(qList, (Nextp)^3);
           counter := counter + 1;
        end if;
     end if;
  end while;

  return(qList);
end function;

SortFindqle := function(N)
//finds all q less than or equal to N
  
  counter := 1;
  cubed := false;
  Nextp := 1;
  qList := [];
   
  while (Nextp lt N) do
     Nextp := NextPrime(Nextp);
     if (((Nextp mod 7 eq 1 mod 7) or (Nextp mod 7 eq 6 mod 7) or (Nextp eq 7)) and (Nextp le N)) then
        qList := Append(qList, Nextp);
        counter := counter + 1;
     else 
        if (Nextp^3 le N) then  
           qList := Append(qList, (Nextp)^3);
           counter := counter + 1;
        end if;
     end if;
  end while;

  return(Sort(qList));
end function;

tr7poly := function(gamma)
// trace 7 polynomial
  return(gamma^3 + gamma^2 - 2*gamma - 1);
end function;

FindABCR := function(q)
// finds A,B,C,R \in SL(2,q) s.t. A has order 4, B has order 3, 
// C has order 7, and R corresponds to the automorphism theta.
// (t-shirt eqns).
// The matrices A, B, and C in SL(2,q) correspond through projectivization 
// to matrices a, b, and c in PSL(2,q) whose orders are 
// |a| = 2, |b| = 3, and |c|=7.  
// 
// This program assumes q is either 7, prime and \pm 1 mod 7,
// or has the form p^3 where p = \pm 2 or \pm 3 mod 7.
//
// returns: 
// ABCList := [ [A1, B1, C1, R1], ... ]  (either 1 or 3 tuples)

  ABCRList := [];

  //counter GenNum for number of A,B,C:
  if (not(q eq 7) and IsPrime(q)) then
     GenCount := 3;
     GenNum := false;
  else 
     GenCount := 1;
     GenNum := true;
  end if;

  for counter in [1..GenCount] do
      Append(~ABCRList, []);
  end for;

  F := FiniteField(q); 
  SL2q := MatrixAlgebra(F,2);

  //choose a conjugation such that 
  A := SL2q![0, 1, -1, 0];

  //choose B
  elementsF := {@ k : k in F @};

  counter := 0;
  for gamma in elementsF do
    if (tr7poly(gamma) eq Zero(F)) then 
       for x in elementsF do
          for y in elementsF do
            B := SL2q![-(x+1), y, y+gamma, x];
            if (Determinant(B) eq One(F)) then
                 if (Order(B) eq 3) then 
                   if (Order(A*B) eq 7) then 
                      C := (A*B)^(-1);

                      for r in elementsF do
                        // Can't this be moved elsewhere? 
			  if (not((2*y+gamma) eq Zero(F))) then
			       r1 := r;	                   	   
			       r2 := (2*x+1)*(2*y+gamma)^(-1)*r1;
			  elif (not ((2*x+1) eq Zero(F))) then
			       r2 := r;
			       r1 := (2*x+1)^-1*(2*y+gamma)*r2;
			  else
			       r1 := Id(SL2q);
			  end if;
			  if not(r1 eq Id(SL2q)) then
                              if (not(-(r1^2 + r2^2) eq Zero(F))) then
                                  R := SL2q![r1, r2, r2, -r1];
                                  if ((R*A eq A^(-1)*R) and (R*B eq B^(-1)*R)) then
                                       ABCRList[counter+1] := [A,B,C,R];
                                       if IsSquare(Determinant(R)) then
                                           break r;
                                       end if;
                                  end if;
			      end if;
                           end if;
                      end for;
                      if ((GenNum) or (counter eq 2)) then 
                         break gamma;
                      else 
                         counter := counter + 1;
                         break x;
                      end if;
                    end if;
                 end if;
              end if;
           end for;
        end for;
     end if;
  end for;
 

   // find R for each generating triple.
   // Let R := SL2q![r1, r2, r3, r4];
   //
   // Define equations for R. 
   //      RA = A^{-1}R  <---> R := [r1, r2, r2, -r1]
   //      RB = B^{-1}R 
   // 
   //  => (2x+1)r1 = (2y+gamma)r2
   // 
   // |R| = 1 so -r1^2 - r2^2 = 1_F
   
   return(ABCRList);   

end function;

PrintABCR := procedure(N, filename) 
// Finds all ABCR tuples for q le N and outputs them in a file.
   Pileofq := SortFindqle(N);
   PileofABCR := [];
   for q in Pileofq do
       ABCR_q := FindABCR(q);
       PrintFile(filename, "******************");
       PrintFile(filename, q);
       PrintFile(filename, "******************");
       PrintFile(filename, ABCR_q);
   end for;
end procedure;

SanityTest:=function(Q)
   results := [];
   Gstar := FindGstar(FindABCR(Q),Q);      
   for num in [1..#Gstar] do
   Gs := Gstar[num][1];                   
   p := Gstar[num][2];
   q := Gstar[num][3];
   r := Gstar[num][4];
   G := sub<Gs | p*q, q*r>;
   if Order(G) eq Order(Gs) then 
      Append(~results,"bad");
   else Append(~results,"good");
   end if;
   end for;
   return results;
end function;

CompileSanityTest := procedure(N)
   Z := IntegerRing();
   index := SortFindq(N);
   for n in index do
	 if IsPrime(n) then 
           p := n;
         else p := Z!(n^(1/3));
         end if;
	 if p mod 7 gt 3 then
	     mo := (p mod 7) - 7;
	 else
	     mo := p mod 7;
         end if;
     	 print "p =",p,":: p =",mo,"mod 7 :: q =",n,"::",SanityTest(n);
   end for;
end procedure;


Gstar := function(Q)
  return FindGstar(FindABCR(Q),Q);
  end function;