/***************************************************
  fns.m -- Functions for ...

  Jennifer Sinnott, 2004
 ***************************************************/

/**********
Here is a collection of functions and procedures that will allow you to 
generate large tables of elliptic curves over a chosen number field with 
class number 1.

To use, make directories "newdata" and "tables" and "isotables".
Then open Magma and load these commands; e.g., 

>load "fns.m";

Then calculate data as follows: 

> m := <irreducible polynomial of your choice>; (e.g., m := x^2+1;)
> B := <positive integer, a "bound" of your choice; (e.g., B := 3;)
> K<a> := NumberField(m);
> calculate(K, B);

Data will be calculated and saved into the directory "newdata".
The "bound" B will be explained later; note that data will be calculated 
up through bound 1 and saved, then from bound 1 to bound 2 and saved, and 
so forth (so if you set B so high that it takes too long to calculate all 
the way up through B, you will still get saved data for lower bounds). 
Then make tables as follows:

> goodtablefile(m, B);;

This will generate a table sorted by norm of the conductor of information 
about the curve, including its torsion subgroup.  This table will be 
stored in "tables". 
One can also generate tables of curves sorted into probable isogeny 
classes by:

>isotablefile(m, B);

This table will be stored in "isotables".
**********/



R<x>:=PolynomialRing(RationalField()); // for convenience.



/**********
I. Useful functions.
**********/



/***
Increment function.  Use for generating strings of coefficients in a given 
numerical range.  Input:  X is a sequence of nonnegative integers.  B is a 
positive integer > 1.  Output:  a sequence got by â€œincrementingâ€� the 
sequence X by 1 in base B.
***/

function inc(X, B)

   assert Type(X) eq SeqEnum;
   assert Type(B) eq RngIntElt;

   if X eq [] then
      return [];
   end if;
   X[#X] +:= 1;
   if X[#X] eq B then
      Y := inc([X[i] : i in [1..#X-1]], B);
      return Y cat [0];
   end if;

   return X;

end function;



/***
Flatten function.  Flattens a sequence of sequences into a single 
sequence.  Input:  x is a sequence of sequences.  Output:  a sequence got 
by concatenating the sequences in x.
***/

function flatten(x)
  
   return &cat x;
  
end function;



/***
Filename maker.  Changes a sequence (e.g. [a, b, c]) to a string 
(â€œa_b_câ This and the function parsename below are used for making and 
understanding filenames when we will later save data.
***/

function makename(seq)
   
   s2 := Sprintf("%o", seq);
   n := #s2; 
   S := "";
   for i in [1..n] do
      if s2[i] eq "[" then  
      elif s2[i] eq "," then
         S cat:= "_";
      elif s2[i] eq " " then
      elif s2[i] eq "]" then
      else 
         S cat:= s2[i];
      end if;
   end for;
   
   return S;

end function;



/***
Filename understand-er.  Takes a file name as generated by makename and 
gives back the integer sequence.
***/

function parsename(str)
   n := #str;
   k := "";  
   for i in [1..n] do
      if str[i] eq "_" then
         k cat:= " ";
      else
         k cat:= str[i];
      end if;
   end for;  
   return StringToIntegerSequence(k);
end function;



/**********
II. Curve listing functions.  These generate elliptic curves with 
coefficients in a certain "bound": we take curves in the form Y^2 = X^3 
+AX + B, where A and B are linear combinations of the power basis ( 1, a, 
a^2, ..., a^(n-1) ) of the field K=Q(a), where a is a root of the 
defining polynomial for K.  Then, curves "bounded" by an integer N will 
be curves whose A and B have coefficients on the power basis in the 
range [-B, B].  We will be generating lists of all curves in a certain 
"bound", and also generating "cleanlists" in which the curves are 
pairwise nonisomorphic over K.
**********/



/***
A basic list function. Generates a list of all elliptic curves with 
coefficients in a certain range Input:  K is a number field.  B is a 
positive integer. Output:  A list of elliptic curves Y^2 = X^3 +aX + b, 
where a and b are linear combinations of the power basis of the field K 
defined by m, with coefficients in the range [-B, B].
Uses: inc function.
***/
   
function list(K, B)
   
   assert Type(K) eq FldNum;
   assert Type(B) eq RngIntElt;

   K<a> := K;
   n := Degree(K);
   L := [];
   A4 := [0 : i in [1..n]];
   A6 := [0 : i in [1..n]];
   s := K![-B : i in [1..n]];
   for i in [1..(2*B+1)^n] do
      a4 := K!(K!A4 + s);  A4 := inc(A4,2*B+1);
      for j in [1..(2*B+1)^n] do
         a6 := K!(K!A6 + s);  A6 := inc(A6,2*B+1);
         if IsEllipticCurve([a4,a6]) then
            Append(~L, EllipticCurve([a4,a6]));
         end if;
      end for;
   end for;
      
   return L;
   
end function;
      
      
   
/***
A basic list extending function.  Generates a list of all elliptic curves 
with coefficients in some range that are not all in a specified smaller 
range. Input:  K is a number field; B & C are integers with B < C.
Output:  A list of elliptic curves Y^2 = X^3 +aX + b, where a and b are 
linear combinations of the integral basis elements with coefficients in 
the range [-C, C] which are not all in the range [-B, B].
Uses: inc function.
***/

function extendlist(K, B, C);

   assert Type(K) eq FldNum;
   assert Type(B) eq RngIntElt;
   assert Type(C) eq RngIntElt;

   K<a> := K;
   L := [];
   n := Degree(K);
   A4 := [0 : i in [1..n]];
   A6 := [0 : i in [1..n]];
   s := K![-C : i in [1..n]];
   for i in [1..(2*C+1)^(n)] do
      f := false;
      for k in [1..n] do
         if A4[k] in [0..C-B-1] or A4[k] in [C+B+1..2*C] then
            f := true;
            break;
         end if;
      end for;
      if f then
         a4 := K!(K!A4 + s);  A4 := inc(A4,2*C+1);
         for j in [1..(2*C+1)^n] do
            a6 := K!(K!A6 + s);  A6 := inc(A6,2*C+1);
            if IsEllipticCurve([a4,a6]) then
               Append(~L, EllipticCurve([a4,a6]));
            end if;
         end for;
      else  
         a4 := K!(K!A4 + s);  A4 := inc(A4,2*C+1);
         for j in [1..(2*C+1)^n] do
            g := false;
            for k in [1..n] do
              if A6[k] in [0..C-B-1] or A6[k] in [C+B+1..2*C] then
                 g := true;
                 break;
              end if;
            end for;
            if g then
               a6 := K!(K!A6 + s);  A6 := inc(A6,2*C+1);
               if IsEllipticCurve([a4,a6]) then
                  Append(~L, EllipticCurve([a4,a6]));
               end if;
            else
               A6 := inc(A6, 2*C+1);
            end if;
         end for;
      end if;
   end for;
   
   return L;
   
end function;
   
   
   
/***
A clean-list function.  Generates a list of all pairwise non-isomorphic 
elliptic curves with coefficients in some range.
Input:  K is a number field.  B is a positive integer.
Output:  A list of elliptic curves Y^2 = X^3 +aX + b, where a and b are
linear combinations of the integral basis elements with coefficients in
the range [-B, B], then throws out a curve if it is isomorphic to one
earlier in the list; returns also a list of their j-invariants.
Uses: list function, inc function.
***/
         
function cleanlist(K, B)
            
   assert Type(K) eq FldNum;
   assert Type(B) eq RngIntElt;
            
   K<a> := K;
   L := list(K, B);
   L2 := [];
   J := [];
   J2 := [];
   for i in [1..#L] do
      j := jInvariant(L[i]);  Append (~J, j);
      t := true;
      if #J gt 1 then  
         for d in [1..#J-1] do
            if J[#J] eq J[d] then
               if IsIsomorphic(L[#J], L[d]) then
                  t:=false;
                  break;
               end if;
            end if;   
         end for;
      end if;
      if t eq true then
         Append(~L2, L[i]);
         Append(~J2, J[i]);
      end if;
   end for;
   
   return <L2, J2>;

end function;
   
   

/***
A clean-list extending function.  Extends a list of nonisomorphic curves 
over a field K up through a bound B and extends it through a bound C.
Input:  A number field K, a list L2 of nonisomorphic curves up 
through B with associated j-invariants C, and extends the list up through 
C.  The resulting list is "clean"; also returns list of their 
j-invariants.  
Uses: extendlist function, inc function.
***/

function extendcleanlist(K, L2, J2, B, C);

   assert Type(K) eq FldNum;
   assert Type(L2) eq SeqEnum;
   assert Type(J2) eq SeqEnum;
   assert Type(B) eq RngIntElt;
   assert Type(C) eq RngIntElt;
   
   K<a> := K;
   S := #L2;
   L3 :=extendlist(K, B, C);
   for i in [1..#L3] do
      j := jInvariant(L3[i]);
      t := true;
      if #J2 ge 1 then
         for d in [1..#J2] do
            if j eq J2[d] then
               if IsIsomorphic(L3[i], L2[d]) then
                  t:=false;
                  break;   
               end if;  
            end if;   
         end for;     
      end if;
      if t eq true then
         Append(~L2, L3[i]);
         Append(~J2, j);   
      end if;
   end for;  
   L4:=[]; 
   J4:=[];
   for i in [S+1..#L2] do
      Append(~L4, L2[i]);
      Append(~J4, J2[i]);
   end for;
   
   return <L4, J4>;

end function;



/**********
III.  Some file saving, loading, etc. procedures.  Things are being saved 
in a 
directory called "newdata"; one can go through and change that name in 
the following four functions so that everything is saved in a different 
directory.
**********/  
   


/***
Save-to-file procedure.  Saves an integer sequence to a file. (Creates or 
overwrites the file given by filename.)
***/
   
procedure save_to_file(seq, filename)
   
   assert Type(filename) eq MonStgElt;
   assert Type(seq) eq SeqEnum;

   file := Open("newdata/"*filename, "w");
   for x in seq do
      fprintf file, "%o ", x;
   end for;

end procedure;



/***
Append-to-file procedure.  This would append a sequence to a file; I 
don't actually use this, but it could be useful so I include it.
***/

procedure append_to_file(seq, filename)
   
   assert Type(filename) eq MonStgElt;
   assert Type(seq) eq SeqEnum;
   
   file := Open("newdata/"*filename, "a");
   for x in seq do
      fprintf file, "%o ", x;
   end for;
   
end procedure;
   


/***
Load-from-file function.  Opens a file and converts it to an integer 
string.  Will get mad if file does not consist of a string of integers.
***/

function load_from_file(filename)
      
   str := Read("newdata/"*filename);
   
   return StringToIntegerSequence(str);
   
end function;
   
   

/***
Checks if a file exists.
***/  

function Exists(path)
   assert Type(path) eq MonStgElt;
   x := Gets(POpen("ls newdata/"*path,"r"));
   return not IsEof(x);
end function;



/**********
IV. Elldata functions.  We are interested in saving "elldata": 

elldata is a 4-tuple < elliptic curve, j-invariant, conductor,
torsion subgroup > where 
  elliptic_curve -- elliptic curve E over a number field K
  j-invariant -- the j-invariant of E
  conductor -- the conductor of E, as an integral ideal of K
  torsion subgp -- an abstract abelian group
   
These functions convert elldata to integer strings so that they can be
saved, or interprete back integer strings into elldata.  Here is where we
will be using that the class number of the field K is 1: Magma's
"Conductor" only works when that assumption is true.
**********/



/***
Elldata to integer sequence conversion.
Input: elldata is a 4-tuple < elliptic curve, j-invariant, conductor,
torsion subgroup > where
Output: an integer list representing the curve, of length 4n+4, where
n is the degree of the field extension.
Uses: flatten function.
***/

function elldata_to_sequence(elldata)

   Z := RingOfIntegers();
   ans := [];
   E, j, N, tor := Explode(elldata);
   K := BaseField(E);
   _,_,_,a,b := Explode(aInvariants(E));
   Append(~ans, Eltseq(a));
   Append(~ans, Eltseq(b));
   Append(~ans, Eltseq(Numerator(j)));
   Append(~ans, [Denominator(j)]);
   truth, gen := IsPrincipal(N);  
   Append(~ans, Eltseq(Numerator(K!gen)));
   Append(~ans, [Denominator(K!gen)]);
   tor := Invariants(tor);
   for i in [1..2-#tor] do
      Append(~tor,[1]);   
   end for;
   Append(~ans, tor);
   f := flatten(ans);
   K := #f;
   f2:= [];
   for i in [1..K] do
      assert Denominator(f[i]) eq 1;
      Append(~f2, Z!f[i]);
   end for;
   return f2;
   
end function;      



/***
Integer sequence to elldata conversion.
Input: K is a number field of degree n.  seq is an integer sequence of
length 4n+4.
Output: elldata 4-tuple.
***/
   
function sequence_to_elldata(K, seq)
      
   assert Type(K) eq FldNum;
   assert Type(seq) eq SeqEnum;
   
   n := Degree(K);
   k := 0; 
   a := K![seq[k+i] : i in [1..n]]; k +:= n;
   b := K![seq[k+i] : i in [1..n]]; k +:= n;
   E := EllipticCurve([a,b]);
   jn := K![seq[k+i] : i in [1..n]]; k +:= n;
   jd := K!seq[k+1]; k +:= 1;
   j := jn/jd;
   cn := K![seq[k+i] : i in [1..n]]; k +:= n;
   cd := K!seq[k+1]; k +:= 1;
   N := (cn/cd)*MaximalOrder(K);
   tor := AbelianGroup([seq[k+i] : i in [1..2]]); k +:= 2;
   assert k eq #seq;

   return <E, j, N, tor>;

end function;

      
         
/***
A data extracting-from-file function.
The calculate function below will take a number field and a "bound" 
and produce files which are lists of curves given in integer 
sequence form (as in the above elldata_to_sequence function); this 
function takes a number field and a file name and converts the 
integer sequence into a list of elldata four-tuples.
Uses: parsename, load_from_file, sequence_to_elldata.
***/
   
function Extract(K, filename)

   R<x> := PolynomialRing(RationalField());
   Z := RingOfIntegers();
   P := parsename(filename);

  assert Degree(K) eq P[1];

   n := Degree(K);

  assert R!DefiningPolynomial(K) eq R!Polynomial([P[i+1] : i in [1..n+1]]);

   List:=[];

   N := 4*n+4;
   datastring := load_from_file(filename);
   l := #datastring/N;
   assert Denominator(l) eq 1;
   l := Z!l;
   for j in [0..l-1] do
      L := [];
      for k in [1..N] do
         Append(~L, datastring[j*N+k]);
      end for;
      s := sequence_to_elldata(K, L);
      Append(~List, s);
   end for;
   return List;

end function;


            
/**********
V.  Sorting procedures: uses quicksort.       
**********/



/***
Procedure underlying quicksort.
***/

procedure QuickSort_recurse(~items, low, high, lessthan)
  // sorts Seqenum using the quick sort algorithm
   if low ge high then
      return; 
   end if;
   if #items le 2 then
      t := items[high];
      items[high] := items[low];
      items[low] := t;
      assert lessthan(items[1],items[2]);
      return;  
   end if;

   lo := low;
   hi := high + 1;
   elem := items[low]; 
   while true do
      while lo lt high and
            lessthan(items[lo+1],elem) do
         lo +:= 1;  
      end while;
      lo +:= 1;
      while lessthan(elem,items[hi-1]) do
         hi -:= 1;
      end while;
      hi -:= 1;
      if lo lt hi then
         t := items[hi];
         items[hi] := items[lo];
         items[lo] := t;
      else
         break;
      end if;
   end while;
   t := items[hi];
   items[hi] := items[low];
   items[low] := t;
   QuickSort_recurse(~items,low,hi-1,lessthan);
   QuickSort_recurse(~items,hi + 1,high,lessthan);
end procedure;  

    

/***
Quicksort procedure that will be called.
***/

procedure QuickSort(~items, lessthan)
   // items    -- a sequence
   // lessthan -- a function that takes two arguments, and returns true
   //             if and only if the first is "less than" the second.  
   QuickSort_recurse(~items,1,#items, lessthan);
end procedure;


               
/***
One lessthan function for use with quicksort.
***/

function lessthan(a, b)
   if Norm(a[3]) lt Norm(b[3]) then
      return true;
   else
      return false;
   end if;   
end function;



/***
Another lessthan function for use with quicksort.
***/

function llessthan(A, B)
   if StringToInteger(A[1]) lt StringToInteger(B[1]) then
      return true;
   end if;   
   if StringToInteger(A[1]) gt StringToInteger(B[1]) then
      return false;
   end if;
   for i in [8..#A] do         // bad is less than everything (but bad)
      if A[i] eq "bad" then   
         if B[i] ne "bad" then
             return true;
         else   
             continue;   
         end if;
      end if;
      a := StringToInteger(A[i]);
      if B[i] eq "bad" then   
         return false;
      end if;
      b := StringToInteger(B[i]);
      if a lt b then
         return true;
      elif a gt b then
         return false;
      end if;
   end for;
   return false;
end function;



/**********
VI.  Calculating and table-making functions and procedures.  
**********/


   
/***
Generates data; checks if it already exists before making it anew; saves 
it to files.  File names are given by, for a field K with defining 
polynomial m:
<n=degree of K> _ <constant term of m> _ <coefficient on x in m> _ 
... _ <coefficient on x^n in m> _ <bound>
where bound is a bound in the sense described in the listing functions, 
and where curves "bounded" by B-1 do not appear in the file saved with 
bound B.
Uses: makename, flatten, Exists, Extract, cleanlist, extendcleanlist, 
elldata_to_sequence, save_to_file.
***/
      
procedure calculate(K, B);
         
   assert Type(K) eq FldNum;
   assert Type(B) eq RngIntElt;

   K<a> := K;
   m := DefiningPolynomial(K);
   n := Degree(K);
   M := [];
   L2 := [];
   J2 := [];

   for i in [1..B] do
      name:= []; 
      Append(~name, [Degree(K)]);
      Append(~name, Coefficients(m));
      Append(~name, [i]);
      Append(~M, makename(flatten(name)));
   end for;
      
   for i in [1..B] do
      if Exists(M[i]) eq true then
         print "found data for size", i;
         Starter := Extract(K, M[i]);
         for k in [1..#Starter] do
            Append(~L2, Starter[k][1]);
            Append(~J2, Starter[k][2]);
         end for;
      else  
         if i eq 1 then
            L := cleanlist(K, 1);
            seq_i := [];  
            Starter := [];
            for j in [1..#L[1]] do
               G:=TorsionSubgroup(L[1][j]);
               C:=Conductor(L[1][j]);
               Info:= <L[1][j], L[2][j], C, G>;
               Append(~Starter, Info);
               seq:=elldata_to_sequence(Info);
               Append(~seq_i, seq);
            end for;
            for k in [1..#Starter] do
               Append(~L2, Starter[k][1]);
               Append(~J2, Starter[k][2]);
            end for;
            finalseq := flatten(seq_i);
            save_to_file(finalseq, M[i]);
            print "done with 1";
         else
            P := extendcleanlist(K, L2, J2, i-1, i);
            seq_i := [];
            Starter := [];
            for j in [1..#P[1]] do
               G := TorsionSubgroup(P[1][j]);
               C := Conductor(P[1][j]);
               Info := <P[1][j], P[2][j], C, G>;
               Append(~Starter, Info);
               seq := elldata_to_sequence(Info);
               Append(~seq_i, seq);
            end for;
            for k in [1..#Starter] do
               Append(~L2, Starter[k][1]);
               Append(~J2, Starter[k][2]);
            end for;
            finalseq := flatten(seq_i);
            save_to_file(finalseq, M[i]);
            print "done with", i;
         end if;
      end if;
   end for;

end procedure;



/***
An initial table-making procedure. Once data has been calculated, this 
accesses it and sorts it.  Data is given as a sequence of <>-sequences, 
each <>-sequence listing information about one curve.
Uses: makename, flatten, Extract, Quicksort with lessthan.
***/

function maketable(K, B);

  assert Type(K) eq FldNum;
  assert Type(B) eq RngIntElt;

   R<x>:=PolynomialRing(Rationals());
   K<a>:=K;
   m:=DefiningPolynomial(K);
   M:=[];

   for i in [1..B] do
      name := [];
      Append(~name, [Degree(K)]); 
      Append(~name, Coefficients(m));  
      Append(~name, [i]);
      Append(~M, makename(flatten(name)));
   end for; 

   for i in [1..B] do
      if Exists(M[i]) eq false then
         return "data has not yet been computed for", i;
      end if;
   end for;

   BigList:=[];
   for i in [1..B] do 
      Append(~BigList, Extract(K, M[i]));
   end for;

   L:=flatten(BigList);
   QuickSort(~L, lessthan);

   LL:=[];
   for i in [1..#L] do
      s0 := Norm(L[i][3]);
      _,_,_,a,b := Explode(aInvariants(L[i][1]));
      s1:=a;
      s2 := b;
      s3 := L[i][2];
      t, s := IsPrincipal(L[i][3]);
      s4 := K!s;
      tor := Invariants(L[i][4]);
      for i in [1..2-#tor] do
         Append(~tor, [1]);
      end for;
      s5 := tor[1];
      s6 := tor[2];
      Append(~LL, <s0, s1, s2, [s5, s6], s5*s6, s3, s4>);
   end for;

   return LL;

end function; 



/***
Makes a table in the style of those on the website.
Uses: maketable.
***/

function goodtable(K, B)
  assert Type(K) eq FldNum;
  assert Type(B) eq RngIntElt;
   K<a>:=K;
   n:=Degree(K);
   m:=DefiningPolynomial(K);
   M:=maketable(K, B);
   list:="";
   length:=12*n-9;
   for i in [1..#M] do 
      f:="";
      s1:=Sprintf("%o", M[i][2]);
      s11:="";
      for j in [1..#s1] do
         if s1[j] eq " " then
         else
            s11 cat:= s1[j];
         end if;
      end for;    
      s2:=Sprintf("%o", M[i][3]);
      s22:="";
      for j in [1..#s2] do
         if s2[j] eq " " then
         else
            s22 cat:= s2[j];
         end if;
      end for;
      length1:=#s11;
      length2:=#s22;
      L:=length1+length2+3;
      for j in [1..#M[1]] do
         s:=Sprintf("%o", M[i][j]);
         if j eq 2 then 
            f cat:= "["*s*",";
         elif j eq 3 then
            gap:=length-L;     
            f cat:=s*"]"*"_"^(gap+1);
         elif j eq 6 then
            if #s le 7 then
               f cat:= s*"\t\t\t";
            elif #s gt 7 and #s le 14 then
               f cat:= s*"\t\t";
            else
               f cat:= s*"\t";
            end if; 
         elif j eq 7 then
            f cat:= "("*s*")"*"\n";
         else
            f cat:= s*"\t";
         end if;
      end for;
      f2:="";
      for j in [1..#f] do
         if f[j] eq " " then
         elif f[j] eq "_" then
            f2 cat:= " ";
         else
            f2 cat:= f[j];
         end if;
      end for;
      list cat:= f2;   
   end for;   
   return list;
end function;



/***
Creates a table in the style of those on the website, and saves it into a 
directory "tables".
Uses: makename, goodtable. 
***/

procedure goodtablefile(m, B)
      R<x>:=PolynomialRing(Rationals());
      m:=R!m;
      K<a>:=NumberField(m);
      n:=Sprintf("%o", Degree(K));
      c:=makename(Coefficients(m));
      b:=Sprintf("%o", B);
      filename := "table_degree_"*n*"_poly_"*c*"_size_"*b*".txt";
      s := goodtable(K,B);
      F := Open("tables/"*filename,"w");
      fprintf F, "%o", s;
end procedure;



/***
This function produces a table in which curves, given by conductor, norm 
of the conductor, and coefficients a and b, are divided up into probable 
isogeny classes, where two curves in the same probable isogeny class both 
with conductor N are listed with the same capital letter (e.g., two 
curves of conductor with norm 16 which are probably isogenous would be 
listed as 16A and 16A; two definitely non-isogenous such curves would be 
listed as 16A and 16B.)
Uses: maketable.
 ***/

function divideup(m, B)

   assert Type(m) eq RngUPolElt;
   assert Type(B) eq RngIntElt;

   Z:=RingOfIntegers;
   K<a>:=NumberField(m);
   n:=Degree(K);
   OK := MaximalOrder(K);
   M:=maketable(K, B);
   M2:=[];
   for i in [1..#M] do
      if M[i][1] le 100000 then
         Append(~M2, M[i]);
      else
         break;
      end if;
   end for;
   primes := [];
   Zprimes := [ p: p in [1..100] | IsPrime(p) eq true];
   for i in [1..#Zprimes] do
      F:=Factorization(Zprimes[i]*OK);
      for j in [1..#F] do
         Append(~primes, F[j][1]);
      end for;
   end for;
   list:=[ ];
   for i in [1..#M2] do
      M3:=[];
      for j in [1..#M2[i]] do
         Append(~M3, Sprintf("%o", M2[i][j]));   
      end for;
      K:=[];
      bad:=BadPlaces(EllipticCurve([M2[i][2], M2[i][3]]));
      for j in [1..#primes] do
         if primes[j] in bad then
            Append(~K, "bad");
         else   
            ER:=Reduction(EllipticCurve([M2[i][2], M2[i][3]]), primes[j]);
            Append(~K, Sprintf("%o", TraceOfFrobenius(ER)));
         end if;
      end for;
      Append(~list, M3 cat K);
   end for;
   QuickSort(~list, llessthan);
   alphabet:=["A", "B", "C", "D", "E", "F", "G", "H", "I", "J", "K", 
"L", "M", "N", "O", "P", "Q", "R", "S", "T", "U", "V", "W", "X", "Y", 
"Z"];
   for i in [1..#list] do
      if i eq 1 then
         k := 1;
         list[i][1] cat:= " "*alphabet[k];
      elif list[i][1]*" "*alphabet[k] eq list[i-1][1] then
         t:=true;
         for j in [8..#list[1]] do
            if list[i][j] eq list[i-1][j] then
            else
               t:=false;
               break;
            end if;
         end for;
         if t eq true then
            list[i][1] cat:= " "*alphabet[k];
         else
            k +:= 1;
            list[i][1] cat:= " "*alphabet[k];
         end if;
      else
         k := 1;
         list[i][1] cat:= " "*alphabet[k];
      end if;
   end for;
   M:=list;
   list:="";
   length:=12*n-9;
   for i in [1..#M] do
      f:="";
      s1:=M[i][2];
      s11:="";
      for j in [1..#s1] do
         if s1[j] eq " " then
         else
            s11 cat:= s1[j];
         end if;
      end for;  
      s2:=M[i][3];
      s22:="";
      for j in [1..#s2] do
         if s2[j] eq " " then
         else
            s22 cat:= s2[j];
         end if;
      end for;
      length1:=#s11;
      length2:=#s22;
      L:=length1+length2+3; 
      for j in [1..#M[1]] do  
         s:=M[i][j];
         if j eq 2 then
            f cat:= "["*s*",";
         elif j eq 3 then
            gap:=length-L;
            f cat:=s*"]"*"_"^(gap+1);
         elif j eq 6 or j eq 5 or j eq 4 then
         elif j eq 7 then
            s1:="";
            for i in [1..#s] do
               if s[i] eq " " then
               else
                  s1 cat:= s[i];
               end if;
            end for;
            if #s1+1 lt 7 then
               f cat:= "("*s*")\t\t";
            else
               f cat:= "("*s*")\t";
            end if;
         elif j eq #M[1] then
            f cat:= s*"\n";
         else
            f cat:= s*"\t";
         end if;
      end for;
      f2:="";
      for j in [1..#f] do
         if f[j] eq " " then
         elif f[j] eq "_" then
            f2 cat:= " ";
         else 
            f2 cat:= f[j];
         end if;
      end for; 
      list cat:= f2;
   end for;
   return list;
end function;



/***
Computes the same table of probable isogeny classes, this time saving to 
a file.
Uses: makename, divideup.
***/

procedure isotablefile(m, B) 
      R<x>:=PolynomialRing(Rationals());
      m:=R!m;   
      K<a>:=NumberField(m);
      n:=Sprintf("%o", Degree(K));
      c:=makename(Coefficients(m));
      b:=Sprintf("%o", B);
      filename := "isotable_degree_"*n*"_poly_"*c*"_size_"*b*".txt";
      s := divideup(m, B);
      F := Open("isotables/"*filename,"w");
      fprintf F, "%o", s;
end procedure;



/**********
Miscellaneous.
**********/


/***
Arbitrary curve generator. Gives you elldata of some curve
over your field (i.e., produces <curve, j-invariant, conductor, torsion
subgroup>).  This is useful to have around to test things out, but is not
used in any later function.
***/

function somecurve(K);

   assert Type(K) eq FldNum;

   K<a> := K;
   n := Degree(K);
   A4 := [Random(-50, 50) : i in [1..n]];
   A6 := [Random(-50, 50) : i in [1..n]];
   a4 := K!(K!A4); a6 := K!(K!A6);
   if IsEllipticCurve([a4, a6]) eq true then
      EC := EllipticCurve([a4, a6]);
      J := jInvariant(EC);
      C := Conductor(EC);
      G := TorsionSubgroup(EC);
      return <EC, J, C, G>;
   else
      return "please try again";
   end if;

end function;
