CoCalc -- Collaborative Calculation in the Cloud
Sharedwww / eck / oldfns.mOpen in CoCalc
/***************************************************
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
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;

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