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 r1 in elementsF do
                         if (not((2*y+gamma) eq Zero(F))) then
                            r2 := (2*x+1)*(2*y+gamma)^(-1)*r1;
                            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
 print "square";
                                        break r1;
                                  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;