/***************************************************
  oldfns.m -- Functions for ...

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

R<x>:=PolynomialRing(RationalField());



/***
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;



/***
Generates a list of all elliptic curves with coefficients in a certain 
range.
Input:  m is a polynomial in R<x> defining 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 of the ring of 
integers 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;
      
      
   
/***
Generates a list of all elliptic curves with coefficients in some range 
that are not all in a specified smaller range.
Input:  m is a polynomial in R<x> defining 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;
   
   
   
/***
Generates a list of all pairwise non-isomorphic elliptic curves with
coefficients in some range.
Input:  m is a polynomial in R<x> defining 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;
   
   

/***
Extends a list of nonisomorphic curves over a field K up through a bound B
and extends it through a bound C.
Input:  A polynomial m in R<x>, 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, increment 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;



/***
Gives you elldata of some curve over your field (i.e., produces <curve, 
j-invariant, conductor, torsion subgroup>).
***/

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;
      
         
         
/***
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;



/***
Input: 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
Output: an integer list representing the curve, of length n*(3+n)+3, where
n is the degree of the field extension.
***/
   
function elldata_to_sequence(elldata)
   
   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)]);
   for b in Basis(N) do
      Append(~ans, Eltseq(K!b));
   end for;
   tor := Invariants(tor);
   for i in [1..2-#tor] do
      Append(~tor,[1]);
   end for;
   Append(~ans, tor);
         
   return flatten(ans);

end function;


      
/***
Input: K is a number field.  seq is an integer sequence.
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;
   B := [];
   for ii in [1..n] do
      b := K![seq[k+i] : i in [1..n]]; k +:= n;
      Append(~B, b);
   end for;
   N := ideal<MaximalOrder(K) | B>;
   tor := AbelianGroup([seq[k+i] : i in [1..2]]); k +:= 2;
   assert k eq #seq;
   
   return <E, j, N, tor>;
   
end function;
   
   
   
   
procedure save_to_file(seq, filename)
   
   assert Type(filename) eq MonStgElt;
   assert Type(seq) eq SeqEnum;

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

end procedure;




procedure append_to_file(seq, filename)
   
   assert Type(filename) eq MonStgElt;
   assert Type(seq) eq SeqEnum;
   
   file := Open("data/"*filename, "a");
   for x in seq do
      fprintf file, "%o ", x;
   end for;
   
end procedure;
   
   
   
   
function load_from_file(filename)
      
   str := Read("data/"*filename);
   
   return StringToIntegerSequence(str);
   
end function;
   
   

/***
Changes a sequence (e.g. [a, b, c]) to a string (“a_b_c”).
***/
      
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;
   
   
   
/***
taking a file name 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;
   


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

function Exists(path)
   assert Type(path) eq MonStgElt;
   x := Gets(POpen("ls data/"*path,"r"));
   return not IsEof(x);
end function;
      
      
         
/***
Input: a polynomial over Q defining a number field, and an integer i
giving the bound we want.
Extracts from a file a list of elliptic curve integral sequence.
***/
   
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 := (n+3)*n+3;
   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;


   
/***
Generates data; checks if it already exists before making it anew.
***/
      
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;
         L2;
         J2;
      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;



/***
a procedure to sort curves in standard form by norm of the conductor

This is bubble sort, which is not very fast when the length of the sequence
is large.  Maybe past in my quicksort, which is quicker.
***/

function sortcurves(seq)
   assert Type(seq) eq SeqEnum;
   for i in [1..#seq] do
      for j in [1..#seq-1] do
         if Norm(seq[j][3]) gt Norm(seq[j+1][3]) then
            Insert(~seq, j, seq[j+1]);
            Remove(~seq, j+2);
         end if;
      end for;
   end for;
return seq;
end function;


/***
quicksort procedures
***/

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

procedure QuickSort_recurse(~items, low, high, lessthan) 
  // sorts Seqenum using the quick sort algorithm
   if low ge high then 
      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;


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;


/***
the ultimate table-making procedure
***/

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, ",please hold 
while it is being computed" ;
         calculate(K, B);
      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; 




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;



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(filename,"w");
      fprintf F, "%o", s;
end procedure;
