///////////////////////////////////////////////////////////////////////////////
// Computing all integral points on elliptic curves over number fields
// Based on Smart & Stephens' paper (Math. Proc. Camb. Phil. Soc. 122 (1997),
// pp. 9-16) with some modification
//
// Written by Thotsaphon Thongjunthug, BSc(Hons)
// Last update: 14 May 2009
//
// This package requires some functions on Silverman's bound implemented in
// John Cremona's "nfhtbound.m", available freely on his webpage
// http://www.warwick.ac.uk/~masgaj/ftp/progs/magma/index.html
///////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////
// Commands:
// IntegralPoints() : Find all integral points on elliptic curves over
// number fields
///////////////////////////////////////////////////////////////////////////////

// Declare printing verbose
// 0 = return result only, no detail will be printed
// 1 = minimal amount of details shown
// 2 = all details (e.g. values of all constants) shown (for debugging only)
declare verbose Detail, 2;

///////////////////////////////////////////////////////////////////////////////
// Constants for an upper bound on linear form in logarithm
// The indices of c are as in Smart & Stephens' paper
///////////////////////////////////////////////////////////////////////////////

// Compute constant c3
// Input: L = a sequence of points on elliptic curves over number field
function c3(L)
M := Eigenvalues(HeightPairingMatrix(L));
M := SetToSequence(M);
M := [m[1] : m in M]; // ignore multiplicity
M := Minimum(M); // least eigenvalue
return M;
end function;

// Compute constant c6
// Input: E = elliptic curve over real/complex numbers
function c6(E)
P<x> := PolynomialRing(ComplexField());
C4, C6 := Explode(cInvariants(E));
f := x^3 - (P! C4/48)*x - (P! C6/864);
R := Roots(f);
R := [r[1] : r in R]; // ignore multiplicity
R := [Abs(r) : r in R];
R := 2 * Maximum(R);
return R;
end function;

// Compute constant c8, and the periods of the lattice L where E(C) = C/L
// This requires complex periods of elliptic curves over C
// Require: complexell.m
// Input: E = elliptic curve over real/complex numbers
function c8(E)
org_w1, org_w2 := Explode(ComplexPeriods(E));
// Apply transformation by SL(2,Z) so that tau = w2/w1 is in the
// fundamental domain
_, A := ReducedTau(org_w2 / org_w1); // require "complexell.m"
w2 := A[1,1]*org_w2 + A[1,2]*org_w1;
w1 := A[2,1]*org_w2 + A[2,2]*org_w1;
L := [Abs(w) : w in [w1, w2, w1+w2]];
L := Maximum(L);
return L, [w1, w2];
end function;

// Compute absolute logarithmic height of an element in n-projective space
// over a number field K
// Input: X = a sequence of elements in a number field K
function AbsLogHeight(X)
// X cannot be zero vector
if #X eq #[x : x in X | x eq 0] then
error "AbsLogHeightExt: X cannot be zero vector";
end if;

// Find all prime ideal that divides some denominators of x_i
I := {};
K := Parent(X[1]); // assume each x_i is in the same field
O := RingOfIntegers(K);
for x in X do
den := O ! Denominator(x);
L := Decomposition(den);
L := {l[1] : l in L}; // ignore multiplicity
I := I join L;
end for;

h := 0;
// Non-archimedean contributions
for p in I do
M := [Rationals()| ];
for x in X do
if x eq 0 then
Append(~M, 0);
else
Append(~M, Norm(p)^(-Valuation(x, p)));
end if;
end for;
M := Maximum(M);
M := Log(M);
h +:= M;
end for;

// Archimedean contributions
s1, s2 := Signature(K);
// M = [log max(|x_1|_v,..., |x_n|_v) : v in M_K]
M := [];
for x in X do
C := Conjugates(x);
newC := [];
// Real embedding contributions
for i := 1 to s1 do
Append(~newC, Abs(C[i]));
end for;
// Complex embedding contributions
for i := 1 to s2 do
Append(~newC, Abs(C[s1+(2*i-1)])^2);
end for;

if #M eq 0 then
M := newC;
else
for i := 1 to (s1+s2) do
if newC[i] gt M[i] then
M[i] := newC[i];
end if;
end for;
end if;
end for;
M := [Log(m) : m in M];

// Overall absolute logarithmic height
h +:= (&+M);
h /:= Degree(K);
return h;
end function;

///////////////////////////////////////////////////////////////////////////////
// Constants for David's lower bound on linear form in logarithm
// Notation used as in Appendix A of Smart's book
// "The Algorithmic Resolution of Diophantine Equations", with c's being
// replaced by d's
///////////////////////////////////////////////////////////////////////////////

// Compute the "height" of elliptic curve
// Input: E = elliptic curve over number field
function h_E(E)
j := jInvariant(E);
C4, C6 := Explode(cInvariants(E));
g2 := C4/12;
g3 := C6/216;
m := Maximum([1, AbsLogHeight([1, g2, g3]), AbsLogHeight([1, j])]);
return m;
end function;

// Compute the list of modified height of a point P in E(K)
// according to embeddings
// Input: E = elliptic curve over a number field K
// P = a point in E(K)
// ElogEmbedP = elliptic logarithm of P when embedded in some
// embedding
// D7 = constant d7 (depending on embedding)
function h_m(E, P, ElogEmbedP, D7)
K := BaseRing(E);
L := [Height(P), h_E(E), D7/Degree(K)*Abs(ElogEmbedP)^2];
L := Maximum(L);
return L;
end function;

// Compute two extra h_m's based on the two periods
// Similar to h_m, but now ElogEmbedP becomes a period of the fundamental
// parallelogram of E(C) for some embedding
// Input: E = elliptic curve over a number field
// Periods = a sequence of two periods of the fundamental
// parallelogram of E(C) for some embbeding
// D7 = constant d7 (depending on embedding)
function Extra_h_m(E, Periods, D7)
D := Degree(BaseRing(E));
h := h_E(E);
tmp1 := Maximum([0, h, D7/D*Abs(Periods[1])^2]);
tmp2 := Maximum([0, h, D7/D*Abs(Periods[2])^2]);
return [tmp1, tmp2];
end function;

// Compute constant d8
// Input: E = elliptic curve over a number field K
// L = a sequence of points in E(K) (e.g. Mordell-Weil basis)
// Elog = a sequence of (pre-computed) elliptic logarithm of
// each point in L at a fixed embedding
// Periods = a sequence of two periods of the fundamental
// parallelogram of E(C) for some embbeding
// D7 = constant d7 (depending on embedding)
function d8(E, L, Elog, Periods, D7)
C := [Exp(1) * h_E(E)];
D := Degree(BaseRing(E));
for i := 1 to #L do
Append(~C, h_m(E, L[i], Elog[i], D7) / D);
end for;
C := C cat [t / D : t in Extra_h_m(E, Periods, D7)];
C := Maximum(C);
return C;
end function;

// Compute constant d9
// Input: E = elliptic curve over a number field
// L = a sequence of points in E(K) (e.g. Mordell-Weil basis)
// Elog = a sequence of (pre-computed) elliptic logarithm of
// each point in L at a fixed embedding
// Periods = a sequence of two periods of the fundamental
// parallelogram of E(C) for some embbeding
// D7 = constant d7 (depending on embedding)
function d9(E, L, Elog, Periods, D7)
D := Degree(BaseRing(E));
C := [];
for i := 1 to #L do
tmp := Exp(1) * Sqrt(D * h_m(E, L[i], Elog[i], D7) / D7);
tmp /:= Abs(Elog[i]);
C[i] := tmp;
end for;

// Take minimum among extra_h_m
Ehm := Extra_h_m(E, Periods, D7);
tmp1 := Exp(1) * Sqrt(D*Ehm[1]/D7) / Abs(Periods[1]);
tmp2 := Exp(1) * Sqrt(D*Ehm[2]/D7) / Abs(Periods[2]);
C := C cat [tmp1, tmp2];
C := Minimum(C);
return C;
end function;

// Compute constant d10
// Input: E = elliptic curve over a number field
// L = a sequence of points in E(K) (e.g. Mordell-Weil basis)
// Elog = a sequence of (pre-computed) elliptic logarithm of
// each point in L at a fixed embedding
// Periods = a sequence of two periods of the fundamental
// parallelogram of E(C) for some embbeding
// D7 = constant d7 (depending on embedding)
function d10(E, L, Elog, Periods, D7)
D := Degree(BaseRing(E));
n := #L+2;
scalar := 2 * 10^(8 + 7*n) * (2/Exp(1))^(2*n^2);
scalar *:= (n+1)^(4*n^2 + 10*n) * D^(2*n + 2);
scalar *:= (Log(d9(E, L, Elog, Periods, D7)))^(-2*n-1);
for i := 1 to #L do
scalar *:= h_m(E, L[i], Elog[i], D7);
end for;
scalar *:= &*(Extra_h_m(E, Periods, D7));
return scalar;
end function;

// Compute the right-hand side of the Principal Inequality
// Input: D = Degree(K), K = number field
// r = rank(E(K))
// C9 = constant c9
// C10 = constant c10
// D9 = constant d9
// D10 = constant d10
// h = h_E(E), E = elliptic curve over K
// Q = initial bound for the coefficients of P_i's, where
// P_i's are in the Mordell-Weil basis
function RHS(D, r, C9, C10, D9, D10, h, Q, expTors)
bound := (Log(Log(Q*r*expTors)) + h + Log(D*D9))^(r+3);
bound *:= D10*(Log(Q*r*expTors) + Log(D*D9));
bound +:= Log(C9*expTors);
bound /:= C10;
return bound;
end function;

// Approximate initial bound on Q = max_{1\le i\le r}{|q_i|}
// Input: D = Degree(K), where K = number field
// r = rank of E(K)
// Q0 = constant Q0 (in S&S paper, this is called K0)
// C9 = constant c9
// C10 = constant c10
// D8 = constant d8 (from function d8())
// D9 = constant d9 (from function d9())
// D10 = constant d10 (from function d10())
// h = h_E(E)
// expTors = exponent of the torsion subgroup of E(K)
// Revised: 5 May 2009
function InitialQ(D, r, Q0, C9, C10, D8, D9, D10, h, expTors)
minQ := Maximum(Q0, Exp(D8));

// Try to approximate Q such that Q^2 = RHS(Q) (i.e. Q makes both sides
// of the Principal Inequality equal)
// Firstly, set a guess for Q, say minQ + 1 (so that Q > minQ)
// For simplicity, let's round Q up to the nearest power of 10
Q := minQ + 1;
x := Ceiling(Log(10, Q)); // x = log_10(Q)

// Check if Q satisfies the Principal Inequality, i.e. if Q^2 < RHS(Q)
// If so, repeat with the larger Q until we find the first Q that
// violates the P.I.
// N.B. This loop will eventually terminate
exp_OK := 0; // the exponent that satisfies P.I.
exp_fail := 0; // the exponent that fails P.I.
while 10^(2*x) lt RHS(D, r, C9, C10, D9, D10, h, 10^x, expTors) do
exp_OK := x; // Principal Inequality satisfied
x *:= 2; // double x, and retry
end while;
exp_fail := x; // x that fails the Principal Inequality

// So now x = log_10(Q) must lie between exp_OK and exp_fail
// Refine x further using "binary search"
repeat
x := Floor((exp_OK + exp_fail)/2);
if 10^(2*x) ge RHS(D, r, C9, C10, D9, D10, h, 10^x, expTors) then
exp_fail := x;
else
exp_OK := x;
end if;
until (exp_fail - exp_OK) le 1;
return exp_fail; // over-estimate
end function;

// Reduce the bound Q by LLL reduction until no further improvement
// is possible. This function initially requires high precision to
// proceed, although this should be done automatically by now
// Require: complexell.m
// Input: EmbedE = model of E embedded in some (complex) embedding
// EmbedL = a sequence of points on EmbedE (e.g. points in the
// Mordell-Weil basis when embedded into EmbedE)
// Elog = a sequence of (pre-computed) elliptic logarithm of
// each point in EmbedL
// C9 = constant c9
// C10 = constant c10
// Periods = sequence of the double periods of the fundamental
// parallelogram associated to EmbedE
// expTors = exponent of the torsion subgroup of E(K), K = number field
// initQ = initial guess for Q to be reduced by LLL
function ReducedQ(EmbedE, EmbedL, Elog, C9, C10, Periods, expTors, initQ)
w1, w2 := Explode(Periods);
r := #EmbedL;
newQ := initQ;

// Repeat LLL reduction until no further reduction is possible
repeat
Q := newQ;
S := r*(Q^2)*(expTors^2);
T := 3*r*Q*expTors/Sqrt(2);

// Create the basis matrix
C := 1;
repeat
C *:= Q^Ceiling((r+2)/2);
L := ZeroMatrix(Integers(), r+2, r+2);

// Elliptic logarithm should have precision "suitable to" C
// e.g. If C = 10^100, then Re(Elog[i]) should be computed
// correctly to at least 100 decimal places
pow10_C := Ceiling(Log(10, C));
if pow10_C gt Precision(Elog[1]) then
vprint Detail, 2:
"Need higher precision, recompute elliptic logarithm ...";
// Re-compute elliptic logarithm to the right precision
Elog := [ComplexEllLog(EmbedE, P : Prec := pow10_C+10) : P in EmbedL];
vprint Detail, 2: "Elliptic logarithm recomputed";
end if;
// Assign all non-zero entries
for i := 1 to r do
L[i,i] := 1;
L[r+1, i] := Truncate(C*Re(Elog[i]));
L[r+2, i] := Truncate(C*Im(Elog[i]));
end for;
L[r+1, r+1] := Truncate(C*Re(w1));
L[r+1, r+2] := Truncate(C*Re(w2));
L[r+2, r+1] := Truncate(C*Im(w1));
L[r+2, r+2] := Truncate(C*Im(w2));
L := Transpose(L); // In Magma, basis is spanned by row vector!

// LLL reduction and constants
L := LLL(L);
b1 := L[1]; // 1st row of reduced basis
// Norm(b1) = square of Euclidean norm
lhs := 2^(-r-1)*Norm(b1) - S;
until (lhs ge 0) and (Sqrt(lhs) gt T);

newQ := ( Log(C*C9*expTors) - Log(Sqrt(lhs) - T) ) / C10;
newQ := Floor(Sqrt(newQ));
pow10 := Floor(Log(10, C));
vprintf Detail, 2: "Choose C = %.4o x 10^%o. ", 1.*C/10^pow10, pow10;
vprintf Detail, 2: "After reduction, Q <= %o\n", newQ;
until ((Q - newQ) le 1) or (newQ le 1);
return newQ;
end function;

///////////////////////////////////////////////////////////////////////////////
// Main intrinsic functions
///////////////////////////////////////////////////////////////////////////////

// Search for all integral points on elliptic curves over number fields
// within a given bound
// Input: E = elliptic curve over a number field K
// L = a sequence of points in the Mordell-Weil basis for E(K)
// Q = a maximum for the absolute bound on all coefficients
// in the linear combination of points in L
// Output: S1 = sequence of all integral points on E(K) modulo [-1]
// S2 = sequence of tuples representing the points as a
// linear combination of points in L
// Option: tol = tolerance used for checking integrality of points.
// (Default = 0 - only exact arithmetic will be performed)
intrinsic IntegralPoints(E::CrvEll, L::[PtEll], Q::RngIntElt : tol := 0) ->
SeqEnum, SeqEnum
{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.}

// Check input validity
require IsNumberField(BaseRing(E)):
"Elliptic curve must be defined over a number field";
require tol ge 0: "Tolerance must be non-negative";

// Find the generators of the torsion subgroup of E(K)
Tors, map := TorsionSubgroup(E);
expTors := Exponent(Tors);
G := Generators(Tors);
Tors := [map(g) : g in G]; // each generator of E(K)_tors
OrdG := [Order(g) : g in G]; // order of each generator

if (#L eq 0) and (#Tors eq 0) then
return [], [];
end if;

// Create all possible r+#Tors-tuples
r := #L; // r = rank of E(K)
C := [0 : i in [1..(r + #Tors)]];
ListC := [];
for i := 0 to Q do
C[1] := i;
Append(~ListC, C);
end for;

for i := 2 to r do
tmp := [];
for j := 1 to Q do
for l in ListC do
tup := l;
tup[i] := j;
Append(~tmp, tup);

// Avoid having its negative in the list
// Only use when all previous entries in tuple are zero
for k := 1 to i-1 do
if tup[k] ne 0 then
tup[i] := -j;
Append(~tmp, tup);
break;
end if;
end for;

end for;
end for;
ListC := ListC cat tmp;
end for;

// Add torsion point, if any
if #Tors ne 0 then
for i := 1 to #Tors do
tmp := [];
for j := 1 to (OrdG[i]-1) do
for l in ListC do
tup := l;
tup[r+i] := j;
Append(~tmp, tup);
end for;
end for;
ListC := ListC cat tmp;
end for;
end if;
Remove(~ListC, 1); // remove point at infinity

L := L cat Tors;
vprint Detail, 2: "Generators = ", L;
PtsList := [];
CoeffList := [];

// Skip the complex arithmetic and only perform exact arithmetic if tol = 0
if tol eq 0 then
vprint Detail : "Exact arithmetic";
for l in ListC do
P := &+[l[i]*L[i] : i in [1..#L]];
if IsIntegral(P[1]) and IsIntegral(P[2]) then
vprintf Detail: "%o ---> %o\n", l, P;
Append(~PtsList, P);
TupList := [ <L[i], l[i]> : i in [1..#L] | l[i] ne 0 ];
Append(~CoeffList, TupList);
end if;
end for;
vprint Detail: "*"^45;
return PtsList, CoeffList;
end if;

// Suggested by John Cremona
// Point search. This is done via arithmetic on complex points on each
// embedding of E. Exact arithmetic will be carried if the resulting
// complex points are "close" to being integral, subject to some tolerance

// Embed each generator of the torsion subgroup
X := [Conjugates(P[1]) : P in (L cat Tors)];
Y := [Conjugates(P[2]) : P in (L cat Tors)];
// Create all embeddings of E
K := BaseRing(E);
s1, s2 := Signature(K);
a1, a2, a3, a4, a6 := Explode(aInvariants(E));
a1 := Conjugates(a1);
a2 := Conjugates(a2);
a3 := Conjugates(a3);
a4 := Conjugates(a4);
a6 := Conjugates(a6);
EmbedEList := [EllipticCurve([a1[j], a2[j], a3[j], a4[j], a6[j]]): j in
[1..(s1+2*s2)]];

// Set tolerance - This should be larger than 10^(-current precision) to
// avoid missing any integral points. Too large tolerance will slow the
// computation, too small tolerance may lead to missing some integral points.
// This is now made as an option.
vprint Detail: "Tolerance = ", tol * 1.;

// Create the matrix containing all embeddings of the integral basis of K
// as its columns
IntBasis := IntegralBasis(K);
B := Matrix([Conjugates(a) : a in IntBasis]);
// Note that B is always invertible, so we can take its inverse
B := B^(-1);

// Embeddings of each point in L
// Format: [[P1 ... P1], [P2 ... P2], ...]
EmbedLList := [];
for i := 1 to (r+#Tors) do
EmbedLi := [];
for j := 1 to (s1+2*s2) do
P_i := Points(EmbedEList[j], X[i][j])[1];
// Choose the right sign for the y-coordinate
if (Y[i][j] ne 0) and (Re(P_i[2]/Y[i][j]) lt 0) then
P_i := -P_i;
end if;
Append(~EmbedLi, P_i);
end for;
Append(~EmbedLList, EmbedLi);
end for;

// Point search
for l in ListC do
// Compute P by complex arithmetic for every embedding
EmbedPList := [E![0,1,0] : E in EmbedEList];
for j := 1 to (s1+2*s2) do
EmbedPList[j] := &+[l[i]*EmbedLList[i][j] : i in [1..(r + #Tors)]];
end for;

// Check if the x-coordinate of P is "close to" being integral
// If so, compute P exactly and check if it is integral; skip P otherwise
XofP := Matrix([[P[1] : P in EmbedPList]]);
// Write x(P) w.r.t. the integral basis of (K)
// Due to the floating arithmetic, some entries in LX may have very tiny
// imaginary part, which can be thought as zero
LX := XofP * B;
LX := [Abs( Re(LX[1,i]) - Round(Re(LX[1,i])) ): i in [1..#IntBasis]];
LX := &and[dx lt tol : dx in LX];
if not LX then
// x-coordinate of P is not integral, skip P
continue;
end if;

// Now check P by exact arithmetic
// Add P and the list of tuples representing P into the list
// if P is integral
P := &+[l[i]*L[i] : i in [1..#L]];
if IsIntegral(P[1]) and IsIntegral(P[2]) then
vprintf Detail: "%o ---> %o\n", l, P;
Append(~PtsList, P);
TupList := [ <L[i], l[i]> : i in [1..#L] | l[i] ne 0 ];
Append(~CoeffList, TupList);
end if;
end for;
vprint Detail: "*"^45;
return PtsList, CoeffList;
end intrinsic;

// Compute all integral points on elliptic curve over a number field
// This function simply computes a suitable bound Q, and return
// IntegralPoints(E, L, Q : tol := ...).
// Input : E = elliptic curve over a number field K
// L = a sequence of points in the Mordell-Weil basis for E(K)
// Output: S1 = sequence of all integral points on E(K) modulo [-1]
// S2 = sequence of tuples representing the points as a
// linear combination of points in L
// Option: tol = tolerance used for checking integrality of points.
// (Default = 0 - only exact arithmetic will be performed)
intrinsic IntegralPoints(E::CrvEll, L::[PtEll] : tol := 0) -> SeqEnum, SeqEnum
{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 := ...)}

// Check input validity
require IsNumberField(BaseRing(E)):
"Elliptic curve must be defined over a number field";
require tol ge 0:
"Tolerance must be non-negative";

if #L eq 0 then
return IntegralPoints(E, [], 0 : tol := tol);
end if;

a1, a2, a3, a4, a6 := Explode(aInvariants(E));

// Embed E into various (real/complex) embeddings
a1 := Conjugates(a1);
a2 := Conjugates(a2);
a3 := Conjugates(a3);
a4 := Conjugates(a4);
a6 := Conjugates(a6);
K := BaseRing(E);
s1, s2 := Signature(K);
b2 := Conjugates(bInvariants(E)[1]);
pi := Pi(RealField());

// Embed generators in the Mordell-Weil basis
X := [Conjugates(P[1]) : P in L];
Y := [Conjugates(P[2]) : P in L];

// Find the generators of the torsion subgroup of E(K)
Tors, map := TorsionSubgroup(E);
expTors := Exponent(Tors);
G := Generators(Tors);
Tors := []; OrdG := [];
for g in G do
Append(~Tors, map(g)); // torsion point
Append(~OrdG, Order(g)); // order of torsion point
end for;

// Global constants (independent to the embedding of E)
C2 := -SilvermanLowerHeightBound(E); // require "nfhtbound.m"
C3 := c3(L);
h := h_E(E);
vprint Detail, 2: "Global constants";
vprintf Detail, 2: "c2 = %.4o\n", C2;
vprintf Detail, 2: "c3 = %.4o\n", C3;
vprintf Detail, 2: "h_E = %.4o\n", h;
vprint Detail, 2: "-"^45;

Q := [];
// Find the most reduced bound on each embedding of E
// But first let's adjust the index
for i := 1 to (s1+s2) do
if i le s1 then
j := i;
nv := 1;
vprintf Detail, 2: "Real embedding #%o\n", j;
else
j := s1 + (2*(i-s1)-1);
nv := 2;
vprintf Detail, 2: "Complex embedding #%o\n", i-s1;
end if;

// Create complex embedding of E
ee := EllipticCurve([a1[j], a2[j], a3[j], a4[j], a6[j]]);

// Embedding of all points in Mordell-Weil basis
// Find complex elliptic logarithm of each embedded point
// Let's set precision initially to, say, 1000
EmbedL := [[X[i][j], Y[i][j]] : i in [1..#L]];
Elog := [ComplexEllLog(ee, P: Prec := 1000) : P in EmbedL];
// Local constants (depending on embedding)
// C9, C10 are used for the upper bound on the linear form in logarithm
C4 := C3 * Degree(K) / (nv*(s1+s2));
vprintf Detail, 2: "c4 = %.4o\n", C4;
C5 := C2 * Degree(K) / (nv*(s1+s2));
vprintf Detail, 2: "c5 = %.4o\n", C5;
C6 := c6(ee);
vprintf Detail, 2: "c6 = %.4o\n", C6;
delta := 1 + (nv-1)*pi;
C8, Periods := c8(ee);
vprintf Detail, 2: "c8 = %.4o\n", C8;
C7 := 8*(delta^2) + (C8^2)*Abs(b2[j])/12;
vprintf Detail, 2: "c7 = %.4o\n", C7;
C9 := C7*Exp(C5/2);
vprintf Detail, 2: "c9 = %.4o\n", C9;
C10 := C4/2;
vprintf Detail, 2: "c10 = %.4o\n", C10;
Q0 := Sqrt( ( Log(C6+Abs(b2[j])/12) + C5) / C4 );
vprintf Detail, 2: "Q0 = %.4o\n", Q0;

// Constants for David's lower bound on the linear form in logarithm
w1, w2 := Explode(Periods); // N.B. periods are already in "standard" form
D7 := 3*pi / ((Abs(w2)^2) * Im(w2/w1));
D8 := d8(E, L, Elog, Periods, D7);
D9 := d9(E, L, Elog, Periods, D7);
D10 := d10(E, L, Elog, Periods, D7);
vprintf Detail, 2: "d7 = %.4o\n", D7;
vprintf Detail, 2: "d8 = %.4o\n", D8;
vprintf Detail, 2: "d9 = %.4o\n", D9;
vprintf Detail, 2: "d10 = %.4o\n", D10;

// Find the reduced bound for the coefficients in the linear
// logarithmic form
initQ := InitialQ(Degree(K), #L, Q0, C9, C10, D8, D9, D10, h, expTors);
vprintf Detail, 2: "Initial Q <= 10^%o\n", initQ;
initQ := 10^initQ;
Append(~Q, ReducedQ(ee, EmbedL, Elog, C9, C10, Periods, expTors, initQ));
vprint Detail, 2: "-"^45;
end for;
Q := Maximum(Q);
vprintf Detail: "Maximum absolute bound on coefficients = %o\n", Q;
return IntegralPoints(E, L, Q : tol := tol);
end intrinsic;