A numerical transcendental method in algebraic geometry

Based on high precision computation of periods and lattice reduction techniques, we compute the Picard group of smooth surfaces. We also study the lattice reduction technique that is employed in order to quantify the possibility of numerical error in terms of an intrinsic measure of complexity of each surface. The method applies more generally to the computation of the lattice generated by Hodge cycles of middle dimension on smooth projective hypersurfaces. We demonstrate the method by a systematic study of thousands of quartic surfaces (K3s) defined by sparse polynomials. As an application, we count the number of rational curves of a given degree lying on each surface. For quartic surfaces we also compute the endomorphism ring of their transcendental lattice.


Introduction.
In "A transcendental method in algebraic geometry" [30], Griffiths emphasized the role of certain multivariate integrals, known as periods, "to construct a continuous invariant of arbitrary smooth projective varieties". Periods often determine the projective variety completely and therefore its algebraic invariants. Translating periods into discrete algebraic invariants is a difficult problem, exemplified by the long standing Hodge conjecture which describes how periods determine the algebraic cycles within a projective variety.
Recent progress in computer algebra makes it possible to compute periods with high precision and put transcendental methods into practice. We focus mainly on algebraic surfaces and give a numerical method to compute Picard groups. As an application, we count smooth rational curves on quartic surfaces using the Picard group. These methods apply more generally to hypersurfaces in a projective space of arbitrary dimension and we show some examples.
Structure of the Picard group and main results. There are many curves in a smooth surface X ⊂ P 3 C , the basic ones are those obtained by intersecting X with another surface S in P 3 C . If S 1 and S 2 are two surfaces of the same degree then the curve C 1 = S 1 ∩ X can be deformed into the curve C 2 = S 2 ∩ X by varying continuously the coefficients of the defining equation of S 1 . The curves C 1 and C 2 are said to be linearly equivalent. The notion of linear equivalence extends to formal Z-linear combinations of curves and the Picard group of X is defined by Pic(X) def =Z algebraic curves in X / linear equivalence relations .
The Picard group is an algebraic invariant that reflects the nature of the algebraic curves lying on X. It is a free abelian group, that is, Pic(X) Z ρ for a positive integer ρ called the Picard number of X. As Zariski wrote, "The evaluation of ρ for a given surface presents in general grave difficulties" [67, p. 110].
There is more to the Picard group than the Picard number. The intersection product, which for any two curves C 1 and C 2 in X associates an integer C 1 · C 2 , induces a bilinear map Pic(X) × Pic(X) → Z. The intersection product is an intrinsic algebraic invariant of X that is finer than the Picard number. There is also an extrinsic invariant in the Picard group, called the polarization, recording much of the geometry of X within P 3 C . The polarization is the linear equivalence class of any curve obtained by intersecting X with a plane in P 3 C . The problem we address is then the following: Given the defining equation of X, compute the Picard number ρ of X, the ρ × ρ matrix of the intersection product and the ρ coordinates of the polarization in some basis of Pic(X) Z ρ .
We approach the problem using transcendental methods, that is, we use the complex geometry of the hypersurfaces and compute multivariate integrals on topological cycles, namely the periods. For surfaces, Lefschetz (1, 1) theorem identifies the Picard group of a surface with the lattice of integer linear relations between periods. The rank, intersection product and polarization of the Picard group can be computed from a high precision computation of the periods [57] and well-established techniques in lattice reduction. We apply these techniques also to the computation of the endomorphism ring of the transcendental lattice in order to compute Charles' gap [10], see below. Counting rational curves of a given degree lying on a surface is an interesting application of the computation of the Picard group with its intersection product and polarization.
The method extends to higher dimensional hypersurfaces in order to compute the group of Hodge cycles. For a hypersurface X in a projective space of odd dimension P 2k+1 C with k > 1 there are two interesting objects to study, replacing the Picard group for surfaces under present discussion: the group of algebraic cycles Alg k (X) generated by the cohomology classes of k-dimensional algebraic subvarieties of X or the group of Hodge cycles Hdg k (X) generated by integral linear relations between periods. The Hodge conjecture states in greater generality that, after tensoring with rational numbers, the two groups Alg k (X) ⊗ Z Q and Hdg k (X) ⊗ Z Q coincide [19]. The resolution of this conjecture is one of the seven Millennium Prize Problems posed by the Clay Institute. Let us point out that this problem is far from being resolved even for hypersurfaces beyond the case of surfaces in P 3 C . Perhaps the current method will allow for experimentation in this direction with the ability to compute Hdg k (X) together with its intersection product and polarization, see §6.
Related work. For surfaces of degree at most three the Picard group does not depend on the defining equations and the main arguments to compute it have been known since the 19th century [20]. Starting with surfaces of degree four, the Picard group is sensitive to the defining equations and poses an entirely different kind of challenge, where a complete solution must be algorithmic in nature. Noether and Lefschetz [41] proved that a very general quartic surface is expected to have Picard number one. However, the first quartic with ρ = 1 defined by a polynomial with integer coefficients appeared in 2007 with van Luijk's seminal paper [63] where he used techniques involving reduction to finite characteristic. Since then the reduction techniques have been going through a phase of rapid development. The original argument of van Luijk was refined by Elsenhans and Jahnel [22,21] but the computational bottle neck persisted: in working with surfaces over finite fields the computation of their Zeta function initially required the expensive process of point counting. This bottleneck has been alleviated using ideas from p-adic cohomology and gave rise to two different approaches: one dealing directly with the surface [37,1,16] and another which deforms the given surface to a simpler one [40,50].
To complement the upper bounds coming from prime reductions, lower bounds on the Picard number can be obtained, at least in theory, by enumerating all the algebraic curves in X. One could compute infinite sequences of lower and upper bounds that may eventually determine the Picard number. However, Charles [10] proved that the upper bounds obtained from finite characteristic could significantly overestimate the Picard number and he expressed the gap in terms of the endomorphism ring of the transcendental lattice, see §2.3. In practice, computing this gap appears to be just as difficult as the computation of the Picard number. However, Charles demonstrated at last that the Picard number of a K3 surface (e.g. a quartic surface) defined with coefficients in a number field is computable. On the theoretical side, effective algorithms have been developed with a broader reach but with low practicability [10,32,53]. There is recent work addressing the issue of practicability [27].
Concerning numerical methods, high precision computation of periods has been successfully applied to many problems concerning algebraic curves e.g. [64,8], even with the possibility of a posteriori symbolic certification [17].
Tools. For the purpose of exposition, we focus mainly on quartic surfaces. Our techniques, as well as our code, work for hypersurfaces of any degree and dimension, given sufficient computational resources. The main computational tool on the one hand is an algorithm to approximate periods [57], based on Picard-Fuchs differential equations [52] and algorithms to compute them [13, 38, 49, 47, 39, 7, e.g.] via numerical analytic continuation [62,61,12,44,45]. On the other hand, we use algorithms to compute integer linear relations between vectors of real numbers [26,42,33,9,25,11].
The reliability of numerical computations. Although the periods vary continuously with the coefficients of the defining equation of a surface, the Picard number is nowhere continuous for surfaces of degree at least 4 [14], behaving like the indicator function of the rational numbers on the real line. This fact suggests that deducing the Picard group from approximate periods must be hopeless, since the numerical computation of the Picard group is based on finding integer relations between real numbers that we know only with finite precision. However, working with sparse polynomials with small rational coefficients, we observe a remarkable tolerance for error, see §2.2 for a typical example. Our computations agree with the literature whenever a check is possible. In particular, we compared our Picard number computations against controlledreduction 4 [15] and Shioda's algorithm for Delsarte surfaces [58].
The possibility of error and its nature is quantified precisely in § §4-5 and compared to an intrinsic measure of complexity of the surface. This measure of complexity is out of reach but its determination would allow certification of the numerical computation.
Outline. In Section 2, we describe the computation of the Picard group of surfaces and the endomorphism ring of the transcendental lattice of quartic surfaces from an approximation of periods. In Section 3 we apply our computations to count smooth rational curves in quartic surfaces. In Section 4 we describe and analyze a standard procedure to recover integer relations between approximate real vectors. In Section 5 we quantify the nature of error in a way that is independent of the methods employed and express it in terms of an intrinsic measure of complexity of the given surface. Section 6 explains the situation for higher dimensional hypersurfaces where the general idea of the method as explained in Section 2 applies verbatim. Section 7 summarizes the experimental results obtained for thousands of quartic surfaces. Section 8 explains how we compute the polarization by fleshing out the argument given in [57].

Periods and Picard group.
2.1. Principles. Following Picard, Lefschetz and Hodge, algebraic curves on a smooth complex surface X can be characterized among all topological 2-dimensional cycles of X in terms of multivariate integrals (for a historical perspective, see [46] and [34]).
Let X ⊂ P 3 be a smooth complex surface. An algebraic curve C ⊂ X is supported on a topological 2-dimensional cycle. Lefschetz proved that two algebraic curves are topologically homologous if and only if they are linearly equivalent [41], see also [46,Chap. 9]. In other words, the Picard group comes with a natural inclusion into the homology group This homology group is a topological invariant that depends only on the degree of X, while Pic(X) is a much finer invariant of X.
Recall that for a 2-dimension cycle γ and any holomorphic differential 2-form ω on X, the integral γ ω is well defined on the homology class of γ. When X has degree 4, H 2 (X, Z) Z 22 and X admits a unique non-zero holomorphic 2-form up to scaling, which we denote ω X . Given a basis γ 1 , . . . , γ 22 of H 2 (X, Z), Theorem 2.1 rewords as The integrals γ i ω X appearing here are called the periods of X. All periods can be expressed as a sum of integrals in the affine chart we can form a thin tube τ ⊂ C 3 \ X around γ so that where f is the degree 4 homogeneous polynomial defining X [29].
In general, when X has degree d ≥ 4, H 2 (X, Z) has rank m = d 3 − 4d 2 + 6d − 2 and the space V of holomorphic 2-forms on X is of dimension r = d−1 3 . Fixing bases H 2 (X, Z) = Z γ 1 , . . . , γ m , V = C ω 1 , . . . , ω r and applying Lefschetz (1,1) theorem we get: In view of (2.2) and (2.4) we can determine Pic(X) by computing the matrix of periods [ γ j ω i ] i,j and then finding integer linear relations between the rows. The algorithm presented in [57] to compute the periods of X takes care of the first step. We may then use lattice reduction algorithms [42, p. 525] to find generators for Pic(X).
We briefly recall how periods of X are computed in [57]. The surface X is put into a single parameter family of surfaces containing the Fermat surface The matrix of periods along the family vary holomorphically in terms of the parameter and these entries satisfy ordinary differential equations which are computed exactly. The value of this one parameter period matrix-as well as its derivatives-at Y are given by closed formulas involving Gamma functions. The differential equations together with the periods on Y expresses the periods of X as the solution to an initial value problem. This initial value problem is solved using Mezzarobba's implementation of numerical analytic continuation [45] to arbitrary precision with rigorous error bounds.
The reconstruction of integer relations between transcendental numbers that are only approximately given is not possible in general. However, when the transcendental numbers are well behaved, this reconstruction may be possible. We devote §4 to the study of this problem.

2.
2. An example. Consider the quartic surface X ⊂ P 3 defined by the polynomial As described above, we may compute a 1 × 22 matrix of periods to arbitrary precision. In this example, the differential equations that arise are of order 5 with polynomial coefficients of degree at most 59. On a laptop, the determination of this differential equation takes about two seconds and it takes 30 seconds to integrate it to 100 digits of accuracy, with rigorous error bounds. All of this can be done with the command PeriodHomotopy([f]) using the package period-suite written for Magma [6] and utilizing ore algebra-analytic [45]. Applying the LLL algorithm to these approximate periods of X gives a basis of integer relations between the approximate periods. More precisely, we consider for the lattice Λ of integer vectors (u, v, a 1 , . . . , a 22 ) ∈ Z 24 satisfying  where [−] denotes the rounding to nearest integer. Equation (2.6) should be compared with (2.2). Short vectors in Pic(X) give rise to short vectors in Λ, and a short vector in Λ is likely to come from a vector in Pic(X), unless a suprising numerical cancellation happens.
Concerning the example, Figure 2.1 shows a matrix whose columns form a LLL-reduced basis for the lattice Λ. We observe an important gap in size, between the 14th and 15th column. We conclude that the Picard number of X is most likely 14 and that the columns of the lower left 22 × 14 submatrix is a basis of Pic(X). The norm of the first dismissed column, about 10 25 , fits precisely the expected situation described in Proposition 4.3.
This numerical approach may fail in two ways: either by missing a relation or by returning a false relation which nevertheless holds up to high precision. Proposition 4.1 quantifies the way in which such a failure may occur: either the computation of Pic(X) is correct; or Pic(X) is not generated by elements of norm < 10 20 ; or there is some (a i ) ∈ Z m with i a 2 i ≤ 4 such that i a i γ i ω X is not zero but smaller than 10 −99 . Section 5 expresses these quantities in terms of an intrinsic norm that is independent from the coordinates that were used to carry out the computations.
The intersection product on Pic(X) is readily computed from the generators, as we now describe. The basis of homology on X is obtained by carrying a basis from the Fermat surface Y by parallel transport. On Y the intersection numbers γ i · γ j are known exactly as well as the polarization, i.e., the coordinates of the homology class of a general place section H ∩ X in the basis As these values remain constant during parallel transport, we know the intersection product on the homology of X as well as the polarization. Computing the intersection product of the 14 generators of Pic(X) in homology, we obtain the intersection product on Pic(X). Since the polarization lies in Pic(X) we express it in terms of these generators of Pic(X). The result of this operation is displayed in Figure 2.2. The command HodgeLattice of period-suite performs all the operations starting from the computation of the periods. Applying standard methods to be discussed in §3, we find from the Picard lattice of X that there are 4 lines, 102 quartic curves and no twisted cubics inside X.

Transcendental lattice and reduction to finite characteristic.
Definition and properties. Let X be a quartic surface. Beyond the Picard group of X, we can compute its transcendental lattice and its endomorphism ring. The transcendental lattice of X is defined as which is a Q-linear space of dimension 22−rk Pic(X). Furthermore, let T C = T ⊗C ⊂ H 2 (X, C) and observe ω X ∈ T C by (2.2). The endomorphism ring E ⊂ End Q (H 2 (X, Q)) is defined as the subring of all linear maps e such that e(ω X ) ∈ C ω X , where e has been extended canonically to H 2 (X, C). The map ϕ : E → C defined by e(ω X ) = ϕ(e)ω X is an injective ring morphism and every element in E is invertible, therefore E is a number field [35, Corollary 3.3.6]. In fact, E is either totally real or a CM-field [66], see also [35]. Charles [10] determined in terms of E the overestimation of reduction methods to compute the Picard number of K3 surfaces. We give here a quick overview, see [35,60] for further results. Although we state these results for quartic surfaces over Q much of it holds for any K3 surface over a number field.
If the quartic X ⊂ P 3 is defined by a polynomial f with integer coefficients, we may consider for all but finitely many prime p the smooth quartic surface X p defined over F p by the reduction of f modulo p. Let ρ and ρ p denote the (geometric) Picard numbers of X and X p respectively. Let ρ red be the minimum of the set {ρ p | p prime and X p smooth}. The starting point of reduction methods is the inequality ρ ≤ ρ red and the relative ease with which the numbers ρ p are computed. A key issue is to determine whether ρ = ρ red .
Although ρ can be either even or odd, ρ p is always even. This issue was partially overcome by van Luijk [63] who gave necessary conditions for ρ = ρ red . He used his argument to give the first example of a K3 surface defined over the rationals with Picard number 1 by exhibiting a surface X with ρ red = 2 that does not satisfy his necessary condition. It was asked by Elsenhans and Jahnel [23] whether ρ = ρ red if ρ is even and ρ = ρ red − 1 if ρ is odd. Charles [10] settled the question in the negative.
Theorem 2.2 (Charles). The equality ρ red = ρ holds unless E is totally real and the dimension of T over E is odd, in which case ρ red = ρ + dim Q E.
Computation. From the numerical computation of periods, we obtain approximations For e ∈ End Q (H 2 (X, Q)) the condition e ∈ E can be rewritten as: Writing A = (a 1 , . . . , a 22 ) t for the coefficient vector of ω X , we can compute the endomorphism ring E via the following formulation: Just as with the computation of Pic(X), the problem of computing E is now a problem of computing integer solutions to linear equations with approximate real coefficients. We approach it once again with lattice reduction algorithms, see §4. Examples are provided in §7.

3.
Smooth rational curves in K3 surfaces. The data of the matrix of the intersection product in some basis of the Picard group of a smooth quartic surface X ⊂ P 3 together with the coordinates of the class of hyperplace section in the same basis (as in Figure 2.2) is enough to count all smooth rational curves of a given degree lying on X. We describe an algorithm here that, at least in broad strokes, seems to be folklore. 5 In principle, smooth rational curves in a surface can be enumerated using purely symbolic methods and for lines this process is routine. However, it is a challenge to enumerate even the quadric curves in quartic surfaces, let alone higher degree curves in higher degree surfaces. The computation of the Picard group offers an indirect solution to this problem.
Fix a smooth quartic X ⊂ P 3 and for each positive integer d let R d be the set of all smooth rational curves of degree d lying in X. In order to compute the cardinality of the set R d we will first observe that a smooth rational curve is completely determined by its linear equivalence class. Recall that we denote by h X ∈ Pic(X) the class of a hyperplane section.
Lemma 3.1. A smooth rational curve in X is isolated in its linear equivalence class. Moreover, the map R d → Pic(X) which maps a rational curve to its linear equivalence class injects Proof. Let C ∈ R d and D = [C] ∈ Pic(X). As C is of degree d, it intersects a general hyperplane in d points so that C · h X = d. Recall that the canonical class K X of the K3 surface X is trivial so that adjunction formula reads Now we show that C is isolated in its linear system. Indeed, if C is a curve linearly equivalent to but different from C, then the intersection number [C] · [C ] must be positive, as this number can be obtained by counting the points in C ∩ C with multiplicity. This leads to a contradiction: Note that when d = 1 there are no constraints and we have N 1 = M 1 .
Proof. For any two distinct irreducible curves C and C we have C · C ≥ 0. Upon taking C ∈ R d and C ∈ R d for d < d we see that R d injects in to N d . Now take any D ∈ N d . From the Riemann-Roch theorem for surfaces [31, V.1.6] we get: so that either D or −D must be effective. Since D · h X > 0, −D can not be effective and therefore D must be. Let us write D as a sum of classes of distinct irreducible curves i n i C i with n i > 0. Since D 2 < 0 there exists an index i such that C i · D < 0. Moreover, C i · C j ≥ 0 for every j = i, so C 2 i < 0. By adjunction formula, we conclude that C i must be a smooth rational curve [31, Ex. IV.1.8]. Furthermore, let d = C i · h X and observe d ≤ d. By definition of N d we must have d = d and therefore D = C i . Therefore, R d surjects onto N d . Proposition 3.2 implies that in order to compute the cardinality of the set R d it suffices to compute the set N d (see Algorithm 3.1). The latter can be easily computed from the sets M d for d ≤ d. We now reduce the computation of M d for each d > 0 to the enumeration of all vectors of a given norm in a lattice with a negative definite quadratic form.
Let Pic 0 (X) = {D ∈ Pic(X) | D · h X = 0}. The intersection product on Pic 0 (X) is negative definite [35,Proposition 1.2.4]. Recalling that h 2 X = 4, we define a map π : The map π maps M d bijectively on to the following set: The inverse map M d → M d is given by E → 1 4 (E + dh X ). In order to compute M d we first find the finitely many elements E ∈ Pic 0 (X) of norm −(32+ 4d 2 ), for example using KFP algorithm [36,28]. Then, among all such E, we select those Algorithm 3.1 Finding rational curve classes on smooth quartic surfaces.

Input.
The Picard group (i.e. matrix of the interction product in some basis and the coordinates of the class of hyperplane section in the same basis) of a smooth quartic surface X ⊂ P 3 ; an integer d > 0. Output. The set {[C] ∈ Pic(X) | C ⊂ X is a smooth rational curve of degree d}.
We find that X has Picard number 18 with the following representation of (Pic(X), h X ): Applying Algorithm 3.1 we see that there are 16 lines, 288 quadrics and 1536 twisted cubics as determined by this lattice of X. The 16 lines, and their incidence correspondence, as we compute from this lattice are in agreement with what we can compute rigorously using symbolic methods.

Numerical reconstruction of integer relations.
In view of (2.2) and (2.4), recovering Pic(X) boils down to finding integer linear relations between the period vectors. With the methods employed here, a finite but high enough precision will successfully recover Pic(X). It seems difficult to decide if a given precision is "high enough". Instead, we will study the process of finding linear relations between approximate vectors of real numbers and quantify the expected behavior of "noise", that is, of relations that are an artifact of the finite approximation. We will thus select relations whose behavior significantly differs from the expected behavior of noise.
The reconstruction of integer relations between real numbers is a well known application of the Lenstra-Lenstra-Lovász lattice basis reduction algorithm [42, p. 525], see also [9,11]. There are many other algorithms for the problem of computing integer relations, in particular the HJLS [33] and PSLQ [25,4,24] families and the first successful algorithm by Ferguson and Forcade [26]. A strong point in favor of the folklore LLL approach is that efficient LLL implementations are available in most computer algebra systems. To the best of our knowledge, existing work do not address the problem of computing the full lattice of integer relations (not just one) between real vectors (not just real numbers) given by approximations (not assuming exact data and exact arithmetic).
In this section, we recall and analyze the LLL approach to solve the following problem: Given a numerical approximation of a real matrix P ∈ R m×p , with p ≤ m, recover a basis of the lattice Λ = {x ∈ Z m | xP = 0}.
In our setting, the coefficients of P are the real and imaginary parts of the periods γ i ω j of the surface X under consideration. For B > 0 and ε > 0 let Λ B,ε be the lattice A rigorous numerical computation of Λ faces two obstacles: the lack of an a priori bound on the norm of generators and the inability to recognize zero among periods. In contrast, the lattice Λ B,ε can be computed exactly for given large B and small ε.
If B is larger than the length of the largest vector in a generating family of Λ, then for all ε > 0 small enough Λ B,ε = Λ. No a priori bounds are known about the values of B and ε that would ensure the desired equality Λ B,ε = Λ. If we choose B too small, we may miss integer relations in Λ. If we choose ε too big, we may compute relations that do not belong to Λ. Yet, meaningful results can be obtained by comparing with the expected situation.
Assume that, for some large β > 0 (typically 10 300 ), we are given the exact value of the m × p integer matrix P β obtained by entry wise rounding to the nearest integer the coefficients of βP , that is We complement the folklore LLL approach with the following heuristic. If β is large enough, Proposition 4.3 suggests that for ρ = rk Λ the norm b ρ is small but the norm b ρ+1 is large and comparable to β

Algorithm 4.1 Computation of the lattice of integer relations between approximate real vectors with a heuristic check.
Input. Q ∈ Z p×m and β > 0. Precondition. Q = βP + E for some P ∈ R p×m and E ∈ [− 1 2 , 1 2 ] p×m . Output.
function IntegerRelationLattice(Q, β) Compute an LLL-reduced basis b 1 , . . . , b m of the lattice spanned by the rows Q I m .
The growth of this function governs the ability to numerically reconstruct Λ.
As above, assume that, for some β > 0, we are given the exact value of the integer m × p matrix P β obtained by entry wise rounding to the nearest integer the coefficients of βP . Having · op denotes the operator norm, and where we used 1 2 √ pm ≤ 1 2 m ≤ m − 1, as m ≥ 2. Let R be the lattice generated by the rows of the integer m × (p + m) matrix M = P β I m and let b 1 , . . . , b m be an LLL-reduced basis of R. We denote B 0 = 0 and B i = b i , for 1 ≤ i ≤ m. In particular B 0 ≤ B 1 ≤ · · · ≤ B m . Gaps in this sequence typically separate the elements of R that come from genuine integer relations from spurious relations coming from the inaccuracy of the approximations. (ii) Λ is not generated by elements of norm ≤ κB i+1 ; (iii) ϕ(β) ≤ B i .
Proof. By Lemma 4.2 below, we have If Λ is generated by elements of norm ≤ κB i+1 then Λ ⊆ Λ κB i+1 ,mβ −1 κB i+1 , and therefore Λ ⊆ Proof. Let Λ i ⊂ Z m be the lattice generated by pr(b 1 ), . . . , pr(b i ) and let R i ⊂ R be the lattice generated by b 1 , . . . , b i . Let R| τ denote the sublattice of R generated by vectors of length at most τ . We Conversely, let x ∈ Z m such that x ≤ B and xP < mBβ −1 . Let r = xP β x . We check that The properties of an LLL-reduced basis [48,Thm. 9] imply that no family of i + 1 vectors in R with norms less than 2 −(m−1)/2 B i+1 is independent. Since b 1 , . . . , b i are independent and of norm ≤ B, it follows that r ∈ QR i . Moreover R is a primitive lattice (that is QR ∩ Z p+m = R) therefore any subset of the basis b 1 , . . . , b m of R generates a primitive lattice, so r ∈ R i . And therefore x = pr(r) ∈ Λ i .
Minkowski's Theorem on linear forms shows that if ε p β m−rk Λ−p ≥ det(P T P ), there is an x ∈ Λ ⊥ ∩ Z m such that xP ≤ pε. Therefore We define the irrationality measure of P , denoted µ(P ) as the infimum of all µ > 0 such that ε(β) = O(β 1−µ ) as β → ∞. As for the usual irrationality of real numbers, we can show with Borel-Cantelli Lemma that µ(P ) = m−rk Λ p for allmost all P ∈ R m×p with a given lattice Λ of integer relations. Generalizing Roth's Theorem on rational approximation of algebraic numbers, Schmidt [55] proved that if P has algebraic coefficients, with some additional hypotheses, then it again holds that µ(P ) = m−rk Λ p . All in all, this leads to Algorithm 4.1. The heuristic check relies on assuming µ(P ) = m−rk Λ p , approximating ϕ(β) β 1/µ(P ) , that is ϕ(β) β p m−rk Λ , and applying Proposition 4.3.

5.
Developing an intrinsic measure of error. Working with a finite approximation of periods, there is the possibility of miscomputing the Picard group. Although there will be no miscomputation if the periods are approximated to sufficient precision, we do not know a priori what constitutes "sufficient precision" for any given example. We can not solve this problem here but we will attempt to facilitate an a posteriori certification scheme.
More precisely, we define two intrinsic quantities B min and τ N related to an algebraic surface X. Then we give an algorithm that computes numerically the Picard group of X alongside with half-a-certificate (B, N, ε) ∈ R 3 >0 such that if B min < B and ε < τ N then the computation of the Picard group is correct. Unfortunately, the quantities B min and τ N are not easily computable.

Summary of the results in this section.
There is a canonical positive definite inner product on H 2 (X, R) = H 2 (X, Z) ⊗ Z R whose exact formulation we recall in Section 5.2. Briefly, it is obtained from the intersection product on cohomology by flipping the sign of the intersection matrix on a certain subspace of H 2 (X, R). We will refer to the norm induced by this inner product as the canonical norm.
Let us point out that the homology group H 2 (X, R) and the cohomology group H 2 (X, R) are canonically identified via Poincaré duality. If {γ i } m i=1 is a basis for homology and the cohomology is equipped with the corresponding dual basis, then the identification H 2 (X, R) ∼ → H 2 (X, R) is given by the intersection product matrix (γ i · γ j ) i,j . For the rest of this section, homology and cohomology will thus be identified.
Let P ∈ C r×m be the period matrix of X as in Section 2. Unlike the norm, the matrix P depends on the choice of a basis for homology and a basis for the space of holomorphic forms in H 2 (X, C). As outlined in Section 2, we can compute a finite approximation of P and then, via lattice reduction, compute a lattice Λ ⊂ H 2 (X, Z) Z m meant to represent the Picard group of X.
We give an explicit construction in Section 5.4 of a real number B associated to the computation of Λ. Let B min be the infimum over real numbers c such that Pic(X) is generated by elements of canonical norm at most c.
Let P R ∈ R 2r×m be the matrix obtained by joining the real and imaginary parts of P vertically. Define U ⊂ H 2 (X, R) R m to be the subspace generated by the rows of P R . With respect to the Hodge decomposition, one would write U = H 2 (X, R) ∩ H 2,0 (X) ⊕ H 0,2 (X) . Let U ⊥ be the orthogonal complement of U with respect to the cap product on cohomology.
In Section 5.4 we give an explicit construction for a pair of positive real numbers (N, ε) associated to our lattice Λ. Denote by τ N the minimum non-zero distance of vectors of canonical norm at most N in H 2 (X, Z) to U ⊥ , that is, Theorem 5.2. If ε < τ N then Λ ⊂ Pic(X).
The two theorems stated above are proven in Section 5.4.

Canonical norm.
Let X ⊂ P 3 be a smooth surface of degree d. There is a canonical positive definite intersection product on H 2 (X, R) [35, Example 3.1.7(ii)]. For the proofs of the statements made here, we refer to loc. cit.
If (γ i · γ j ) i,j is the intersection matrix on homology H 2 (X, R) with respect to a basis {γ i } m i=1 , then the inverse matrix (γ i · γ j ) −1 i,j represents the matrix for the cup product (γ * i ∪ γ * j ) i,j on cohomology H 2 (X, R) with respect to the dual basis {γ * i } m i=1 . Note that the cup product on cohomology is available to us.
By abuse of notation, we will denote by h X the image of the polarization with respect to the identification H 2 (X, R) ∼ → H 2 (X, R). Recall that U has been defined as the subspace of cohomology generated by the real and imaginary parts of a period matrix P of X. Let U ⊥ ⊂ H 2 (X, R) be the space orthogonal to U . If we denote by W ⊂ U ⊥ the subspace in U ⊥ orthogonal to h X then we get an orthogonal decomposition of cohomology with respect to the cap product: The cup product on cohomology is positive definite on U and R h X but negative definite on W . The canonical positive definite metric on H 2 (X, R) is obtained by flipping the sign of the cup product on W . In particular, for v ∈ H 2 (X, R) if we write v = v U + v W + v h by respecting the decomposition (5.2) then the canonical norm is given by:

It is easy to see that
Substituting these two terms into (5.3) we get the following expression: Let us emphasize that every term in this expression, with the exception of v U , is available to us. We now turn to the problem of approximating v U .
Algorithm 5.1 Computing a rectangular neighbourhood of the projection operator.

Input.
A rectangular neighbourhood of P R , that is, (Q R , E). Output. An operator P U which takes v ∈ H 2 (X, Z) and gives a rectangular neighbourhood containing v U .
1. Using ball arithmetic, orthonormalize the matrix of intervals with center Q R and radius E. Let B denote the resulting matrix of intervals. 2. The orthonormalization P on R of P R is contained in B. Therefore, for any v the operator P U defined by v → B · I · v is what we want.

Rigorous bounds for the canonical norm.
In order to use (5.4), we want to find the orthogonal projection operator P U : H 2 (X, R) → U ; v → v U . This operator can be constructed from the period matrix P of X. We note here that (5.4) varies continuously in v U , that is, small errors in expressing P U will miscompute v by a small constant of proportionality.
Let us denote by I ∈ Z m×m the matrix representing the cup product on cohomology, with respect to the trivialization H 2 (X, Z) Z m used to compute P. Recall that P R is the vertical join of the real and imaginary parts of P. As the cup product is positive definite on U , we can define P on R ∈ R 2r×m to be the matrix whose rows are obtained by the Gramm-Schmidt orthonormalization process from the rows of P R with respect to I. In coordinates, the projection operator P U : H 2 (X, R) → U is given by the matrix P on R · I. In practice, we only have a rational approximation Q R of P R and a matrix E such that the absolute value of the (i, j)-th entry of the difference P R − Q R is bounded from above by the (i, j)-th entry of E. In other words, we have a rectangular neighbourhood containing P R . This allows us to compute an approximation P U of the projection operator P U using Algorithm 5.1 such that P U (v) returns a rectangular neighbourhood v U of v U . Substituting v U for v U in (5.4) we get an interval which contains v .
In particular, for Theorem 5.2 we need the ability to bound the distance to U ⊥ . For any v, this distance is dist U ⊥ (v) = P U (v) . We can therefore bound dist U ⊥ (v) by computing v U and taking the maximum of the interval v U ∪ v U .

5.4.
Computing half of the certificate. Let H 2 (X, Z) be identified with Z m by the choice of a basis γ * 1 , . . . , γ * m in which the period matrix P is computed.
Recall that our approximation Q R of P R comes with rigorous error bounds which we can collect into the matrix E to form the rectangular neighbourhood (Q R , E).
Let Λ ⊂ Z m be the candidate Picard group of X computed using the approximation Q R of P R as in Section 2. The lattice Λ is computed according to an exact procedure detailed in Section 4. In order to claim that Λ contains Pic(X) we need a careful study the lattice reduction algorithm by which Λ is produced. This study takes place in Section 4 and that section allows us to compute Λ together with a real number B Λ . This number B Λ has the property that a vector v ∈ Pic(X) ⊂ Z m is contained in Λ if the coordinate norm |v| is less than B Λ . We will now unshackle this bound from the coordinates chosen and phrase it in terms of the canonical norm. Lemma 5.3. For any v ∈ H 2 (X, Z) the canonical norm is bounded by the coordinate norm as follows: |v| Proof. Use the triangle inequality and Cauchy-Schwartz on · for the upper bound. The lower bound follows from the observation |v i | γ * i ≤ v , as · is induced from an inner product.
In particular, upper and lower bounds of the canonical norms of γ * i can be computed as in Section 5.3. Using these bounds and Lemma 5.3 we can compute ξ 1 , ξ 2 ∈ R >0 such that for every v ∈ H 2 (X, R) we have ξ 1 |v| ≤ v ≤ ξ 2 |v|. We now define B to be ξ 2 B Λ .
Proof of Theorem 5.1. By construction of B Λ , if v ∈ Pic(X) satisfies |v| < B Λ then v ∈ Λ. Using Lemma 5.3 and by definition of B we have v ≤ ξ 2 |v| < ξ 2 B λ = B. By hypothesis, Pic(X) is generated by its elements of canonical norm at most B min < B. Therefore, a generating set for Pic(X) is contained in Λ.
Let v 1 , . . . , v ρ ∈ Λ be a basis. Using Section 5.3 compute for each i the interval v i containing the canonical norm v i and compute the interval containing the distance dist U ⊥ (v i ). Using these intervals we define N, ε ∈ R >0 such that v i < N and dist U ⊥ (v i ) < ε.
Proof of Theorem 5.2. An element v ∈ H 2 (X, Z) is in Pic(X) if and only if dist U ⊥ = 0 by Lefschetz (1,1) theorem. If any v of H 2 (X, Z) having canonical norm at most N is either in Pic(X) or satisfies dist U ⊥ (v) > ε then, by definition of N and ε, Λ ⊂ Pic(X). The hypothesis ε < τ N ensures precisely this condition.
6. Hypersurfaces of arbitrary even dimension. Let k be a positive integer and let X ⊂ P 2k+1 be a smooth hypersurface. Using Lefschetz hyperplane theorem and Poincaré duality we see that the cohomology groups H i (X, Z) are either trivial or Z except when i = 2k. The Hodge decomposition on de Rham cohomology gives H p,q (X, C).
Algebraic cycles of dimension k in X give cohomology classes in As a generalization of Theorem 2.1, Hodge conjecture predicts that the vector space Hdg k (X)⊗ Z Q is spanned by algebraic cycles [65]. The Hodge group Hdg k (X) comes with an intersection pairing obtained by restricting the cup product on cohomology H 2k (X, C). Furthermore, there is a polarization h k X ∈ Hdg k (X) where h X is the class of a generic hyperplane section of X. The tools we used to tackle the computation of Picard groups apply to the following following problem: Given the defining equation of X ⊂ P 2k+1 , compute the rank ρ of Hdg k (X), the ρ × ρ matrix of the intersection product and the ρ coordinates of the polarization h k X in some basis of Hdg k (X) Z ρ .
Suppose now that γ 1 , . . . , γ m ∈ H 2k (X, Z) is a basis for the middle homology group of X. We can then identify the cohomology H 2k (X, C) = Hom(H 2k (X, Z), C) with C m via the dual basis of {γ i } m i=1 . Let us write F 2k, (X, C) = j=0 H 2k−j,j (X, C) for the corresponding Hodge filtration.
Let ω 1 , . . . , ω s ∈ F 2k,k−1 (X, C) be a basis for the (k − 1)-th part of the Hodge filtration. Suppose that we have the coordinates of ω i with respect to the identification H n (X, C) C m , that is, suppose that for each i = 1, . . . , s the following integrals are known These periods of X are listed as the columns of the following matrix: The matrix P induces the linear map P Z : Z m → C s by acting on the integral vectors from the right.
Lemma 6.1. We have a natural isomorphism Hdg k (X) ker P Z .
Proof. The kernel of P Z computes in H 2k (X, Z) the classes annihilated by F 2k,k−1 (X, C). Any integral (or real) class annihilated by F 2k,k−1 will also be annihilated by its complex conjugate F 2k,k−1 . We now use the equality H k,k (X, C) = F 2k,k−1 ⊕ F 2k,k−1 ⊥ and the definition of Poincaré duality.
The kernel of P Z sits most naturally in homology H 2k (X, Z) and is denoted by Hdg k (X). Using Sertöz [57] we can approximate the matrix P to the desired degree of accuracy for an automatically generated basis of F 2k,k−1 (X, C) and some implicit basis of H 2k (X, Z). The basis of H 2k (X, Z) comes with an intersection pairing as well as the coordinates of the polarization h k X in this basis.
In light of Lemma 6.1, we need to compute integral linear relations between the columns of P, that we compute numerically and we recover Hdg k (X) with lattice reduction algorithms, see §4.
To be more precise, identify H 2k (X, Z) with Z m by choosing a basis. In turn, Hdg k (X) is identified with a sublattice Λ ⊂ Z m . What we can compute is the following sublattice for some B, ε ∈ R >0 : (6.5) Λ B,ε = {a ∈ Z m | P Z (a) < ε, a < B}.
Example 6.2. With minimal effort, the period computations explained in [57] can be extended to mildly singular hypersurfaces in P 2k+1 for k > 1. We computed the periods of the Delsarte surface cut out by the quintic polynomial (6.6) We see that the Picard number of this surface is 13, in agreement with [56, §6].
The study of the Hodge groups of cubic fourfolds is an active area of research [54,2]. Although generic cubic fourfolds provide a computational challenge, we can quickly compute the Hodge rank of sparse cubic fourfolds if most of the monomial terms are cubes of a single variable.
7. Experimental results. We documented the Picard groups of 2790 quartic surfaces defined by sparse polynomials. For these computations, setting up the initial value problem for the periods was not the limiting factor but rather the numerical solution of these initial value problems took the greatest amount of time. With our current methods, the periods of a quartic surface defined by dense polynomials can not be computed in a reasonable amount of time, say, less than a day.
The database presented here relies on a systematic exploration of quartics that are defined by a sum of at most six monomials in x, y, z, w with coefficients 0 or 1. We built a graph whose vertices store the defining polynomials and an edge between two polynomial is constructed if the difference of the two polynomials is supported on at most two monomials (this is done to ease the computations). Then, for each edge, we setup and attempt to solve the initial value problem defining the transition matrix from one set of periods onto the other, using 300 decimal digits of precision, see [57] for details. Computation is stopped if it takes longer than an hour and the edge deleted. Having explicit formulas for the periods of Fermat surface x 4 + y 4 + z 4 + w 4 = 0 , we can compute the periods of any vertex in the connected component of x 4 + y 4 + z 4 + w 4 in the resulting graph by simply multiplying the transition matrices of each edge along a path.
For each of the 2790 polynomials in our database, we computed the Picard group, the polarization, the intersection product, the endomorphism ring and the number of smooth  Example 7.2. {x 3 y + x 3 z + y 3 z + yz 3 + z 4 + xw 3 = 0}. This surface has Picard number 10. It contains 13 lines that generate the Picard group. The endomorphism ring is Q(exp( 2πi 18 )), a cyclotomic extension of Q of degree 6. Up to degree 10, the smooth rational curves inside X count as follows. Example 7.4. {x 3 y + z 4 + y 3 w + zw 3 = 0}. This surface has Picard number 4. It contains exactly 4 lines that generate the Picard group and no other smooth rational curve of degree < 100. The endomorphism ring is Q(exp( 2πi 54 )), a cyclotomic extension of Q of degree 18.
8. Computing the coordinates of the polarization. Let X = Z(f X ) ⊂ P n+1 be a smooth hypersurface of degree d and assume n is even. We compute a basis for the middle integral homology H n (X, Z) by carrying over a basis from a hypersurface of Fermat type [57, §1.3]. If h X = [X ∩ H] denotes the hyperplane class in X, then h n/2 X ∈ H n (X, Z) is the polarization. The orthogonal complement of h n/2 X is the primitive homology, denoted P H n (X, Z). In order to compute the periods of X, it is sufficient to work only with the primitive homology P H n (X, Z) as is done in [57]. In §4.5 of loc. cit. there is a sketch on how to complete the given basis for the primitive homology to a basis of homology. In this section we flesh out the details as the particular choices we make in completing the basis determine the coordinates of the polarization. The problem that must be addressed is that h n/2 X and P H n (X, Z) do not generate H n (X, Z) but a full rank sublattice.
In [57] the Fermat surface Y = Z(x d 0 + · · · + x d n − x d n+1 ) was used for the construction of a basis of primitive homology. This basis is formed using a Pham cycle and the Pham cycle itself is formed by gluing translates of the following simplex: D = {[s 0 : s 1 : · · · : s n : 1] | s i ∈ [0, 1], s d 0 + · · · + s d n = 1} ⊂ Y.
For β = (β 0 , . . . , β n+1 ) ∈ Z n+2 we define the translations t β : P n+1 → P n+1 by the action on the coordinates x i → ξ β i x i . Then the Pham cycle S is defined by: where summation is union and negation is change of orientation [51]. It is possible to compute a subset B ⊂ Z n+2 for which the set {t β S} β∈B is a basis for the primitive homology P H n (Y, Z), for instance, use Corollary 4.8 [57]. We will now add one more cycle to complete {t β S} β∈B to a basis of homology.
With d being the degree of X and Y , we denote the d-th root of −1 by η := exp( π √ −1 d ) and the d-th root of 1 by ξ := exp( 2π √ −1 d ). Let P n/2 be a projective space with coordinate functions µ 0 , . . . , µ n/2 and consider the linear map P n/2 → P n+1 defined by x 2k = µ k , x 2k+1 = ηµ k k = 0, . . . , n 2 − 1, (8.1) x n = µ n 2 , x n+1 = µ n 2 . The image of this map is a linear space L which is evidently contained in Y . Let [L] be the homology class of L and let γ β be the homology class of t β S. The set {[L]} ∪ {γ β } β∈B is a basis for the integral homology H n (X, Z).
As Y is deformed into X, the homology class of L will typically deform into a class which no longer supports an algebraic subvariety and therefore this class will typically have non-zero periods. Nevertheless, we can deduce the periods of L as it deforms based on the following two observations: The polarization h n/2 Y deforms in to h n/2 X and will always remain algebraic throughout the deformation. We will know the periods of the Pham basis {t β S} β∈B as it deforms.
The homology with rational coefficients H n (Y, Q) splits into the direct sum P H n (Y, Q) ⊕ Q h n/2 Y so that we may write: i=0 . The Pham cycles t β S and t β S intersect as follows: For a proof of Proposition 8.1 see any one of [3,46,43]. We reformulated the statement here for the choices that were made in [57] and in the style that was first communicated to us by Degtyarev  With respect to this basis, and the ordering presented above, we find that the vector a B,L of (8.3) is given by