## Attachment 'intptsanf.m'

Download

```   1 ///////////////////////////////////////////////////////////////////////////////
2 // Computing all integral points on elliptic curves over number fields
3 // Based on Smart & Stephens' paper (Math. Proc. Camb. Phil. Soc. 122 (1997),
4 // pp. 9-16) with some modification
5 //
6 // Written by Thotsaphon Thongjunthug, BSc(Hons)
7 // Last update: 14 May 2009
8 //
9 // This package requires some functions on Silverman's bound implemented in
10 // John Cremona's "nfhtbound.m", available freely on his webpage
11 // http://www.warwick.ac.uk/~masgaj/ftp/progs/magma/index.html
12 ///////////////////////////////////////////////////////////////////////////////
13
14
15 ///////////////////////////////////////////////////////////////////////////////
16 // Commands:
17 // IntegralPoints() : Find all integral points on elliptic curves over
18 //                    number fields
19 ///////////////////////////////////////////////////////////////////////////////
20
21
22 // Declare printing verbose
23 // 0 = return result only, no detail will be printed
24 // 1 = minimal amount of details shown
25 // 2 = all details (e.g. values of all constants) shown (for debugging only)
26 declare verbose Detail, 2;
27
28
29 ///////////////////////////////////////////////////////////////////////////////
30 // Constants for an upper bound on linear form in logarithm
31 // The indices of c are as in Smart & Stephens' paper
32 ///////////////////////////////////////////////////////////////////////////////
33
34 // Compute constant c3
35 // Input: L = a sequence of points on elliptic curves over number field
36 function c3(L)
37   M := Eigenvalues(HeightPairingMatrix(L));
38   M := SetToSequence(M);
39   M := [m : m in M]; // ignore multiplicity
40   M := Minimum(M); // least eigenvalue
41   return M;
42 end function;
43
44 // Compute constant c6
45 // Input: E = elliptic curve over real/complex numbers
46 function c6(E)
47   P<x> := PolynomialRing(ComplexField());
48   C4, C6 := Explode(cInvariants(E));
49   f := x^3 - (P! C4/48)*x - (P! C6/864);
50   R := Roots(f);
51   R := [r : r in R]; // ignore multiplicity
52   R := [Abs(r) : r in R];
53   R := 2 * Maximum(R);
54   return R;
55 end function;
56
57 // Compute constant c8, and the periods of the lattice L where E(C) = C/L
58 // This requires complex periods of elliptic curves over C
59 // Require: complexell.m
60 // Input: E = elliptic curve over real/complex numbers
61 function c8(E)
62   org_w1, org_w2 := Explode(ComplexPeriods(E));
63   // Apply transformation by SL(2,Z) so that tau = w2/w1 is in the
64   // fundamental domain
65   _, A := ReducedTau(org_w2 / org_w1); // require "complexell.m"
66   w2 := A[1,1]*org_w2 + A[1,2]*org_w1;
67   w1 := A[2,1]*org_w2 + A[2,2]*org_w1;
68   L := [Abs(w) : w in [w1, w2, w1+w2]];
69   L := Maximum(L);
70   return L, [w1, w2];
71 end function;
72
73 // Compute absolute logarithmic height of an element in n-projective space
74 // over a number field K
75 // Input: X = a sequence of elements in a number field K
76 function AbsLogHeight(X)
77   // X cannot be zero vector
78   if #X eq #[x : x in X | x eq 0] then
79     error "AbsLogHeightExt: X cannot be zero vector";
80   end if;
81
82   // Find all prime ideal that divides some denominators of x_i
83   I := {};
84   K := Parent(X); // assume each x_i is in the same field
85   O := RingOfIntegers(K);
86   for x in X do
87     den := O ! Denominator(x);
88     L := Decomposition(den);
89     L := {l : l in L}; // ignore multiplicity
90     I := I join L;
91   end for;
92
93   h := 0;
94   // Non-archimedean contributions
95   for p in I do
96     M := [Rationals()| ];
97     for x in X do
98       if x eq 0 then
99         Append(~M, 0);
100       else
101         Append(~M, Norm(p)^(-Valuation(x, p)));
102       end if;
103     end for;
104     M := Maximum(M);
105     M := Log(M);
106     h +:= M;
107   end for;
108
109   // Archimedean contributions
110   s1, s2 := Signature(K);
111   // M = [log max(|x_1|_v,..., |x_n|_v) : v in M_K]
112   M := [];
113   for x in X do
114     C := Conjugates(x);
115     newC := [];
116     // Real embedding contributions
117     for i := 1 to s1 do
118       Append(~newC, Abs(C[i]));
119     end for;
120     // Complex embedding contributions
121     for i := 1 to s2 do
122       Append(~newC, Abs(C[s1+(2*i-1)])^2);
123     end for;
124
125     if #M eq 0 then
126       M := newC;
127     else
128       for i := 1 to (s1+s2) do
129         if newC[i] gt M[i] then
130           M[i] := newC[i];
131         end if;
132       end for;
133     end if;
134   end for;
135   M := [Log(m) : m in M];
136
137   // Overall absolute logarithmic height
138   h +:= (&+M);
139   h /:= Degree(K);
140   return h;
141 end function;
142
143
144 ///////////////////////////////////////////////////////////////////////////////
145 // Constants for David's lower bound on linear form in logarithm
146 // Notation used as in Appendix A of Smart's book
147 // "The Algorithmic Resolution of Diophantine Equations", with c's being
148 // replaced by d's
149 ///////////////////////////////////////////////////////////////////////////////
150
151 // Compute the "height" of elliptic curve
152 // Input: E = elliptic curve over number field
153 function h_E(E)
154   j := jInvariant(E);
155   C4, C6 := Explode(cInvariants(E));
156   g2 := C4/12;
157   g3 := C6/216;
158   m := Maximum([1, AbsLogHeight([1, g2, g3]), AbsLogHeight([1, j])]);
159   return m;
160 end function;
161
162 // Compute the list of modified height of a point P in E(K)
163 // according to embeddings
164 // Input: E = elliptic curve over a number field K
165 //        P = a point in E(K)
166 //        ElogEmbedP = elliptic logarithm of P when embedded in some
167 //                     embedding
168 //        D7 = constant d7 (depending on embedding)
169 function h_m(E, P, ElogEmbedP, D7)
170   K := BaseRing(E);
171   L := [Height(P), h_E(E), D7/Degree(K)*Abs(ElogEmbedP)^2];
172   L := Maximum(L);
173   return L;
174 end function;
175
176 // Compute two extra h_m's based on the two periods
177 // Similar to h_m, but now ElogEmbedP becomes a period of the fundamental
178 // parallelogram of E(C) for some embedding
179 // Input: E = elliptic curve over a number field
180 //        Periods = a sequence of two periods of the fundamental
181 //                  parallelogram of E(C) for some embbeding
182 //        D7 = constant d7 (depending on embedding)
183 function Extra_h_m(E, Periods, D7)
184   D := Degree(BaseRing(E));
185   h := h_E(E);
186   tmp1 := Maximum([0, h, D7/D*Abs(Periods)^2]);
187   tmp2 := Maximum([0, h, D7/D*Abs(Periods)^2]);
188   return [tmp1, tmp2];
189 end function;
190
191 // Compute constant d8
192 // Input: E = elliptic curve over a number field K
193 //        L = a sequence of points in E(K) (e.g. Mordell-Weil basis)
194 //        Elog = a sequence of (pre-computed) elliptic logarithm of
195 //               each point in L at a fixed embedding
196 //        Periods = a sequence of two periods of the fundamental
197 //                  parallelogram of E(C) for some embbeding
198 //        D7 = constant d7 (depending on embedding)
199 function d8(E, L, Elog, Periods, D7)
200   C := [Exp(1) * h_E(E)];
201   D := Degree(BaseRing(E));
202   for i := 1 to #L do
203     Append(~C, h_m(E, L[i], Elog[i], D7) / D);
204   end for;
205   C := C cat [t / D : t in Extra_h_m(E, Periods, D7)];
206   C := Maximum(C);
207   return C;
208 end function;
209
210 // Compute constant d9
211 // Input: E = elliptic curve over a number field
212 //        L = a sequence of points in E(K) (e.g. Mordell-Weil basis)
213 //        Elog = a sequence of (pre-computed) elliptic logarithm of
214 //               each point in L at a fixed embedding
215 //        Periods = a sequence of two periods of the fundamental
216 //                  parallelogram of E(C) for some embbeding
217 //        D7 = constant d7 (depending on embedding)
218 function d9(E, L, Elog, Periods, D7)
219   D := Degree(BaseRing(E));
220   C := [];
221   for i := 1 to #L do
222     tmp := Exp(1) * Sqrt(D * h_m(E, L[i], Elog[i], D7) / D7);
223     tmp /:= Abs(Elog[i]);
224     C[i] := tmp;
225   end for;
226
227   // Take minimum among extra_h_m
228   Ehm := Extra_h_m(E, Periods, D7);
229   tmp1 := Exp(1) * Sqrt(D*Ehm/D7) / Abs(Periods);
230   tmp2 := Exp(1) * Sqrt(D*Ehm/D7) / Abs(Periods);
231   C := C cat [tmp1, tmp2];
232   C := Minimum(C);
233   return C;
234 end function;
235
236 // Compute constant d10
237 // Input: E = elliptic curve over a number field
238 //        L = a sequence of points in E(K) (e.g. Mordell-Weil basis)
239 //        Elog = a sequence of (pre-computed) elliptic logarithm of
240 //               each point in L at a fixed embedding
241 //        Periods = a sequence of two periods of the fundamental
242 //                  parallelogram of E(C) for some embbeding
243 //        D7 = constant d7 (depending on embedding)
244 function d10(E, L, Elog, Periods, D7)
245   D := Degree(BaseRing(E));
246   n := #L+2;
247   scalar := 2 * 10^(8 + 7*n) * (2/Exp(1))^(2*n^2);
248   scalar *:= (n+1)^(4*n^2 + 10*n) * D^(2*n + 2);
249   scalar *:= (Log(d9(E, L, Elog, Periods, D7)))^(-2*n-1);
250   for i := 1 to #L do
251     scalar *:= h_m(E, L[i], Elog[i], D7);
252   end for;
253   scalar *:= &*(Extra_h_m(E, Periods, D7));
254   return scalar;
255 end function;
256
257 // Compute the right-hand side of the Principal Inequality
258 // Input: D   = Degree(K), K = number field
259 //        r   = rank(E(K))
260 //        C9  = constant c9
261 //        C10 = constant c10
262 //        D9  = constant d9
263 //        D10 = constant d10
264 //        h   = h_E(E), E = elliptic curve over K
265 //        Q   = initial bound for the coefficients of P_i's, where
266 //              P_i's are in the Mordell-Weil basis
267 function RHS(D, r, C9, C10, D9, D10, h, Q, expTors)
268   bound := (Log(Log(Q*r*expTors)) + h + Log(D*D9))^(r+3);
269   bound *:= D10*(Log(Q*r*expTors) + Log(D*D9));
270   bound +:= Log(C9*expTors);
271   bound /:= C10;
272   return bound;
273 end function;
274
275 // Approximate initial bound on Q = max_{1\le i\le r}{|q_i|}
276 // Input: D   = Degree(K), where K = number field
277 //        r   = rank of E(K)
278 //        Q0  = constant Q0 (in S&S paper, this is called K0)
279 //        C9  = constant c9
280 //        C10 = constant c10
281 //        D8  = constant d8 (from function d8())
282 //        D9  = constant d9 (from function d9())
283 //        D10 = constant d10 (from function d10())
284 //        h   = h_E(E)
285 //        expTors = exponent of the torsion subgroup of E(K)
286 // Revised: 5 May 2009
287 function InitialQ(D, r, Q0, C9, C10, D8, D9, D10, h, expTors)
288   minQ := Maximum(Q0, Exp(D8));
289
290   // Try to approximate Q such that Q^2 = RHS(Q) (i.e. Q makes both sides
291   // of the Principal Inequality equal)
292   // Firstly, set a guess for Q, say minQ + 1 (so that Q > minQ)
293   // For simplicity, let's round Q up to the nearest power of 10
294   Q := minQ + 1;
295   x := Ceiling(Log(10, Q)); // x = log_10(Q)
296
297   // Check if Q satisfies the Principal Inequality, i.e. if Q^2 < RHS(Q)
298   // If so, repeat with the larger Q until we find the first Q that
299   // violates the P.I.
300   // N.B. This loop will eventually terminate
301   exp_OK := 0;   // the exponent that satisfies P.I.
302   exp_fail := 0; // the exponent that fails P.I.
303   while 10^(2*x) lt RHS(D, r, C9, C10, D9, D10, h, 10^x, expTors) do
304     exp_OK := x; // Principal Inequality satisfied
305     x *:= 2;     // double x, and retry
306   end while;
307   exp_fail := x; // x that fails the Principal Inequality
308
309   // So now x = log_10(Q) must lie between exp_OK and exp_fail
310   // Refine x further using "binary search"
311   repeat
312     x := Floor((exp_OK + exp_fail)/2);
313     if 10^(2*x) ge RHS(D, r, C9, C10, D9, D10, h, 10^x, expTors) then
314       exp_fail := x;
315     else
316       exp_OK := x;
317     end if;
318   until (exp_fail - exp_OK) le 1;
319   return exp_fail; // over-estimate
320 end function;
321
322 // Reduce the bound Q by LLL reduction until no further improvement
323 // is possible. This function initially requires high precision to
324 // proceed, although this should be done automatically by now
325 // Require: complexell.m
326 // Input: EmbedE  = model of E embedded in some (complex) embedding
327 //        EmbedL  = a sequence of points on EmbedE (e.g. points in the
328 //                  Mordell-Weil basis when embedded into EmbedE)
329 //        Elog    = a sequence of (pre-computed) elliptic logarithm of
330 //                  each point in EmbedL
331 //        C9      = constant c9
332 //        C10     = constant c10
333 //        Periods = sequence of the double periods of the fundamental
334 //                  parallelogram associated to EmbedE
335 //        expTors = exponent of the torsion subgroup of E(K), K = number field
336 //        initQ   = initial guess for Q to be reduced by LLL
337 function ReducedQ(EmbedE, EmbedL, Elog, C9, C10, Periods, expTors, initQ)
338   w1, w2 := Explode(Periods);
339   r := #EmbedL;
340   newQ := initQ;
341
342   // Repeat LLL reduction until no further reduction is possible
343   repeat
344     Q := newQ;
345     S := r*(Q^2)*(expTors^2);
346     T := 3*r*Q*expTors/Sqrt(2);
347
348     // Create the basis matrix
349     C := 1;
350     repeat
351       C *:= Q^Ceiling((r+2)/2);
352       L := ZeroMatrix(Integers(), r+2, r+2);
353
354       // Elliptic logarithm should have precision "suitable to" C
355       // e.g. If C = 10^100, then Re(Elog[i]) should be computed
356       // correctly to at least 100 decimal places
357       pow10_C := Ceiling(Log(10, C));
358       if pow10_C gt Precision(Elog) then
359 	vprint Detail, 2:
360 	  "Need higher precision, recompute elliptic logarithm ...";
361         // Re-compute elliptic logarithm to the right precision
362         Elog := [ComplexEllLog(EmbedE, P : Prec := pow10_C+10) : P in EmbedL];
363         vprint Detail, 2: "Elliptic logarithm recomputed";
364       end if;
365       // Assign all non-zero entries
366       for i := 1 to r do
367         L[i,i] := 1;
368         L[r+1, i] := Truncate(C*Re(Elog[i]));
369         L[r+2, i] := Truncate(C*Im(Elog[i]));
370       end for;
371       L[r+1, r+1] := Truncate(C*Re(w1));
372       L[r+1, r+2] := Truncate(C*Re(w2));
373       L[r+2, r+1] := Truncate(C*Im(w1));
374       L[r+2, r+2] := Truncate(C*Im(w2));
375       L := Transpose(L); // In Magma, basis is spanned by row vector!
376
377       // LLL reduction and constants
378       L := LLL(L);
379       b1 := L; // 1st row of reduced basis
380       // Norm(b1) = square of Euclidean norm
381       lhs := 2^(-r-1)*Norm(b1) - S;
382     until (lhs ge 0) and (Sqrt(lhs) gt T);
383
384     newQ := ( Log(C*C9*expTors) - Log(Sqrt(lhs) - T) ) / C10;
385     newQ := Floor(Sqrt(newQ));
386     pow10 := Floor(Log(10, C));
387     vprintf Detail, 2: "Choose C = %.4o x 10^%o. ", 1.*C/10^pow10, pow10;
388     vprintf Detail, 2: "After reduction, Q <= %o\n", newQ;
389   until ((Q - newQ) le 1) or (newQ le 1);
390   return newQ;
391 end function;
392
393
394 ///////////////////////////////////////////////////////////////////////////////
395 // Main intrinsic functions
396 ///////////////////////////////////////////////////////////////////////////////
397
398 // Search for all integral points on elliptic curves over number fields
399 // within a given bound
400 // Input:    E = elliptic curve over a number field K
401 //           L = a sequence of points in the Mordell-Weil basis for E(K)
402 //           Q = a maximum for the absolute bound on all coefficients
403 //               in the linear combination of points in L
404 // Output:  S1 = sequence of all integral points on E(K) modulo [-1]
405 //          S2 = sequence of tuples representing the points as a
406 //               linear combination of points in L
407 // Option: tol = tolerance used for checking integrality of points.
408 //               (Default = 0 - only exact arithmetic will be performed)
409 intrinsic IntegralPoints(E::CrvEll, L::[PtEll], Q::RngIntElt : tol := 0) ->
410   SeqEnum, SeqEnum
411 {Given an elliptic curve E over a number field, its Mordell-Weil basis L, and a positive integer Q, return the sequence of all integral points modulo [-1] of the form P = q_1*L + ... + q_r*L[r] + T with some torsion point T and |q_i| <= Q, followed by a sequence of tuple sequences representing the points as a linear combination of points. An optional tolerance may be given to speed up the computation when checking integrality of points.}
412
413   // Check input validity
414   require IsNumberField(BaseRing(E)):
415     "Elliptic curve must be defined over a number field";
416   require tol ge 0: "Tolerance must be non-negative";
417
418   // Find the generators of the torsion subgroup of E(K)
419   Tors, map := TorsionSubgroup(E);
420   expTors := Exponent(Tors);
421   G := Generators(Tors);
422   Tors := [map(g) : g in G];   // each generator of E(K)_tors
423   OrdG := [Order(g) : g in G]; // order of each generator
424
425   if (#L eq 0) and (#Tors eq 0) then
426     return [], [];
427   end if;
428
429   // Create all possible r+#Tors-tuples
430   r := #L; // r = rank of E(K)
431   C := [0 : i in [1..(r + #Tors)]];
432   ListC := [];
433   for i := 0 to Q do
434     C := i;
435     Append(~ListC, C);
436   end for;
437
438   for i := 2 to r do
439     tmp := [];
440     for j := 1 to Q do
441       for l in ListC do
442         tup := l;
443         tup[i] := j;
444         Append(~tmp, tup);
445
446         // Avoid having its negative in the list
447         // Only use when all previous entries in tuple are zero
448         for k := 1 to i-1 do
449           if tup[k] ne 0 then
450             tup[i] := -j;
451             Append(~tmp, tup);
452             break;
453           end if;
454         end for;
455
456       end for;
457     end for;
458     ListC := ListC cat tmp;
459   end for;
460
461   // Add torsion point, if any
462   if #Tors ne 0 then
463     for i := 1 to #Tors do
464       tmp := [];
465       for j := 1 to (OrdG[i]-1) do
466         for l in ListC do
467           tup := l;
468           tup[r+i] := j;
469           Append(~tmp, tup);
470         end for;
471       end for;
472       ListC := ListC cat tmp;
473     end for;
474   end if;
475   Remove(~ListC, 1); // remove point at infinity
476
477   L := L cat Tors;
478   vprint Detail, 2: "Generators = ", L;
479   PtsList := [];
480   CoeffList := [];
481
482   // Skip the complex arithmetic and only perform exact arithmetic if tol = 0
483   if tol eq 0 then
484     vprint Detail : "Exact arithmetic";
485     for l in ListC do
486       P := &+[l[i]*L[i] : i in [1..#L]];
487       if IsIntegral(P) and IsIntegral(P) then
488         vprintf Detail: "%o ---> %o\n", l, P;
489         Append(~PtsList, P);
490         TupList := [ <L[i], l[i]> : i in [1..#L] | l[i] ne 0 ];
491         Append(~CoeffList, TupList);
492       end if;
493     end for;
494     vprint Detail: "*"^45;
495     return PtsList, CoeffList;
496   end if;
497
498   // Suggested by John Cremona
499   // Point search. This is done via arithmetic on complex points on each
500   // embedding of E. Exact arithmetic will be carried if the resulting
501   // complex points are "close" to being integral, subject to some tolerance
502
503   // Embed each generator of the torsion subgroup
504   X := [Conjugates(P) : P in (L cat Tors)];
505   Y := [Conjugates(P) : P in (L cat Tors)];
506   // Create all embeddings of E
507   K := BaseRing(E);
508   s1, s2 := Signature(K);
509   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
510   a1 := Conjugates(a1);
511   a2 := Conjugates(a2);
512   a3 := Conjugates(a3);
513   a4 := Conjugates(a4);
514   a6 := Conjugates(a6);
515   EmbedEList := [EllipticCurve([a1[j], a2[j], a3[j], a4[j], a6[j]]): j in
516     [1..(s1+2*s2)]];
517
518   // Set tolerance - This should be larger than 10^(-current precision) to
519   // avoid missing any integral points. Too large tolerance will slow the
520   // computation, too small tolerance may lead to missing some integral points.
521   // This is now made as an option.
522   vprint Detail: "Tolerance = ", tol * 1.;
523
524   // Create the matrix containing all embeddings of the integral basis of K
525   // as its columns
526   IntBasis := IntegralBasis(K);
527   B := Matrix([Conjugates(a) : a in IntBasis]);
528   // Note that B is always invertible, so we can take its inverse
529   B := B^(-1);
530
531   // Embeddings of each point in L
532   // Format: [[P1 ... P1], [P2 ... P2], ...]
533   EmbedLList := [];
534   for i := 1 to (r+#Tors) do
535     EmbedLi := [];
536     for j := 1 to (s1+2*s2) do
537       P_i := Points(EmbedEList[j], X[i][j]);
538       // Choose the right sign for the y-coordinate
539       if (Y[i][j] ne 0) and (Re(P_i/Y[i][j]) lt 0) then
540         P_i := -P_i;
541       end if;
542       Append(~EmbedLi, P_i);
543     end for;
544     Append(~EmbedLList, EmbedLi);
545   end for;
546
547   // Point search
548   for l in ListC do
549     // Compute P by complex arithmetic for every embedding
550     EmbedPList := [E![0,1,0] : E in EmbedEList];
551     for j := 1 to (s1+2*s2) do
552       EmbedPList[j] := &+[l[i]*EmbedLList[i][j] : i in [1..(r + #Tors)]];
553     end for;
554
555     // Check if the x-coordinate of P is "close to" being integral
556     // If so, compute P exactly and check if it is integral; skip P otherwise
557     XofP := Matrix([[P : P in EmbedPList]]);
558     // Write x(P) w.r.t. the integral basis of (K)
559     // Due to the floating arithmetic, some entries in LX may have very tiny
560     // imaginary part, which can be thought as zero
561     LX := XofP * B;
562     LX := [Abs( Re(LX[1,i]) - Round(Re(LX[1,i])) ): i in [1..#IntBasis]];
563     LX := &and[dx lt tol : dx in LX];
564     if not LX then
565       // x-coordinate of P is not integral, skip P
566       continue;
567     end if;
568
569     // Now check P by exact arithmetic
570     // Add P and the list of tuples representing P into the list
571     // if P is integral
572     P := &+[l[i]*L[i] : i in [1..#L]];
573     if IsIntegral(P) and IsIntegral(P) then
574       vprintf Detail: "%o ---> %o\n", l, P;
575       Append(~PtsList, P);
576       TupList := [ <L[i], l[i]> : i in [1..#L] | l[i] ne 0 ];
577       Append(~CoeffList, TupList);
578     end if;
579   end for;
580   vprint Detail: "*"^45;
581   return PtsList, CoeffList;
582 end intrinsic;
583
584 // Compute all integral points on elliptic curve over a number field
585 // This function simply computes a suitable bound Q, and return
586 // IntegralPoints(E, L, Q : tol := ...).
587 // Input :   E = elliptic curve over a number field K
588 //           L = a sequence of points in the Mordell-Weil basis for E(K)
589 // Output:  S1 = sequence of all integral points on E(K) modulo [-1]
590 //          S2 = sequence of tuples representing the points as a
591 //               linear combination of points in L
592 // Option: tol = tolerance used for checking integrality of points.
593 //               (Default = 0 - only exact arithmetic will be performed)
594 intrinsic IntegralPoints(E::CrvEll, L::[PtEll] : tol := 0) -> SeqEnum, SeqEnum
595 {Given an elliptic curve over a number field and its Mordell-Weil basis, return the sequence of all integral points modulo [-1], followed by a sequence of tuple sequences representing the points as a linear combination of points. An optional tolerance may be given to speed up the computation when checking integrality of points. (This function simply computes a suitable Q and call IntegralPoints(E, L, Q: tol := ...)}
596
597   // Check input validity
598   require IsNumberField(BaseRing(E)):
599     "Elliptic curve must be defined over a number field";
600   require tol ge 0:
601     "Tolerance must be non-negative";
602
603   if #L eq 0 then
604     return IntegralPoints(E, [], 0 : tol := tol);
605   end if;
606
607   a1, a2, a3, a4, a6 := Explode(aInvariants(E));
608
609   // Embed E into various (real/complex) embeddings
610   a1 := Conjugates(a1);
611   a2 := Conjugates(a2);
612   a3 := Conjugates(a3);
613   a4 := Conjugates(a4);
614   a6 := Conjugates(a6);
615   K := BaseRing(E);
616   s1, s2 := Signature(K);
617   b2 := Conjugates(bInvariants(E));
618   pi := Pi(RealField());
619
620   // Embed generators in the Mordell-Weil basis
621   X := [Conjugates(P) : P in L];
622   Y := [Conjugates(P) : P in L];
623
624   // Find the generators of the torsion subgroup of E(K)
625   Tors, map := TorsionSubgroup(E);
626   expTors := Exponent(Tors);
627   G := Generators(Tors);
628   Tors := []; OrdG := [];
629   for g in G do
630     Append(~Tors, map(g));   // torsion point
631     Append(~OrdG, Order(g)); // order of torsion point
632   end for;
633
634   // Global constants (independent to the embedding of E)
635   C2 := -SilvermanLowerHeightBound(E); // require "nfhtbound.m"
636   C3 := c3(L);
637   h := h_E(E);
638   vprint Detail, 2: "Global constants";
639   vprintf Detail, 2: "c2 = %.4o\n", C2;
640   vprintf Detail, 2: "c3 = %.4o\n", C3;
641   vprintf Detail, 2: "h_E = %.4o\n", h;
642   vprint Detail, 2: "-"^45;
643
644   Q := [];
645   // Find the most reduced bound on each embedding of E
646   // But first let's adjust the index
647   for i := 1 to (s1+s2) do
648     if i le s1 then
649       j := i;
650       nv := 1;
651       vprintf Detail, 2: "Real embedding #%o\n", j;
652     else
653       j := s1 + (2*(i-s1)-1);
654       nv := 2;
655       vprintf Detail, 2: "Complex embedding #%o\n", i-s1;
656     end if;
657
658     // Create complex embedding of E
659     ee := EllipticCurve([a1[j], a2[j], a3[j], a4[j], a6[j]]);
660
661     // Embedding of all points in Mordell-Weil basis
662     // Find complex elliptic logarithm of each embedded point
663     // Let's set precision initially to, say, 1000
664     EmbedL := [[X[i][j], Y[i][j]] : i in [1..#L]];
665     Elog := [ComplexEllLog(ee, P: Prec := 1000) : P in EmbedL];
666     // Local constants (depending on embedding)
667     // C9, C10 are used for the upper bound on the linear form in logarithm
668     C4 := C3 * Degree(K) / (nv*(s1+s2));
669     vprintf Detail, 2: "c4 =  %.4o\n", C4;
670     C5 := C2 * Degree(K) / (nv*(s1+s2));
671     vprintf Detail, 2: "c5 = %.4o\n", C5;
672     C6 := c6(ee);
673     vprintf Detail, 2: "c6 = %.4o\n", C6;
674     delta := 1 + (nv-1)*pi;
675     C8, Periods := c8(ee);
676     vprintf Detail, 2: "c8 = %.4o\n", C8;
677     C7 := 8*(delta^2) + (C8^2)*Abs(b2[j])/12;
678     vprintf Detail, 2: "c7 = %.4o\n", C7;
679     C9 := C7*Exp(C5/2);
680     vprintf Detail, 2: "c9 = %.4o\n", C9;
681     C10 := C4/2;
682     vprintf Detail, 2: "c10 = %.4o\n", C10;
683     Q0 := Sqrt( ( Log(C6+Abs(b2[j])/12) + C5) / C4 );
684     vprintf Detail, 2: "Q0 = %.4o\n", Q0;
685
686     // Constants for David's lower bound on the linear form in logarithm
687     w1, w2 := Explode(Periods); // N.B. periods are already in "standard" form
688     D7 := 3*pi / ((Abs(w2)^2) * Im(w2/w1));
689     D8 := d8(E, L, Elog, Periods, D7);
690     D9 := d9(E, L, Elog, Periods, D7);
691     D10 := d10(E, L, Elog, Periods, D7);
692     vprintf Detail, 2: "d7 =  %.4o\n", D7;
693     vprintf Detail, 2: "d8 = %.4o\n", D8;
694     vprintf Detail, 2: "d9 = %.4o\n", D9;
695     vprintf Detail, 2: "d10 = %.4o\n", D10;
696
697     // Find the reduced bound for the coefficients in the linear
698     // logarithmic form
699     initQ := InitialQ(Degree(K), #L, Q0, C9, C10, D8, D9, D10, h, expTors);
700     vprintf Detail, 2: "Initial Q <= 10^%o\n", initQ;
701     initQ := 10^initQ;
702     Append(~Q, ReducedQ(ee, EmbedL, Elog, C9, C10, Periods, expTors, initQ));
703     vprint Detail, 2: "-"^45;
704   end for;
705   Q := Maximum(Q);
706   vprintf Detail: "Maximum absolute bound on coefficients = %o\n", Q;
707   return IntegralPoints(E, L, Q : tol := tol);
708 end intrinsic;
```

## Attached Files

To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.
• [get | view] (2010-06-09 19:20:06, 25.7 KB) [[attachment:intptsanf.m]]
All files | Selected Files: delete move to page copy to page

You are not allowed to attach a file to this page.