attachment:intptsanf.m of days22/cremona

   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[1] : 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[1] : 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[1]); // 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[1] : 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[1])^2]);
 187   tmp2 := Maximum([0, h, D7/D*Abs(Periods[2])^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[1]/D7) / Abs(Periods[1]);
 230   tmp2 := Exp(1) * Sqrt(D*Ehm[2]/D7) / Abs(Periods[2]);
 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[1]) 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[1]; // 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[1] + ... + 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[1] := 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[1]) and IsIntegral(P[2]) 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[1]) : P in (L cat Tors)];
 505   Y := [Conjugates(P[2]) : 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])[1];
 538       // Choose the right sign for the y-coordinate
 539       if (Y[i][j] ne 0) and (Re(P_i[2]/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[1] : 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[1]) and IsIntegral(P[2]) 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)[1]);
 618   pi := Pi(RealField());
 619 
 620   // Embed generators in the Mordell-Weil basis
 621   X := [Conjugates(P[1]) : P in L];
 622   Y := [Conjugates(P[2]) : 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;