Open in CoCalc
1
/***************************************************
2
oldfns.m -- Functions for ...
3
4
Jennifer Sinnott, 2004
5
***************************************************/
6
7
R<x>:=PolynomialRing(RationalField());
8
9
10
11
/***
12
Increment function. Use for generating strings of coefficients in a given
13
numerical range.
14
Input: X is a sequence of nonnegative integers. B is a positive integer
15
> 1.
16
Output: a sequence got by “incrementing” the sequence X by 1 in base
17
B.
18
***/
19
20
function inc(X, B)
21
22
assert Type(X) eq SeqEnum;
23
assert Type(B) eq RngIntElt;
24
25
if X eq [] then
26
return [];
27
end if;
28
X[#X] +:= 1;
29
if X[#X] eq B then
30
Y := inc([X[i] : i in [1..#X-1]], B);
31
return Y cat [0];
32
end if;
33
34
return X;
35
36
end function;
37
38
39
40
/***
41
Generates a list of all elliptic curves with coefficients in a certain
42
range.
43
Input: m is a polynomial in R<x> defining a number field. B is a
44
positive integer.
45
Output: A list of elliptic curves Y^2 = X^3 +aX + b, where a and b are
46
linear combinations of the integral basis elements of the ring of
47
integers of the field K defined by m, with coefficients in
48
the range [-B, B].
49
Uses: inc function.
50
***/
51
52
function list(K, B)
53
54
assert Type(K) eq FldNum;
55
assert Type(B) eq RngIntElt;
56
57
K<a> := K;
58
n := Degree(K);
59
L := [];
60
A4 := [0 : i in [1..n]];
61
A6 := [0 : i in [1..n]];
62
s := K![-B : i in [1..n]];
63
for i in [1..(2*B+1)^n] do
64
a4 := K!(K!A4 + s); A4 := inc(A4,2*B+1);
65
for j in [1..(2*B+1)^n] do
66
a6 := K!(K!A6 + s); A6 := inc(A6,2*B+1);
67
if IsEllipticCurve([a4,a6]) then
68
Append(~L, EllipticCurve([a4,a6]));
69
end if;
70
end for;
71
end for;
72
73
return L;
74
75
end function;
76
77
78
79
/***
80
Generates a list of all elliptic curves with coefficients in some range
81
that are not all in a specified smaller range.
82
Input: m is a polynomial in R<x> defining a number field; B & C are
83
integers with B < C.
84
Output: A list of elliptic curves Y^2 = X^3 +aX + b, where a and b are
85
linear combinations of the integral basis elements with coefficients in
86
the range [-C, C] which are not all in the range [-B, B].
87
Uses: inc function.
88
***/
89
90
function extendlist(K, B, C);
91
92
assert Type(K) eq FldNum;
93
assert Type(B) eq RngIntElt;
94
assert Type(C) eq RngIntElt;
95
96
K<a> := K;
97
L := [];
98
n := Degree(K);
99
A4 := [0 : i in [1..n]];
100
A6 := [0 : i in [1..n]];
101
s := K![-C : i in [1..n]];
102
for i in [1..(2*C+1)^(n)] do
103
f := false;
104
for k in [1..n] do
105
if A4[k] in [0..C-B-1] or A4[k] in [C+B+1..2*C] then
106
f := true;
107
break;
108
end if;
109
end for;
110
if f then
111
a4 := K!(K!A4 + s); A4 := inc(A4,2*C+1);
112
for j in [1..(2*C+1)^n] do
113
a6 := K!(K!A6 + s); A6 := inc(A6,2*C+1);
114
if IsEllipticCurve([a4,a6]) then
115
Append(~L, EllipticCurve([a4,a6]));
116
end if;
117
end for;
118
else
119
a4 := K!(K!A4 + s); A4 := inc(A4,2*C+1);
120
for j in [1..(2*C+1)^n] do
121
g := false;
122
for k in [1..n] do
123
if A6[k] in [0..C-B-1] or A6[k] in [C+B+1..2*C] then
124
g := true;
125
break;
126
end if;
127
end for;
128
if g then
129
a6 := K!(K!A6 + s); A6 := inc(A6,2*C+1);
130
if IsEllipticCurve([a4,a6]) then
131
Append(~L, EllipticCurve([a4,a6]));
132
end if;
133
else
134
A6 := inc(A6, 2*C+1);
135
end if;
136
end for;
137
end if;
138
end for;
139
140
return L;
141
142
end function;
143
144
145
146
/***
147
Generates a list of all pairwise non-isomorphic elliptic curves with
148
coefficients in some range.
149
Input: m is a polynomial in R<x> defining a number field. B is a
150
positive integer.
151
Output: A list of elliptic curves Y^2 = X^3 +aX + b, where a and b are
152
linear combinations of the integral basis elements with coefficients in
153
the range [-B, B], then throws out a curve if it is isomorphic to one
154
earlier in the list; returns also a list of their j-invariants.
155
Uses: list function, inc function.
156
***/
157
158
function cleanlist(K, B)
159
160
assert Type(K) eq FldNum;
161
assert Type(B) eq RngIntElt;
162
163
K<a> := K;
164
L := list(K, B);
165
L2 := [];
166
J := [];
167
J2 := [];
168
for i in [1..#L] do
169
j := jInvariant(L[i]); Append (~J, j);
170
t := true;
171
if #J gt 1 then
172
for d in [1..#J-1] do
173
if J[#J] eq J[d] then
174
if IsIsomorphic(L[#J], L[d]) then
175
t:=false;
176
break;
177
end if;
178
end if;
179
end for;
180
end if;
181
if t eq true then
182
Append(~L2, L[i]);
183
Append(~J2, J[i]);
184
end if;
185
end for;
186
187
return <L2, J2>;
188
189
end function;
190
191
192
193
/***
194
Extends a list of nonisomorphic curves over a field K up through a bound B
195
and extends it through a bound C.
196
Input: A polynomial m in R<x>, a list L2 of nonisomorphic curves up
197
through B
198
with associated j-invariants C, and extends the list up through C. The
199
resulting list is "clean"; also returns list of their j-invariants.
200
Uses: extendlist function, increment function.
201
***/
202
203
function extendcleanlist(K, L2, J2, B, C);
204
205
assert Type(K) eq FldNum;
206
assert Type(L2) eq SeqEnum;
207
assert Type(J2) eq SeqEnum;
208
assert Type(B) eq RngIntElt;
209
assert Type(C) eq RngIntElt;
210
211
K<a> := K;
212
S := #L2;
213
L3 :=extendlist(K, B, C);
214
for i in [1..#L3] do
215
j := jInvariant(L3[i]);
216
t := true;
217
if #J2 ge 1 then
218
for d in [1..#J2] do
219
if j eq J2[d] then
220
if IsIsomorphic(L3[i], L2[d]) then
221
t:=false;
222
break;
223
end if;
224
end if;
225
end for;
226
end if;
227
if t eq true then
228
Append(~L2, L3[i]);
229
Append(~J2, j);
230
end if;
231
end for;
232
L4:=[];
233
J4:=[];
234
for i in [S+1..#L2] do
235
Append(~L4, L2[i]);
236
Append(~J4, J2[i]);
237
end for;
238
239
return <L4, J4>;
240
241
end function;
242
243
244
245
/***
246
Gives you elldata of some curve over your field (i.e., produces <curve,
247
j-invariant, conductor, torsion subgroup>).
248
***/
249
250
function somecurve(K);
251
252
assert Type(K) eq FldNum;
253
254
K<a> := K;
255
n := Degree(K);
256
A4 := [Random(-50, 50) : i in [1..n]];
257
A6 := [Random(-50, 50) : i in [1..n]];
258
a4 := K!(K!A4); a6 := K!(K!A6);
259
if IsEllipticCurve([a4, a6]) eq true then
260
EC := EllipticCurve([a4, a6]);
261
J := jInvariant(EC);
262
C := Conductor(EC);
263
G := TorsionSubgroup(EC);
264
return <EC, J, C, G>;
265
else
266
return "please try again";
267
end if;
268
269
end function;
270
271
272
273
/***
274
Flattens a sequence of sequences into a single sequence.
275
Input: x is a sequence of sequences.
276
Output: a sequence got by concatenating the sequences in x.
277
***/
278
279
function flatten(x)
280
281
return &cat x;
282
283
end function;
284
285
286
287
/***
288
Input: elldata is a 4-tuple < elliptic curve, j-invariant, conductor,
289
torsion subgroup > where
290
elliptic_curve -- elliptic curve E over a number field K
291
j-invariant -- the j-invariant of E
292
conductor -- the conductor of E, as an integral ideal of K
293
torsion subgp -- an abstract abelian group
294
Output: an integer list representing the curve, of length n*(3+n)+3, where
295
n is the degree of the field extension.
296
***/
297
298
function elldata_to_sequence(elldata)
299
300
ans := [];
301
E, j, N, tor := Explode(elldata);
302
K := BaseField(E);
303
_,_,_,a,b := Explode(aInvariants(E));
304
Append(~ans, Eltseq(a));
305
Append(~ans, Eltseq(b));
306
Append(~ans, Eltseq(Numerator(j)));
307
Append(~ans, [Denominator(j)]);
308
for b in Basis(N) do
309
Append(~ans, Eltseq(K!b));
310
end for;
311
tor := Invariants(tor);
312
for i in [1..2-#tor] do
313
Append(~tor,[1]);
314
end for;
315
Append(~ans, tor);
316
317
return flatten(ans);
318
319
end function;
320
321
322
323
/***
324
Input: K is a number field. seq is an integer sequence.
325
Output: elldata 4-tuple.
326
***/
327
328
function sequence_to_elldata(K, seq)
329
330
assert Type(K) eq FldNum;
331
assert Type(seq) eq SeqEnum;
332
333
n := Degree(K);
334
k := 0;
335
a := K![seq[k+i] : i in [1..n]]; k +:= n;
336
b := K![seq[k+i] : i in [1..n]]; k +:= n;
337
E := EllipticCurve([a,b]);
338
jn := K![seq[k+i] : i in [1..n]]; k +:= n;
339
jd := K!seq[k+1]; k +:= 1;
340
j := jn/jd;
341
B := [];
342
for ii in [1..n] do
343
b := K![seq[k+i] : i in [1..n]]; k +:= n;
344
Append(~B, b);
345
end for;
346
N := ideal<MaximalOrder(K) | B>;
347
tor := AbelianGroup([seq[k+i] : i in [1..2]]); k +:= 2;
348
assert k eq #seq;
349
350
return <E, j, N, tor>;
351
352
end function;
353
354
355
356
357
procedure save_to_file(seq, filename)
358
359
assert Type(filename) eq MonStgElt;
360
assert Type(seq) eq SeqEnum;
361
362
file := Open("data/"*filename, "w");
363
for x in seq do
364
fprintf file, "%o ", x;
365
end for;
366
367
end procedure;
368
369
370
371
372
procedure append_to_file(seq, filename)
373
374
assert Type(filename) eq MonStgElt;
375
assert Type(seq) eq SeqEnum;
376
377
file := Open("data/"*filename, "a");
378
for x in seq do
379
fprintf file, "%o ", x;
380
end for;
381
382
end procedure;
383
384
385
386
387
function load_from_file(filename)
388
389
str := Read("data/"*filename);
390
391
return StringToIntegerSequence(str);
392
393
end function;
394
395
396
397
/***
398
Changes a sequence (e.g. [a, b, c]) to a string (“a_b_c”).
399
***/
400
401
function makename(seq)
402
403
s2 := Sprintf("%o", seq);
404
n := #s2;
405
S := "";
406
for i in [1..n] do
407
if s2[i] eq "[" then
408
elif s2[i] eq "," then
409
S cat:= "_";
410
elif s2[i] eq " " then
411
elif s2[i] eq "]" then
412
else
413
S cat:= s2[i];
414
end if;
415
end for;
416
417
return S;
418
419
end function;
420
421
422
423
/***
424
taking a file name gives back the integer sequence
425
***/
426
427
function parsename(str)
428
n := #str;
429
k := "";
430
for i in [1..n] do
431
if str[i] eq "_" then
432
k cat:= " ";
433
else
434
k cat:= str[i];
435
end if;
436
end for;
437
return StringToIntegerSequence(k);
438
end function;
439
440
441
442
/***
443
Checks if a file exists.
444
***/
445
446
function Exists(path)
447
assert Type(path) eq MonStgElt;
448
x := Gets(POpen("ls data/"*path,"r"));
449
return not IsEof(x);
450
end function;
451
452
453
454
/***
455
Input: a polynomial over Q defining a number field, and an integer i
456
giving the bound we want.
457
Extracts from a file a list of elliptic curve integral sequence.
458
***/
459
460
function Extract(K, filename)
461
462
R<x> := PolynomialRing(RationalField());
463
Z := RingOfIntegers();
464
P := parsename(filename);
465
466
assert Degree(K) eq P[1];
467
468
n := Degree(K);
469
470
assert R!DefiningPolynomial(K) eq R!Polynomial([P[i+1] : i in [1..n+1]]);
471
472
List:=[];
473
474
N := (n+3)*n+3;
475
datastring := load_from_file(filename);
476
l := #datastring/N;
477
assert Denominator(l) eq 1;
478
l := Z!l;
479
for j in [0..l-1] do
480
L := [];
481
for k in [1..N] do
482
Append(~L, datastring[j*N+k]);
483
end for;
484
s := sequence_to_elldata(K, L);
485
Append(~List, s);
486
end for;
487
return List;
488
489
end function;
490
491
492
493
/***
494
Generates data; checks if it already exists before making it anew.
495
***/
496
497
procedure calculate(K, B);
498
499
assert Type(K) eq FldNum;
500
assert Type(B) eq RngIntElt;
501
502
K<a> := K;
503
m := DefiningPolynomial(K);
504
n := Degree(K);
505
M := [];
506
L2 := [];
507
J2 := [];
508
509
for i in [1..B] do
510
name:= [];
511
Append(~name, [Degree(K)]);
512
Append(~name, Coefficients(m));
513
Append(~name, [i]);
514
Append(~M, makename(flatten(name)));
515
end for;
516
517
for i in [1..B] do
518
if Exists(M[i]) eq true then
519
print "found data for size", i;
520
Starter := Extract(K, M[i]);
521
for k in [1..#Starter] do
522
Append(~L2, Starter[k][1]);
523
Append(~J2, Starter[k][2]);
524
end for;
525
L2;
526
J2;
527
else
528
if i eq 1 then
529
L := cleanlist(K, 1);
530
seq_i := [];
531
Starter := [];
532
for j in [1..#L[1]] do
533
G:=TorsionSubgroup(L[1][j]);
534
C:=Conductor(L[1][j]);
535
Info:= <L[1][j], L[2][j], C, G>;
536
Append(~Starter, Info);
537
seq:=elldata_to_sequence(Info);
538
Append(~seq_i, seq);
539
end for;
540
for k in [1..#Starter] do
541
Append(~L2, Starter[k][1]);
542
Append(~J2, Starter[k][2]);
543
end for;
544
finalseq := flatten(seq_i);
545
save_to_file(finalseq, M[i]);
546
print "done with 1";
547
else
548
P := extendcleanlist(K, L2, J2, i-1, i);
549
seq_i := [];
550
Starter := [];
551
for j in [1..#P[1]] do
552
G := TorsionSubgroup(P[1][j]);
553
C := Conductor(P[1][j]);
554
Info := <P[1][j], P[2][j], C, G>;
555
Append(~Starter, Info);
556
seq := elldata_to_sequence(Info);
557
Append(~seq_i, seq);
558
end for;
559
for k in [1..#Starter] do
560
Append(~L2, Starter[k][1]);
561
Append(~J2, Starter[k][2]);
562
end for;
563
finalseq := flatten(seq_i);
564
save_to_file(finalseq, M[i]);
565
print "done with", i;
566
end if;
567
end if;
568
end for;
569
end procedure;
570
571
572
573
/***
574
a procedure to sort curves in standard form by norm of the conductor
575
576
This is bubble sort, which is not very fast when the length of the sequence
577
is large. Maybe past in my quicksort, which is quicker.
578
***/
579
580
function sortcurves(seq)
581
assert Type(seq) eq SeqEnum;
582
for i in [1..#seq] do
583
for j in [1..#seq-1] do
584
if Norm(seq[j][3]) gt Norm(seq[j+1][3]) then
585
Insert(~seq, j, seq[j+1]);
586
Remove(~seq, j+2);
587
end if;
588
end for;
589
end for;
590
return seq;
591
end function;
592
593
594
/***
595
quicksort procedures
596
***/
597
598
function lessthan(a, b)
599
if Norm(a[3]) lt Norm(b[3]) then
600
return true;
601
else
602
return false;
603
end if;
604
end function;
605
606
procedure QuickSort_recurse(~items, low, high, lessthan)
607
// sorts Seqenum using the quick sort algorithm
608
if low ge high then
609
return;
610
end if;
611
lo := low;
612
hi := high + 1;
613
elem := items[low];
614
while true do
615
while lo lt high and
616
lessthan(items[lo+1],elem) do
617
lo +:= 1;
618
end while;
619
lo +:= 1;
620
while lessthan(elem,items[hi-1]) do
621
hi -:= 1;
622
end while;
623
hi -:= 1;
624
if lo lt hi then
625
t := items[hi];
626
items[hi] := items[lo];
627
items[lo] := t;
628
else
629
break;
630
end if;
631
end while;
632
t := items[hi];
633
items[hi] := items[low];
634
items[low] := t;
635
QuickSort_recurse(~items,low,hi-1,lessthan);
636
QuickSort_recurse(~items,hi + 1,high,lessthan);
637
end procedure;
638
639
640
procedure QuickSort(~items, lessthan)
641
// items -- a sequence
642
// lessthan -- a function that takes two arguments, and returns true
643
// if and only if the first is "less than" the second.
644
QuickSort_recurse(~items,1,#items, lessthan);
645
end procedure;
646
647
648
/***
649
the ultimate table-making procedure
650
***/
651
652
function maketable(K, B);
653
654
assert Type(K) eq FldNum;
655
assert Type(B) eq RngIntElt;
656
657
R<x>:=PolynomialRing(Rationals());
658
K<a>:=K;
659
m:=DefiningPolynomial(K);
660
M:=[];
661
for i in [1..B] do
662
name := [];
663
Append(~name, [Degree(K)]);
664
Append(~name, Coefficients(m));
665
Append(~name, [i]);
666
Append(~M, makename(flatten(name)));
667
end for;
668
for i in [1..B] do
669
if Exists(M[i]) eq false then
670
return "data has not yet been computed for", i, ",please hold
671
while it is being computed" ;
672
calculate(K, B);
673
end if;
674
end for;
675
BigList:=[];
676
for i in [1..B] do
677
Append(~BigList, Extract(K, M[i]));
678
end for;
679
L:=flatten(BigList);
680
QuickSort(~L, lessthan);
681
LL:=[];
682
for i in [1..#L] do
683
s0 := Norm(L[i][3]);
684
_,_,_,a,b := Explode(aInvariants(L[i][1]));
685
s1:=a;
686
s2 := b;
687
s3 := L[i][2];
688
t, s := IsPrincipal(L[i][3]);
689
s4 := K!s;
690
tor := Invariants(L[i][4]);
691
for i in [1..2-#tor] do
692
Append(~tor, [1]);
693
end for;
694
s5 := tor[1];
695
s6 := tor[2];
696
Append(~LL, <s0, s1, s2, [s5, s6], s5*s6, s3, s4>);
697
end for;
698
return LL;
699
end function;
700
701
702
703
704
function goodtable(K, B)
705
assert Type(K) eq FldNum;
706
assert Type(B) eq RngIntElt;
707
708
K<a>:=K;
709
n := Degree(K);
710
m := DefiningPolynomial(K);
711
M:=maketable(K, B);
712
list:="";
713
length := 12*n-9;
714
for i in [1..#M] do
715
f:="";
716
s1:=Sprintf("%o", M[i][2]);
717
s11:="";
718
for j in [1..#s1] do
719
if s1[j] eq " " then
720
else
721
s11 cat:= s1[j];
722
end if;
723
end for;
724
s2:=Sprintf("%o", M[i][3]);
725
s22:="";
726
for j in [1..#s2] do
727
if s2[j] eq " " then
728
else
729
s22 cat:= s2[j];
730
end if;
731
end for;
732
length1:=#s11;
733
length2:=#s22;
734
L:=length1+length2+3;
735
for j in [1..#M[1]] do
736
s:=Sprintf("%o", M[i][j]);
737
if j eq 2 then
738
f cat:= "["*s*",";
739
elif j eq 3 then
740
gap:=length-L;
741
f cat:=s*"]"*"_"^(gap+1);
742
elif j eq 6 then
743
if #s le 7 then
744
f cat:= s*"\t\t\t";
745
elif #s gt 7 and #s le 14 then
746
f cat:= s*"\t\t";
747
else
748
f cat:= s*"\t";
749
end if;
750
elif j eq 7 then
751
f cat:= "("*s*")"*"\n";
752
else
753
f cat:= s*"\t";
754
end if;
755
end for;
756
f2:="";
757
for j in [1..#f] do
758
if f[j] eq " " then
759
elif f[j] eq "_" then
760
f2 cat:= " ";
761
else
762
f2 cat:= f[j];
763
end if;
764
end for;
765
list cat:= f2;
766
end for;
767
return list;
768
end function;
769
770
771
772
procedure goodtablefile(m, B)
773
R<x>:=PolynomialRing(Rationals());
774
m:=R!m;
775
K<a>:=NumberField(m);
776
n:=Sprintf("%o", Degree(K));
777
c:=makename(Coefficients(m));
778
b:=Sprintf("%o", B);
779
filename := "table_degree_"*n*"_poly_"*c*"_size_"*b*".txt";
780
s := goodtable(K,B);
781
F := Open(filename,"w");
782
fprintf F, "%o", s;
783
end procedure;
784