Paolo Cattaneo
email: paolo.cattaneo@cea.fr

Roland Lenain

Elsa Merle

Cyril Patricot

Didier Schneider

Numerical optimization of a multiphysics calculation scheme

Keywords: MULTIPHYSICS, PARTIAL-CONVERGENCE, APOLLO3 ®, FLICA4

This work concerns the numerical optimization of a multiphysics calculation scheme. The considered application is a 5x5 Pressurized Water Reactor (PWR) assemblies mini-core surrounded by radial and axial reflectors. The scenario adopted for the analysis is steady-state nominal conditions and fission products set to the equilibrium concentration. The neutronics is modelled at the pin-cell scale and the thermal-hydraulics at the subchannel level. Depending on the scenario, the damped fixed-point algorithm might not be sufficiently robust or efficient enough. For this reason, a technique based on the partial convergence of the solvers is tested. In every multiphysic iteration, a maximum number of iterations is imposed for both the neutronics and the thermal-hydraulics solvers. In combination with that, the solver restarts from the results of the last calculation. In this way, if the method is convergent, the initialization progresses towards the fixed-point solution. The results show that such a technique improves both the robustness and the speed of the algorithm. Within this approach, the range of relaxation factors that makes the algorithm converge is significantly broadened and the importance of this parameter on the global performance is reduced. The computing time also decreases by a factor between 10 and 20. Furthermore, this gain does not strongly depend on the exact imposed maximum number of iterations. Some preliminary observations are also reported in respect with the application of such a technique to the Anderson acceleration method.

INTRODUCTION

The aim of this work is to improve the simulations of pressurized water reactor during start-up and normal operation. Nuclear reactor modelling is generally a coupled multiphysics problem, in the sense that merely solving the neutronics problem would lead to satisfying results only in a limited set of scenarios, mostly at zero power. In fact, the heat generation affects the temperature and density 3D distributions and in turn, P. Cattaneo et al., Numerical optimization of a multiphysics calculation scheme Proceedings of the PHYSOR 2020, Cambridge, United Kingdom those fields influence neutronic macroscopic cross-sections, mainly through Doppler and moderator effects. The considered calculation scheme aims to predict variables at a fine scale: the neutronics is discretized at the pin-cell level, the thermal-hydraulics is of the subchannel type and 1D radial heat conduction is modelled in every fuel rod. Such a refinement implies the use of complex models and contributes to the increase of the problem's dimension. Therefore, algorithm optimization becomes even more crucial in order to obtain a sufficiently low computing time.

The development of a dedicated multiphysics solver able to simulate all the involved physics might be prohibitive, so the common tendency is to use operators splitting: reproducing the entire model through the combination of already existing solvers dedicated to the resolution of a subset of the physics. Under this approach, if the solvers are not customized, normally only a limited amount of information is shared among the solvers, not every variable might be easy to access. Hence, not every algorithm can be implemented, resulting in a potential limit in the range of possible applications. Furthermore, naïve implementations of this method might eventually cause problems of stability and robustness. Another point to be considered is that, for the considered applications, data manipulation in the supervisor takes a small amount of time as compared to the full resolution of the neutronics or the thermal-hydraulics. Therefore, it might be very effective to implement advanced numerical methods, in order to exploit, to the greatest extent, the available information; hence to optimize the numerical scheme. For a more general and extended explanation on this type of problem it is possible to refer to [START_REF] Keyes | Multiphysics simulations: Challenges and opportunities[END_REF]. Currently, in reactor physics, many works focus on the numerical optimization in multiphysics calculations. The most recurrent numerical algorithms are damped fixed-point, Anderson acceleration and Jacobian-Free Newton-Krylov (JFNK) [START_REF] Hamilton | An assessment of coupling algorithms for nuclear reactor core physics simulations[END_REF][START_REF] Zhang | The comparison between nonlinear and linear preconditioning JFNK method for transient neutronics/thermal-hydraulics coupling problem[END_REF][START_REF] Toth | Analysis of Anderson acceleration on a simplified neutronics/thermal hydraulics system[END_REF]. Here the attention is concentrated on the partial convergence strategy for the first method. Preliminary observations on the application of such an approach to the Anderson methods are also drawn.

CASE STUDY

The considered scenario is a steady state, nominal conditions and fission products at the equilibrium concentrations. The chosen domain is a mini-core made of 5x5 PWR assemblies surrounded by a radial, an upper and a lower reflector. The fuel assemblies are loaded with 4% enriched urania and the 2D burnup distribution follows a chessboard-like configuration with values 0, 15 or 30 MWd/kg (see Fig. 1). The average linear power is set to 160 W/cm and the mass flux is 3900 kg/m 2 /s. The size of the core is determined by the trade-off between computing time and representativeness of the simulation. The current choice is suitable for systematic studies and it implies a power peaking factor slightly bigger than the one of standard PWRs. Furthermore, such a form factor implies a higher maximum power, hence qualitatively a stronger local coupling, which might be desirable for numerical tests.

THEORETICAL BACKGROUND

This section formalizes the problem into equations and introduces the numerical challenges associated to it. Particular focus is set on some general observations on the fixed-point's convergence. For a more detailed

Physics of Reactors Transition to a Scalable Nuclear Future

Proceedings of the PHYSOR 2020, Cambridge, United Kingdom analysis on the fixed-point and Anderson's methods, it is possible to refer to [START_REF] Quarteroni | Numerical mathematics[END_REF][START_REF] Patricot | Couplages multi-physiques : évaluation des impacts méthodologiques lors de simulations de couplages neutronique/thermique/mécanique[END_REF]. The details on the modelling and the practical applications are left for the following sections.

Problem formalization

A general formulation of the problem is given in the system of Eq. 1.

{ (𝝆 𝑤 , 𝑻 𝑓 ) = 𝑇𝐻(𝒒 𝑓 , 𝒒 𝑤 ) (𝝓, 𝒒 𝑓 , 𝒒 𝑤 ) = 𝑁(𝝆 𝑤 , 𝑻 𝑓 , 𝑪 𝐹𝑃 )

(𝑪 𝐹𝑃 ) = 𝐹𝑃(𝝓) (1)
Where, the TH operator is the thermal-hydraulic one, it determines the water density (𝝆 𝑤 ) and the fuel temperature (𝑻 𝑓 ) fields for a given power distribution. N represents the neutronics operator, its main goal is to compute the scalar neutron flux (𝝓), from which are derived the power fields in the fuel (𝒒 𝑓 ) and in the water (𝒒 𝑤 ). The neutron cross-sections depend on the temperature and density fields and the isotopic concentration of fission products (CFP). FP corresponds to the operator that finds the equilibrium concentration of the fission products for a given neutron flux. In terms of equations, N solves the keigenvalue problem associated to the Boltzmann equation for the neutron flux and associates a power distribution to the latter. TH solves for both the convection of heat in the water and the heat conduction in the fuel. In particular, this operator finds the steady-state solution as asymptotic result of a transient problem.

FP calculates the equilibrium concentration of fission products following the Bateman equations. The global problem is a nonlinear system of equations.

Observations on the Fixed-Point Method

The fixed-point is one of the simplest and most diffused numerical method for the resolution of large nonlinear problems. However, many variants might be implemented, following different strategies. Its versatility derives from the fact that the algorithm only requires one evaluation of the function per iteration. The fixed-point formalism of the problem is defined as it follows. Let 𝐆 ∶ 𝐷 ⊂ ℝ 𝑛 → ℝ 𝑛 be the global function of the system of Eq. 1 and 𝐱 * its fixed-point solution, so that 𝐆(𝐱 * ) = 𝐱 * . 𝐱 is the vector of all the variables (power density, water density, fuel temperature, etc.). The fixed-point sequence is given in Eq. 2.

𝐱 𝑘+1 = 𝐆(𝐱 𝑘 ) (2)
Noting the Jacobian matrix of 𝐆 as 𝐉 𝐆 = 𝝏𝐆 𝝏𝐱 and ρ(𝐉 𝐆 (𝐱 * )) its spectral radius (i.e. the maximum in module of the eigenvalues). If 𝐆 has a fixed-point in 𝐷 , 𝐆 ∈ C 1 in its proximity and ρ(𝐉 𝐆 (𝐱 * )) < 1 , then there exists a neighborhood of 𝐱 * in which, for every initial guess belonging to the latter, the fixed-point sequence is entirely contained and converges to the solution. In practical applications, ρ(𝐉 𝐆 (𝐱 * )) is never calculated since it requires to know the solution. Unless all the derivatives are equal to zero in 𝐱 * , the fixed point converges linearly and rate of convergence is equal to the spectral radius.

A simple and effective way to accelerate the convergence of the fixed-point is to introduce a relaxation factor, which can be seen as a preconditioning method, forming the sequence in Eq. 3. The relaxed function 𝐆 𝐑 has eigenvalues linearly related to the ones of 𝐆, as specified in Eq. 4. Through relaxation, from the non-converging sequence associated to 𝐆, it is possible to build a converging one with 𝐆 𝐑 : if max (λ i ) < 1 ∃ α > 0 ∶ |λ R,i | < 1 ∀ i, with α being the relaxation factor and λ i the i-eigenvalue of 𝐆 and λ R,i the i-eigenvalue of 𝐆 𝐑 . The optimum convergence speed is obtained when the eigenvalues tend to zero. Since the relaxation factor affects every eigenvalue, the aim is to minimize the spectral radius. The more the maximum and minimum eigenvalues are distant, the less there is room for acceleration. An optimum relaxation can be found exploiting the fact that the relaxed minimum and maximum eigenvalues are equal in module when the relaxation factor is the best. The relation is presented in Eq. 5.

𝐱 𝑘+1 = 𝐆 𝐑 (𝐱 𝑘 ) = α * 𝐆(𝐱 𝑘 ) + (1 -α) * 𝐱 𝑘 (3) λ R,i (α) = α * (λ i -1) + 1 (4)
α opt = 1 1- λ max +𝜆 𝑚𝑖𝑛 2
(5)

METHODOLOGY

The multiphysics tool used in this work is CORPUS [START_REF] Le Pallec | Neutronics/Fuel Thermomecanics coupling in the Framework of a REA (Rod Ejection Accident) Transient Scenario Calculation[END_REF], which is derived from the SALOME platform [START_REF] Ribes | Salomé platform component model for numerical simulation[END_REF].

CORPUS has the role of gathering several independent codes and providing the user with tools to help data exchanges among them. The here considered neutronics solver is APOLLO3 ® [START_REF] Schneider | APOLLO3 ® : CEA/DEN deterministic multi-purpose code for reactor physics analysis[END_REF], which is capable of performing both lattice and core calculations with a wide set of models. In particular, in this work it is adopted a two-steps scheme: lattice calculation with the Method Of Characteristics [START_REF] Santandrea | Accelerated polynomial axial expansions for full 3D neutron transport MOC in the APOLLO3 ® code system as applied to the ASTRID fast breeder reactor[END_REF], SuPer-Homogenization equivalence technique [START_REF] Kavenoky | The SPH homogenization method[END_REF] and Simplified Spherical Harmonics [START_REF] Baudron | MINOS : A simplified Pn solver for core calculation[END_REF]. The homogenization is at the pin-cell scale and the cross-sections are parametrized in fuel Doppler temperature, moderator density, boron concentration and burnup. The core model uses Simplified Spherical Harmonics of order three and eight energy groups, this choice resulted from a comparative analysis available at [START_REF] Cattaneo | Development of a multiphysics Best-Estimate approach for LWR reference calculation[END_REF].

The depletion calculations are performed by MENDEL-solver library [START_REF] Lahaye | First verification and validation steps of MENDEL release 1.0 cycle code system[END_REF]. Subchannel thermal-hydraulics (4-equations modelling in porous media) and radial 1D heat conduction in every fuel rod are simulated by FLICA4 [START_REF] Toumi | FLICA4: a three dimensional two-phase flow computer code with advanced numerical methods for nuclear applications[END_REF]. The scale of the three radial meshes is given in Fig. 2. Before entering in the details of the algorithms, it is important to observe that the power iteration method adopted in the neutronics operator naturally allows the possibility to impose a maximum number of iterations, in order to obtain partial convergence. The general effect is similar to an under-relaxation (α < 1), but it may bring to significant savings in computational cost. On the other hand, such an expedient should be used carefully, because the output of the operator then depends on the solver's initialization and on the maximum iteration number, unless convergence is fully met. One strategy consists in keeping the conditions at the end of the last calculation as initialization for the following one. In this case, for stable algorithms, the initial conditions progressively approach the final result and a significant number of iterations can be saved. The drawback of such an approach is that the operator is only consistent, not anymore strongly consistent. In other words, the real operator 𝑁 is substituted by an approximated 𝑁 𝑛 (where n is the multiphysics iteration number) which is equal to the real one only at convergence. The same

Physics of Reactors Transition to a Scalable Nuclear Future

Proceedings of the PHYSOR 2020, Cambridge, United Kingdom reasoning applies to imposing a limit on the number of the time steps for the thermal-hydraulics. The fixedpoint is directly implemented in the calculation scheme, while Anderson is simply imported from the opensource Scipy-Optimize package [START_REF] Jones | SciPy: Open Source Scientific Tools for Python[END_REF]. The only informatic development required was the design of a convertor that changes the format from dictionary of MED fields (SALOME standards) to one Numpyarray (Python standard). The convergence criteria have been determined for the L2-norm of the absolute residuals. The convergence of every variable is required as well as the internal convergence of every solver.

Fixed-Point with Relaxation Method

The relaxation factor is applied on the neutronics operator only. The method is presented in Algorithm 1 and in the flowchart of Fig. 3. The algorithm starts assuming a flat power distribution. At every iteration, the solver is initialized with the results of the last calculation.

Residual formalism

For Anderson one unique function of the type 𝐅(𝐱 * ) = 𝟎 is written, it is available in Eq. 6. The algorithms are not reported here since they are fully open-source [START_REF] Jones | SciPy: Open Source Scientific Tools for Python[END_REF] and about their numerical formulations it is possible to refer to [START_REF] Quarteroni | Numerical mathematics[END_REF][START_REF] Patricot | Couplages multi-physiques : évaluation des impacts méthodologiques lors de simulations de couplages neutronique/thermique/mécanique[END_REF]. In this case, all the variables are normalized by their convergence criterion, in order to give an equal weight to the minimization of each residual.

𝐅(𝐱

FIXED-POINT EXPLORATION

The results deal with a detailed analysis of the stability and efficiency of the fixed-point method for different maximum numbers of iterations within a single operator call. The conditions of all the calculations here presented are always the same, the one described in Section 2 (Case Study). An estimator of the performance of the algorithm is defined according to Eq. 7. I eqv is a computing time estimator, here referred to as number of equivalent iterations, I N is the number of power iterations, I TH of thermal-hydraulic time steps and I FP of calls to the FP operator and they are weighted on their respective computing time. In the following results, I eqv is normalized, so it is expressed as number of Relative Equivalent Iterations.

I eqv = I N * t 1N + I TH * t 1TH + I FP * t FP (7)
The first calculation executed is simple fixed-point iterations with no relaxation and no limitations on the maximum number of single solver iterations. This simple algorithm does not converge for the considered conditions. The power peak oscillates axially and diverges. Given the diverging oscillations in the module of the residual, it is expected that the spectral radius is larger than the unity, specifically the minimum of the eigenvalues smaller than minus one. For this reason, different under-relaxation factors are tested, the results are available in Fig. 4 on the left (Inf.-Inf. curve). The number of equivalent iterations is normalized to the best result in order to easily read the difference in computing time. As expected, there is a range of relaxation factors that makes the scheme converge, in this case it is within zero and 0.4, endpoints excluded. Nevertheless, for some values smaller than 0.1, the convergence is obtained, but the computational cost becomes prohibitive. For α ≥0.4 the calculation does not converge. Considering these results, the fixedpoint might be considered as not suitable for such a strongly coupled problem due to the very narrow range of acceptable relaxation factors. It is therefore tried to limit to different extents the number of iterations of the neutronics and thermal-hydraulics operator per call. The global results appear in Fig. 4 on the left; in the legend it appears the number of maximum iterations per solver under the format 𝑁 𝑁 -𝑁 𝑇𝐻 , where 𝑁 𝑁 and 𝑁 𝑇𝐻 are respectively the limits on neutronic and thermal-hydraulic iterations. It is found that limiting the number of iterations brings two main advantages: significant saving on the total number of iterations and a decreasing sensibility to the relaxation factor. Moreover, even if the α opt may vary significantly among the different operator settings, the performance does not seem to be very sensitive to the exact number of maximum iterations per solver. It is important to notice that I eqv does not account for the data manipulation time. So if this time becomes important compared to the solver times, small limits on the iteration numbers would become much less competitive. Another test that has been done is to asymmetrically reduce the number of iterations for the two solvers, the results are available in Fig. 4 on the right. Limiting the neutronic iterations seems to be more stabilizing than the opposite. The 100-40 case shows a more abrupt increase of sensitivity for relaxation factor higher than the optimum than the 20-150.

Physics of Reactors Transition to a Scalable Nuclear Future

Proceedings of the PHYSOR 2020, Cambridge, United Kingdom

PRELIMINARY RESULTS FOR ANDERSON

The Anderson algorithm may be seen as a technique to accelerate the fixed-point [START_REF] Anderson | Iterative procedures for nonlinear integral equations[END_REF]. Briefly, it consists in combining the M last iterations to find the new one that minimizes the residual. It is important to notice that Scipy-Anderson offers more than the simple Anderson algorithm. It provides the possibility to perform Armijo line search and to introduce a regularization parameter [START_REF] Eyert | A comparative study on methods for convergence acceleration of iterative vector sequences[END_REF]. Firstly, the robustness of Anderson's method is assessed by running the same case without single solver iteration limits, without relaxation, but with Armijo line search and a regularization parameter of 0.01. The algorithm converges for any M larger than zero. The residuals are rescaled so that they are equal to one when convergence is met. The results are available in Fig. 5 on the left, only the convergence of 𝒒 𝑤 (power in the water) is reported, the other variables follow a similar trend. Afterwards, it is tested to impose different single solver iteration limits, with constant M = 5. The results are available in Fig. 5 on the right. The convergence may still be reached, but limiting the number of iterations has a negative effect on the stability of the method and it even leads to a non-convergence in the case 20-40. This can be intuitively explained with the solver initialization strategy that may imply inconsistency within groups of iterations, especially among the firsts. Reducing M alleviates the problem, but to completely solve it, a different initialization strategy should be considered. Similar problems are expected to arise for any quasi-Newton method, in particular for the JFNK.

CONCLUSIONS

A calculation scheme allowing to perform steady-state calculation of the coupled problem of neutronics, thermal-hydraulics, thermal-mechanics and fission products evolution is implemented in CORPUS. Different neutronics and thermal-hydraulics models can be used. The present work focuses on pin-cell and subchannel simulation. A 5x5 PWR assemblies mini-core is modelled in nominal conditions and fission products at equilibrium concentration. Already at nominal power, the coupling has shown to be too strong for a simple fixed-point algorithm. Adding the relaxation method is sufficient to make the scheme converge. However, the range of acceptable relaxation factors is narrow: 0.1< α <0.38. A variant of the fixed-point algorithm is tested combining two actions. The first is to impose a limit on the number of power iterations for the neutronics operator and on the number of time steps for the thermal-hydraulics operator. The second is to keep the solver status at the last iteration as initialization for the following one. In this way, in the case of converging iterations, the solver initialization approaches more and more the final solution. The computing time is significantly reduced, up to twenty times faster and the range of acceptable relaxation factors is broadened to 0.1 < α < 1.1. Moreover, the computing time becomes less sensitive to the relaxation. The Anderson algorithm has proven to converge for any tested M in the interval between one and twenty, even with no relaxation and no limits on the iterations on the solvers. The strategy of imposing limits on the number of iterations of the solvers and letting the solver initialization progress towards the final solution is negatively affecting the stability of the method. The cause may be found in the corresponding inconsistency between the operators at the M last iterations. Resolving such a problem with a different initialization strategy would lead to promising results.

Figure 1 :

 1 Figure 1: Burnup map at beginning of cycle, values relative to the fuel assembly in MWd/kg.

Figure 2 :

 2 Figure 2: Radial mesh refinement (north-east quarter section of the core, reflector included for 𝒒 𝑓 ): pin-cell neutronics, subchannel thermal-hydraulics and 1D heat conduction for every rod.

Algorithm 1 :Figure 3 :

 13 Figure 3: Flowchart of two consecutive iterations of Algorithm 1.

Figure 4 :

 4 Figure 4: Effect of the damping on the convergence of the algorithm for different limitations on the iteration numbers per call. The labels follow the definition given in the current paragraph. DIV stands for divergent behavior, NC for non-convergent trend, SC for slow convergence.

Figure 5 :

 5 Figure 5: Convergence of 𝒒 𝑤 . On the left, it is reported for the Anderson algorithm for several values of the 𝐌 parameter and for the fixed-point iterations (fp), both with no limits on the iteration numbers. On the right, same algorithm with 𝐌=5 and different limits on the number of iterations.

Proceedings of the PHYSOR 2020, Cambridge, United Kingdom


Pierre Lairez

Emre Can Sertöz

A numerical transcendental method in algebraic geometry: Computation of Picard groups and related invariants

Based on high precision computation of periods and lattice reduction techniques, we compute the Picard group of smooth surfaces in P 3 . 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. 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 (K3 surfaces) defined by sparse polynomials. The results are only supported by strong numerical evidence; yet, the possibility of error is quantified in intrinsic terms, like the degree of curves generating the Picard group.

1. Introduction. Griffiths emphasized the role of certain multivariate integrals, known as periods, "to construct a continuous invariant of arbitrary smooth projective varieties" [START_REF] Griffiths | A transcendental method in algebraic geometry[END_REF]. Periods often determine the projective variety completely and therefore its algebraic invariants. Translating periods into discrete algebraic invariants raises difficult questions, 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 [START_REF] Sertöz | Computing periods of hypersurfaces[END_REF] and put transcendental methods into practice. This opens new possibilities in numerical algebraic geometry. While path tracking methods from numerical algebraic geometry are successful at computing topological invariants-irreducible decomposition [START_REF] Sommese | Numerical Decomposition of the Solution Sets of Polynomial Systems into Irreducible Components[END_REF], monodromy groups [START_REF] Hauenstein | Numerical Computation of Galois Groups[END_REF] or Chern numbers [START_REF] Di Rocco | Chern numbers of smooth varieties via homotopy continuation and intersection theory[END_REF]-the computation of periods reveal finer algebraic invariants. For example, they are used to compute the endomorphism ring of the Jacobian variety of algebraic curves [START_REF] Van Wamelen | Examples of genus two CM curves defined over the rationals[END_REF][START_REF] Booker | A database of genus-2 curves over the rational numbers[END_REF][START_REF] Costa | Rigorous computation of the endomorphism ring of a Jacobian[END_REF][START_REF] Bruin | Numerical computation of endomorphism rings of Jacobians[END_REF]. We extend this transcendental method to higher dimensions and focus mainly on algebraic surfaces. We describe and implement a numerical method to compute Picard groups of smooth surfaces. As an application, we count smooth rational curves on quartic surfaces using the Picard group. We demonstrate that these methods apply more generally to smooth hypersurfaces in a projective space of arbitrary dimension. Throughout the paper we work over the complex field but the implementation is for hypersurfaces defined over the rationals.

Structure of the Picard group and main results. There are many curves in a smooth surface X ⊂ P 3 . The basic ones are those obtained by intersecting X with another surface S in P 3 . 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 [40, p. 143] 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" [83, 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 . The polarization is the linear equivalence class of any curve obtained by intersecting X with a plane in P 3 . 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, the 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 [START_REF] Sertöz | Computing periods of hypersurfaces[END_REF] 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's gap [START_REF] Charles | On the Picard number of K3 surfaces over number fields[END_REF], 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. We study rigorously the reliability of using approximate periods to make deductions about integral relations between them. In addition, we built a database containing 180, 000 quartic surfaces to demonstrate the feasibility of the approach.

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 [START_REF] Deligne | The Hodge conjecture[END_REF]. The resolution of this conjecture is one of the seven Millennium Prize Problems posed by the Clay Institute. 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 §5.

Applications. Let X be the space of solutions to a polynomial system of equations, as might occur in integer programming, algebraic statistics or geometric modeling. It is natural to consider subspaces of solutions that can be freely parametrized algebraically, such as linear spaces inside X. When X ⊂ P 3 is a smooth quartic surface, there is no dominant rational map P 2 -→ • X but there are often non-constant maps P1 → X. One of the applications of having the Picard group of a quartic surface is the immediate ability to count smooth rational curves, i.e. embeddings P 1 → X, of any degree, see §3. Symbolic methods can tackle this problem too but to the best of our knowledge, the computational complexity is too high to count anything else that curves of degree one or two. Attempts at counting the number of maps P 1 → X led to deep connections between mathematics, string theory and mirror symmetry [START_REF] Hori | Mirror symmetry[END_REF].

Periods of curves play a crucial role in finding algebro-geometric solutions to soliton equations, e.g. integrable models in the KdV hierarchy. There is a large amount of work, theoretical and applied, building upon this connection [START_REF] Shiota | Characterization of Jacobian varieties in terms of soliton equations[END_REF][START_REF] Gesztesy | Soliton equations and their algebro-geometric solutions[END_REF][START_REF] Deconinck | Computing with Plane Algebraic Curves and Riemann Surfaces: The Algorithms of the Maple Package "Algcurves[END_REF]. It is conceivable that periods of higher dimensional varieties, now that they are easily accessible, will find new applications.

Numerical period computations are used in constructing higher terms in Feynman integrals appearing in perturbative quantum field theory [START_REF] Vanhove | The physics and the mixed Hodge structure of Feynman integrals[END_REF][START_REF] Broadhurst | Knots and numbers in ϕ 4 theory to 7 loops and beyond[END_REF][START_REF] Todorov | Number theory meets high energy physics[END_REF]. Interpretations of the numerical periods as rational linear sums of some expected transcendental numbers play an important role in understanding the structure of these numbers [START_REF] Broadhurst | Knots and numbers in ϕ 4 theory to 7 loops and beyond[END_REF]. A better understanding of the reasons underlying this expected decomposition led to a surprising connection between high energy physics and number theory [START_REF] Todorov | Number theory meets high energy physics[END_REF]. Of direct relevance to our work is the paper [START_REF] Festi | Bhabha scattering and a special pencil of K3 surfaces[END_REF] A concrete application of computing Picard group of K3 surfaces is the proof of the impossibility of resolving a square root in the parametrization of a two-loop Feynman integral [START_REF] Festi | Bhabha scattering and a special pencil of K3 surfaces[END_REF].

Beyond physical interpretations, the behavior of periods is described in the global framework of motives developed in the past decade [START_REF] Brown | Mixed tate motives over Z[END_REF]. The tools we develop open new ways of experimentation. The computation of the Picard number of K3 surfaces in particular play a crucial role in understanding the dimension of the mirror family associated to a given K3 surface [START_REF] Dolgachev | Mirror symmetry for lattice polarized K3 surfaces[END_REF]: computing the lattice structure of the Picard group gives an exact determination for the base of this mirror family.

To further develop number theory as well as its applications to physics and cryptography the LMFDB (The L-functions and Modular Forms Database) [START_REF]The l-functions and modular forms database[END_REF] was constructed (see also [START_REF]Elliptic curves, modular forms and cryptography[END_REF]), partially relying on period computations of curves [START_REF] Lombardo | Computing the geometric endomorphism ring of a genus-2 Jacobian[END_REF][START_REF] Costa | Rigorous computation of the endomorphism ring of a Jacobian[END_REF]. The numerical transcendental method introduced in our work applies well enough to surfaces to allow generation of databases of similar nature 1 .

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 [START_REF] Dolgachev | Luigi Cremona and cubic surfaces[END_REF]. Starting with surfaces of degree four, the Picard group is sensitive to the defining equations and its computation poses an entirely different kind of challenge. Noether and Lefschetz [START_REF] Lefschetz | L'analysis situs et la géométrie algébrique[END_REF] proved that a very general quartic surface has 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 [START_REF] Van Luijk | K3 surfaces with Picard number one and infinitely many rational points[END_REF] 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 [START_REF] Elsenhans | The picard group of a k3 surface and its reduction modulo p[END_REF][START_REF] Elsenhans | On the computation of the Picard group for K3 surfaces[END_REF] but the computational bottleneck 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 [START_REF] Kedlaya | Computing zeta functions via p-adic cohomology[END_REF][START_REF] Abbott | Bounding picard numbers of surfaces using p-adic cohomology[END_REF][START_REF] Costa | Zeta functions of nondegenerate hypersurfaces in toric varieties via controlled reduction in p-adic cohomology[END_REF] and another which deforms the given surface to a simpler one [START_REF] Lauder | Counting solutions to equations in many variables over finite fields[END_REF][START_REF] Pancratz | Improvements to the deformation method for counting points on smooth projective hypersurfaces[END_REF].

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 [START_REF] Charles | On the Picard number of K3 surfaces over number fields[END_REF] 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 [START_REF] Charles | On the Picard number of K3 surfaces over number fields[END_REF][START_REF] Hassett | Effective computation of Picard groups and Brauer-Manin obstructions of degree two K3 surfaces over number fields[END_REF][START_REF] Poonen | Computing Néron-Severi groups and cycle class groups[END_REF]. There is recent work addressing the issue of practicability [START_REF] Festi | A practical algorithm to compute the geometric Picard lattice of K3 surfaces of degree 2[END_REF].

Concerning numerical methods, high precision computation of periods has been successfully applied to many problems concerning algebraic curves [START_REF] Van Wamelen | Examples of genus two CM curves defined over the rationals[END_REF][START_REF] Bruin | Numerical computation of endomorphism rings of Jacobians[END_REF], even with the possibility of a posteriori symbolic certification [START_REF] Costa | Rigorous computation of the endomorphism ring of a Jacobian[END_REF].

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 [START_REF] Sertöz | Computing periods of hypersurfaces[END_REF], based on Picard-Fuchs differential equations [START_REF] Picard | Sur les périodes des intégrales doubles et sur une classe d'équations différentielles linéaires[END_REF], algorithms to compute them, [START_REF] Chyzak | An extension of Zeilberger's fast algorithm to general holonomic functions[END_REF][START_REF] Koutschan | A fast approach to creative telescoping[END_REF][START_REF] Pancratz | Computing Gauss-Manin connections for families of projectives hypersurfaces[END_REF][START_REF] Movasati | Calculation of mixed Hodge structures, Gauss-Manin connections and Picard-Fuchs equations[END_REF][START_REF] Lairez | Computing periods of rational integrals[END_REF][START_REF] Bostan | Creative telescoping for rational functions using the Griffiths-Dwork method[END_REF], and numerical analytic continuation [START_REF] Van Der Hoeven | On effective analytic continuation[END_REF][START_REF] Van Der Hoeven | Fast evaluation of holonomic functions near and in regular singularities[END_REF][START_REF] Chudnovsky | Computer algebra in the service of mathematical physics and number theory[END_REF][START_REF] Mezzarobba | NumGFun: A package for numerical and analytic computation with D-finite functions[END_REF][START_REF] Mezzarobba | Rigorous multiple-precision evaluation of D-finite functions in Sagemath[END_REF]. On the other hand, we use algorithms to compute integer linear relations between vectors of real numbers [START_REF] Ferguson | Generalization of the Euclidean algorithm for real numbers to all dimensions higher than two[END_REF][START_REF] Lenstra | Factoring polynomials with rational coefficients[END_REF][START_REF] Hastad | Polynomial time algorithms for finding integer relations among real numbers[END_REF][START_REF] Buchmann | Computing a lattice basis from a system of generating vectors[END_REF][START_REF] Ferguson | Analysis of PSLQ, an integer relation finding algorithm[END_REF][START_REF] Chen | A New View on HJLS and PSLQ: Sums and Projections of Lattices[END_REF]. We used the computer algebra system SageMath [START_REF]SageMath, the Sage Mathematics Software System[END_REF] with ore algebra-analytic2 [START_REF] Mezzarobba | Rigorous multiple-precision evaluation of D-finite functions in Sagemath[END_REF] for performing numerical analytic continuation and our code base is bundled as numperiods. 3Early experiments have been performed in Magma [START_REF] Bosma | The Magma algebra system. I. The user language[END_REF], with PeriodSuite4 [START_REF] Sertöz | Computing periods of hypersurfaces[END_REF] and periods5 [START_REF] Lairez | Computing periods of rational integrals[END_REF].

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 [START_REF] Ciliberto | General components of the Noether-Lefschetz locus and their density in the space of all surfaces[END_REF] and behaves like the indicator function of the rational numbers on the real line. This suggests that finite precision approximations of the periods of a surface are not enough to formally determine the Picard group. Although sufficiently high precision will compute the Picard group correctly, we do not know what precision is sufficient.

Instead, given the precision of the computation we determine a number B such that the computed lattice contains the sublattice of the Picard group generated by all elements of norm at most B. For example, using 600 decimal digits we have typically B ∼ 10 100 . In terms of degree, it means that the group we compute contains all curves of degree at most B. With over 180,000 examples, we never observed a generator of norm > 100. On the other extreme, small relations may yield periods that are 0 to 600 decimal places but are non-zero. Such a numerical coincidence would invalidate the computation by introducing false generators into the lattice. We quantify the possibility of these errors precisely and relate them to quantities intrinsic to the surface, see §4. If these intrinsics of the surface could be computed they would allow for the certification of the numerical computations, however this computation appears out of reach.

We checked our results against the literature whenever possible. In particular, we compared our Picard number computations against controlledreduction 6 [19] and Shioda's algorithm for Delsarte surfaces defined by a sum of four monomials [START_REF] Shioda | An explicit algorithm for computing the Picard number of certain algebraic surfaces[END_REF]. Whenever a check was possible, our computations gave the correct result.

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. 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 5 explains the situation for higher dimensional hypersurfaces where the general idea of the method as explained in Section 2 applies verbatim. Section 6 summarizes the experimental results obtained for over thousands of quartic surfaces. Section 7 explains how we compute the polarization by fleshing out the argument given in [START_REF] Sertöz | Computing periods of hypersurfaces[END_REF].

Acknowledgments. We would like to thank Bernd Sturmfels for putting us together, for his guidance and constant support. We also would like to thank Alex Degtyarev for explaining to us how to count rational curves in a K3 lattice. We are also grateful for insightful conversations with Simon Brandhorst, John Cannon, Edgar Costa, Stephan Elsenhans, Jon Hauenstein, Marc Mezzarobba, Mateusz Micha lek, Matthias Schütt and Don Zagier, and for the valuable comments of the referees.

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 [START_REF] Movasati | A course in Hodge theory, with emphasis in multiple integrals[END_REF][START_REF] Houzel | La Géométrie Algébrique[END_REF]).

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 [START_REF] Lefschetz | L'analysis situs et la géométrie algébrique[END_REF], see also [START_REF] Movasati | A course in Hodge theory, with emphasis in multiple integrals[END_REF]Chap. 9]. In other words, the Picard group comes with a natural inclusion into the homology group

(2.1) Pic(X) → H 2 (X, Z).
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 γ.

Theorem 2.1 (The Lefschetz (1,1) theorem). A homology class

γ ∈ H 2 (X, Z) is in Pic(X)
if and only if γ ω = 0 for every holomorphic 2-form ω on X.

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 becomes

(2.2) Pic(X) = (a 1 , . . . , a 22 ) ∈ Z 22 22 i=1 a i γ i ω X = 0 .
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

C 3 = {w = 1} ⊂ P 3 . For a 2-cycle γ ⊂ X ∩ {w = 1} we can form a thin tube τ ⊂ C 3 \ X around γ so that (2.3) γ ω X = 1 2π √ -1 τ dx dy dz f (x, y, z, 1)
,

where f is the degree 4 homogeneous polynomial defining X [START_REF] Griffiths | On the periods of certain rational integrals. I, II[END_REF].

In general, when

X has degree d ≥ 4, H 2 (X, Z) has rank m = d 3 -4d 2 + 6d -2 and the space Ω 2 (X) of holomorphic 2-forms on X is of dimension r = d-1 3 . Fixing bases H 2 (X, Z) = Z γ 1 , . . . , γ m , Ω 2 (X) = C ω 1 , .
. . , ω r and applying Lefschetz (1,1) theorem we get:

(2.4) Pic(X) = (a 1 , . . . , a m ) ∈ Z m ∀1 ≤ j ≤ r, m i=1 a i γ i ω j = 0 .
In view of (2.2) and (2.4) we can determine Pic(X) by computing the matrix of periods P = ( γ j ω i ) i,j and then finding integer linear relations between the rows, that is, by computing the integral right kernel of P. The algorithm presented in [START_REF] Sertöz | Computing periods of hypersurfaces[END_REF] to compute the periods of X takes care of the first step. We may then use lattice reduction algorithms [53, p. 525] to find generators for Pic(X).

We briefly recall how periods of X are computed in [START_REF] Sertöz | Computing periods of hypersurfaces[END_REF]. The surface X is put into a single parameter family of surfaces containing the Fermat surface Y = {x d + y d + z d + w d = 0}. 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 Picard-Fuchs differential equations. The value of this one parameter period matrix-as well as its derivatives-at Y are given by 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 [START_REF] Mezzarobba | Rigorous multiple-precision evaluation of D-finite functions in Sagemath[END_REF] to arbitrary precision with rigorous error bounds.

The reconstruction of integer relations between numbers that are only approximately given is not possible in general. However, we can guess generators of the integer relations, and prove that these putative generators are genuine generators when the precision is high enough. Section 4 is devoted to the study of this problem.

                                          
0 0 0 0 0 0 0 0 0 0 0 0 0 -1669083212117905913652734 0 1937019641160560221317687 . . . 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1669083212117905913652734 1937019641160560221317687 . . . 1 0 -1 0 0 0 1 1 0 0 0 0 0 -146511829901195443671789 84478429044587822467823 -365980228690630104919296 . . . 0 0 0 1 0 0 0 0 0 0 0 0 0 -337167720252678310258177 224110151973403946221421 -743116955936487279910552 . . . 0 0 0 0 0 0 0 0 0 0 0 1 -1 357031479253522311483650 768066337666351099432748 940525994719391079998435 . . . 0 0 0 0 1 0 0 1 0 1 0 0 0 -552756671828854153114905 -126018248279583585486071 535095811953165917210863 . . . 0 -1 0 0 0 0 0 1 0 0 -1 0 0 104335431129908645825133 -231616284585318363570849 502730408585962411025306 . . . 0 0 0 0 0 0 0 0 0 0 0 0 -1 -649159586430203173692632 770784867967071100945665 -2152014469737999315531272 . . . 0 0 0 0 0 0 0 0 1 1 0 0 0 277747983934797690835205 -28625739873061372966384 -638732179408358479990097 . . . 1 0 0 0 0 0 0 0 0 0 1 0 0 146511829901195443671790 -84478429044587822467823 365980228690630104919296 . . . 0 0 0 0 0 0 0 0 0 0 -1 1 1 250899146775406645936761 575615030011256031395007 -114830012426104078247291 . . . 0 1 0 0 0 0 1 0 0 -1 0 0 0 104335431129908645825133 -231616284585318363570849 502730408585962411025307 . . . 0 0 0 0 0 -1 0 0 0 0 0 1 -1 -140644950443454586919439 -393058206212350140614235 429933080833930208291557 . . . 0 0 0 0 0 0 0 1 0 0 0 0 0 594933070600140950961561 273156103820314126589096 -671845991848498223316874 . . . 0 0 0 1 0 0 -1 0 0 0 0 0 0 337167720252678310258177 -224110151973403946221421 743116955936487279910552 . . . 0 0 0 0 0 0 0 0 0 0 0 0 1 -824317154838996681984621 177119763197465887754938 -236792300924643740702432 . . . 0 0 0 0 0 0 1 0 0 1 0 0 0 379344119023965108104833 -76972296432673405118395 606366776041154973804541 . . . 0 0 0 0 1 0 0 0 0 0 0 0 0 552756671828854153114905 126018248279583585486070 -535095811953165917210864 . . . 0 0 0 0 0 1 0 0 0 0 0 0 -1 -140644950443454586919440 -393058206212350140614234 429933080833930208291557 . . . 0 0 0 0 0 0 0 0 0 0 0 0 0 -104335431129908645825133 231616284585318363570849 -502730408585962411025307 . . . 0 0 0 0 0 0 0 0 0 0 0 1 0 -467285675585474370500971 -950623161465256990213520 -1255629063127217210042702 . . . 0 0 1 0 0 0 0 0 0 0 0 0 0 -146511829901195443671790 84478429044587822467823 -365980228690630104919296 . . . 0 0 0 0 0 0 0 0 1 0 -1 0 0 -277747983934797690835206 28625739873061372966384 638732179408358479990097 . . . 0 0 0 0 0 0 0 0 0 0 1 0 0 -69025235930677842745100 457102914343586863258366 660652346877586707848817 . . .

                                           Figure 2.1.
Lattice of integer relations between approximate periods. The last 5 columns are omitted.

2.2. An example. Consider the quartic surface X ⊂ P 3 defined by the polynomial

(2.5) f = 3x 3 z -2x 2 y 2 + xz 3 -8y 4 -8w 4 .
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.

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 (2.6)

22 i=1 a i 10 100 γ i ω X = u + v √ -1,
where [-] denotes the rounding of real and imaginary part 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 surprising 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 infer that the Picard number of X is 14 and that the columns of the lower left 22 × 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

                        -4 0 0 2 2 2 0 -3 -1 -2 0 -1 1 -1 0 -4 2 0 -1 1 2 -1 2 0 4 -2 0 0 0 2 -4 0 2 0 -2 -1 -1 0 0 2 1 -1 2 0 0 -4 -1 -1 0 3 2 0 0 0 0 0 2 -1 2 -1 -4 0 1 3 1 2 1 0 -1 -1 2 1 0 -1 0 -4 -1 1 -1 2 -3 -1 1 3 0 2 -2 0 1 -1 -4 1 -2 0 -2 0 2 2 -3 -1 -1 3 3 1 1 -6 -1 -1 1 0 1 -1 -1 2 -1 2 1 -1 -2 -1 -4 1 -2 0 2 0 -2 0 0 0 2 2 0 -1 1 -4 0 1 -1 1 0 4 0 0 1 -3 -2 1 -2 0 -10 -1 0 3 -1 -2 2 0 0 -1 0 0 0 1 -1 -6 2 3 1 0 1 0 -1 1 2 1 2 -1 0 2 -4 0 -1 0 -1 0 -1 3 2 -1 0 1 3 3 0 -10                                                 -4 -5 0 -2 4 3 1 3 -1 6 -2 4 0 2                         Figure 2.2.
Matrix of the intersection product and the coordinates of the hyperplane section in Pic(X).

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 4.2 expresses these quantities independently of any choice occurring in the computation.

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 {γ i } 22 i=1 , see §7. 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.

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 given by (2.7)

T def = ω ∈ H 2 (X, Z) ∀γ ∈ Pic(X), ω • γ = 0 .
We will denote the associated spaces by

T Q = T ⊗ Q and T C = T ⊗ C. Observe that ω X ∈ T C ⊂ H 2 (X, C) by (2.
2). The endomorphism ring E is defined as the subring of all Qlinear maps e : T Q → T Q preserving the intersection product and satisfying e(ω X ) ∈ C ω X , for the canonical extension of e to H 2 (X, C). The map ϕ : E → C defined by e(ω X ) = ϕ(e)ω X is an injective ring morphism and every nonzero element in E is invertible, therefore E is a number field [START_REF] Huybrechts | Lectures on K3 surfaces[END_REF]Corollary 3.3.6]. In fact, E is either totally real or a CM-field [START_REF] Zarhin | Hodge groups of K3 surfaces[END_REF], see also [START_REF] Huybrechts | Lectures on K3 surfaces[END_REF]. Charles [START_REF] Charles | On the Picard number of K3 surfaces over number fields[END_REF] 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 [START_REF] Huybrechts | Lectures on K3 surfaces[END_REF][START_REF] Tretkoff | Periods and special functions in transcendence[END_REF] 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 > 5 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 [START_REF] Van Luijk | K3 surfaces with Picard number one and infinitely many rational points[END_REF] 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 [START_REF] Elsenhans | Kummer surfaces and the computation of the Picard group[END_REF] whether ρ = ρ red if ρ is even and ρ = ρ red -1 if ρ is odd. Charles [START_REF] Charles | On the Picard number of K3 surfaces over number fields[END_REF] settled the question in the negative.

Theorem 2.2 (Charles). The equality ρ red = ρ holds unless E is totally real and the dimension of

T Q over E is odd, in which case ρ red = ρ + dim Q E.
Computation. From the numerical computation of periods, we obtain approximations of a i def = γ i ω X ∈ C for some basis γ 1 , . . . , γ 22 of H 2 (X, Z). The cohomology group H 2 (X, C) is endowed with the dual basis γ * 1 , . . . , γ * 22 so that ω X = i a i γ * i . Once a basis u 1 , . . . , u ρ of the Picard group Pic(X) is computed, a basis v 1 , . . . , v ρ of T = Pic(X) ⊥ ⊆ H 2 (X, Z) is found readily. Let us remark that during our computations, the intersection products γ * i • γ * j are known exactly and therefore the intersection product on T is immediately deduced.

For e ∈ End Q (H 2 (X, Q)) the condition e ∈ E can be rewritten as:

(2.8) ∃λ ∈ C : e(ω X ) = λω X ⇔ ω X , e(ω X ) ω X = ω X , ω X e(ω X ).

Writing A = (a 1 , . . . , a 22 ) t for the coefficient vector of ω X , we can compute the endomorphism ring E via the following formulation:

(2.9)

E = Q • M ∈ Z ρ ×ρ ( Āt M A)A = ( Āt A)M A .
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 §6.

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 hyperplane 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. 7 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. For d > 0 we define the set

M d = D ∈ Pic(X) D 2 = -2, D • h X = d .
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

R d into M d . 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 D 2 = D • (K X + D) = 2g(P 1 ) -2 = -2 [6, §II.11][40, Ex. V.1.3]. This proves that the image of R d lies in M d .
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:

-2 = [C] 2 = [C] • [C ] > 0.
dim H 0 (X, O X (D)) + dim H 0 (X, O X (-D)) ≥ 1 2 D 2 + 2 = 1,
so that either D or -D must be linearly equivalent to an effective divisor. Since D • h X > 0, -D can not be so 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 is the class of a smooth rational curve The Picard group (i.e. matrix of the intersection 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}. function RationalCurves(Pic(X), d) Compute a basis of Pic 0 (X) = {D | D • h X = 0} ⊂ Pic(X) Z ρ Compute a basis of 4 Pic(X) + Zh X ⊂ Pic(X) Compute a basis of Λ = Pic 0 (X) ∩ (4 Pic(X) + Zh X ) ⊂ Pic(X) S ← D ∈ Λ -D 2 = 32 + 4d 2 e.g. with KFP algorithm [47, 36] M d ← 1 4 (D + dh X ) D ∈ S ∩ Z ρ return {D ∈ M d | ∀d < d, ∀D ∈ RationalCurves(Pic(X), d ), D • D ≥ 0} end function
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 neg- ative definite [46, Proposition 1.2.4]. Recalling that h 2 X = 4, we define a map π : Pic(X) → Pic 0 (X) with π(D) = 4D -(D • h X )h X .
The map π maps M d bijectively on to the following set:

(3.2) M d = E ∈ Pic 0 (X) E 2 = -(32 + 4d 2 ) and E + dh X ∈ 4 Pic(X) .
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 [START_REF] Kannan | Improved algorithms for integer programming and related lattice problems[END_REF][START_REF] Fincke | A procedure for determining algebraic integers of given norm[END_REF]. Then, among all such E, we select those where 1 4 (E + dh X ) has integer coordinates to obtain M d . In practice, it is sufficient and more efficient to enumerate the elements of length -(32 + 4d 2 ) in the sublattice π(Pic(X)) = Pic 0 (X) ∩ (4 Pic(X) + Zh X ).

Example 3.3. Take f X = 14x 4 -85x 3 z-2xz 3 +83y 4 -17y 3 w-96w 4 and let X = Z(f X ) ⊂ P 3 . We find that X has Picard number 18 with the following representation of (Pic(X), h X ):

               -4 0 -1 0 -1 2 0 0 0 -2 0 -2 0 1 -1 -1 0 -2 0 -6 -3 -3 -3 3 0 -2 -2 -1 3 -1 -1 3 -1 0 1 0 -1 -3 -4 -2 -2 2 0 0 -1 0 2 0 -2 3 -2 0 1 0 0 -3 -2 -4 -2 2 1 -1 0 -1 1 -1 -1 2 -1 0 1 0 -1 -3 -2 -2 -4 2 0 0 -1 -1 2 0 0 2 -2 0 1 0 2 3 2 2 2 -4 1 2 2 2 -1 2 0 -2 2 0 -1 2 0 0 0 1 0 1 -4 -2 -2 1 1 1 0 0 -1 2 1 1 0 -2 0 -1 0 2 -2 -4 -2 0 1 -1 0 1 0 1 1 0 0 -2 -1 0 -1 2 -2 -2 -4 0 2 0 0 1 -1 1 2 0 -2 -1 0 -1 -1 2 1 0 0 -4 0 -2 1 0 0 -2 -1 -2 0 3 2 1 2 -1 1 1 2 0 -4 0 2 -2 1 -1 -2 -2 -2 -1 0 -1 0 2 1 -1 0 -2 0 -4 0 0 0 0 0 -2 0 -1 -2 -1 0 0 0 0 0 1 2 0 -4 3 -1 1 2 2 1 3 3 2 2 -2 0 1 1 0 -2 0 3 -6 3 0 -2 0 -1 -1 -2 -1 -2 2 -1 0 -1 0 1 0 -1 3 -4 1 2 -1 -1 0 0 0 0 0 2 1 1 -2 -1 0 1 0 1 -4 -2 -1 0 1 1 1 1 -1 1 1 2 -1 -2 0 2 -2 2 -2 -4 -1 -2 0 0 0 0 2 1 0 0 -2 -2 -2 2 0 -1 -1 -1 -4                ,                -2 -1 1 1 3 4 4 0 0 0 0 2 2 0 -2 2 -1 4               
.

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.

4. 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 (the 15th to 22nd columns in Figure 2.1). 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 [53, p. 525], see also [START_REF] Buchmann | Computing a lattice basis from a system of generating vectors[END_REF][START_REF] Chen | A New View on HJLS and PSLQ: Sums and Projections of Lattices[END_REF]. There are many other algorithms for the problem of computing integer relations, in particular Ferguson and Forcade's algorithm [START_REF] Ferguson | Generalization of the Euclidean algorithm for real numbers to all dimensions higher than two[END_REF] and the HJLS [START_REF] Hastad | Polynomial time algorithms for finding integer relations among real numbers[END_REF] and PSLQ [START_REF] Ferguson | Analysis of PSLQ, an integer relation finding algorithm[END_REF][START_REF] Bailey | Integer relation detection[END_REF][START_REF] Feng | The PSLQ algorithm for empirical data[END_REF] families. A strong point in favor of the folklore LLL approach is that efficient LLL implementations are available in most computer algebra systems.

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. 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. However, we can compute a candidate lattice Λ that satisfy a triple alternative: either Λ = Λ; or Λ is not generated by elements of norm less that some explicit number B; or some unexpected numerical cancellation happens (Proposition 4.1).

The computation of a canditate lattice with LLL proceeds as follows. 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 (4.1)

P β = βP + E, with E ∈ [-1 2 , 1 2 ] p×m .
Then, we build the m × (p + m) integer matrix M = P β I m and compute an LLL-reduced basis b 1 , . . . , b m of the lattice spanned by the rows of M . 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 β p m-ρ . In this case, Λ = pr(b 1 ), . . . , pr(b ρ ) , where pr : Z p+m → Z m is the projection on to the last m coordinates.

B ← 1 m 2 -m+1 2 b ρ+1 N ← u ρ ε ← mβ -1 N Return ( Λ, B, N, ε) end function
Equivalently, ε(B) is the largest real number such that Λ B,ε(B) ⊆ Λ. Since ε(B) is nonincreasing as a function of B, the quotient B/ε(B) is strictly increasing as function of B. In particular, for s ≥ 0 we may define a non-decreasing function ϕ with

(4.4) ϕ(s) = max{B ≥ 0 | mB/ε(B) ≤ s}.
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 coefficients in [- 1 2 , 1 2 ], the error matrix

E = P β -βP satisfies E op ≤ 1 2 √ pm ≤ m -1, where
• 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.

Proposition 4.1. Let κ = m -1 2 -m+1 2 .
For any i ∈ {0, . . . , m -1} such that B i ≤ κB i+1 , at least one of the following propositions holds:

(i) {pr(b 1 ), . . . , pr(b i )} is a basis of Λ; (ii) Λ is not generated by elements of norm ≤ κB i+1 ; (iii) ϕ(β) ≤ B i .
Proof. By Lemma 4.2 below, we have

Λ B i ,mβ -1 B i = Λ κB i+1 ,mβ -1 κB i+1 = pr(b 1 ), . . . , pr(b i ) .
If Λ is generated by elements of norm ≤ κB i+1 then Λ ⊆ Λ κB i+1 ,mβ -1 κB i+1 , and therefore Λ ⊆

Λ B i ,mβ -1 B i . If moreover ϕ(β) > B i , then mβ -1 B i ≤ ε(B i )
, by definition of ϕ, and this implies that Λ B i ,mβ -1 B i ⊆ Λ. 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 first show

Λ i ⊆ Λ B,mBβ -1 . Let x = pr(b j ), with j ≤ i. We have x ≤ b j ≤ b i ≤ B. Moreover b j = xP β x , so xP β < b j ≤ B. Since xP = β -1 (xP β -xE), we obtain (4.5) xP ≤ β -1 ( xP β + (m -1) x ) < mBβ -1 . Conversely, let x ∈ Z m such that x ≤ B and xP < mBβ -1 . Let r = xP β x . We check that r ≤ xP β + x ≤ β xP + xE + x (4.6) ≤ 2mB < 2 -(m-1)/2 B i+1 .
The properties of an LLL-reduced basis [60, 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 .

The size of the gap between B rk Λ and B rk Λ+1 can be described more precisely in terms of ϕ(β). Proposition 4.3. Let ρ = rk Λ and let C be the smallest real number such that Λ is generated by elements of norm at most C. For any β > 0:

(i) B ρ ≤ 1 2 κ -1 C; (ii) ϕ(β) ≤ B ρ+1 ; Moreover, if C ≤ 2ϕ(β), then (iii) κB ρ+1 ≤ ϕ(β); (iv) pr(b 1 ), . . . , pr(b ρ ) is a basis of Λ. Proof. For x ∈ Z m let r(x) = xP β x ∈ R. If x ∈ Λ, that is xP = 0, (4.7) r(x) ≤ xP β + x ≤ β xP + xE + x ≤ m x , using P β = βP + E.
In particular, R contains ρ independent elements of norm at most mC.

This implies that B ρ ≤ m2 m-1 2 C = 1 2 κ -1 C; this is (i). For (ii), since Λ has rank ρ, at least one of the pr(b 1 ), . . . , pr(b ρ+1 ) is not in Λ, say pr(b i ), denoted x. Since x = pr(b i ) ∈ Λ, xP = 0 and xP ≥ ε(B ρ+1 ). Moreover B ρ+1 ≥ b i ≥ xP β , because b i = r(x). It follows (4.8) B ρ+1 ≥ β xP -xE ≥ βε(B ρ+1 ) -(m -1)B ρ+1 , which implies ϕ(β) ≤ B ρ+1 . To check (iii), let x ∈ Z m such that x ≤ ϕ(β) and xP = ε(ϕ(β)). By construction, x ∈ Λ. The element r(x) ∈ R satisfies (4.9) r(x) ≤ β xP + xE + x ≤ βε(ϕ(β)) + mϕ(β).
By definition of ϕ, βε(ϕ(β)) = mϕ(β) and therefore r(x) ≤ 2mϕ(β). As shown above, R contains ρ independent elements of norm ≤ mC that project to elements of Λ. The vector r(x) ∈ R does not project on Λ, so R contains ρ + 1 independent elements of norm ≤ m max(C, 2ϕ(β)) = 2ϕ(β). It follows that κB ρ+1 ≤ ϕ(β).

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 (4.10) ε(β) = O β 1-m-rk Λ p and ϕ(s) = O s p m-rk Λ .
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 almost all P ∈ R m×p with a given lattice Λ of integer relations. Generalizing Roth's Theorem on rational approximation of algebraic numbers, Schmidt [START_REF] Schmidt | Linear forms with algebraic coefficients. I[END_REF] 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.

4.2.

Intrinsic error bounds for the computation of the Picard group. Let X ⊂ P 3 be a smooth surface of degree d, (γ i ) i a basis of H 2 (X, Z) and (ω j ) j a basis of the space of holomorphic 2-forms Ω 2 (X) on X. Given a numerical approximation of the period matrix γ i ω j i,j , one can compute with Algorithm 4.1 a sublattice Λ ⊂ H 2 (X, Z) and three positive numbers B, N and ε such that one of the following hypotheses must hold:

1. Λ = Pic(X); or 2. Pic(X) is not generated by the set

{ i a i γ i ∈ Pic(X) | i a 2 i ≤ B 2 }; or 3. there is some i a i γ i ∈ H 2 (X, Z) such that i a 2 i ≤ N 2 and 0 < j i a i γ i ω j 2 ≤ ε 2 .
For B large enough and ε small enough the latter two possibilities can not be realized. As we have no rigorous criterion to exclude Hypotheses 2 or 3, it is important to be able to interpret them as well as possible. Currently, these conditions are not expressed in intrinsic terms and the choice of bases may affect the relevance of the triple (B, N, ε). For example, if γ 1 is replaced by γ 1 + Bγ 2 , Hypothesis 2 becomes plausible. First, we note that Ω 2 (X) is endowed with a natural Hermitian structure, defined by

ω 2 def = X ω ∧ ω.
If we choose a basis (ω j ) j of Ω 2 (X) orthonormal with respect to this Hermitian structure, then the statement of Hypothesis 3 will be independent of the choice of the orthonormal basis. Second, following [46, Example 3.1.7(ii)], we define a canonical norm on H 2 (X, R). Let α j , β j ∈ H 2 (X, R) be the unique homology classes such that for any γ ∈ H 2 (X, R),

(4.11) γ • (α j + √ -1β j ) = γ ω j .
Let U be the R-linear subspace of H 2 (X, R) generated by h X , α j and β j . Let U ⊥ be the orthogonal complement of U with respect to the intersection product so that H 2 (X, R) = U ⊕U ⊥ . The intersection product is positive definite on U and negative definite on U ⊥ . Therefore, the quadratic form q defined on H 2 (X, R) by (4.12)

q(u + v) def = u • u -v • v,
for any u ∈ U and v ∈ U ⊥ , is positive definite. The matrix of q can be computed numerically from the intersection matrix on H 2 (X, R) and bases of U and U ⊥ . We can compute bounds for σ min (q) and σ max (q), the smallest and largest eigenvalues of q respectively. For any v = i a i γ i ∈ H 2 (X, R), we have (4.13) σ min (q)q(v) ≤ i a 2 i ≤ σ max (q)q(v).

Therefore, with N def = σ max (q) 1/2 N and B def = σ min (q) 1/2 B, we obtain the following alternative hypotheses which are independent of the choice of bases:

1. Λ = Pic(X); or 2'. Pic(X) is not generated by the set {v ∈ Pic(X) | q(v) ≤ (B ) 2 }; or 3'. there is some v ∈ H 2 (X, Z) such that q(v) ≤ N 2 and 0 <

j v ω j 2 ≤ ε 2 .
We obtain rigorous bounds on B , N and ε using interval arithmetic. Lemma 4.4 allows for another interpretation for any surface. For simplicity, let us assume X has degree four. If Hypothesis 2' fails for a surface X of degree d, then Pic(X) is not generated by classes of irreducible curves of degree less than B . In our computations, B is typically larger than 10 100 which proves that the sublattice of H 2 (X, R) that we compute contains the classes of all irreducible algebraic curves in X of degree at most 10 100 . In particular, our count of smooth rational curves (up to degree 10 100 ) gives rigorous upper bounds for the correct numbers.

Lemma 4.4. If X is a smooth surface of degree d and D ∈ Pic(X) the class of an irreducible curve of degree δ, then q(D) ≤ 2 + (d -4)δ + 2 d δ 2 . Proof. Any D ∈ Pic(X) decomposes as

D = 1 d (D • h X )h X + D -1 d (D • h X )h X ,
with the first and second terms respectively in U and U ⊥ . Therefore

q(D) = 1 d 2 (D • h X ) 2 (h X • h X ) -D -1 d (D • h X )h X 2 = -D 2 + 2 d (D • h X ) 2 .
In particular, if D is the class of an irreducible curve C ⊂ X of degree δ and arithmetic genus p a , then

-D • (D + K X ) = 2 -2p a ≤ 2 and D • h X = δ. Moreover K X = (d -4)h X [40, Ex. 8.20.3], therefore q(D) = -D • (D + K X ) + (d -4)D • h X + 2 d (D • h X ) 2 ≤ 2 + (d -4)δ + 2 d δ 2 .
5. 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 (5.1)

H 2k dR (X, C) = p+q=2k H p,q (X, C).
Algebraic cycles of dimension k in X give cohomology classes in

(5.2) Hdg k (X) def = H k,k (X, C) ∩ H 2k (X, Z).
As a generalization of Theorem 2.1, the Hodge conjecture predicts that the vector space Hdg k (X) ⊗ Z Q is spanned by algebraic cycles [START_REF] Voisin | The Hodge conjecture[END_REF]. 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 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 (5.3)

γ 1 ω i , . . . , γm ω i ∈ C m .
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 5.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). 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) [START_REF] Sertöz | Computing periods of hypersurfaces[END_REF]. 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 5.1, it remains to compute integral linear relations between the columns of P to recover Hdg k (X). That is the problem studied in §4.

The study of the Hodge groups of cubic fourfolds is an active area of research [START_REF] Ranestad | Variety of power sums and divisors in the moduli space of cubic fourfolds[END_REF][START_REF] Addington | Some non-special cubic fourfolds[END_REF]. 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.

Example 5.2. Let X be the cubic fourfold in P 5 cut out by the equation

(5.5) 6x 3 0 + 10x 0 x 2 x 4 + 9x 0 x 2 x 5 + 4x 3 1 + 2x 1 x 2 2 + 4x 3 2 + 3x 3 3 + 4x 3 4 + 9x 3 5 .
We find that Hdg 2 (X) is of rank 13.

6. Experimental results. We put our methods into practice and computed the Picard lattices of thousands of smooth 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 may take weeks to compute. Further work will address the problem.

As a first experiment, we performed 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 [START_REF] Sertöz | Computing periods of hypersurfaces[END_REF] 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 this graph, we computed the Picard group, the polarization, the intersection product, the endomorphism ring and the number of smooth Defining polynomial Picard number [START_REF] Costa | Effective computations of Hasse-Weil zeta functions[END_REF] and symbolic line counting. Afterwards, we extended our collection of quartics by following a random walk in the space of homogeneous polynomials of degree 4 in 4 variables, where each step adds or substracts a monomial. We computed the periods of more than 180,000 quartics as well as the related data: Picard lattice, endomorphism ring of the transcendental lattice, number of smooth rational curves of degree at most 4. 8We give a few interesting examples below. Of course, all the assertions were obtained through numerical computations, so they are not proved. Example 6.1. {x 4 + y 3 z + xyzw + z 3 w + yw 3 = 0}. This surface has Picard number 19. It contains no smooth rational curves of degree < 4 but has 133056 smooth rational curves of degree 4. These generate the Picard group. Example 6.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 number of smooth rational curves inside X are as follows. It contains no smooth rational curves of degree ≤ 21 (the space complexity of the search makes it difficult to search for higher degree curves). The Picard lattice is generated by smooth quartic curves of genus 1. It would be interesting to settle the existence of rational curves on this surface. Example 6.5. {wx 3 + w 3 y + y 4 + xz 3 + z 4 = 0}. This surface has Picard number 1. The first example of such a quartic surface defined over Q was given recently by van Luijk [START_REF] Van Luijk | K3 surfaces with Picard number one and infinitely many rational points[END_REF]. As Picard number 1 is the smallest possible value, implementations of finite characteristic methods, e.g. controlledreduction [START_REF] Costa | Effective computations of Hasse-Weil zeta functions[END_REF], can readily prove that this is the correct number.

wx 3 + w 3 y + y 4 + xz 3 + z 4 1 x 3 y + z 4 + y 3 w + zw 3 4 x 3 y + y 4 + z 3 w + yw 3 + zw 3 6 w 3 x + x 4 + wx 2 z + x 3 z + xy 2 z -y 3 z + wxz 2 + x 2 z 2 -xz 3 + z 4 7 x 3 y + z 4 + y 3 w + xw 3 + w 4 8 w 4 + wx 2 y + y 4 + x 3 z -xy 2 z + z 4 9 x 3 y + z 4 + y 3 w + w 4 10 w 4 + x 4 + x 2 y 2 + y 4 -w 3 z -2 xy 2 z + x 2 z 2 + z 4 11 x 3 y + y 4 + z 3 w + x 2 w 2 +

7.

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 [68, §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 [START_REF] Sertöz | Computing periods of hypersurfaces[END_REF]. 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 [START_REF] Sertöz | Computing periods of hypersurfaces[END_REF] 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:

S = (1 -t -1 0 ) • • • (1 -t -1 n )D,
where summation is union and negation is change of orientation [START_REF] Pham | Formules de Picard-Lefschetz généralisées et ramification des intégrales[END_REF]. 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 [START_REF] Sertöz | Computing periods of hypersurfaces[END_REF]. We will now add one more cycle to complete ). 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, (7.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:

1 = M 1 . 3 . 2 .

 1132 Typically, the inclusion R d → M d is strict. We now demonstrate that with knowledge of M d for each d ≤ d one can compute the image of R d in M d . For each d > 0 define inductively a subset N d ⊂ M d as follows: (3.1) N d = D ∈ M d ∀d < d, ∀D ∈ N d , D • D ≥ 0 . Note that when d = 1 there are no constraints and we have N Proposition For d > 0, the image of the inclusion R d → M d is N d . 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 [40, V.1.6] we get:

[ 40 ,Algorithm 3 . 1

 4031 Ex. IV.1.8]. Furthermore, let d = C i • h X and observe d ≤ d. By definition of N d we 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 Finding rational curve classes on smooth quartic surfaces. Input.

4. 1 .Algorithm 4 . 1

 141 Quantitative results. For B > 0 and ε > 0 let Λ B,ε be the lattice(4.2) Λ B,ε = x ∈ Z m | x ≤ B and xP < ε . For B ≥ 1, let (4.3) ε(B) = min { uP | u ∈ Z m , u ≤ B and uP = 0} .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.Fail or return a sublattice Λ ⊂ Z m and B, N, ε > 0.Postcondition. Either Λ = {x ∈ Z m | xP = 0}; or {x ∈ Z m | xP = 0} isnot generated by vectors of norm at most B; or ∃x ∈ Z m : x ≤ N , xP ≤ ε and xP = 0. function IntegerRelationLattice(Q, β) Compute an LLL-reduced basis b 1 , . . . , b m of the lattice spanned by the rows Q I m . Find ρ such that b ρ ≤ 2 -m b ρ+1 and β p m-ρ ≈ b ρ+1 . Fail if there is none. Λ ← pr(b 1 ), . . . , pr(b ρ ) , where pr : Z p+m → Z m takes the last m coordinates.

Lemma 4 . 2 .

 42 For any i ∈ {0, . . . , m -1} and any B ∈ [B i , κB i+1 ] Λ B,mBβ -1 = pr(b 1 ), . . . , pr(b i ) .

Figure 6 . 2 .Example 6 . 3 .

 6263 Figure 6.2. Occurence of Picard numbers in the graph of 4-nomials and 5-nomials

  {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

  The coefficients {a β } β∈B ⊂ Q of this relation remain constant as we carry the basis {[L]} ∪ {γ β } β∈B to a basis of H n (X, Z). The problem of computing the periods of L as it deforms is therefore reduced to computing the coefficients {a β } β∈B . Put an ordering on B and let(7.3) a B,L = (a β ) β∈B ∈ Q #B

  w 4 12 w 4 + 3 x 4 + wy 3 + y 2 z 2 + wz 3 + 2 xz 3 13 x 3 y + y 4 + z 3 w + yw 3 + w 4 14 x 3 y + y 3 z + z 4 + xy 2 w + zw 3 15 x 3 y + y 4 + z 3 w + xyw 2 + y 2 w 2 + w 4 16 x 3 y + y 4 + z 4 + x 2 w 2 + zw 3 17 x 3 y + x 3 z + y 3 z + yz 3 + w 4 18 x 3 y + z 4 + y 3 w + xyzw + xw 3 19 x 3 y + z 4 + y 3 w + xw 3 20

Figure 6.1. Specimen polynomials for each Picard number found rational curves of degree up to 4. We found quartic surfaces with Picard number 4, 6, 8, 10, 12, 14, 15, 16, 17, 18, 19 and 20, see Figures 6.1 and 6.2. When possible, we checked that our results were consistent with Shioda's formula for 4-nomial quartic surfaces [69], reduction methods with Costa's implementation

A proof-of-concept database containing more that 180, 000 quartic surfaces is available at https://pierre. lairez.fr/quarticdb.

http://marc.mezzarobba.net/code/ore algebra-analytic

https://gitlab.inria.fr/lairez/numperiods

https://github.com/emresertoz/PeriodSuite

https://github.com/lairez/periods

https://github.com/edgarcosta/controlledreduction

We are indebted to Alex Degtyarev for sharing his understanding with us.

Results are compiled at https://pierre.lairez.fr/quarticdb.

denote the row vector of coefficients defined in (7.2).

Let b B,L = ([L] • γ β ) β∈B ∈ Q #B be the intersection numbers of L with the Pham basis and let M B = (γ β • γ β ) β,β ∈B be the matrix of intersections of the Pham basis. We see that a B,L = M -1 B b B,L so it remains to compute M B and b B,L . Fix d ≥ 2 and define a function χ : Z → {-1, 0, 1} as follows:

. The Pham cycles t β S and t β S intersect as follows:

For a proof of Proposition 7.1 see any one of [START_REF] Arnold | Singularities of differentiable maps[END_REF][START_REF] Movasati | A course in Hodge theory, with emphasis in multiple integrals[END_REF][START_REF] Looijenga | Fermat varieties and the periods of some hypersurfaces[END_REF]. We reformulated the statement here for the choices that were made in [START_REF] Sertöz | Computing periods of hypersurfaces[END_REF] and in the style that was first communicated to us by Degtyarev and Shimada.

Define the function τ d : Z → {-1, 0, 1} where:

The intersection pairing of the linear space L with the translates of the Pham cycle S can be expressed as follows:

Lemma 7.2 is proven by a straightforward application of Theorem 2.2 in [START_REF] Degtyarev | On the topology of projective subspaces in complex Fermat varieties[END_REF].

Example 7.3. Let us consider quartic surfaces in P 3 , that is d = 4 and n = 2. Using Corollary 4.8 [START_REF] Sertöz | Computing periods of hypersurfaces[END_REF] we find B = {(0, 0, 0, 0), (0, 0, 1, 0), (0, 1, 0, 0), (1, 0, 0, 0), (0, 0, 2, 0), (0, 1, 1, 0), (1, 0, 1, 0), (0, 2, 0, 0), (1, 1, 0, 0), (2, 0, 0, 0), (0, 1, 2, 0), (1, 0, 2, 0), (0, 2, 1, 0), (1, 1, 1, 0), (2, 0, 1, 0), (1, 2, 0, 0), (2, 1, 0, 0), (0, 2, 2, 0), (1, 1, 2, 0), (2, 0, 2, 0), (1, 2, 1, 0)}.

With respect to this basis, and the ordering presented above, we find that the vector a B,L of (7.3) is given by

). The set {γ β } β∈B is completed to a basis with the addition of [L]. In this basis, Equation (7.2) gives us the coordinates of the polarization: (7.5) h X = (0, 4, -2, 0, 0, -2, 4, 0, -2, 0, -3, -1, 1, 2, 1, 3, -1, 0, 0, -2, 2, 4).