Florent Bréhard
email: florent.brehard@univ-lille.fr

Nicolas Brisebarre
email: nicolas.brisebarre@ens-lyon.fr

Efficient and Validated Numerical Evaluation of Abelian Integrals

Keywords: CCS Concepts:, Mathematics of computing → Quadrature, Ordinary differential equations, Approximation, Interval arithmetic, Interpolation, Arbitrary-precision arithmetic Hilbert's 16th problem, Abelian integral, limit cycles, rigorous numerics, trigonometric polynomial interpolation, Newton-like operator

Abelian integrals play a key role in the infinitesimal version of Hilbert's 16th problem. Being able to evaluate such integrals -with guaranteed error bounds -is a fundamental step in computer-aided proofs aimed at this problem. Using interpolation by trigonometric polynomials and quasi-Newton-Kantorovitch validation, we develop a validated numerics method for computing Abelian integrals in a quasi-linear number of arithmetic operations. Our approach is both effective, as exemplified on two practical perturbed integrable systems, and amenable to an implementation in a formal proof assistant, which is key to provide fully reliable computer-aided proofs.

INTRODUCTION

After more than a century since the early stages of automated reasoning and mechanized theorem proving, notably with Hilbert's program in the beginning of the 20th century, it can be safely said that machines have not replaced mathematicians so far. Yet, for an increasing number of problems involving intense calculations, algorithms running on computers have already proven to be much more efficient than pen-and-paper work. Rather than fully machine-generated, the resulting proofs are computer-assisted. The field of dynamical systems, in particular, has benefited from automated techniques in computer algebra, numerical analysis and rigorous numerics over the past decades. Famous proofs highlighting these achievements are, for instance, the universality of the Feigenbaum constant [START_REF] Lanford | A computer-assisted proof of the Feigenbaum conjectures[END_REF]], the existence of chaos in the Kuramoto-Sivashinsky equations [START_REF] Wilczak | Chaos in the Kuramoto-Sivashinsky equations-a computer-assisted proof[END_REF]], and the (almost) finiteness of relative equilibria for the 5-body problem of celestial mechanics [START_REF] Albouy | Finiteness of central configurations of five bodies in the plane[END_REF].

The proliferation of computer-assisted proofs raises several central questions concerning their acceptability and utility by the mathematical community. 1. How efficient can algorithms be made to tackle hard problems from a computational point of view? 2. Can we trust the algorithms used to compute the solutions to these problems? 3. Can we moreover trust the implementations of these algorithms, written in practical programming languages and run on real world computers?

Having those questions in mind, this article deals with a computational problem originating from the infinitesimal Hilbert's 16th problem in dynamical systems, namely the evaluation of so-called Abelian integrals. The challenge is the following, with further details postponed to the rest of the introduction.

Problem 1.1. For a polynomial potential function 𝐻 : R 2 → R, a rational rescaling factor 𝜇 : R 2 → R, a polynomial perturbation (𝑃, 𝑄) : R 2 → R 2 and a regular1 level value ℎ, compute the Abelian integral along the oval Γ(ℎ) in the plane:

ℑ(ℎ) = ∫ Γ (ℎ)
𝑃 (𝑥, 𝑦)d𝑦 -𝑄 (𝑥, 𝑦)d𝑥 𝜇 (𝑥, 𝑦)

, (1.1)
where an oval Γ(ℎ) denotes a bounded connected component of the level set 𝐻 (𝑥, 𝑦) = ℎ for a regular value ℎ2.

We propose a validated numerical algorithm, which, given polynomial or rational functions 𝐻, 𝜇, 𝑃, 𝑄, together with the level ℎ, computes an interval that contains the exact value of ℑ(ℎ).

It is intended to be part of computer-assisted proofs around infinitesimal Hilbert's 16th problem, as in the examples at the end of this article. The method has to be both efficient -since a very large number of accurate digits may be necessary to return a tight enough enclosure of ℑ(ℎ) -and strongly reliable -since the (sometimes quite ill-conditioned) numerical calculations are actual parts of the computer-assisted proof. To overcome this, we combine higher-order techniques based on Fourier series with fixed-point a posteriori validation to rigorously approximate the integration path and perform the integration. We obtain an algorithm with quasi-linear complexity in terms of arithmetic operations (see main Theorem 1.4), for which we provide an efficient implementation3 in Julia4 [START_REF] Bezanson | Julia: a fresh approach to numerical computing[END_REF]]. The method is sufficiently self-contained and is tailored to the specific needs of a formal proof. This will allow us to implement the algorithm within a proof assistant like Coq [START_REF] Bertot | Interactive Theorem Proving and Program Development: Coq'Art: The Calculus of Inductive Constructions[END_REF] in the near future. The end goal is to provide an efficient and fully certified calculator for the infinitesimal Hilbert's 16th problem.

Hilbert's 16th problem

In 1900, at the International Congress of Mathematics held in Paris, David Hilbert presented ten open problems in mathematics, and later published a more comprehensive list of 23 problems [START_REF] Hilbert | Mathematische Probleme. Vortrag, gehalten auf dem internationalen Mathematiker-Congress zu Paris[END_REF]] aimed at challenging the mathematical community. Today, most of the Hilbert problems have been resolved (two of them were deemed to be unresolvable), but a few ones still remain open: one of these is Hilbert's 16th problem.

Hilbert's 16th problem has two distinct parts: one in real algebraic geometry, and one in dynamical systems. We shall address the latter which asks for H (𝑛) -the maximal number of limit cycles (i.e., isolated periodic orbits [Christopher and Li 2007, §II.1.1]) the family of two-dimensional polynomial vector fields of (total) degree at most 𝑛 can display. Note that the bound H (𝑛) should be uniform, that is, it should not depend on the particular polynomial vector field, only on its degree 𝑛. As of today, this question is not resolved even in the simplest case 𝑛 = 2. Even finding non-trivial lower bounds for H (𝑛) appears to be very hard.

1It means that ∇𝐻 (𝑥, 𝑦) ≠ 0 for all (𝑥, 𝑦) s.t. 𝐻 (𝑥, 𝑦) = ℎ. 2This makes Γ (ℎ) a Jordan curve, 𝐶 ∞ -diffeomorphic to the unit circle. 3The repositories for our code are available from https://gitlab.inria.fr/abintvalid. 4https://julialang.org/ In light of the lack of progress regarding bounds for H (𝑛), in the mid-seventies, V.I. Arnold [START_REF] Arnold | Loss of stability of self-induced oscillations near resonance, and versal deformations of equivariant vector fields[END_REF][START_REF] Arnold | Ten problems[END_REF] proposed to study a restricted version of the original problem, now known as the infinitesimal (or weak, or tangential) Hilbert's 16th problem. Rather than considering the class of all polynomial vector fields of a certain degree, Arnold suggested that only small perturbations of Hamiltonian polynomial vector fields be considered. Thus, the corresponding question can be asked: Problem 1.2. Consider the differential system in R 2 : where 𝐻 (𝑥, 𝑦) is a polynomial of degree at most 𝑛 + 1. The polynomials 𝑃 and 𝑄 of degree at most 𝑛, and the positive number 𝜀 > 0 define the small perturbation of the Hamiltonian system.

         𝑥 = -
Is there a bound Z(𝑛) on the number of limit cycles the system (1.2) can have (for small 𝜀), that only depends on the degree 𝑛?

For the infinitesimal problem, significant progress has been made, notably the proof of the finiteness of Z(𝑛) for all 𝑛 [START_REF] Binyamini | On the number of zeros of Abelian integrals[END_REF], and the uniform bound for the quadratic case Z(2) = 2 [START_REF] Chen | A unified proof on the weak Hilbert 16th problem for 𝑛 = 2[END_REF]

The Poincaré-Pontryagin theorem and Abelian integrals

The study of perturbed Hamiltonian (or even integrable) systems heavily depends on a theorem by Poincaré and Pontryagin, that makes a strong connection between the existence of a limit cycle and a zero of the Abelian integral (1.1).

The Poincaré return map is the key tool to understand this connection. Consider the unperturbed Hamiltonian system in a region made of a continuous foliation of periodic orbits. Take a transversal Σ, that is an open interval of a line crossing all orbits it encounters non-tangentially. This interval may be parameterized by the parameter ℎ (the energy level) of the unperturbed system for some domain ℎ -< ℎ < ℎ + . By continuity, for every ℎ ∈ (ℎ -, ℎ + ) and sufficiently small 𝜀, the trajectory originating from the point Σ(ℎ) on the transversal in the perturbed Hamiltonian system will cross Σ again. The Poincaré return map Π associates to ℎ and 𝜀 the parameter Π(ℎ, 𝜀) corresponding to the first point of return to Σ. We call 𝑑 (ℎ, 𝜀) = Π(ℎ, 𝜀)ℎ the displacement function. Clearly, the point Σ(ℎ) belongs to a periodic orbit of the perturbed system if and only if 𝑑 (ℎ, 𝜀) = 0, and this is a limit cycle if and only if the zero of 𝑑 (•, 𝜀) is isolated.

The Poincaré-Pontryagin theorem roughly states that, for a Hamiltonian (or more generally, integrable) system perturbed by 𝜀-small terms, the Abelian integral -which is nothing but the integral of the perturbation along a non-perturbed periodic orbit -is the first-order approximation in 𝜀 of the Poincaré return map. We now state the following adaptation of this result (see [START_REF] Christopher | Limit cycles of differential equations[END_REF] for a proof) that will prove useful in the sequel. . Let 𝐻 : R 2 → R be a real analytic potential function, 𝑃, 𝑄 : R 2 → R real analytic functions, 𝜀 > 0, and 𝜇 : R 2 → R an analytic rescaling factor. Consider the perturbed integrable system: (1.3) Let Γ(ℎ) be an oval of 𝐻 of level ℎ over which 𝜇 does not vanish. Then the displacement function 𝑑 (ℎ, 𝜀) is approximated as 𝑑 (ℎ, 𝜀) = 𝜀ℑ(ℎ) + 𝑂 (𝜀 2 ), as 𝜀 → 0,

         𝑥 = -𝜇 (𝑥,
where the Abelian integral ℑ(ℎ) was defined in (1.1).

In particular, 𝑑 (ℎ, 𝜀) and ℑ(ℎ) have the same sign for small 𝜀 > 0. Hence, the number of isolated zeros of ℑ(ℎ) where a change of sign occurs (in particular, simple zeros of ℑ(ℎ)) provides a lower bound for the number of limit cycles of (1.3) that exist for small 𝜀 > 0. Considering zeros of ℑ(ℎ) of higher multiplicity it is possible to get an upper bound on the number of limit cycles that can bifurcate from the unperturbed periodic orbit(s) Γ(ℎ), see [START_REF] Christopher | Limit cycles of differential equations[END_REF].

A pessimistic, yet constructive, upper bound for Z(𝑛) was obtained by the authors of [START_REF] Binyamini | On the number of zeros of Abelian integrals[END_REF]] by bounding the number of zeros of ℑ(ℎ) in terms of the degree 𝑛 only, using the Picard-Fuchs differential equations satisfied by ℎ ↦ → ℑ(ℎ). On the other hand, counting the sign alternations of the Abelian integrals for well chosen integrable systems and perturbations can provide lower bounds for some H (𝑛) [START_REF] Li | A cubic system with thirteen limit cycles[END_REF] or Z(𝑛) [START_REF] Johnson | An improved lower bound on the number of limit cycles bifurcating from a quintic Hamiltonian planar vector field under quintic perturbation[END_REF].

Rigorous computation of Abelian integrals: challenge and related works

In general, we lack closed forms for the Abelian integrals, and except for specific families of systems where the analytic behavior of these integrals was investigated by pure -but involved! -pen-and-paper techniques (e.g., [START_REF] Li | A cubic system with thirteen limit cycles[END_REF]] for H (3) ⩾ 13), the numerical evaluation of ℑ(ℎ) requires a delicate strategy to certify a lower bound on the number of sign changes. Since these computations are a part of the proof, the results must come with strong guarantees, typically validated bounds. This is the field of rigorous (or validated) numerics [START_REF] Moore | Interval analysis[END_REF][START_REF] Tucker | Validated numerics: a short introduction to rigorous computations[END_REF], where interval and higher order methods are used to enclose the actual value (number or function) in a guaranteed set-valued representation. Moreover, high precision is often necessary, even though only the sign of ℑ(ℎ) needs to be validated in the end. This is because Abelian integrals with many zeros are often expressed as linear combinations of simpler integrals. These linear combinations are often subject to high cancellations. (see the example of Section 7.1 with Figure 2).

A first work to mention in this regard is the rigorous computation of Poincaré maps [START_REF] Kapela | Recent advances in rigorous computation of Poincaré maps[END_REF]]5 using the celebrated CAPD6 [START_REF] Kapela | CAPD::DynSys: a flexible C++ toolbox for rigorous numerical analysis of dynamical systems[END_REF]] library for validated numerics in connection with dynamical systems. Although CAPD's efficiency is well-established, we prefer not using rigorous (nonlinear) ODE integrators. Our aim is to keep the method as minimal as possible, in order to ease its implementation in a formal proof assistant.

Several works deal with the rigorous approximation of implicitly defined curves. The closest related one [START_REF] Johnson | On a computer-aided approach to the computation of Abelian integrals[END_REF], by one of the authors of this article, computes the Abelian integral by subdividing the interior of the oval Γ(ℎ) and enclosing its border using parallelotopes. A similar approach can be found in [START_REF] Martin | Certified parallelotope continuation for one-manifolds[END_REF] for the rigorous continuation of one-dimensional varieties. Other techniques we can mention for rigorous path continuation, yet for a different problem, namely rigorous homotopy tracking, are [START_REF] Beltrán | Robust certified numerical homotopy tracking[END_REF][START_REF] Xu | An Approach for Certifying Homotopy Continuation Paths: Univariate Case[END_REF]. A common feature of these methods is that the curve is represented by piecewise low-order approximations, with the accuracy controlled by the step size.

Since we target high precision without compromising the efficiency, we promote the use of global, high-order approximations for a suitable parameterization of the oval Γ(ℎ). Parameterizing algebraic curves is a deep topic in algebraic geometry and geometric modelling for computer-aided design. However, existence theorems for exact (polynomial, rational or trigonometric) parameterizations are mostly restricted to curves of genus 0 [START_REF] Abhyankar | Computations with Algebraic Curves[END_REF][START_REF] Hong | Algorithms for trigonometric curves (simplification, implicitization, parameterization)[END_REF], while approximate methods (piecewise rational, splines, etc.) developed for more general curves most of the time come without guaranteed error bounds (see e.g., [START_REF] Bajaj | Piecewise rational approximations of real algebraic curves[END_REF][START_REF] Gao | Rational quadratic approximation to real algebraic curves[END_REF][START_REF] Yang | An evolution-based approach for approximate parameterization of implicitly defined curves by polynomial parametric spline curves[END_REF]). These methods usually target the standard binary64 precision (i.e. double precision), which is not sufficient for our problem. A notable exception is the computation of certain elliptic integrals [START_REF] Bost | Moyenne arithmético-géométrique et périodes des courbes de genre 1 et 2[END_REF], with quadratic convergence due to the arithmetic-geometric mean.

Our approach and contributions

We propose a higher-order method where the oval Γ(ℎ) is approximated not by parallelotopes, but by interpolation trigonometric polynomials. This allows us to take advantage of the excellent (typically, exponentially converging) approximation properties in Fourier analysis, while keeping a minimalist framework for rigorous computations, namely trigonometric polynomials, which we develop in Section 2. The result is a fully automated and rigorous Algorithm AbIntValid stated in Section 5, together with the following result, that we shall state more precisely as Theorem 5.2. Theorem 1.4. Algorithm AbIntValid computes an interval enclosure 𝑰 for ℑ(ℎ) (cf. (1.1)) in 𝑂 (𝑁 log 𝑁 ) arithmetic operations, where 𝑁 is the degree used for trigonometric polynomials to approximate the curve Γ(ℎ). As 𝑁 → ∞, the diameter of 𝑰 tends exponentially fast to 0.

The theorem says that in order to compute ℑ(ℎ) rigorously to 𝑘 digits of accuracy, one has to use a degree 𝑁 = 𝑂 (𝑘), where the constant hidden in this notation depends on the Hamiltonian system and ℎ, but not on 𝑘.

Our approach is based on an a posteriori validation scheme, where we first approximate the oval Γ(ℎ) with a trigonometric parameterization, and afterwards validate a tube around it, over which the Abelian integral is eventually computed rigorously. The main algorithm AbIntValid is decomposed as the following subroutines:

(1) First, the oval Γ(ℎ) is approximated by the two routines OvalApproxInit and OvalAp-proxRefine (both described in Section 3). The result is a pair of degree 𝑁 trigonometric polynomials ( x, y), such that 𝑡 ∈ [0, 2𝜋] ↦ → ( x (𝑡), y (𝑡)) is a smooth approximate parameterization of Γ(ℎ).

(2) Next, Algorithm OvalValid (described in Section 4) computes an a posteriori error bound for the distance between ( x, y) and Γ(ℎ), thus defining a tube. This relies on the Newton-Kantorovitch validation principle summarized in Section 2.3. (3) After that, Algorithm AbIntValidQuad (described in Section 5) rigorously computes an interval containing ℑ(ℎ) by integrating the perturbation over the tube. Therefore, the smaller the obtained error bound is, the tighter the final enclosure is; hence the need for accurate validation techniques. Finally, in Section 6, we present our implementation and compare it with alternative approaches, and we assess its efficiency on two practical examples borrowed from the literature in Section 7.

APPROXIMATION TOOLS FOR FOURIER ANALYSIS

In this section, we collect useful results for the sequel of the text. After recalling some basic facts about interpolation by trigonometric polynomials, we introduce rigorous trigonometric approximations, a tool for performing validated computations, and an arithmetic on these objects.

Interpolation of periodic functions

Throughout this paper, we shall need to approximate parameterizations of the ovals Γ(ℎ). Seeing that these parameterizations are continuous and 2𝜋-periodic, a natural choice for approximating them is via trigonometric polynomials [START_REF] Zygmund | Trigonometric series[END_REF]].

Definition 2.1 (Trigonometric Polynomials).

A degree 𝑁 Trigonometric Polynomial (TP) f is a real-valued trigonometric polynomial over [0, 2𝜋]:

f (𝑡) = 𝑎 0 + 𝑁 ∑ 𝑘=1 𝑎 𝑘 cos(𝑘𝑡) + 𝑏 𝑘 sin(𝑘𝑡) ,
where 𝑎 𝑁 or 𝑏 𝑁 is nonzero. We represent f in terms of its 2𝑁 + 1 coefficients.

Let TP 𝑁 (R) denote the vector space of trigonometric polynomials with real coefficients and degree at most 𝑁 , and let (𝐶 0 2𝜋 , ∥ • ∥ ∞ ) denote the vector space of real continuous 2𝜋-periodic functions, equipped with the supremum norm:

∥𝑓 ∥ ∞ = max 0⩽𝑥⩽2𝜋 |𝑓 (𝑥)|.
A usual way to obtain quasi-optimal trigonometric polynomial approximants is to consider the Fourier series expansion of the function 𝑓 ∈ 𝐶 0 2𝜋 truncated to order 𝑁 . And yet, from a computational point of view, it is often more interesting to consider the degree (at most) 𝑁 trigonometric polynomial that interpolates the function under consideration at the equispaced points 2𝑘𝜋 2𝑁 +1 2𝑁 𝑘=0 [START_REF] Wright | Extension of Chebfun to periodic functions[END_REF]. We then define the interpolation operator

I 𝑁 : 𝐶 0 2𝜋 → TP 𝑁 (R) 𝑓 ↦ → f s.t. f 2𝑘𝜋 2𝑁 + 1 = 𝑓 2𝑘𝜋 2𝑁 + 1
for 𝑘 = 0, . . . , 2𝑁 .

(2.1)

The Lebesgue constant Λ 𝑁 [START_REF] Powell | Approximation theory and methods[END_REF]] associated to the trigonometric interpolation operator I 𝑁 is the operator norm of I 𝑁 , that is to say

Λ 𝑁 = sup 𝑓 ∈𝐶 0 2𝜋 ,∥ 𝑓 ∥ ∞ =1 ∥I 𝑁 (𝑓 ) ∥ ∞ .
A first nice feature of these interpolation trigonometric polynomials is that Λ 𝑁 = 𝑂 (log 𝑁 ) [START_REF] Ehlich | Auswertung der Normen von Interpolationsoperatoren[END_REF], which makes them quasi-optimal approximations with respect to ∥ • ∥ ∞ .

A second nice feature is that we can compute the coefficients of the interpolation polynomial by the Fast Fourier Transform, hence in 𝑂 (𝑁 log(𝑁 )) arithmetic operations. Indeed, if f (𝑡) = 𝑎 0 + 𝑁 𝑘=1 (𝑎 𝑘 cos(𝑘𝑡) + 𝑏 𝑘 sin(𝑘𝑡)) and 𝑓 𝑗 = f 2𝑗𝜋 2𝑁 +1 for 𝑗 = 0, . . . , 2𝑁 , we have, for 𝑘 = 1, . . . , 𝑁 ,

𝑎 0 = 1 2𝑁 + 1 2𝑁 ∑ 𝑗=0 𝑓 𝑗 , 𝑎 𝑘 = 1 2𝑁 + 1 2𝑁 ∑ 𝑗=0 𝑓 𝑗 cos 2𝑗𝑘𝜋 2𝑁 + 1 , 𝑏 𝑘 = 1 2𝑁 + 1 2𝑁 ∑ 𝑗=0 𝑓 𝑗 sin 2𝑗𝑘𝜋 2𝑁 + 1 .
Having this in mind, we define

FFT 𝑁 : R 2𝑁 +1 → TP 𝑁 (R) (𝑓 𝑗 ) 2𝑁 𝑗=0 ↦ → f s.t. f (𝑡 𝑗 ) = 𝑓 𝑗 for 0 ⩽ 𝑗 ⩽ 2𝑁 ,
and its inverse transform

IFFT 𝑁 : TP 𝑁 (R) → R 2𝑁 +1 f ↦ → ( f (𝑡 𝑗 )) 2𝑁 𝑗=0 where 𝑡 𝑗 = 2𝜋 𝑗 2𝑁 +1 for 0 ⩽ 𝑗 ⩽ 2𝑁 ,
which allows us to interpolate a function and to evaluate a TP on the equispaced grid in only 𝑂 (𝑁 log(𝑁 )) arithmetic operations. By combining both routines, we can compute the degree 2𝑁 product of two degree 𝑁 TPs with the same asymptotic complexity 𝑂 (𝑁 log(𝑁 )). In practice, this outperforms the "naive" multiplication algorithm in 𝑂 (𝑁 2 ) arithmetic operations for degrees 𝑁 of the order of magnitude of 100 or more. We assume this 𝑂 (𝑁 log(𝑁 )) asymptotic complexity for the product of degree 𝑁 TPs in what follows.

When working with real trigonometric approximations, a convenient way to bound the supremum norm of functions is to use the ∥ • ∥ ℓ 1 norm directly defined from the Fourier coefficients 𝑎 𝑘 , 𝑏 𝑘 of 𝑓 :

∥ 𝑓 ∥ ℓ 1 = |𝑎 0 | + ∞ ∑ 𝑘=1 |𝑎 𝑘 | + |𝑏 𝑘 | .
This norm is well-defined as long as the sum converges, and then satisfies ∥ 𝑓 ∥ ℓ 1 ⩾ ∥𝑓 ∥ ∞ . Note that the sum does not necessarily converge for 𝑓 ∈ 𝐶 0 2𝜋 . However, this becomes true under very mild regularity assumptions, and a fortiori7 for analytic 𝑓 as assumed throughout this article. Finally, we have that

∥ 𝑓 𝑔∥ ℓ 1 ⩽ ∥𝑓 ∥ ℓ 1 ∥𝑔∥ ℓ 1 .

Rigorous trigonometric approximations

By analogy with intervals used to represent real numbers rigorously on computers [START_REF] Moore | Interval analysis[END_REF]], (generalized) rigorous polynomial approximations are used in rigorous numerics to enclose mathematical functions in set-valued representations. Taking its roots in the so-called ultra-arithmetics developed in the early 80s [Epstein et al. 1982a,b], this concept has been extended and implemented in various settings, notably with Taylor [Makino and Berz 2003a,b] or Chebyshev [START_REF] Brisebarre | Chebyshev interpolation polynomial-based tools for rigorous computing[END_REF] approximations. Since the functions considered throughout this article are univariate, real-valued and periodic, the trigonometric polynomials expressed in the Fourier trigonometric basis are the natural candidates to build rigorous, set-valued representations. For short, we call them Rigorous Trigonometric Approximations (RTA). Some previous works already made use of quite similar concepts (see, e.g., [START_REF] Figueras | Numerical computations and computer assisted proofs of periodic orbits of the Kuramoto-Sivashinsky equation[END_REF][START_REF] Hungria | Rigorous numerics for analytic solutions of differential equations: the radii polynomial approach[END_REF]). However, for the sake of completeness, and in order to provide sound bases for the complexity analysis of the computation method for Abelian integrals presented in this article, we provide some elementary definitions, basic routines and lemmas concerning RTAs.

Definition 2.2 (Rigorous Trigonometric Approximations).

A degree 𝑁 Rigorous Trigonometric Approximation (RTA) is a pair 𝒇 = ( f , 𝜀), with f a degree 𝑁 TP and 𝜀 ⩾ 0 representing the closed 𝜀-ball around f in the Banach space 𝐶 0 2𝜋 :

𝐵( f , 𝜀) = 𝑓 ∈ 𝐶 0 2𝜋 ∥𝑓 -f ∥ ∞ ⩽ 𝜀 .
Therefore, we say that a degree 𝑁 RTA 𝒇 = ( f , 𝜀) represents a function 𝑓 ∈ 𝐶 0 2𝜋 (or 𝒇 contains 𝑓 , or simply 𝑓 ∈ 𝒇 ) if there exists a function 𝑠 ∈ 𝐶 0 2𝜋 with ∥𝑠 ∥ ∞ ⩽ 𝜀 such that 𝑓 (𝑡) = f (𝑡) + 𝑠 (𝑡) for all 𝑡 ∈ R.

Linear operations are trivially defined on RTAs to match the (Banach) linear space structure of (𝐶 0 2𝜋 , ∥ • ∥ ∞ ): for all degree 𝑁 RTAs ( f , 𝜀), ( ǧ, 𝜂) and 𝜆 ∈ R, we have

( f , 𝜀) + ( ǧ, 𝜂) = ( f + ǧ, 𝜀 + 𝜂), ( f , 𝜀) -( ǧ, 𝜂) = ( f -ǧ, 𝜀 + 𝜂), 𝜆( f , 𝜀) = (𝜆 f , |𝜆|𝜀).
Also, Bound(𝒇 ), defined as:

Bound(( f , 𝜀)) = 𝑁 ∑ 𝑘=0 |𝑎 𝑘 | + 𝑁 ∑ 𝑘=1 |𝑏 𝑘 | + 𝜀.
7This is a straightforward corollary of the exponential decay of Fourier coefficients for analytic functions, see e.g. [Wright et al. 2015, Thm. 4.1].

satisfies Bound(𝒇

) ⩾ ∥ f ∥ ℓ 1 + 𝜀 ⩾ ∥ f ∥ ∞ + 𝜀, so that ∥ 𝑓 ∥ ∞ ⩽ Bound(𝒇 ) for all 𝑓 ∈ 𝒇 .
The product of two degree 𝑁 RTAs is the degree 2𝑁 RTA defined by

( f , 𝜀) × ( ǧ, 𝜂) = f ǧ, Bound( f )𝜂 + Bound( ǧ)𝜀 + 𝜀𝜂 ,
which can be computed in 𝑂 (𝑁 log(𝑁 )) arithmetic operations using the fast multiplication algorithm for TPs mentioned previously. We shall also define the integral of 𝒇 = ( f , 𝜀) with f = 𝑎 0 + 𝑁 𝑘=1 (𝑎 𝑘 cos(𝑘𝑡) + 𝑏 𝑘 sin(𝑘𝑡)) over one period:

∫ 2𝜋 0 𝒇 (𝑡)d𝑡 = 2𝜋 (𝑎 0 + [-𝜀, 𝜀]) .
Lemma 2.3. If the RTAs 𝒇 = ( f , 𝜀) and 𝒈 = ( ǧ, 𝜂) represent functions 𝑓 and 𝑔, then the RTA 𝜆𝒇 represents 𝜆𝑓 for any 𝜆 ∈ R, and the RTA 𝒇 ★ 𝒈 represents 𝑓 ★ 𝑔 for ★ ∈ {+, -, ×}. Also,

∫ 2𝜋 0 𝑓 (𝑡)d𝑡 ∈ ∫ 2𝜋 0 𝒇 (𝑡)d𝑡 for all 𝑓 ∈ 𝒇 .
Proof. We prove the least trivial case only, namely the multiplication. By hypothesis, there exist functions 𝑢 and 𝑣 in 𝐶 0 2𝜋 with ∥𝑢 ∥ ∞ ⩽ 𝜀 and ∥𝑣 ∥ ∞ ⩽ 𝜂, such that 𝑓 = f + 𝑢 and 𝑔 = ǧ + 𝑣. Therefore,

∥ 𝑓 𝑔 -f ǧ∥ ∞ ⩽ ∥ f 𝑣 + ǧ𝑢 + 𝑢𝑣 ∥ ∞ ⩽ ∥ f ∥ ∞ ∥𝑣 ∥ ∞ + ∥ ǧ∥ ∞ ∥𝑢 ∥ ∞ + ∥𝑢 ∥ ∞ ∥𝑣 ∥ ∞ ⩽ Bound( f )𝜂 + Bound( ǧ)𝜀 + 𝜀𝜂,
which by definition implies 𝑓 𝑔 ∈ 𝒇 × 𝒈. □ Remark 2.4. By analogy with a globally set precision for floating-point and interval arithmetic, we can choose to fix a global degree 𝑁 used for TP / RTAs computations. Since some operations (e.g. multiplication) may increase the degree, the obtained result is truncated. In the case of RTAs, the higher-order terms must be added to the error component:

Trunc(𝒇, 𝑁 ) = 𝑎 0 + 𝑁 ∑ 𝑘=1 𝑎 𝑘 cos(𝑘𝑡) + 𝑏 𝑘 sin(𝑘𝑡) , 𝑚 ∑ 𝑘=𝑁 +1 |𝑎 𝑘 | + |𝑏 𝑘 | + 𝜀 , for 𝒇 = 𝑎 0 + 𝑚 𝑘=1 𝑎 𝑘 cos(𝑘𝑡) +𝑏 𝑘 sin(𝑘𝑡) , 𝜀 when 𝑚 > 𝑁 . Clearly, if 𝑓 ∈ 𝒇 , then 𝑓 ∈ Trunc(𝒇, 𝑁 ).
These elementary operations on RTAs are the building blocks of more complex ones. For instance, if 𝑝 (𝑥) is a polynomial and 𝒇 an RTA for 𝑓 , then substituting 𝑥 by 𝒇 and applying an evaluation scheme of choice with arithmetic operations on RTAs gives an RTA 𝑝 (𝒇 ) for 𝑝 (𝑓 ).

Finally, we also need a routine IsPositive to determine whether, for an RTA 𝒇 , 𝑓 > 0 holds for every 𝑓 ∈ 𝒇 . Such a routine is discussed in Appendix A.

A posteriori validation using Newton-like operators

Many operations defined on functions do not preserve the ring of trigonometric polynomials. Examples include division, the square root, or the more complex operation performed by Algorithm OvalValid in Section 4 used to rigorously approximate a parameterization of the oval Γ(ℎ). More generally, the target function 𝜑 * can be viewed as the solution of a functional equation F (𝜑) = 0 for a suitable operator F mapping a function space to another one. For such cases, the principle of a posteriori validation is widely used in validated numerics:

• First, a TP φ approximating the solution 𝜑 * is computed using nonrigorous algorithms, e.g. by trigonometric interpolation as in Section 2.1.

• Next, an upper bound on the approximation error 𝜀 ⩾ ∥ φ -𝜑 * ∥ ∞ is rigorously computed by expressing 𝜑 * as a fixed point of a contracting operator T , and applying the Banach fixed-point theorem (see Theorem 2.5 below). The result is an RTA 𝝋 = ( φ, 𝜀) such that 𝜑 * ∈ 𝝋.

Theorem 2.5. (see, e.g., [Yamamoto 1998, Thm. 2.1]) Let (𝐸, ∥ • ∥) be a Banach space, φ ∈ 𝐸, and T : 𝐸 → 𝐸. If we can find 𝑟, 𝜆 ⩾ 0 such that:

• T is 𝜆-contracting over the closed ball 𝐵( φ, 𝑟 ), i.e., 𝜆 < 1 and:

∥T (𝜑 1 ) -T (𝜑 2 ) ∥ ⩽ 𝜆∥𝜑 1 -𝜑 2 ∥ for all 𝜑 1 , 𝜑 2 ∈ 𝐵( φ, 𝑟 ); • We have 𝑑 + 𝜆𝑟 ⩽ 𝑟, where 𝑑 := ∥T ( φ) -φ ∥ (called the defect);
Then T admits a unique fixed point 𝜑 * in the ball 𝐵( φ, 𝑟 ), and we have the following enclosure for the approximation error:

𝑑 1 + 𝜆 ⩽ ∥ φ -𝜑 * ∥ ⩽ 𝑑 1 -𝜆 .
Assume the functional equation is of the form F (𝜑) = 0 with F : 𝐶 0 2𝜋 → 𝐶 0 2𝜋 being of class 𝐶 1 , and let φ be an approximate zero of F . A classical method in rigorous numerics (see e.g., [START_REF] Hungria | Rigorous numerics for analytic solutions of differential equations: the radii polynomial approach[END_REF][START_REF] Yamamoto | A numerical verification method for solutions of boundary value problems with local uniqueness by Banach's fixed-point theorem[END_REF]) to obtain an equivalent fixed-point equation T (𝜑) = 𝜑 is to construct T : 𝐶 0 2𝜋 → 𝐶 0 2𝜋 as a Newton-like operator:

T (𝜑) = 𝜑 -A (F (𝜑)).
Here, A : 𝐶 0 2𝜋 → 𝐶 0 2𝜋 is an invertible bounded linear operator approximating (DF ( φ)) -1 , the inverse of the differential of F at φ, making T a local contraction in a neighborhood of φ.

Newton-like validation in the case of polynomial equations.

For the purpose of validating approximate parameterizations of planar algebraic curves, we apply the above strategy to polynomial equations of the form

𝑐 𝑟 (𝑡)𝜑 (𝑡) 𝑟 + 𝑐 𝑟 -1 (𝑡)𝜑 (𝑡) 𝑟 -1 + • • • + 𝑐 1 (𝑡)𝜑 (𝑡) + 𝑐 0 (𝑡) = 0, 𝑡 ∈ [0, 2𝜋],
with 𝑐 0 , . . . , 𝑐 𝑟 ∈ 𝐶 0 2𝜋 , that is to say, F ∈ 𝐶 0 2𝜋 [𝑋 ] is a polynomial with coefficients in 𝐶 0 2𝜋 . In this case, the differential DF ( φ) coincides with the multiplication by F ′ ( φ) ∈ 𝐶 0 2𝜋 , where the prime symbol denotes the usual differentiation in a ring of polynomials, here 𝐶 0 2𝜋 [𝑋 ]. Therefore, the linear operator A is set to be the multiplication by a TP ǎ constructed by interpolation to approximate 1/F ′ ( φ):

A (𝜑) = ǎ𝜑, with ǎ(𝑡) ≈ 1 F ′ ( φ) (𝑡)
.

The following lemma provides a simple way to bound the Lipschitz constant of the resulting Newton-like operator T .

Lemma 2.6.

If F ∈ 𝐶 0 2𝜋 [𝑋 ], A (𝜑) (𝑡) = ǎ(𝑡)𝜑 (𝑡), and
T (𝜑) = 𝜑 -A (F (𝜑)), then T is 𝜆(𝑟 )- Lipschitz over 𝐵( φ, 𝑟 ) with 𝜆(𝑟 ) = 𝜆 0 + 𝛼𝑟𝜆 1 (𝑟 )
, obtained from the (rigorously computed) bounds:

𝜆 0 ⩾ ∥1 -ǎ F ′ ( φ) ∥ ∞ , 𝛼 ⩾ ∥ ǎ∥ ∞ , 𝜆 1 (𝑟 ) ⩾ sup 𝜑 ∈𝐵 ( φ,𝑟 ) ∥F ′′ (𝜑) ∥ ∞ .
Proof. To get an upper bound for the Lipschitz constant for

T ∈ 𝐶 0 2𝜋 [𝑋 ] over 𝐵( φ, 𝑟 ), we simply bound ∥T ′ (𝜑) ∥ ∞ for ∥𝜑 -φ ∥ ∞ ⩽ 𝑟 by the triangle inequality. First, ∥T ′ ( φ) ∥ ∞ = ∥1-ǎ F ′ ( φ) ∥ ∞ ⩽ 𝜆 0 by definition of 𝜆 0 . Then, ∥T ′ (𝜑) -T ′ ( φ) ∥ ∞ = ∥ ǎ (F ′ (𝜑) -F ′ ( φ)) ∥ ∞ ⩽ 𝛼𝑟𝜆 1 (𝑟 )
using the mean value theorem, by definition of 𝛼 and 𝜆 1 , and since 𝜑 ∈ 𝐵( φ, 𝑟 ). Adding both bounds, we obtain the expected upper bound for the Lipschitz constant 𝜆(𝑟 ). □

Note that A being invertible is a byproduct of the second assumption of Theorem 2.5: except in the trivial case 𝑑 = 0, the inequality implies 𝜆 0 < 1, so ǎ never vanishes.

Example 2.7 (Division of RTAs). For two functions 𝑔, ℎ ∈ 𝐶 0 2𝜋 with ℎ(𝑡) ≠ 0 for all 𝑡, represented by RTAs 𝒈 and 𝒉, Algorithm RTADiv(𝒈, 𝒉, 𝑁 ) computes a degree 𝑁 TP φ approximating the quotient 𝜑 * = 𝑔/ℎ by trigonometric interpolation, and validates it using Lemma 2.6, since 𝜑 * is the unique zero of F ∈ 𝐶 0 2𝜋 [𝑋 ] defined by F (𝜑) = ℎ𝜑 -𝑔. When a global degree 𝑁 for the RTA computations is fixed, we may use the notation 𝒈/𝒉 for RTADiv(𝒈, 𝒉, 𝑁 ).

Algorithm 1 RTADiv(𝒈, 𝒉, 𝑁 ) Input: RTAs 𝒈 = ( ǧ, 𝛿) and 𝒉 = ( ȟ, 𝜂), and approximation degree 𝑁 ∈ N Output: Degree 𝑁 RTA 𝝋 representing the quotient 𝒈/𝒉 ⊲ Build degree 𝑁 candidate TP φ by interpolation

1: (𝑔 𝑗 ) 2𝑁 𝑗=0 ← IFFT 𝑁 ( ǧ) and (ℎ 𝑗 ) 2𝑁 𝑗=0 ← IFFT 𝑁 ( ȟ) 2: φ ← FFT 𝑁 𝑔 𝑗 ℎ 𝑗 2𝑁 𝑗=0
⊲ Build TP ǎ defining the Newton-like operator and bound the Lipschitz constant

3: ǎ ← FFT 𝑁 1 ℎ 𝑗 2𝑁 𝑗=0 4: 𝜆 ← Bound(1 -ǎ𝒉)
⊲ Compute a posteriori error bound, if possible

5: if 𝜆 < 1 then 6: 𝑑 ← Bound ǎ (𝒉 φ -𝒈) 7: 𝜀 ← 𝑑 1-𝜆 8:
return 𝝋 = ( φ, 𝜀) 9: else 10:

return 𝐹𝐴𝐼𝐿 11: end if

We now prove that the division operator, as implemented in Algorithm1, is correct, convergent, and has a good complexity.

Lemma 2.8. Let 𝒈 = ( ǧ, 𝛿) and 𝒉 = ( ȟ, 𝜂) be RTAs. Then, (𝑖) For all 𝑁 ∈ N, if RTADiv(𝒈, 𝒉, 𝑁 ) returns an RTA 𝝋, then for all 𝑔 ∈ 𝒈 and ℎ ∈ 𝒉, ℎ(𝑡) ≠ 0 for all 𝑡 and 𝑔/ℎ ∈ 𝝋.

(𝑖𝑖) If ȟ(𝑡) ≠ 0 for all 𝑡 and 𝜂 < ∥ ȟ-1 ∥ -1
ℓ 1 , then there exists an 𝑁 0 ∈ N such that for all 𝑁 ⩾ 𝑁 0 , Algorithm RTADiv(𝒈, 𝒉, 𝑁 ) does not fail, and the remainder of 𝝋 = ( φ, 𝜀) satisfies

𝜀 ⩽ ∥ ȟ-1 ∥ ℓ 1 𝛿 + ∥ ǧ/ ȟ2 ∥ ℓ 1 𝜂 1 -∥ ȟ-1 ∥ ℓ 1 𝜂 + 𝑂 (𝜅 -𝑁 ), as 𝑁 → ∞,
for some 𝜅 > 1 depending on ǧ and ȟ. (𝑖𝑖𝑖) RTADiv(𝒈, 𝒉, 𝑁 ) runs in 𝑂 (𝑁 ′ log(𝑁 ′ )) arithmetic operations, where

𝑁 ′ = max(deg 𝒈, deg 𝒉, 𝑁 ).
Proof. (𝑖) Soundness. The Newton-like operator associated to F is defined by T (𝜑) = 𝜑ǎ (ℎ𝜑 -𝑔) with ǎ a degree 𝑁 TP approximating 1/F ′ ( φ) = 1/ℎ(𝑡) (line 3). Since F (and hence T ) has degree 1, F ′′ = 0. Hence the bound 𝜆 1 (𝑟 ) in Lemma 2.6 is identically 0 , and we have 𝜆(𝑟 ) = 𝜆 0 = ∥1 -ǎℎ∥ ∞ , rigorously upper-bounded by 𝜆 in line 4. Hence, as soon as 𝜆 < 1, we can set 𝑟 arbitrary large in Theorem 2.5 and return the bound implemented in lines 6 and 7:

∥ φ -𝜑 * ∥ ∞ ⩽ ∥T ( φ) -φ ∥ ∞ 1 -𝜆 ⩽ ∥ ǎ (ℎ φ -𝑔) ∥ ∞ 1 -𝜆 ⩽ Bound ǎ (𝒉 φ -𝒈) 1 -𝜆 .
Therefore, the degree 𝑁 RTA 𝝋 = ( φ, 𝜀) returned in line 8 contains 𝑔/ℎ for all 𝑔 ∈ 𝒈 and ℎ ∈ 𝒉. Note moreover that the nonvanishing of all ℎ ∈ 𝒉 is a byproduct of the inequality ∥1 -ǎℎ∥ ∞ ⩽ 𝜆 < 1.

(𝑖𝑖) Convergence. Suppose that ȟ(𝑡) ≠ 0 for all 𝑡, and that [Wright et al. 2015, Thm. 4.2]8).

𝜂 < ∥ ȟ-1 ∥ -1 ℓ 1 . In particular, 𝜂 < ∥ ȟ-1 ∥ -1 ∞ = min 𝑡 ∈R | ȟ(𝑡)|. Since ȟ-1 (resp. ǧ/ ȟ) is analytic over R, then ǎ = I 𝑁 ( ȟ-1 ) (resp. φ = I 𝑁 ( ǧ/ ȟ)) converges to ȟ-1 (resp. ǧ/ ȟ) w.r.t. ∥ • ∥ ℓ 1 exponentially fast as 𝑁 → ∞ (see e.g.
Now for the bound 𝜆 computed in line 4, we have for all ℎ ∈ 𝒉,

∥1 -ǎℎ∥ ℓ 1 ⩽ ∥1 -ǎ ȟ∥ ℓ 1 + ∥ ǎ∥ ℓ 1 ∥ℎ -ȟ∥ ℓ 1 ⩽ ∥ ǎ -ȟ-1 ∥ ℓ 1 ∥ ȟ∥ ℓ 1 + ∥ ǎ∥ ℓ 1 𝜂 ⩽ ∥ ȟ-1 ∥ ℓ 1 𝜂 + 𝑂 (𝜅 -𝑁 ),
for some 𝜅 > 1, as 𝑁 → ∞. In particular, Algorithm RTADiv always terminates for sufficiently large 𝑁 . Similarly, for the bound 𝑑 computed in line 6, for all 𝑔 ∈ 𝒈 and ℎ ∈ 𝒉,

∥ ǎ(ℎ φ -𝑔) ∥ ℓ 1 ⩽ ∥ ǎ( ȟ φ -ǧ) ∥ ℓ 1 → 0 + ∥ ǎ∥ ℓ 1 → ∥ ȟ-1 ∥ ℓ 1 𝛿 + ∥ ǎ φ ∥ ℓ 1 → ∥ ǧ/ ȟ2 ∥ ℓ 1 𝜂 ⩽ ∥ ȟ-1 ∥ ℓ 1 𝛿 + ∥ ǧ/ ȟ2 ∥ ℓ 1 𝜂 + 𝑂 (𝜅 -𝑁 ),
for some 𝜅 > 1, as 𝑁 → ∞. This finally yields the expected estimate for 𝜀 computed in line 7.

(𝑖𝑖𝑖) Complexity. The asymptotic complexity is determined by the FFT 𝑁 /IFFT 𝑁 routines used for evaluation on the equispaced grid, interpolation and TP/RTA multiplications. □ Remark 2.9. In the case where the polynomial operator F has degree greater than 1 (which will be the case for Algorithm OvalValid in Section 4 as long as the potential function 𝐻 (𝑥, 𝑦) has degree at least 2), an explicit radius 𝑟 satisfying the hypothesis of Theorem 2.5 must be rigorously computed. Given the function 𝑟 ↦ → 𝜆(𝑟 ), this amounts to locating the smallest positive zero of (1 -𝜆(𝑟 ))𝑟 -𝑑. A simple bisection method NewtonBall is proposed in Appendix C for this purpose.

APPROXIMATION OF THE OVAL Γ(ℎ) AND NUMERICAL INTEGRATION

This section is devoted to the computation of a smooth approximate parameterization 𝑡 ∈ [0, 2𝜋] → ( x (𝑡), y (𝑡)) of the oval Γ(ℎ) with x, y TPs of degree 𝑁 .

First, an implicit parameterization for Γ(ℎ) is given by the unscaled, unperturbed Hamiltonian system itself (i.e., System (1.2) with 𝜀 = 0):

         𝑥 = - 𝜕𝐻 𝜕𝑦 (𝑥, 𝑦), 𝑦 = 𝜕𝐻 𝜕𝑥 (𝑥, 𝑦). (3.1)
For a level ℎ on the portion of Σ of interest, the trajectory of (3.1) starting at (𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 ) = Σ ∩ 𝐻 -1 (ℎ) is periodic since it exactly follows the oval Γ(ℎ). Call 𝑡 end the time of first return onto Σ. By a linear change of the independent variable 𝑡, we may assume that 𝑡 end = 2𝜋, giving a periodic analytic parameterization 𝑡 ∈ [0, 2𝜋] ↦ → (𝑥 * (𝑡), 𝑦 * (𝑡)) of Γ(ℎ), (3.2) with a winding number equal to ±1 with respect to any point inside Γ(ℎ). The approximation of Γ(ℎ) consists of two steps:

• First, Algorithm OvalApproxInit presented in Section 3.1 computes an initial approximation (𝑥 • , 𝑦 • ) ∈ (𝐶 0 2𝜋 ) 2 of (𝑥 * , 𝑦 * ), with TPs 𝑥 • , 𝑦 • of a low degree 𝑁 0 . The accuracy of this initial guess (and therefore 𝑁 0 ) should ensure the convergence of the Newton iterations in the next routine OvalApproxRefine. It only depends on the system under consideration (i.e., 𝐻 , Σ and ℎ) and not on the (possibly very high) final precision that we target. Hence, the uniform convergence of the output of Algorithm OvalApproxInit (without asymptotic estimate) is sufficient for our complexity analysis.

• Next, Algorithm OvalApproxRefine refines this initial guess by computing TPs x, y of higher degree 𝑁 . As explained in Section 3.2, this is done by applying Newton iterations to the initial guess so as to "project" it onto Γ(ℎ). We provide a complexity analysis and an asymptotic convergence estimate of ( x, y) to Γ(ℎ) in terms of the degree 𝑁 . We emphasize the fact that the convergence results in this section are asymptotic only. Indeed, our objective is to prove the quasi-linear arithmetic complexity of this approximation procedure with respect to the number of correct digits. On the other hand, computing effective error bounds is the role of the algorithms presented in the Sections 4 and 5.

Algorithm OvalApproxInit: initial approximation of the curve

Algorithm OvalApproxInit is given an initial point (𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 ) ∈ Γ(ℎ) ∩ Σ on the transversal. It first calls a numerical ODE solver to approximately solve the differential system (3.1) until the first return onto Σ. It provides us with a continuous function 𝑡 ↦ → ( x (𝑡), ỹ (𝑡)), which we assume to be rescaled over the time interval [0, 2𝜋]. Finally, in order to obtain smooth approximations to (𝑥 * , 𝑦 * ), cf. (3.2), we interpolate ( x, ỹ) at evenly spaced points and obtain low degree TPs 𝑥 • , 𝑦 • (lines 2-5) representing our analytic initial guess for Γ(ℎ). The proof of the following Proposition is postponed to the end of this section. Proposition 3.1. Under Assumption 3.2, Algorithm OvalApproxInit(𝐻, Σ, 𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 , 𝑁 0 ) computes degree 𝑁 0 TPs 𝑥 • and 𝑦 • that converge uniformly to the rescaled parameterization

(𝑥 * , 𝑦 * ) ∈ (𝐶 0 2𝜋 ) 2 of Γ(ℎ): ∥𝑥 • -𝑥 * ∥ ∞ , ∥𝑦 • -𝑦 * ∥ ∞ → 0 as 𝑁 0 → ∞.
Algorithm 2 OvalApproxInit(𝐻, Σ, 𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 , 𝑁 0 )

Input: Potential function 𝐻 , transversal Σ, initial point (𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 ), and approximation degree 3.1.1 The numerical solver. We need a numerical routine DSolve(𝐻, Σ, 𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 , 𝛿) that integrates the ODE system (3.1) with some precision parameter 𝛿, starting from the initial point (𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 ) ∈ Γ(ℎ) ∩ Σ. This routine has to detect numerically the first return onto Σ and stop at the corresponding time 𝑡 end . Moreover, we assume that this routine rescales the time interval [0, 𝑡 end ] to [0, 2𝜋], in order to work on a fixed domain of definition. Most scientific programming languages offer libraries with sophisticated explicit or implicit iterative schemes to compute a time-discretized solution to (3.1), with excellent timings for the moderate precision, e.g., the standard double precision in most cases, needed at this stage. Here, the precision parameter 𝛿 may represent a maximal time step. Moreover, the values between two consecutive time steps are usually automatically interpolated, so as to define a continuous function 𝑡 ↦ → ( x (𝑡), ỹ (𝑡)).

𝑁 0 ∈ N Output: Degree 𝑁 0 TPs 𝑥 • , 𝑦 • approximating Γ(ℎ) ⊲ Approximation of Γ(ℎ), computed with precision parameter 𝑁 -1 0 1: ( x, ỹ) ← DSolve(𝐻, Σ, 𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 , 𝑁 -1 0 ) 2: for 𝑗 = 0 to 2𝑁 0 do 3: 𝑥 • 𝑗 ← x 2𝑗𝜋 2𝑁 0 + 1 , 𝑦 • 𝑗 ← ỹ 2𝑗𝜋 2𝑁 0 + 1 4: end for ⊲ Interpolate the points 5: 𝑥 • ← FFT 𝑁 0 (𝑥 • 𝑗 ) 2𝑁 0 𝑗=0 , 𝑦 • ← FFT 𝑁 0 (𝑦 • 𝑗 ) 2𝑁 0

Detection of the first return onto Σ.

Another feature often proposed by numerical libraries is the so-called event detection: the solver stops when a user-defined continuous function of the dependent and independent variables vanishes (the "event"). In our case, intersecting the transversal line Σ is trivially recast as a linear equation in the Cartesian coordinates. Note, however, that the event is already satisfied at 𝑡 = 0 (up to numerical errors), since (𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 ) ∈ Σ, which may cause the solver to stop before detecting the actual first return. A more robust solution is to use a polar parameterization with respect to a reference point (𝑥 𝑟 , 𝑦 𝑟 ) of choice inside Γ(ℎ), that is the ODE system:

𝑥 (𝑡) = 𝑥 𝑟 + 𝜌 (𝑡) cos 𝜃 (𝑡), 𝑦 (𝑡) = 𝑦 𝑟 + 𝜌 (𝑡) sin 𝜃 (𝑡), with              𝜌 = sin 𝜃 𝜕𝐻 𝜕𝑥 (𝑥, 𝑦) -cos 𝜃 𝜕𝐻 𝜕𝑦 (𝑥, 𝑦), 𝜃 = cos 𝜃 𝜕𝐻 𝜕𝑥 (𝑥, 𝑦) + sin 𝜃 𝜕𝐻 𝜕𝑦 (𝑥, 𝑦) 𝜌 .
The stopping condition is when 𝜃 (𝑡) reaches 𝜃 (0) + 𝜍2𝜋, with 𝜍 = ±1 the winding number of (𝑥 * , 𝑦 * ) with respect to (𝑥 𝑟 , 𝑦 𝑟 ).

Convergence of the ODE solver and a proof of Proposition 3.1.

To prove Proposition 3.1, we need the following assumption on the numerical routine DSolve to ensure that the numerical solution ( x, ỹ), computed with a precision parameter 𝛿 = 1/𝑁 0 , converges fast enough to (𝑥 * , 𝑦 * ) to balance the overestimation factor due to the degree 𝑁 0 trigonometric interpolation. In practice, such a mild convergence assumption in 𝑂 (𝛿) is guaranteed by most iterative schemes, even a simple Euler scheme with time step 𝛿.

Assumption 3.2. The numerical routine DSolve(𝐻, Σ, 𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 , 𝛿) computes a continuous approx- imation 𝑡 ∈ [0, 2𝜋] ↦ → ( x (𝑡), ỹ (𝑡)) of the rescaled parameterization (𝑥 * , 𝑦 * ), cf. (3.2), of Γ(ℎ), such that: ∥ x -𝑥 * ∥ ∞ , ∥ ỹ -𝑦 * ∥ ∞ = 𝑂 (𝛿) as 𝛿 → 0.
We can now prove the convergence of Algorithm OvalApproxInit.

Proof of Proposition 3.1. Let the symbol 𝑧 denote either 𝑥 or 𝑦. Let 𝜂 = z -𝑧 * be the error between 𝑧 * and the approximate solution z computed by DSolve in line 1. Thanks to the linearity of the operator I 𝑁 0 (see (2.1)), the error between the interpolant 𝑧 • computed in lines 2-5 and 𝑧 * can be bounded from above by the error of interpolation for 𝑧 * , plus the interpolation of the error 𝜂 :

∥𝑧 • -𝑧 * ∥ ∞ = ∥I 𝑁 0 ( z) -𝑧 * ∥ ∞ ⩽ ∥I 𝑁 0 (𝑧 * ) -𝑧 * ∥ ∞ + ∥I 𝑁 0 (𝜂) ∥ ∞ .
• Error of interpolation. Since 𝑧 * is 2𝜋-periodic and continuously differentiable, the trigonometric interpolants I 𝑁 0 (𝑧 * ) converge uniformly to 𝑧 * [Zygmund 2002, Chap. X].

• Interpolation of error. By Assumption 3.2, 𝜂 is a continuous function over [0, 2𝜋] and ∥𝜂 ∥ ∞ converges uniformly to 0 in 𝑂 (𝑁 -1 0 ) as the precision parameter 𝛿 = 𝑁 -1 0 tends to 0. From Section 2.1, we have

∥I 𝑁 0 (𝜂) ∥ ∞ ⩽ Λ 𝑁 0 ∥𝜂 ∥ ∞ = 𝑂 (log(𝑁 0 )/𝑁 0 ) → 0 as 𝑁 0 → ∞. □

Algorithm OvalApproxRefine: refining the curve by projection

As its name suggests, the role of Algorithm OvalApproxRefine is to refine the initial guess (𝑥 • , 𝑦 • ) with TPs ( x, y) of higher degree 𝑁 ⩾ 𝑁 0 , yielding an analytic approximate parameterization of Γ(ℎ) of very high accuracy.

By solving algebraic equations only, it is in general not possible to recover the parameterization (𝑥 * , 𝑦 * ) ∈ (𝐶 0 2𝜋 ) 2 satisfying the differential system (3.1). However, the value of ℑ(ℎ) is independent of the parameterization used for Γ(ℎ). Therefore, we define another one, (𝑥 # , 𝑦 # ) ∈ (𝐶 0 2𝜋 ) 2 , called the projected parameterization, which Algorithm OvalApproxRefine approximates using Newton's method on points sampled on the initial curve (𝑥 • , 𝑦 • ).

Definition 3.3. Let (𝑥 • , 𝑦 • ) ∈ (𝐶 0 2𝜋 ) 2 be an approximate parameterization of Γ(ℎ), (𝑢 • , 𝑣 • ) ∈ (𝐶 0 2𝜋 ) 2 , and 𝜀 > 0. If the equation 𝐻 (𝑥 • (𝑡) + 𝑠 (𝑡) 𝑢 • (𝑡), 𝑦 • (𝑡) + 𝑠 (𝑡) 𝑣 • (𝑡)) = ℎ, 𝑡 ∈ [0, 2𝜋], (3.3) has a unique solution 𝑠 # ∈ 𝐶 0 2𝜋 with ∥𝑠 # ∥ ∞ ⩽ 𝜀, called the shift, then the curve (𝑥 # , 𝑦 # ) ∈ (𝐶 0 2𝜋 ) 2 , 𝑥 # (𝑡) = 𝑥 • (𝑡) + 𝑠 # (𝑡) 𝑢 • (𝑡), 𝑦 # (𝑡) = 𝑦 • (𝑡) + 𝑠 # (𝑡) 𝑣 • (𝑡), (3.4)
is called the projection of (𝑥 • , 𝑦 • ) onto Γ(ℎ) w.r.t. the direction (𝑢 • , 𝑣 • ).

The following Lemma, whose proof is postponed to Appendix B, establishes the good properties of the projected parameterization when the initial guess (𝑥 • , 𝑦 • ) returned by OvalApproxInit is sufficiently close to Γ(ℎ).

Lemma 3.4. There exists 𝜂 > 0 (depending on 𝐻 and ℎ) such that if (𝑥 • , 𝑦 • ) and (𝑢 • , 𝑣 • ) are analytic and 𝜂-close9 to (𝑥 * , 𝑦 * ) and

∇𝐻 (𝑥 * ,𝑦 * ) ∥ ∇𝐻 (𝑥 * ,𝑦 * ) ∥ 2 , respectively, in (𝐶 0 2𝜋 ) 2 , then: (𝑖) (𝑥 # , 𝑦 # ) ∈ (𝐶 0 2𝜋
) 2 is well-defined and it is an analytic parameterization of Γ(ℎ) with same orientation as (𝑥 * , 𝑦 * ). (𝑖𝑖) For all 𝑡 ∈ [0, 2𝜋], letting (𝑥, 𝑦) = (𝑥 • (𝑡), 𝑦 • (𝑡)) and (𝑢, 𝑣) = (𝑢 • (𝑡), 𝑣 • (𝑡)), the Newton operator associated to Equation (3.3) at 𝑡 is well-defined:

N 𝑥,𝑦,𝑢,𝑣 (𝑠) = 𝑠 - 𝐻 (𝑥 + 𝑠𝑢, 𝑦 + 𝑠𝑣) -ℎ 𝑢 𝜕𝐻 𝜕𝑥 (𝑥 + 𝑠𝑢, 𝑦 + 𝑠𝑣) + 𝑣 𝜕𝐻 𝜕𝑦 (𝑥 + 𝑠𝑢, 𝑦 + 𝑠𝑣) , (3.5)
and the Newton iterations 𝑠 (𝑘 ) = N 𝑘 𝑥,𝑦,𝑢,𝑣 (0) converge quadratically fast to 𝑠 # (𝑡). Based on those properties, Algorithm OvalApproxRefine discretizes the initial guess (𝑥 • , 𝑦 • ) over the equispaced grid (𝑡 𝑗 ) of size 2𝑁 + 1 (line 1), then applies several Newton iterations (3.5), implemented in lines 6-8, and finally interpolates the obtained values 𝑠 𝑗 in order to construct a degree 𝑁 TP š approximating the shift 𝑠 # . This defines the approximate parameterization ( x, y) of Γ(ℎ):

x

(𝑡) = 𝑥 • (𝑡) + š (𝑡) 𝑢 • (𝑡), y (𝑡) = 𝑦 • (𝑡) + š (𝑡) 𝑣 • (𝑡), 𝑡 ∈ [0, 2𝜋]. (3.6) 9It means that max( ∥𝑥 • -𝑥 * ∥ ∞ , ∥𝑦 • -𝑦 * ∥ ∞ ) ⩽ 𝜂.
, Vol. 1, No. 1, Article . Publication date: December 2023.

Algorithm 3 OvalApproxRefine(𝐻, ℎ, 𝑥

⊲ Sample initial points 1: (𝑥 𝑗 ) 2𝑁 𝑗=0 ← IFFT 𝑁 (𝑥 • ), (𝑦 𝑗 ) 2𝑁 𝑗=0 ← IFFT 𝑁 (𝑦 • )
⊲ Apply ⌈log 𝑁 ⌉ Newton iterations 2: for 𝑗 = 0 to 2𝑁 do 3:

(𝑢 𝑗 , 𝑣 𝑗 ) ← ∇𝐻 (𝑥 𝑗 ,𝑦 𝑗 ) ∥ ∇𝐻 (𝑥 𝑗 ,𝑦 𝑗 ) ∥ 2 4: 𝑠 𝑗 ← 0 5: for 𝑘 = 1 to ⌈log 𝑁 ⌉ do 6: 𝛿 ← ℎ -𝐻 (𝑥 𝑗 , 𝑦 𝑗 ) 𝑢 𝑗 𝜕𝐻 𝜕𝑥 (𝑥 𝑗 , 𝑦 𝑗 ) + 𝑣 𝑗 𝜕𝐻 𝜕𝑦 (𝑥 𝑗 , 𝑦 𝑗 ) 7: 𝑠 𝑗 ← 𝑠 𝑗 + 𝛿 8: (𝑥 𝑗 , 𝑦 𝑗 ) ← (𝑥 𝑗 , 𝑦 𝑗 ) + 𝛿 (𝑢 𝑗 , 𝑣 𝑗 ) 9:
end for 10: end for

⊲ Interpolate 11: 𝑢 • ← FFT 𝑁 (𝑢 𝑗 ) 2𝑁 𝑗=0 , 𝑣 • ← FFT 𝑁 (𝑣 𝑗 ) 2𝑁 𝑗=0 12: š ← FFT 𝑁 (𝑠 𝑗 ) 2𝑁 𝑗=0
13: return (𝑢 • , 𝑣 • ), š

The following proposition establishes the exponential convergence of this method w.r.t. the approximation degree 𝑁 and its quasi-linear arithmetic complexity. Proposition 3.5. There exists 𝜂 > 0 (depending on 𝐻 and ℎ) such that if the analytic initial guess (𝑥 • , 𝑦 • ) ∈ (𝐶 0 2𝜋 ) 2 is 𝜂-close to (𝑥 * , 𝑦 * ), then there is 𝜅 > 1 such that for any 𝑁 ∈ N, Algorithm OvalApproxRefine(𝐻, ℎ, 𝑥 • , 𝑦 • , 𝑁 ) computes a degree 𝑁 TP š in 𝑂 (𝑁 log 𝑁 ) arithmetic operations, and the corresponding curve ( x, y), cf. (3.6), converges exponentially fast to the projected parameterization

(𝑥 # , 𝑦 # ) of Γ(ℎ): ∥ š -𝑠 # ∥ ∞ , ∥ x -𝑥 # ∥ ∞ , ∥ y -𝑦 # ∥ ∞ = 𝑂 (𝜅 -𝑁 ) as 𝑁 → ∞.
Proof. Let 𝜂 ′ > 0 denote the 𝜂 of Lemma 3.4. Clearly, one can chose an 𝜂 ∈ (0,

𝜂 ′ ] such that if (𝑥 • , 𝑦 • ) is 𝜂-close to (𝑥 * , 𝑦 * ), then (𝑢 • , 𝑣 • ) = I 𝑁 ∇𝐻 (𝑥 • ,𝑦 • ) ∥ ∇𝐻 (𝑥 • ,𝑦 • ) ∥ 2 (lines 3 and 11) is 𝜂 ′ -close to ∇𝐻 (𝑥 * ,𝑦 * )
∥ ∇𝐻 (𝑥 * ,𝑦 * ) ∥ 2 , since by assumption ∇𝐻 does not vanish over Γ(ℎ). Then the analytic shift 𝑠 # and the projected parameterization (𝑥 # , 𝑦 # ) are well-defined. As in the proof of Proposition 3.1, we bound the total error with the sum of the error of interpolation and the interpolation of the error:

∥ š -𝑠 # ∥ ∞ ⩽ ∥I 𝑁 (𝑠 # ) -𝑠 # ∥ ∞ + Λ 𝑁 sup 0⩽𝑗⩽2𝑁 | š (𝑡 𝑗 ) -𝑠 # (𝑡 𝑗 )|,
where 𝑡 𝑗 = 2𝑗𝜋 2𝑁 +1 are the equispaced interpolation points. • Error of interpolation. Since 𝑠 # ∈ 𝐶 0 2𝜋 is 2𝜋-periodic and analytic over R by Lemma 3.4, there is a 𝜌 > 0 such that the open horizontal strip R + (-𝜌, 𝜌)𝑖 avoids all possible singularities of 𝑠 # in the complex plane. Then the interpolants I 𝑁 (𝑠 # ) converge uniformly to 𝑠 # over R in 𝑂 (𝜅 -𝑁 ) for any 𝜅 ∈ (1, 𝑒 𝜌 ) (see [Wright et al. 2015, Thm. 4.2]).

• Interpolation of the error. Since the interpolation in line 12 uses approximations 𝑠 𝑗 of the exact points 𝑠 # (𝑡 𝑗 ), the resulting "interpolation of the error" depends on the accuracy reached after 𝑘 = ⌈log 𝑁 ⌉ Newton iterations (lines 2-10). By Eq. (B.1) in the proof of Lemma 3.4, we have:

|𝑠 𝑗 -𝑠 # (𝑡 𝑗 )| ⩽ |N 𝑘 𝑥 𝑗 ,𝑦 𝑗 ,𝑢 𝑗 ,𝑣 𝑗 (0) -𝑠 # (𝑡 𝑗 )| ⩽ |𝑠 # (𝑡 𝑗 )| 2 2 𝑘 -1 ⩽ 2 -𝑁 ∥𝑠 # ∥ ∞ 2
.

• Combining both errors, we conclude that ∥𝑠 -𝑠 # ∥ ∞ = 𝑂 (𝜅 -𝑁 ) for any 𝜅 ∈ (1, min(𝑒 𝜌 , 2)). It is also clear from Eq. (3.4) and (3.6) that the same convergence happens for ( x, y) to (𝑥 # , 𝑦 # ).

• Complexity. The sampling of initial points and reinterpolation through FFT 𝑁 /IFFT 𝑁 requires 𝑂 (𝑁 log 𝑁 ) arithmetic operations. For the projection step (lines 2-10), each of the 2𝑁 + 1 points is refined through ⌈log 𝑁 ⌉ Newton iterations, where the cost of a single iteration is 𝑂 (1). This justifies the overall complexity in 𝑂 (𝑁 log 𝑁 ) arithmetic operations. □

A POSTERIORI VALIDATION OF THE OVAL Γ(ℎ)

The purpose of Algorithm OvalValid presented in this section is to enclose the oval Γ(ℎ) in a validated tube 𝑇 ( x, y, 𝑢 • , 𝑣 • , 𝜀) around the approximate parameterization ( x, y) (3.6) containing 2𝜋-periodic curves:

𝑇 ( x, y, 𝑢 • , 𝑣 • , 𝜀) = (𝑥, 𝑦) ∈ (𝐶 0 2𝜋 ) 2 ∃𝑠 ∈ 𝐶 0 2𝜋 , ∥𝑠 ∥ ∞ ⩽ 𝜀 and 𝑥 (𝑡) = x (𝑡) + 𝑠 (𝑡)𝑢 • (𝑡) 𝑦 (𝑡) = y (𝑡) + 𝑠 (𝑡)𝑣 • (𝑡)
.

The width 𝜀 of such a tube is rigorously computed by using the fixed-point based a posteriori approach presented in Section 2.3 to bound the distance between š computed by OvalApproxRefine and the exact shift 𝑠 # defining the parameterization (𝑥 # , 𝑦 # ) of Γ(ℎ) (see Definition 3.3). Under a few additional conditions also checked by OvalValid (winding number, etc.10), this ensures that (𝑥 # , 𝑦 # ) exists, satisfies Lemma 3.4, and lies in the tube 𝑇 ( x, y, 𝑢 • , 𝑣 • , 𝜀).

Using Equation (3.3), we construct a polynomial Newton-like operator N ∈ 𝐶 0 2𝜋 [𝑋 ] acting on a function 𝑠 ∈ 𝐶 0 2𝜋 :

N (𝑠) = 𝑠 -ǎ (𝐻 (𝑥 • + 𝑠 𝑢 • , 𝑦 • + 𝑠 𝑣 • ) -ℎ) , (4.1
) for which 𝑠 # is the unique fixed point in some small neighborhood of š. Here,

ǎ(𝑡) ≈ 1 𝑢 • (𝑡) 𝜕𝐻 𝜕𝑥 ( x (𝑡), y (𝑡)) + 𝑣 • (𝑡) 𝜕𝐻 𝜕𝑦 ( x (𝑡), y (𝑡))
, is a TP built by interpolation to make the derivative of N ,

N ′ (𝑠) = 1 -ǎ 𝑢 • 𝜕𝐻 𝜕𝑥 (𝑥 • + 𝑠 𝑢 • , 𝑦 • + 𝑠 𝑣 • ) + 𝑣 • 𝜕𝐻 𝜕𝑦 (𝑥 • + 𝑠 𝑢 • , 𝑦 • + 𝑠 𝑣 • ) ,
small in the neighborhood of 𝑠 = š (that is, geometrically, in the neighborhood of the approximate curve ( x, y)).

Contrary to the a posteriori validation for the division presented in Section 2.3.1, the Newton-like operator N in (4.1) is not of degree 1 in 𝑠 unless 𝐻 itself has degree 1. As a consequence, the value of ∥N ′ (𝑠) ∥ ∞ depends on 𝑠, so that the conditions of Theorem 2.5 must be satisfied with a finite radius 𝑟 > 0. We propose a generic subroutine NewtonBall (not specific to the oval validation problem considered here) in Appendix C to automate this task. Specifically, given:

10Although such kinds of properties may not always be checked by other existing approaches, it turns out to be necessary in order to move towards formally checked computed-assisted proofs.

• a procedure computing an upper bound 𝜆(𝑟 ) for sup ∥𝑠 -š ∥ ∞ ⩽𝑟 ∥N ′ (𝑠) ∥ ∞ for any 𝑟 ⩾ 0,

• an upper bound 𝑑 for the defect ∥N ( š)š ∥ ∞ , NewtonBall(𝑑, 𝜆) either computes a radius 𝑟 and a Lipschitz constant 𝜆 = 𝜆(𝑟 ) satisfying the hypotheses of Theorem 2.5, or returns an error.

Algorithm 4 OvalValid(𝐻, ℎ, 𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 , 𝜀 0 , 𝑥 𝑟 , 𝑦 𝑟 , 𝜍, 𝑥 • , 𝑦 • , 𝑢 • , 𝑣 • , š, 𝑁 ) Input: 𝐻, ℎ, 𝜍 = ±1, (𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 ) ∈ R 2 , 𝜀 0 > 0, (𝑥 𝑟 , 𝑦 𝑟 ) ∈ R 2 , and degree 𝑁 TPs 𝑥 • , 𝑦 • , 𝑢 • , 𝑣 • , š Output: a bound 𝜀 > 0 defining a tube 𝑇 ( x, y, 𝑢 • , 𝑣 • , 𝜀) around ( x, y) containing (𝑥 # , 𝑦 # ) ⊲ Build TP ǎ defining the Newton-like operator N 1: (𝑥 𝑗 ) 2𝑁 𝑗=0 ← IFFT 𝑁 ( x) and (𝑦 𝑗 ) 2𝑁 𝑗=0 ← IFFT 𝑁 ( y) 2: (𝑢 𝑗 ) 2𝑁 𝑗=0 ← IFFT 𝑁 (𝑢 • ) and (𝑣 𝑗 ) 2𝑁 𝑗=0 ← IFFT 𝑁 (𝑣 • ) 3: ǎ ← FFT 𝑁 1 𝑢 𝑗 𝜕𝐻 𝜕𝑥 (𝑥 𝑗 , 𝑦 𝑗 ) + 𝑣 𝑗 𝜕𝐻 𝜕𝑦 (𝑥 𝑗 , 𝑦 𝑗 ) 2𝑁 𝑗=0
⊲ Bound operator norm of N , using rigorous operations with intervals and RTAs

4: 𝒙 ← 𝑥 • + š 𝑢 • and 𝒚 ← 𝑦 • + š 𝑣 • ⊲ (𝒙, 𝒚) RTAs for ( x, y) 5: 𝒙 ′ ← 𝑥 •′ + š′ 𝑢 • + š 𝑢 •′ and 𝒚 ′ ← 𝑦 •′ + š′ 𝑣 • + š 𝑣 •′ 6: 𝑑 ← Bound(𝒔) where 𝒔 ← ǎ (ℎ -𝐻 (𝒙, 𝒚)) ⊲ defect 7: 𝜆 0 ← Bound(ℓ 0 ) where ℓ 0 ← 1 -ǎ 𝑢 • 𝜕𝐻 𝜕𝑥 (𝒙, 𝒚) + 𝑣 • 𝜕𝐻 𝜕𝑦 (𝒙, 𝒚) 8: 𝛼 ← Bound( ǎ) 9: 𝑋 ← Bound(𝒙) • [-1, 1] and 𝑌 ← Bound(𝒚) • [-1, 1] 10: 𝑈 ← Bound(𝑢 • ) • [-1, 1] and 𝑉 ← Bound(𝑣 • ) • [-1,
⊲ Find stable ball ⊲ Validated error bound 𝜀 17: 𝜀 ← 𝑑/(1 -𝜆 + ) ⊲ Check initial conditions 18: if (𝒙 (0) -𝜀𝑢 • (0) -𝑥 𝑖𝑛𝑖 ) 2 + (𝒚(0) -𝜀𝑣 • (0) -𝑦 𝑖𝑛𝑖 ) 2 > 𝜀 2 0 or 19: (𝒙 (0) + 𝜀𝑢 • (0) -𝑥 𝑖𝑛𝑖 ) 2 + (𝒚(0) + 𝜀𝑣 • (0) -𝑦 𝑖𝑛𝑖 ) 2 > 𝜀 2 0 , then return 𝐹𝐴𝐼𝐿 end if
⊲ Check the winding number of the curve

20: 𝒓2 ← (𝒙 -𝑥 𝑟 ) 2 + (𝒚 -𝑦 𝑟 ) 2 21: if ¬ IsPositive 𝒓2 -𝜀 2 Bound(𝑢 •2 + 𝑣 •2 ) then return 𝐹𝐴𝐼𝐿 end if 22: 𝝇 ← 1 2𝜋 ∫ 2𝜋 0 (𝒙 (𝑡) -𝑥 𝑟 )𝒚 ′ -(𝒚(𝑡) -𝑦 𝑟 )𝒙 ′ 𝒓2(𝑡) d𝑡 23: if 𝝇 ∩ Z ≠ {𝜍 }
š -𝑠 # ∥ ∞ ⩽ 𝜀.
Moreover, for a sufficiently large degree 𝑁 , the process never fails and returns a bound 𝜀 that converges exponentially fast to 0 w.r.t. 𝑁 , in 𝑂 (𝑁 log(𝑁 )) arithmetic operations.

Proof. We apply Theorem 2.5 over the space 𝐶 0 2𝜋 to compute an upper bound 𝜀 for ∥𝑠 # ∥ ∞ , where 𝑠 # ∈ 𝐶 0 2𝜋 is solution to (3.3). First (lines 1-3), the Newton-like operator N is built as in Equation ( 4.1) by computing a TP ǎ interpolating the function which represents the inverse differential. Next, the algorithm computes a rigorous bound 𝑑 on the defect ∥N ( š)š ∥ ∞ (line 6), and produces a procedure 𝜆 s.t. 𝜆(𝑟 ) is a Lipschitz constant for N over 𝐵( š, 𝑟 ) in 𝐶 0 2𝜋 w.r.t. to the ∥ • ∥ ∞ norm (lines 7-15). The formula for 𝜆(𝑟 ) in line 14 comes from Lemma 2.6, with 𝑈 , 𝑉 enclosures for the range of 𝑢 • , 𝑣 • , and 𝑋 𝑟 , 𝑌 𝑟 enclosures for the range of 𝑥 • + 𝑠𝑢

• = x + (𝑠 -š)𝑢 • and 𝑦 • + 𝑠𝑣 • = y + (𝑠 -š)𝑣 • ,
for all 𝑠 ∈ 𝐵( š, 𝑟 ). After that, the auxiliary routine NewtonBall (line 16) described in Appendix C rigorously computes a suitable radius 𝑟 + satisfying the hypotheses of Theorem 2.5, if possible. We thus deduce the existence and uniqueness of 𝑠 # ∈ 𝐶 0 2𝜋 (and hence (𝑥 # , 𝑦 # )) together with a validated bound 𝜀 (line 17).

• Initial conditions. By hypothesis, the ball 𝐵((𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 ), 𝜀 0 ) in R 2 intersects 𝐻 -1 (ℎ) in its component Γ(ℎ) only. Hence, by continuity, checking that (𝑥 # (0), 𝑦 # (0)) lies in that ball is sufficient to ensure that (𝑥 # (𝑡), 𝑦 # (𝑡)) belongs to the component Γ(ℎ) for all 𝑡 ∈ [0, 2𝜋]. Line 19 does this by checking that both endpoints of the transversal segment of the tube at time 0 belong to that ball.

• Winding number. First, lines 20-21 ensure that the trajectories in the tube 𝑇 ( x, y, 𝑢 • , 𝑣 • , 𝜀) all avoid the point (𝑥 𝑟 , 𝑦 𝑟 ) inside Γ(ℎ). Then, the winding number of ( x, y) is rigorously computed in line 22 and compared to the reference value 𝜍 at line 23. Since the tube avoids (𝑥 𝑟 , 𝑦 𝑟 ), the winding number remains the same for all curves in it (so in particular (𝑥 # , 𝑦 # )).

• Analyticity. The Banach fixed-point based argument used previously only guarantees the continuity of (𝑥 # , 𝑦 # ). To gain more regularity, we need to consider a larger function space, namely the Banach space of continuous 2𝜋-periodic functions defined over a horizontal strip S 𝜌 := R + [-𝜌, 𝜌]𝑖 of C for some small 𝜌 > 0, equipped with the norm: ∥𝜑 ∥ S 𝜌 = sup 𝑧 ∈ S 𝜌 |𝜑 (𝑧)|11. The a posteriori validation argument based on Theorem 2.5 applies the same way, except that we need to replace ∥ • ∥ ∞ with ∥ • ∥ S 𝜌 in all the bounds defining 𝑑 and 𝑟 ↦ → 𝜆(𝑟 ). If we choose 𝜌 small enough, the inequalities 𝜆(𝑟 + ) < 1 and 𝑑 + 𝑟 + 𝜆(𝑟 + ) < 𝑟 + guaranteed by NewtonBall still hold for the same 𝑟 + , which means that Theorem 2.5 now proves the existence of a unique continuous periodic solution to (3.3) extending 𝑠 # over S 𝜌 . Moreover, the iterations N 𝑘 (0) are all trigonometric polynomials (hence analytic) and they converge uniformly to 𝑠 # over S 𝜌 as 𝑘 → ∞. As a consequence, the limit 𝑠 # is analytic over R.

• Convergence. Here we outline the main arguments and leave the technical details to the reader. First, we know from Proposition 3.5 that the input š computed by OvalApproxRefine tends exponentially fast to 𝑠 # . This holds for the ∥ • ∥ ∞ norm, but also for the ∥ • ∥ ℓ 1 norm since they are TPs of degree 𝑂 (𝑁 ). As a consequence, the computed bounds 𝑑 and 𝜆 0 tend to 0 exponentially fast, while the other bounds 𝛼, 𝑋 , 𝑌 , 𝑈 , 𝑉 and 𝜆 1 (for given 𝑟 ) tend to the mathematical values obtained by replacing the approximations š, x, y by 𝑠 # , 𝑥 # , 𝑦 # . Therefore, in the routine NewtonBall, for 𝑁 sufficiently large, the initial radius 𝑟 max computed in line 1 is bounded away from 0 and satisfies 𝑑 + 𝜆(𝑟 max )𝑟 max < 𝑟 max . Hence, NewtonBall necessarily returns an 𝑟 satisfying 𝑑 + 𝜆(𝑟 )𝑟 < 𝑟 , and then OvalValid computes 𝜀 = 𝑑/(1 -𝜆 + ) ⩽ 2𝑑, which tends exponentially fast to 0. To conclude, similar arguments show that the remaining tests (initial conditions and winding number) never fail for 𝑁 sufficiently large.

11Note that this argument is purely theoretical: there is no need to reimplement Algorithm OvalValid to actually compute with this norm ∥ • ∥ S𝜌 . , Vol. 1, No. 1, Article . Publication date: December 2023.

• Complexity. The asymptotic complexity w.r.t. the degree 𝑁 is determined by the FFT 𝑁 /IFFT 𝑁 routines used for evaluation, interpolation, RTA multiplications, the RTA positivity test in line 21 and the RTA division in line 22. Note that the execution time of NewtonBall does not depend on 𝑁 , since the function 𝑟 ↦ → 𝜆(𝑟 ) (lines 11-15) involves intervals rather than RTAs. □

RIGOROUS ENCLOSURE OF THE ABELIAN INTEGRAL ℑ(ℎ)

According to Proposition 4.1, the tube 𝑇 ( x, y, 𝑢 • , 𝑣 • , 𝜀) computed by OvalValid contains a unique analytic parameterization (𝑥 # , 𝑦 # ) of Γ(ℎ) with correct winding number, so that we have:

ℑ(ℎ) = ∫ 2𝜋 0 𝑃 (𝑥 # (𝑡), 𝑦 # (𝑡))𝑦 # ′ (𝑡) -𝑄 (𝑥 # (𝑡), 𝑦 # (𝑡))𝑥 # ′ (𝑡) 𝜇 (𝑥 # (𝑡), 𝑦 # (𝑡))
d𝑡 .

(5.1)

In order to provide an interval enclosure for ℑ(ℎ), we must bound the difference between the exact value given by Equation (5.1) and the corresponding integral taken along the approximate curve ( x, y) defined by š:

𝐼 š = ∫ 2𝜋 0 𝑃 ( x (𝑡), y (𝑡)) y′ (𝑡) -𝑄 ( x (𝑡), y (𝑡)) x′ (𝑡) 𝜇 ( x (𝑡), y (𝑡)) d𝑡, (5.2)
which can be rigorously computed using operations on RTAs. To do so, we make use of Stokes' theorem and bound the resulting term:

𝑐 (𝑥, 𝑦) d𝑥 ∧ d𝑦 = 𝜕(𝑃/𝜇) 𝜕𝑥 (𝑥, 𝑦) + 𝜕(𝑄/𝜇) 𝜕𝑦 (𝑥, 𝑦) d𝑥 ∧ d𝑦, (5.3)
over the region located between two curves ( x, y) and (𝑥 # , 𝑦 # ), which is safely overapproximated by the tube 𝑇 ( x, y, 𝑢 • , 𝑣 • , 𝜀).

Algorithm 5 AbIntValidQuad(𝑃, 𝑄, 𝜇, 𝑥 ⊲ Compute rigorous integral along approximate curve using RTAs

1: 𝒙 ← 𝑥 • + š 𝑢 • and 𝒚 ← 𝑦 • + š 𝑣 • ⊲ (𝒙, 𝒚) RTAs for ( x, y) 2: 𝒙 ′ ← 𝑥 •′ + š′ 𝑢 • + š 𝑢 •′ and 𝒚 ′ ← 𝑦 •′ + š′ 𝑣 • + š 𝑣 •′ 3: 𝑷 ← 𝑃 (𝒙, 𝒚), 𝑸 ← 𝑄 (𝒙, 𝒚), and 𝝁 ← 𝜇 (𝒙, 𝒚) 4: 𝑬 ← (𝑷 𝒚 ′ -𝑸 𝒙 ′ ) / 𝝁 5: 𝑰 š ← ∫ 2𝜋 0 𝑬 (𝑡) d𝑡
⊲ Compute integration error in the tube using Stokes' theorem

6: 𝒙 𝜀 ← (𝑥 • + š 𝑢 • , 𝜀 Bound(𝑢 • ))) and 𝒚 𝜀 ← (𝑦 • + š 𝑣 • , 𝜀 Bound(𝑣 • ))) 7: 𝑐 ← 𝜕 (𝑃 /𝜇 ) 𝜕𝑥 + 𝜕 (𝑄/𝜇 ) 𝜕𝑦
⊲ symbolically as a rational fraction

8: 𝐶 ← Bound(𝑐 (𝒙 𝜀 , 𝒚 𝜀 )) 9: 𝐵 1 ← Bound(𝒙 ′ 𝑣 • -𝒚 ′ 𝑢 • ) 10: 𝐵 2 ← Bound(𝑢 •′ 𝑣 • -𝑣 •′ 𝑢 • ) 11: 𝛿 ← 2𝜋𝜀 (2𝐵 1 + 𝜀𝐵 2 )𝐶 ⊲ Return interval enclosure 12: return 𝑰 ← 𝑰 š + [-𝛿, 𝛿]
Proposition 5.1. Let 𝑃, 𝑄, 𝜇 be as in Problem 1.1. If AbIntValidQuad(𝑃, 𝑄, 𝜇, 𝑥 • , 𝑦 • , 𝑢 • , 𝑣 • , š, 𝜀, 𝑁 ) with 𝜀 > 0 computed by OvalValid does not fail, then it returns an interval enclosure 𝑰 of ℑ(ℎ).

Moreover, if 𝑁 is sufficiently large, then it does never fail and computes an interval 𝑰 whose diameter tends exponentially fast to 0 w.r.t. 𝑁 , in 𝑂 (𝑁 log(𝑁 )) arithmetic operations.

Proof. First, an interval enclosure 𝑰 š of the integral 𝐼 š (5.2) along the approximate curve ( x, y) is computed rigorously by using RTA representations and operations for the integrand (lines 3 and 4), and then integrating it over [0, 2𝜋] (line 5). Next, lines 6-11 compute a bound 𝛿 for the error between 𝐼 š and the exact value ℑ(ℎ) taken along (𝑥 # , 𝑦 # ) of Γ(ℎ), so that AbIntValidQuad returns a validated enclosure 𝑰 = 𝑰 š + [-𝛿, 𝛿] for ℑ(ℎ). The strategy to obtain 𝛿 is detailed below.

• Error bound 𝛿 with Stokes' theorem. Let (𝑥, 𝑦) denote the standard coordinates of the plane R 2 , and add a third one, 𝑧, living in R/2𝜋Z, corresponding to the time 𝑡. The approximate parameterization ( x, y) traces a smooth simple closed curve

𝑡 ∈ R/2𝜋Z ↦ → (𝑥 • (𝑡) + š (𝑡)𝑢 • (𝑡), 𝑦 • (𝑡) + š (𝑡)𝑣 • (𝑡), 𝑡) in R 2 × R/2𝜋Z. Similarly, (𝑥 # , 𝑦 # ) traces 𝑡 ∈ R/2𝜋Z ↦ → (𝑥 • (𝑡) + 𝑠 # (𝑡)𝑢 • (𝑡), 𝑦 • (𝑡) + 𝑠 # (𝑡)𝑣 • (𝑡), 𝑡).
Since 𝑠 # is analytic, it only has a finite number of zeros in [0, 2𝜋]. Hence, since (𝑢 • , 𝑣 • ) never vanishes, the two curves only intersect each other a finite number of times. The "strip" S between them,

S = 𝑥 • (𝑡) + 𝑠𝑢 • (𝑡), 𝑦 • (𝑡) + 𝑠𝑣 • (𝑡), 𝑡 , 𝑡 ∈ R/2𝜋Z, 𝑠 ∈ [š (𝑡), 𝑠 # (𝑡)] or [𝑠 # (𝑡), š (𝑡)] ,
is a finite collection of two-dimensional cells (with two singularities each). Choosing the appropriate orientation for each cell, the difference between ℑ(ℎ) and 𝐼 š is equal to the integral of the 1-form (𝑃/𝜇)d𝑦 -(𝑄/𝜇)d𝑥 along the oriented boundary 𝜕S of S:

ℑ(ℎ) -𝐼 š = ∫ 𝜕S 𝑃 (𝑥, 𝑦)d𝑦 -𝑄 (𝑥, 𝑦) d𝑥 𝜇 (𝑥, 𝑦)
.

Using Stokes' theorem [Lee 2013, Thm. 14.9], and since 𝑃 and 𝑄 do not depend on 𝑧, we obtain:

|ℑ(ℎ) -𝐼 š | ⩽ ∫ S |𝑐 (𝑥, 𝑦)| d𝑥 ∧ d𝑦,
with 𝑐 (𝑥, 𝑦) defined as in Equation (5.3).

A bound 𝐶 for the integrand |𝑐 (𝑥, 𝑦)| is computed in lines 6-8. Moreover, the integration domain S can be safely overapproximated by the geometric realization 𝑇 𝜀 of the tube 𝑇 ( x, y, 𝑢 • , 𝑣 • , 𝜀) in R 2 × R/2𝜋Z, which is the two-dimensional manifold parameterized by:

𝜎 : R/2𝜋Z × [-𝜀, 𝜀] -→ 𝑇 𝜀 (𝑡, 𝑠) ↦ -→ x (𝑡) + 𝑠 𝑢 • (𝑡) y (𝑡) + 𝑠 𝑣 • (𝑡) 𝑡 .
Hence, the difference ℑ(ℎ) -𝐼 š is bounded by:

|ℑ(ℎ) -𝐼 š | ⩽ 𝐶 ∫ T 𝜀 |d𝑥 ∧ d𝑦| = 𝐶 ∫ 2𝜋 0 ∫ 𝜀 -𝜀 x′ (𝑡) 𝑢 • (𝑡) y′ (𝑡) 𝑣 • (𝑡) + 𝑠 𝑢 •′ (𝑡) 𝑢 • (𝑡) 𝑣 •′ (𝑡) 𝑣 • (𝑡) d𝑠 d𝑡 ⩽ 2𝜋𝜀 (2𝐵 1 + 𝜀𝐵 2 ) 𝐶 = 𝛿, with 𝐵 1 ⩾ sup 0⩽𝑡 ⩽2𝜋 x′ (𝑡) 𝑢 • (𝑡) y′ (𝑡) 𝑣 • (𝑡) and 𝐵 2 ⩾ sup 0⩽𝑡 ⩽2𝜋 𝑢 •′ (𝑡) 𝑢 • (𝑡) 𝑣 •′ (𝑡) 𝑣 • (𝑡)
, computed in lines 9 and 10. Hence,

ℑ(ℎ) = 𝐼 š + (ℑ(ℎ) -𝐼 š ) ∈ 𝑰 š + [-𝛿, 𝛿] = 𝑰 .
• Convergence. When 𝑁 is large enough, ( x, y) is sufficiently close to Γ(ℎ), and since 𝜇 does not vanish in a neighborhood of this oval, all the functions involved in the computation of 𝑰 š (lines 3-5) are analytic; hence the diameter of this interval tends exponentially fast to 0 w.r.t. 𝑁 . The same arguments hold for the computation of Stokes' error bound 𝛿 in lines 6-11. Note that 𝐵 1 , 𝐵 2 and 𝐶 are bounded since they converge to the corresponding ∥ • ∥ ℓ 1 norms with ( x, y) and (𝑢 • , 𝑣 • ) replaced by (𝑥 # , 𝑦 # ) and

∇𝐻 (𝑥 • ,𝑦 • ) ∥ ∇𝐻 (𝑥 • ,𝑦 • ) ∥ 2 ,
respectively. Since 𝜀 tends to 0 exponentially fast w.r.t. 𝑁 according to Proposition 4.1, so does 𝛿 (line 11), and consequently the total diameter of 𝑰 in line 12.

• Complexity. The asymptotic complexity w.r.t. the globally set RTA degree 𝑁 is determined by the FFT 𝑁 /IFFT 𝑁 routines used for TP/RTA multiplications and divisions. □

Now we state our main algorithm AbIntValid and our main theorem, which is a straightforward consequence of Propositions 3.1, 3.5, 4.1 and 5.1.

Algorithm 6 AbIntValid(𝐻, ℎ, Σ, 𝑃, 𝑄, 𝜇, 𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 , 𝜀 0 , 𝑥 𝑟 , 𝑦 𝑟 , 𝜍, 𝑁 0 , 𝑁 )

Input: Potential 𝐻 , level ℎ, transversal Σ, 𝑃, 𝑄, 𝜇 as in Problem 1.1, (𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 ) ∈ R 2 , 𝜀 0 > 0, (𝑥 𝑟 , 𝑦 𝑟 ) ∈ R 2 , 𝜍 = ±1, 𝑁 0 ∈ N and 𝑁 ⩾ 𝑁 0 Output: interval enclosure 𝑰 of ℑ(ℎ) 1: (𝑥 • , 𝑦 • ) ← OvalApproxInit(𝐻, Σ, 𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 , 𝑁 0 ) 2: (𝑢 • , 𝑣 • ), š ← OvalApproxRefine(𝐻, ℎ, 𝑥 • , 𝑦 • , 𝑁 ) 3: 𝜀 ← OvalValid(𝐻, ℎ, 𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 , 𝜀 0 , 𝑥 𝑟 , 𝑦 𝑟 , 𝜍, 𝑥 • , 𝑦 • , 𝑢 • , 𝑣 • , š, 𝑁 ) 4: return 𝑰 ← AbIntValidQuad(𝑃, 𝑄, 𝜇, 𝑥 • , 𝑦 • , 𝑢 • , 𝑣 • , š, 𝜀, 𝑁 )
Theorem 5.2. Let 𝐻, 𝑃, 𝑄, 𝜇 be as in Problem 1.1, ℎ be a regular level value on a transversal Σ, the ball of center (𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 ) and radius 𝜀 0 contain a unique point in 𝐻 -1 (ℎ) ∩ Σ, the point (𝑥 𝑟 , 𝑦 𝑟 ) be in the interior of the oval Γ(ℎ), and 𝜍 ∈ {-1, 1} be the winding number of Γ(ℎ) w.r.t. (𝑥 𝑟 , 𝑦 𝑟 ) in the orientation prescribed by the Hamiltonian. Then for a sufficiently large degree 𝑁 0 (depending on the preceding input data but not on the target accuracy), as the approximation degree 𝑁 → ∞, Algorithm AbIntValid(𝐻, ℎ, Σ, 𝑃, 𝑄, 𝜇, 𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 , 𝜀 0 , 𝑥 𝑟 , 𝑦 𝑟 , 𝜍, 𝑁 0 , 𝑁 ) never fails and returns an interval enclosure 𝑰 of ℑ(ℎ) in 𝑂 (𝑁 log 𝑁 ) arithmetic operations. The diameter of 𝑰 tends to 0 exponentially fast, i.e., there exists 𝜅 > 1 depending on 𝐻, ℎ, Σ, 𝑃, 𝑄, 𝜇, 𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 , 𝜀 0 , 𝑥 𝑟 , 𝑦 𝑟 , 𝜍, 𝑁 0 such that diam(𝑰 ) = 𝑂 (𝜅 -𝑁 ).

IMPLEMENTATION DETAILS AND CONCURRENT APPROACHES

We now present the Julia implementation of our method and review concurrent approaches. The repositories for our code are available from https://gitlab.inria.fr/abintvalid.

A brief overview of our Julia implementation

Two main reasons motivating our choice of the Julia language for the prototype implementation accompanying this article are: the convenience of the interactive mode for experimenting, with performances close to compiled languages such as C; and the existing libraries for numerical computing, notably approximations with Fourier series. We developed two Julia packages: RigorousFourier.jl to provide RTA data structures and operations used throughout this article, and AbelianIntegral.jl to implement Algorithm AbIntValid and its subroutines.

RigorousFourier.jl (v.0.0.1). We built this package on top of the existing ApproxFun.jl package12 which provides a very modular framework to approximate and manipulate functions in 12https://juliaapproximation.github.io/ApproxFun.jl/latest/ various bases, including Fourier: A degree 𝑁 TP is stored as the length 2𝑁 + 1 vector of its coefficients represented by floating-point numbers of arbitrary precision (the BigFloat type in Julia is a wrapper for the standard GNU MPFR C library13). Since this data structure is parameterized by the type of the coefficients, TPs can also be defined with interval coefficients14. Therefore, in RigorousFourier.jl, RTAs are implemented as pairs made up of a TP with interval coefficients together with a floatingpoint error bound. Operations on them, described in Section 2, are implemented using interval operations on the interval coefficients, and floating-point operations with the suitable rounding mode activated for the error component. Finally, a global approximation degree 𝑁 is set by the user so that TPs and RTAs are automatically truncated to degree 𝑁 after each operation. Note that the degree 𝑁 ′ = 2 ⌈log 2 (𝑁 ) ⌉ ⩾ 𝑁 is used internally for FFT based operations (such as multiplication) for the sake of efficiency. Remark 6.1. Using floating-point coefficients for TPs f adds rounding errors to the approximation error w.r.t. to the represented functions 𝑓 . Similarly, using floating-point interval coefficients for RTAs 𝒇 enlarges the set of functions contained in 𝒇 , due to the growth of the interval widths. However, it can be easily proved that with 𝑝 ≈ 𝑁 bits for the underlying floating-point arithmetic, the total contribution of rounding errors in the evaluation of a fixed expression involving TPs or RTAs and basic operations (addition, subtraction, coefficientwise multiplication, integration, etc.) of Section 2.2 is bounded by 𝑂 (𝑁 𝑐 2 -𝑁 ) for some 𝑐 > 0. Therefore, these errors will not affect the convergence result of Theorem 5.2.

The next remark concerns FFT-based operations: interpolation, fast multiplication, division, composition, etc. Remark 6.2. Directly using an FFT scheme on intervals rather than floating-point numbers may lead to gross overestimations. Some works addressed the use of FFT for convolutions in a rigorous setting with a focus on several related issues (e.g., aliasing errors [START_REF] Cyranka | Efficient and generic algorithm for rigorous integration forward in time of dPDEs: Part I[END_REF][START_REF] Lessard | Computing discrete convolutions with verified accuracy via Banach algebras and the FFT[END_REF], overestimation of tail coefficients [START_REF] Lessard | Computing discrete convolutions with verified accuracy via Banach algebras and the FFT[END_REF]], wrapping effect when input coefficients are given as large intervals [START_REF] Liu | Fast convolution and fast Fourier transform under interval and fuzzy uncertainty[END_REF]). However, the intrinsic growth of the intervals during the execution of the FFT algorithm remains an open and insufficiently documented problem in the literature (see e.g., [Cyranka 2014, §3.4]). We currently propose two approaches to alleviate this issue when performing fast and rigorous RTA multiplication:

(𝑎) An efficient and numerically stable way is to first "flatten" RTAs (i.e., interval coefficients are shrunk to their floating-point midpoint, with the error component updated accordingly), then to apply a floating-point FFT scheme, and finally to update the error component with known rigorous estimates on the resulting rounding errors (see e.g., [START_REF] Brisebarre | Error Analysis of Some Operations Involved in the Cooley-Tukey Fast Fourier Transform[END_REF]).

Using 𝑝 ≈ 𝑁 bits of precision, the claimed bounds are in 𝑂 (𝑁 2 2 -𝑁 ), which, again, does not affect the convergence result of Theorem 5.2. In our current prototype implementation, though, this last step is still missing, because adapting these error estimates to the precise FFT scheme used in ApproxFun.jl requires some additional work. (𝑏) A more direct approach is to apply an interval FFT scheme on the flattened coefficients.

Although quite simple, this method may generate larger error bounds and leads to an overhead factor of roughly 5 to 8 in timings compared to method (𝑎). We however also provide the timings for this method (𝑏) since it is the only fully rigorous one in our prototype code. AbelianIntegral.jl (v.0.0.1). In this package, we implemented all routines related to the evaluation of Abelian integrals, described in Sections 3 to 5. Its main dependencies are:

• the TypedPolynomials.jl package15 for multivariate polynomials used in the data structures representing integrable systems, oval families and perturbations; • the DifferentialEquations.jl package16 that provides a numerical ODE solver to instantiate the DSolve routine used by Algorithm OvalApproxInit to compute the initial guess (𝑥 • , 𝑦 • ). As of today, it uses by default a Tsitouras 5/4 Runge-Kutta method with 4th order interpolant [START_REF] Tsitouras | Runge-Kutta pairs of order 5(4) satisfying only the first column simplifying assumption[END_REF]]. This package also offers event detection features, which we use to end the integration when the first return onto the transversal Σ is detected. The resulting implementation of DSolve is thus compliant with Assumption 3.2.

Remark 6.3. Since Assumption 3.2 is satisfied in our implementation, and considering the rounding error estimates given in Remarks 6.1 and 6.2 (assuming option (𝑎)), the convergence results of the preceding sections in 𝑂 (𝜅 -𝑁 ) for some 𝜅 ∈ (1, 2) still hold for the floating-point implementation of AbIntValid. Combining this with the bit complexity of 𝑂 (𝑁 log 𝑁 ) for arithmetic operations on 𝑁 -bit floating-point numbers [START_REF] Harvey | Integer multiplication in time 𝑂 (𝑛log 𝑛)[END_REF], a refined version of Theorem 5.2 says that Algorithm AbIntValid, using degree 𝑁 TPs/RTAs and floating-point arithmetic with 𝑁 bits of precision (neglecting considerations about the exponent range), computes a rigorous enclosure of ℑ(ℎ) with 𝑂 (𝑁 ) bits of accuracy in 𝑂 (𝑁 2 log 2 𝑁 ) bit operations.

Concurrent approaches and implementations

Despite the constraint of keeping the mathematical background of the validation routine as minimal as possible to ease the future formal proof implementation, our algorithm remains sufficiently efficient for the rigorous high precision evaluation of Abelian integrals. Timings are compared with the following available software:

• The subdivision algorithm [START_REF] Johnson | On a computer-aided approach to the computation of Abelian integrals[END_REF], by one of the authors of this article, is, to our knowledge, the only fully algorithmic work dedicated to the rigorous evaluation of Abelian integrals. The integration domain (the interior of the oval, by using the Green formula) is safely approximated using boxes. Such a strategy would be a natural candidate for a formal proof implementation. However, contrary to the exponential convergence of our method, this routine has only finite order, which makes it unsuitable for high precision evaluation. • The CAPD library, developed in C++, is well established among the community of computer assisted proofs in dynamical systems, both for its efficiency and reliability. Trajectories are rigorously computed using various higher order methods, including Taylor forms [START_REF] Neumaier | Taylor forms-use and limits[END_REF]]. Our strategy to evaluate an Abelian integral using CAPD is to augment the plane (𝑥, 𝑦) with a third dimension 𝐼 (the integral), and then add a third component to the vector field, representing 𝐼 , which is equal to the integrand of (1.1), where d𝑥 and d𝑦 are replaced by the values of 𝑥 and 𝑦 given by the first two components. However, implementing such a method in Coq would require significant work to formalize some theory of differential equations (e.g., the Picard-Lindelöf theorem) and a rigorous ODE solver. • The Arb library17 [START_REF] Johansson | Arb: efficient arbitrary-precision midpoint-radius interval arithmetic[END_REF]], developed in C, is a general purpose rigorous numerics library providing extremely efficient basic routines, including one for arbitrary-precision numerical integration with rigorous error bounds [START_REF] Johansson | Numerical integration in arbitrary-precision ball arithmetic[END_REF]]. Although it cannot rigorously approximate general algebraic functions yet (in order to approximate the oval), 15https://github.com/JuliaAlgebra/TypedPolynomials.jl. 16https://diffeq.sciml.ai/v2.0/ 17https://arblib.org/ we used it, as a SageMath18 external package, for the first example of Section 7 where explicit parameterizations can be used. In this case, timings easily outperform our method and CAPD. We also wish to mention [START_REF] Lairez | Computing the Volume of Compact Semi-Algebraic Sets[END_REF]. This approach should be quite effective to address Problem 1.1 but we did not include it in our comparisons due to the following two issues: as of today, there is no ready-to-use implementation and, overall, the algorithms rely on tools that are currently out of reach for the proof assistant Coq.

WORKED EXAMPLES, DISCUSSION AND CONCLUSION

We now apply our algorithm on two problems regarding the isolation of limit cycles in near-integrable polynomial planar vector fields:

• The first example (Section 7.1) is an integrable quartic system introduced in [START_REF] Johnson | A quartic system with twenty-six limit cycles[END_REF] to claim H (4) ⩾ 26. Unfortunately, an implementation error led to incorrect evaluation of the Abelian integrals. An analysis led by the authors of this article (see [Bréhard 2011, Chap. 6]) showed that 24 limit cycles could still be obtained from this example. The particular form of the potential function 𝐻 makes it possible to divide the ovals into arcs 𝑥 ↦ → 𝑦 (𝑥) and 𝑦 ↦ → 𝑥 (𝑦) with simple formulas, so that ℑ(ℎ) can easily be computed as a classical integral over a segment of R with an explicit integrand. • The second example (Section 7.2) is a cubic Hamiltonian system used in [START_REF] Li | A cubic system with thirteen limit cycles[END_REF]] to demonstrate Z(3) ⩾ 13. Contrary to our first example, here the ovals of the unperturbed system cannot be easily parameterized using explicit formulas. Therefore, the evaluation of the Abelian integrals necessarily requires a method working on the implicit representation of the ovals. The tests were executed on an Intel i7-6600U 2.60GHz CPU with a 64-bit Linux-based system.

Johnson's symmetric quartic system revisited

In [START_REF] Johnson | A quartic system with twenty-six limit cycles[END_REF]], T. Johnson constructed a perturbed quartic pseudo-Hamiltonian vector field of the form:

𝑥 = -4𝑦 2 (𝑦 2 -𝑌 0 ), 𝑦 = 4𝑥𝑦 (𝑥 2 -𝑋 0 ) + 𝜀𝑔(𝑥, 𝑦), (7.1)
where 𝑋 0 = 9 10 , 𝑌 0 = 11 10 , and 𝑔 is a degree 4 perturbation, 𝑔(𝑥, 𝑦)

= 𝛼 00 + 𝛼 20 𝑥 2 + 𝛼 22 𝑥 2 𝑦 2 + 𝛼 40 𝑥 4 + 𝛼 04 𝑦 4 , (7.2)
with well-chosen coefficients 𝛼 𝑖 𝑗 . Using a rigorous validation integration routine, he claimed to prove the existence of 26 limit cycles, thus surpassing the previously known record H (4) ⩾ 22 [START_REF] Christopher | Limit cycles of differential equations[END_REF]. Unfortunately, a bug in his implementation led him to observe more zeros in the Abelian integral than what actually exists. In [Bréhard 2011, Chap. 6], the authors of this article however showed that by using different values for the coefficients 𝛼 𝑖 𝑗 , one can prove the following Theorem establishing H (4) ⩾ 2419.

Theorem 7.1. The quartic system (7.1) with 𝑋 0 = 9 10 , 𝑌 0 = 11 10 , and the degree 4 perturbation 𝑔 defined as in (7.2) with coefficients:

𝛼 00 = -0.78622148667854837664, 𝛼 20 = 0.87723523612653436051, 𝛼 22 = 1, 𝛼 40 = 0.23742713894293038223, 𝛼 04 = -0.21823846173078863753, 18https://www.sagemath.org/ 19In the meantime, Prohens and Torregrosa [START_REF] Prohens | New lower bounds for the Hilbert numbers using reversible centers[END_REF] showed H (4) ⩾ 28 by a different kind of computer-aided proof (and for a different system).

, Vol. 1, No. 1, Article . Publication date: December 2023.

has at least 24 limit cycles.

A particularity of the potential function 𝐻 associated to this system is that the ovals can be parameterized explicitly using arcs 𝑥 (𝑦) and 𝑦 (𝑥) (see Section 7.1.1 below). This was the strategy adopted in [START_REF] Johnson | A quartic system with twenty-six limit cycles[END_REF]] and [Bréhard 2011, Chap. 6] for the rigorous evaluation of the Abelian integral. In Section 7.1.2, we compare our approach to the other methods described above. Finally, in Section 7.1.3, we redo the computations using the algorithm presented in this article (which does not require explicit parameterizations), which gives a rigorous proof of Theorem 7.1.

Potential function, ovals and symmetries.

The quartic system (7.1) admits the following potential function:

𝐻 (𝑥, 𝑦) = (𝑥 2 -𝑋 0 ) 2 + (𝑦 2 -𝑌 0 ) 2 , (7.3)
together with the rescaling factor 𝜇 (𝑥, 𝑦) = 𝑦.

The level set associated to the parameter ℎ ⩾ 0, represented in the (𝑥 2 , 𝑦 2 ) plane, is the portion of the circle of center (𝑋 0 , 𝑌 0 ) and radius

√

ℎ located in the positive quadrant (see Figure 1a). In the (𝑥, 𝑦) plane, this corresponds to the ovals depicted in Figure 1b. Only those not intersecting the 𝑥-axis (over which 𝜇 vanishes) must be considered. We call them small and big ovals, respectively:

• When ℎ ∈ (0, 𝑋 2 0 ), the circle in the (𝑥 2 , 𝑦 2 )-plane entirely lies in the positive quadrant. In the (𝑥, 𝑦) plane, this results into four symmetric small ovals:

Γ ++ (ℎ) = 𝐻 -1 (ℎ) ∩ R ⩾0 × R ⩾0 , Γ +-(ℎ) = 𝐻 -1 (ℎ) ∩ R ⩾0 × R ⩽0 , Γ -+ (ℎ) = 𝐻 -1 (ℎ) ∩ R ⩽0 × R ⩾0 , Γ --(ℎ) = 𝐻 -1 (ℎ) ∩ R ⩽0 × R ⩽0 .
• When ℎ ∈ (𝑋 2 0 , 𝑌 2 0 ), a portion of this circle in the (𝑥 2 , 𝑦 2 )-plane crosses the 𝑦-axis, yielding two symmetric big ovals:

Γ + (ℎ) = 𝐻 -1 (ℎ) ∩ R × R ⩾0 , Γ -(ℎ) = 𝐻 -1 (ℎ) ∩ R × R ⩽0 .
In what follows, the notation Γ(ℎ) stands for the small oval Γ ++ (ℎ) when ℎ ∈ (0, 𝑋 2 0 ), and for the big oval Γ + (ℎ) when ℎ ∈ (𝑋 2 0 , 𝑌 2 0 ). The corresponding Abelian integral is:

ℑ(ℎ) = - ∫ Γ (ℎ)
𝑔(𝑥, 𝑦) 𝑦 d𝑥 .

By symmetry and the fact that 𝑔(𝑥, 𝑦) contains only even powers of 𝑥 and 𝑦, it is clear that the Abelian integrals taken along the other ovals of 𝐻 -1 (ℎ) are equal, up to the sign. Due to the very specific form of the potential function 𝐻 in (7.3), we can divide the ovals into four arcs 𝑥 ↦ → 𝑦 up (𝑥), 𝑥 ↦ → 𝑦 down (𝑥), 𝑦 ↦ → 𝑥 left (𝑦) and 𝑦 ↦ → 𝑥 right (𝑦), parameterized by explicit algebraic functions involving square roots only (see [Bréhard 2011, §6.2.3]). Therefore, explicit rigorous quadrature routines, e.g., the one provided by Arb, can be used to address this example.

Timings.

In Table 1, we present the execution times of our implementation and its three alternatives. The timings obtained with our method are comparable to those of the reference CAPD implementation, which is very encouraging for a future formal proof implementation in Coq.

Concerning the implementations of [START_REF] Johnson | A quartic system with twenty-six limit cycles[END_REF]] and in Arb/Sage, they both rely on an explicit parameterization of the ovals, which is possible only due to the very specific form of the potential function. While Arb/Sage exhibits remarkable timings due to the particularly efficient implementation of this library, the code of [START_REF] Johnson | A quartic system with twenty-six limit cycles[END_REF]] is limited by the binary64 (double precision) format and the finite order quadrature used to compute ℑ(ℎ). 1. Timings in seconds to rigorously evaluate ℑ(0.25) associated to (7.1) with relative error at most 10 -𝑝 , using our method with degree 𝑁 RTAs and floating-point FFT (𝑎) or interval FFT (𝑏) for RTA multiplication (see Remark 6.2), and other software: the CAPD library, the original code of [START_REF] Johnson | A quartic system with twenty-six limit cycles[END_REF]] and the Arb library.

0 1 2 0 1 2 X 0 , Y 0 X = x 2 Y = y 2 (a) Circles in the (𝑥 2 , 𝑦 2 ) plane -1 0 1 -1 0 Γ -+ Γ -- Γ + Γ - x y (b) Ovals in the (𝑥, 𝑦) plane --ℎ 0 -- 0 < ℎ < 𝑋 2 0 -- ℎ = 𝑋 2 0 --𝑋 2 0 < ℎ < 𝑌 2 0 --ℎ = 𝑌 2 0 --𝑌 2 0 < ℎ < 𝑋 2 0 + 𝑌 2 0 --ℎ = 𝑋 2 0 + 𝑌 2 0 --ℎ > 𝑋 2 0 + 𝑌 2 0 Fig.
7.1.3 Computer-assisted proof of H (4) ⩾ 24. In Table 2, we compute rigorous interval enclosures 𝑰 𝑁 for ℑ(ℎ), for specific values of ℎ, using our algorithm with RTAs of degree 𝑁 as small as possible, as long as the resulting interval guarantees the sign of ℑ(ℎ). The existence of sufficiently many simple zeros needed to prove Theorem 7.1 results from the obtained sign alternation. Figures 2 and3 provide a graphical representation of the sign alternations on small and big ovals. [-5.1962e-7, -5.2103e-8] [-2.8586e-7, -2.8585e-7] big 0.901 0.811801 197 [-1.2620e-5, -3.4042e-7] [-6.4798e-6, -6.4797e Proof of Theorem 7.1. Let ℎ 𝑖 (1 ⩽ 𝑖 ⩽ 9) denote the 9 values of ℎ in Table 2, taken in increasing order. The rigorous intervals enclosures given in this table are sufficient to guarantee the sign of ℑ(ℎ 𝑖 ). According to Theorem 1.3, there exists for each ℎ 𝑖 , an 𝜀 𝑖 > 0 such that 𝑑 (ℎ 𝑖 , 𝜀) and ℑ(ℎ 𝑖 ) share the same (strict) sign whenever 0 < 𝜀 ⩽ 𝜀 𝑖 . Hence, with 𝜀 * = min 0⩽𝑖⩽9 𝜀 𝑖 > 0, we have that the displacement function ℎ ↦ → 𝑑 (ℎ, 𝜀) alternates sign at least 5 times on (0, 𝑋 2 0 ) and at least 2 times on (𝑋 2 0 , 𝑌 2 0 ), for each fixed 0 < 𝜀 ⩽ 𝜀 * , giving respectively at least 5 and 2 isolated zeros in these intervals. Using the symmetries on the four small ovals and the two big ovals, we deduce the existence of at least 5 × 4 + 2 × 2 = 24 limit cycles in the quartic system (7.1) whenever 0 < 𝜀 ⩽ 𝜀 * . □ 7.2 Li, Liu and Yang's cubic system with 13 limit cycles In [START_REF] Li | A cubic system with thirteen limit cycles[END_REF]], C. Li, C. Liu and J. Yang showed that Z(3) ⩾ 13 using the following perturbed cubic Hamiltonian system: with well-chosen20 0 < 𝜆 < 1, 𝜅 > 1, and coefficients 𝛼 1 , 𝛼 2 , 𝛼 3 , 𝛼 4 ∈ R for the 𝜀-small perturbation. The critical points of the corresponding first integral,

𝑥 = -𝑦 3 + 𝜅 2 𝑦, 𝑦 = 𝑥 3 + (1 -𝜆)𝑥 2 -𝜆𝑥 + 𝜀𝑦 (𝛼 1 + 𝛼 2 𝑥 + 𝛼 3 𝑥 2 + 𝛼 4 𝑦 2 ),
𝐻 (𝑥, 𝑦) = 𝑥 4 4 + 1 -𝜆 3 𝑥 3 - 𝜆 2 𝑥 2 𝐹 (𝑥 ) + 𝑦 4 4 - 𝜅 2 𝑦 2 2 , (7.5)
depicted in Figure 4, are ranked according to their level value:

𝐻 (𝑃 3 ) < 𝐻 (𝑃 1 ) < 𝐻 (𝑃 2 ) < 𝐻 (𝑃 6 ) < 𝐻 (𝑃 4 ) < 𝐻 (𝑃 5 ). = = = = 𝐻 (𝑃 9 ) 𝐻 (𝑃 7 ) 𝐻 (𝑃 8 ) 0
The resulting families of ovals can be considered as a (horizontally) dissymmeterized version of the ovals of [START_REF] Johnson | A quartic system with twenty-six limit cycles[END_REF]] discussed in the previous section: Γ 2 (and symmetrically Γ 5 ) when ℎ ∈ (𝐻 (𝑃 3 ), 𝐻 (𝑃 2 )), Γ 1 (and sym. Γ 4 ) when ℎ ∈ (𝐻 (𝑃 1 ), 𝐻 (𝑃 2 )), Γ 3 (and sym. Γ 6 ) when ℎ ∈ (𝐻 (𝑃 2 ), 𝐻 (𝑃 6 )), Γ 8 when ℎ ∈ (𝐻 (𝑃 6 ), 𝐻 (𝑃 4 )), Γ 9 when ℎ ∈ (𝐻 (𝑃 4 ), 𝐻 (𝑃 5 )), and Γ 7 when ℎ > 𝐻 (𝑃 4 ).

The Abelian integrals under consideration are:

ℑ 𝑗 (ℎ) = ∫ Γ 𝑗 (ℎ) 𝑦 (𝛼 1 + 𝛼 2 𝑥 + 𝛼 3 𝑥 2 + 𝛼 4 𝑦 2 ) d𝑥, for 1 ⩽ 𝑗 ⩽ 9.
In Section 7.2.1, we first benchmark our algorithm on the more complex geometry of these ovals and compare the results with competing software. After that, we show in Section 7.2.2 how our approach based on numerical experiments and a posteriori validation provides an easy-to-check instance of Z(3) ⩾ 13 with effective values for the coefficients 𝜆, 𝜅, 𝛼 1 , 𝛼 2 , 𝛼 3 , 𝛼 4 in System (7.4). Then, over their respective interval of definition:

• ℑ 2 (and symmetrically ℑ 5 ) has at least one simple zero;

• ℑ 3 (and symmetrically ℑ 6 ) has at least 5 simple zeros;

• ℑ 7 has at least one simple zero.

As a result, the cubic Hamiltonian system (7.4) has at least 13 limit cycles.

Such an approach complements the detailed but harder-to-check existence proof for such coefficients in [START_REF] Li | A cubic system with thirteen limit cycles[END_REF]], based on the careful analysis of limit cycles for Liénard systems by F. Dumortier and C. Li [START_REF] Dumortier | Perturbation from an elliptic Hamiltonian of degree four. IV. Figure eight-loop[END_REF].

7.2.1 Timings. We benchmark our implementation on ovals Γ 3 (non symmetric, non convex) and Γ 8 (since it is not even star-shaped). We set 𝜆 = 7 10 and 𝜅 = 11 10 (the same values were used to plot Figure 4), and 𝛼 1 = 𝛼 2 = 𝛼 3 = 𝛼 4 = 1. Let 𝑄 1 (resp. 𝑄 2 ) denote the topmost intersection of the homoclinic connection of 𝑃 6 (resp. 𝑃 4 ) with the 𝑦-axis. Then, for Γ 3 , we consider the periodic orbit starting from 1 10 𝑃 2 + 9 10 𝑄 1 , and for Γ 8 , we take as initial point 1 2 𝑄 1 + 1 2 𝑄 2 (this corresponds exactly to the ovals Γ 3 and Γ 8 depicted in Figure 4).

The execution times in terms of the target relative accuracy for the Abelian integral are given in Table 3. They show that our method is particularly well-suited for high-precision evaluation of Abelian integrals, no matter the geometry of the oval. Note however that CAPD still performs better 20Note that it is assumed in [START_REF] Li | A cubic system with thirteen limit cycles[END_REF]] that 𝜅 > 10.

for low precision (up to standard double precision), probably due to its time interval subdivision strategy that we have not implemented yet in our method.

Finally, we also tried this example with the algorithm of [START_REF] Johnson | On a computer-aided approach to the computation of Abelian integrals[END_REF]. The resulting timings were not included in Table 3 since they already exceed 5 minutes for a targeted 4 digit accuracy. The reason is that the underlying oval subdivision strategy converges too slowly to provide more than a few digits in a reasonable time. Table 3. Timings in seconds to rigorously evaluate an Abelian integral along ovals from the Γ 3 and Γ 8 family with relative error at most 10 -𝑝 , using our implementation with degree 𝑁 RTAs and floating-point FFT (𝑎) or interval FFT (𝑏) for RTA multiplication (see Remark 6.2), and CAPD.

(a) Ovals of potential function 𝐻 (7.6).

(b) Wronskian 𝑊 3 of system { ℑ 3𝑖 } 1⩽𝑖⩽4 (7.7) and a linear combination ℑ 3 with 𝜆 = 0.17 and appropriate 𝛽 𝑖 realizing 5 zeros { r𝑖 } 1⩽𝑖⩽5 .

Fig. 5. Obtaining numerical 𝛽 𝑖 realizing 5 zeros for ℑ 3 in Dumortier and Li's system [START_REF] Dumortier | Perturbation from an elliptic Hamiltonian of degree four. IV. Figure eight-loop[END_REF].

7.2.2

Computer-assisted proof of Z(3) ⩾ 13. Before proving Theorem 7.2, we present a heuristic process, directly inspired by [START_REF] Li | A cubic system with thirteen limit cycles[END_REF], that yields potential relevant parameters. To do so, we focus on ℑ 2 and ℑ 3 (once the parameters fixed, the result on ℑ 7 will be straightforward to establish 4. Rigorous evaluations of ℑ 2 (ℎ), ℑ 3 (ℎ) and ℑ 7 (ℎ) with the parameters of Theorem 7.2 using AbIntValid. The ovals Γ 2 (ℎ), Γ 3 (ℎ) and Γ 7 (ℎ) are parameterized with a variable 𝑟 denoting a position on a transversal, with ℎ = 𝐻 (-1 + 𝑟, 𝜅), ℎ = 𝐻 (0, 𝜅 + 𝑟 ), and ℎ = 𝐻 (𝑟, 0), respectively. 𝑁 is the degree used for RTAs, and 𝛿 is the diameter of the computed interval enclosure for ℑ(ℎ).

with our approach). In this case, [START_REF] Li | A cubic system with thirteen limit cycles[END_REF]] uses the change of variable ỹ = 𝑦 2 -𝜅 2 2 to obtain a simpler potential function, 𝐻 (𝑥, ỹ) = ỹ2 + 𝐹 (𝑥), (7.6) whose families of ovals Γ 1 , Γ 2 and Γ 3 are depicted in Figure 5a. The corresponding Abelian integrals ℑ 𝑗 for 1 ⩽ 𝑗 ⩽ 3 are:

ℑ 𝑗 (ℎ) = ∫ Γ 𝑗 (ℎ) ỹ (𝛽 1 + 𝛽 2 𝑥 + 𝛽 3 𝑥 2 + 𝛽 4 ỹ2 ) d𝑥 = 𝛽 1 ℑ 𝑗1 (ℎ) + 𝛽 2 ℑ 𝑗2 (ℎ) + 𝛽 3 ℑ 𝑗3 (ℎ) + 𝛽 4 ℑ 𝑗4 (ℎ). (7.7)
The authors of [START_REF] Li | A cubic system with thirteen limit cycles[END_REF]] establish, as soon as 𝜅 is large enough, an equivalence between the numbers of zeros of ℑ 2 and ℑ 2 , resp. ℑ 3 and ℑ 3 .

Zeros of ℑ 2 and ℑ 3 . By combining the main existence result of [START_REF] Dumortier | Perturbation from an elliptic Hamiltonian of degree four. IV. Figure eight-loop[END_REF] with some asymptotic analysis, C. Li, C. Liu and J. Yang prove [START_REF] Li | A cubic system with thirteen limit cycles[END_REF], Lemma 1] that if 𝜆 is close to a certain critical value 𝜆 * , then there exist coefficients 𝛽 1 , 𝛽 2 , 𝛽 3 , 𝛽 4 such that ℑ 2 has at least one simple zero and ℑ 3 has at least 5 simple zeros. The reason why 5 zeros can be obtained from the 4-term linear combination ℑ 3 is that for this critical value 𝜆 * , the Wronskian 𝑊 3

(ℎ) = det ℑ (𝑖 -1) 3𝑗 (ℎ) 1⩽𝑖,𝑗⩽4
has a double zero h0 = r 2 0 (see Figure 5b). Then, an appropriate perturbation of 𝜆 around 𝜆 * creates two simple zeros for 𝑊 (ℎ), from which one hopes to deduce the existence of values for the 𝛽 𝑖 realizing 5 zeros21.

In contrast with the non-constructive existence proof in [START_REF] Li | A cubic system with thirteen limit cycles[END_REF]], we confirm this heuristic in a constructive way, by determining candidate values for 𝜆 and the 𝛽 𝑖 purely numerically. To do so, we compute Chebyshev approximations of the Wronskian 𝑟 ↦ → 𝑊 3 (𝑟 2 ) 22 for several values of 𝜆 and count the number of zeros. This is done by interpolating the Abelian integrals on a large number of Chebyshev nodes, with very high precision for each value23. The value 𝜆 = 0.17 is found to fulfill those requirements, that is 2 zeros h1 = r 2 1 < r 2 2 = h2 for 𝑊 3 (at least numerically, see Figure 5b), with r1 and r2 not too close for better numerical stability. Now, to obtain values for the 𝛽 𝑖 , we fix 𝛽 1 = 1, and numerically solve for 𝛽 2 , 𝛽 3 and 𝛽 4 the equations ℑ 3 ( r 2 𝑖 ) = 0 for r1 , r2 and the midpoint r3 = r1 + r2 2 . As desired, the numerical plot of the resulting combination 𝑟 ↦ → ℑ 3 (𝑟 ) in Figure 5b exhibits 2 additional zeros r4 and r5 , making a total of 5 zeros for ℑ 3 . Finally, a numerical plot suggests that ℑ 2 also has one simple zero.

Zeros of ℑ 2 , ℑ 3 and ℑ 7 . As claimed by [Li et al. 2009, Lemma 2], by selecting a sufficiently large parameter 𝜅 > 1 in (7.4) and defining the coefficients 𝛼 𝑖 from the 𝛽 𝑖 using the following formulas:

𝛼 1 = 𝛽 1 𝜅 5 + 3𝛽 4 2𝜅 , 𝛼 2 = 𝛽 2 𝜅 5 , 𝛼 3 = 𝛽 3 𝜅 5 , 𝛼 4 = - 𝛽 4 2𝜅 3 ,
the 5 zeros of ℑ 3 and the single zero of ℑ 2 are recovered in ℑ 3 and ℑ 2 . Moreover, with these coefficients, an additional zero can be found for ℑ 7 , according to [Li et al. 2009, Lemma 3].

Fixing 𝜅 = 20 (based on some experimenting) leads to the desired number of zeros (see plots in Figure 6), while being not too large, to avoid numerical issues. For the final values of the 𝛼 𝑖 , given in Theorem 7.2 and used for the computer-assisted proof, no less than 13 digits are necessary in the decimal truncation to keep the expected number of zeros.

We now conclude with a rigorous numerics based proof of the desired sign alternations of the Abelian integrals, hence of Z(3) ⩾ 13.

Proof of Theorem 7.2. We perform rigorous evaluations of Abelian integrals ℑ 2 , ℑ 3 and ℑ 7 on the points given in Table 4 using Algorithm AbIntValid, with the values of 𝜆, 𝜅 and the 𝛼 𝑖 given in the statement of the theorem. By continuity of these three functions, it follows that ℑ 2 , ℑ 3 and ℑ 7 have at least 1, 5 and 1 simple zeros, respectively. Moreover, by the symmetry (𝑥, 𝑦) ↦ → (𝑥, -𝑦), ℑ 5 and ℑ 6 have at least 1 and 5 simple zeros, respectively. Finally, by the same argument used in the proof of Theorem 7.1 and based on Theorem 1.3 (Poincaré-Pontryagin), we obtain that for sufficiently small 𝜀 > 0, System (7.4) has at least 13 limit cycles. □ From a performance point of view, it is noteworthy that rather large degrees for RTAs are necessary for ℑ 2 and ℑ 3 (see Table 4), mainly due to the winding number check in OvalValid since the denominator in line 22 becomes small when the curve gets close to the point (𝑥 𝑟 , 𝑦 𝑟 ). Alternative methods for this should be investigated in the future to avoid this phenomenon.

21Indeed, the case of a single (simple) zero for the Wronskian (together with additional non-vanishing conditions on lowerdimensional Wronskians of the system) would imply that ℑ 3 cannot exhibit more than 4 zeros (see [Novaes and Torregrosa 2017, Cor. 1.4]). 22For convenience, we use 𝑟 = √ ℎ instead of ℎ as the independent variable. 23For this we use a plain floating-point evaluation scheme, adapted from AbIntValid by removing all the rigorous numerics for more efficiency, since error bounds are not necessary at this stage.

, Vol. 1, No. 1, Article . Publication date: December 2023.

Conclusion

We have presented a set of algorithms that makes it possible to evaluate Abelian integrals in a fast and rigorous way. Next steps are the development of a rigorous FFT following [START_REF] Brisebarre | Error Analysis of Some Operations Involved in the Cooley-Tukey Fast Fourier Transform[END_REF], cf. (a) in Remark 6.2, and the formalization in Coq of our method. Our hope is that, in particular, it will be used in the dynamical system community as a basic brick in the quest for new records for the H (𝑛) and Z(𝑛) bounds.

A ALGORITHM ISPOSITIVE: CHECKING POSITIVITY OF RTAS

In rigorous numerics, it is common to determine if a function 𝑓 is positive over a given interval using interval subdivision and branch and bound techniques [Tucker 2011, Chap. 5.1]. In our case, however, evaluating RTAs of high degree over non-thin intervals yields extremely large overapproximations, making this process ineffective for our purpose.

The following routine IsPositive is based on an alternative approach. The RTA 𝒇 = ( f , 𝜀) is preconditioned by an approximate inverse ǧ ≈ 1 f . If 𝒇 does not vanish, then the product ǧ𝒇 is likely to be close to 1, which is measured by the quantity ∥1 -ǧ𝒇 ∥ ℓ 1 . (𝑖𝑖) If moreover f > 0 and 𝜀 < ∥ f -1 ∥ -1 ℓ 1 (which implies 𝑓 > 0 for all 𝑓 ∈ 𝒇 ), then there is an 𝑁 0 such that IsPositive(𝒇, 𝑁 ) returns true for all 𝑁 ⩾ 𝑁 0 . (𝑖𝑖𝑖) IsPositive(𝒇, 𝑁 ) runs in 𝑂 (𝑁 ′ log(𝑁 ′ )) arithmetic operations, where 𝑁 ′ = max(deg 𝒇, 𝑁 ).

Proof. (𝑖) Soundness. Assume IsPositive(𝒇, 𝑁 ) returns true. A first possibility is that the condition in line 4 is true, implying that 𝑎 0 > 0 and that any function 𝑓 ∈ 𝒇 can be written as 𝑎 0 + ℎ with ∥ℎ∥ ∞ < 𝑎 0 , and is therefore positive. The other possibility (line 8) is that we obtained a ǧ such that for any 𝑓 ∈ 𝒇 , ∥1 -ǧ𝑓 ∥ ∞ < 1, implying that (𝑡) ≠ 0 for all 𝑡. By continuity, 𝑓 has constant sign over R. In particular, this is true for f . Since the condition in line 2 was false, i.e., 𝑎 0 = ∫ 2𝜋 0 f d𝑡 > 0,

we get f > 0. Finally, since the ball in 𝐶 0 2𝜋 denoted by 𝒇 is convex, the same statement 𝑓 > 0 is true for all 𝑓 ∈ 𝒇 . (𝑖𝑖) Completeness. Suppose conversely that f > 0 and that 𝜀 < ∥ f -1 ∥ -1 ℓ 1 . In particular, 𝜀 < ∥ f -1 ∥ -1 ∞ = min 𝑡 ∈R f (𝑡), implying that 𝑓 > 0 for all 𝑓 ∈ 𝒇 = ( f , 𝜀). We prove that IsPositive(𝒇, 𝑁 ) returns true for 𝑁 sufficiently large.

First, the condition in line 2 is false since 𝑎 0 = ∫ 2𝜋 0 f d𝑡 > 0. Then, the algorithm either returns true if the condition in line 4 is true, or it computes ǧ ← I 𝑁 ( f -1 ) in line 7. Since 1/ f is analytic, ∥ ǧ -f -1 ∥ ℓ 1 converges to 0 as the approximation degree 𝑁 tends to infinity. Hence,

Bound(1 -ǧ𝒇 ) = ∥1 -ǧ f ∥ ℓ 1 + 𝜀Bound( ǧ) → 0 + 𝜀 ∥ f -1 ∥ ℓ 1 < 1 as 𝑁 → ∞.
Therefore, there is an 𝑁 0 such that for all 𝑁 ⩾ 𝑁 0 , the degree 𝑁 TP ǧ computed by trigonometric interpolation makes the condition in line 8 true, so that IsPositive(𝒇, 𝑁 ) returns true. that serves as a shift to project any point (𝑥, 𝑦) onto Γ(ℎ) in the direction given by (𝑢, 𝑣). Then, we prove that over a possibly smaller neighborhood 𝑈 ′ ⊆ 𝑈 of Γ(ℎ), the Newton iterations N 𝑘 𝑥,𝑦,𝑢,𝑣 (0) (see Eq. satisfies 𝜕𝐹 𝜕𝑠 (𝑥, 𝑦, 𝑢, 𝑣, 0) = ∥∇𝐻 (𝑥, 𝑦) ∥ 2 > 0 for all (𝑥, 𝑦, 𝑢, 𝑣) ∈ Γ(ℎ), since ℎ is a regular value of 𝐻 . By the analytic implicit function theorem [Cartan 1995, Chap. IV, Prop. 6.1], there exists a neighborhood 𝑈 ⊆ R 4 of Γ(ℎ), a neighborhood 𝑉 ⊆ R 5 of Γ(ℎ) × {0}, and a unique analytic function ŝ : 𝑈 → R such that:

((𝑥, 𝑦, 𝑢, 𝑣, 𝑠) ∈ 𝑉 ∧ 𝐹 (𝑥, 𝑦, 𝑢, 𝑣, 𝑠) = 0) ⇔ ((𝑥, 𝑦, 𝑢, 𝑣) ∈ 𝑈 ∧ 𝑠 = ŝ (𝑥, 𝑦, 𝑢, 𝑣)) .

In particular, ŝ (𝑥, 𝑦, 𝑢, 𝑣) = 0 if and only if (𝑥, 𝑦) ∈ Γ(ℎ). • Convergence of Newton iterations. Since Γ(ℎ) is compact, we may assume (up to restricting the neighborhoods 𝑈 and 𝑉 ) that 𝜕𝐹 𝜕𝑠 (𝑥, 𝑦, 𝑢, 𝑣, 𝑠) ⩾ 𝑚 1 for some 𝑚 1 > 0, and | 𝜕 2 𝐹 𝜕𝑠 2 (𝑥, 𝑦, 𝑢, 𝑣, 𝑠)| ⩽ 𝑀 2 for some 𝑀 2 > 0, for all (𝑥, 𝑦, 𝑢, 𝑣, 𝑠) ∈ 𝑉 . Now for any (𝑥, 𝑦, 𝑢, 𝑣, 𝑠) ∈ 𝑉 , the classical error analysis of Newton's method [Tucker 2011, §5.1.2] . By the property of composition of analytic functions, the map 𝑡 ↦ → 𝑠 # (𝑡) = ŝ (𝑥 • (𝑡), 𝑦 • (𝑡), 𝑢 • (𝑡), 𝑣 • (𝑡)) is analytic, and so is (𝑥 # , 𝑦 # ) by Eq. (3.4). Also, by choosing 𝜂 > 0 sufficiently small, we can make (𝑥 # , 𝑦 # ) sufficiently close to (𝑥 * , 𝑦 * ), which ensures that the winding number with respect to a given point inside Γ(ℎ) remains the same (±1 depending on the orientation). This concludes the proof that (𝑥 # , 𝑦 # ) is an analytic parameterization of Γ(ℎ) with same orientation as (𝑥 * , 𝑦 * ), although not equal to (𝑥 * , 𝑦 * ) in general. Finally, the quadratic convergence of the Newton iteration scheme for any point (𝑥 • (𝑡), 𝑦 • (𝑡)) w.r.t. direction (𝑢 • (𝑡), 𝑣 • (𝑡)) follows from the preceding analysis. □

C ALGORITHM NEWTONBALL

In order to make effective the a posteriori validation process presented in Section 2.3, we must be able to automatically determine a radius 𝑟 and a Lipschitz constant 𝜆 satisfying the hypotheses of Theorem 2.5. To this aim, given a function 𝑟 ↦ → 𝜆(𝑟 ) providing a guaranteed upper bound for the Lipschitz constant for the fixed-point operator T over 𝐵( φ, 𝑟 ), and an upper bound for the defect 𝑑 ⩾ ∥T ( φ)φ ∥, Algorithm NewtonBall implements a simple bisection method to identify the smallest zero of 𝑟 ↦ → 𝑓 (𝑟 ) := (1 -𝜆(𝑟 ))𝑟 -𝑑. Specifically, it maintains a stack of subintervals [𝑟 -, 𝑟 + ] to be investigated, each of them satisying 𝑓 (𝑟 -) ⩽ 0 . Under the reasonable assumption that 𝜆(𝑟 ) is an increasing function of 𝑟 , the range of 𝑓 over [𝑟 -, 𝑟 + ] can be overapproximated by [(1 -𝜆 + )𝑟 --𝑑, (1 -𝜆 -)𝑟 + -𝑑] where 𝜆 ± := 𝜆(𝑟 ± ), thus giving a simple test (line 7) to detect the possible existence of a zero of 𝑓 in [𝑟 -, 𝑟 + ]. If moreover 𝑓 (𝑟 + ) ⩾ 0, then [𝑟 -, 𝑟 + ] necessarily contains a zero, and the stack of remaining intervals is cleared (line 8). When the process successfully terminates, it returns a pair (𝑟 + , 𝜆 + ) satisfying the hypotheses of Theorem 2.5.

Remark C.1. Algorithm NewtonBall ensures a slightly stronger property than Theorem 2.5, namely the strict inequality 𝑑 + 𝜆 + 𝑟 + < 𝑟 + . By doing so, we make sure that this inequality remains valid under small perturbations of 𝜆(𝑟 ), which turns out to be essential in the proof of Proposition 4.1 to guarantee that the unique fixed point is not only continuous, but also analytic.

  𝑦) 𝜕𝐻 𝜕𝑦 (𝑥, 𝑦) + 𝜀𝑃 (𝑥, 𝑦), 𝑦 = 𝜇 (𝑥, 𝑦) 𝜕𝐻 𝜕𝑥 (𝑥, 𝑦) + 𝜀𝑄 (𝑥, 𝑦).

  8This theorem provides convergence bounds w.r.t. the ∥ • ∥ ∞ norm, but the same bounds actually hold for the ∥ • ∥ ℓ 1 norm since the proof consists in bounding the Fourier coefficients., Vol. 1, No. 1, Article . Publication date: December 2023.

𝑗=0 6 :

 6 return (𝑥 • , 𝑦 • ) , Vol. 1, No. 1, Article . Publication date: December 2023.

Fig. 2 .

 2 Fig.2. Rigorously computed interval enclosures of ℑ(ℎ) on small ovals, for ℎ ∈ {0.25, 0.6084, 0.7744, 0.7921, 0.801025, 0.80766169}, proving the existence of 5 zeros. The right plot is a zoom of the blue box in the left plot.

Fig. 3 .

 3 Fig.3. Rigorously computed interval enclosures of ℑ(ℎ) on big ovals, for ℎ ∈ {0.811801, 0.8649, 0.9025}, proving the existence of 2 zeros. The right plot is a zoom of the green box in the left plot.

Fig. 4 .

 4 Fig. 4. Ovals of potential function 𝐻 (7.5), with 𝜆 = 0.7 and 𝜅 = 1.1.

  No. 1, Article . Publication date: December 2023.

  Theorem 7.2. Consider System (7.4) with parameters: 𝜆 = 0.17, 𝛼 1 = 8.808855593098, 𝛼 3 = -2.599597169555 • 10 -5 , 𝜅 = 20, 𝛼 2 = -2.078279433211 • 10 -5 , 𝛼 4 = -7.340712733831 • 10 -3 .

Algorithm 7 4 :

 74 IsPositive(𝒇, 𝑁 )Input: RTA 𝒇 = ( f , 𝜀) and approximation degree 𝑁 Output: Boolean true or false, with true guaranteeing that 𝑓 > 0 for all 𝑓 ∈ 𝒇1: 𝑎 0 ← constant coefficient of f 2: if 𝑎 0 ⩽ 0 then else if Bound(𝒇 -𝑎 0 ) < 𝑎 0 then 𝑁 (( 1 𝑓 𝑗 ) 2𝑁 𝑗=0 ), where (𝑓 𝑗 ) 2𝑁 𝑗=0 ← IFFT 𝑁 ( f ) 8:if ǧ is well defined and Bound(1 -ǧ𝒇 )Lemma A.1. Let 𝒇 = ( f , 𝜀) be an RTA. Then, (𝑖) For any approximation degree 𝑁 , if IsPositive(𝒇, 𝑁 ) returns true, then 𝑓 > 0 for all 𝑓 ∈ 𝐶 0 2𝜋 represented by 𝒇 .

  (𝑖𝑖𝑖) Complexity. The asymptotic complexity is determined by the FFT 𝑁 /IFFT 𝑁 routines used for evaluation on the equispaced grid, interpolation and TP/RTA multiplication. □ B PROOF OF LEMMA 3.4 FOR THE PROJECTED PARAMETERIZATION Proof. The proof consists of three steps. First, we prove the existence of an analytic function ŝ : 𝑈 → R over a neighborhood 𝑈 ⊆ R 4 of Γ(ℎ) = (𝑥, 𝑦, 𝑢, 𝑣) ∈ R 4 s.t. (𝑥, 𝑦) ∈ Γ(ℎ) and (𝑢, 𝑣) = ∇𝐻 (𝑥, 𝑦) ∥∇𝐻 (𝑥, 𝑦) ∥ 2 ,

  (3.5)) converge quadratically fast to ŝ (𝑥, 𝑦, 𝑢, 𝑣) for all (𝑥, 𝑦, 𝑢, 𝑣) ∈ 𝑈 ′ . Finally, we compose by 𝑡 ∈ [0, 2𝜋] → (𝑥 • (𝑡), 𝑦 • (𝑡)) to complete the proof. • Analytic projection. Consider the equation: 𝐹 (𝑥, 𝑦, 𝑢, 𝑣, 𝑠) = 𝐻 (𝑥 + 𝑠𝑢, 𝑦 + 𝑠𝑣)ℎ = 0. By definition, 𝐹 (𝑥, 𝑦, 𝑢, 𝑣, 0) = 0 for all (𝑥, 𝑦, 𝑢, 𝑣) ∈ Γ(ℎ). Moreover, the derivative with respect to 𝑠, 𝜕𝐹 𝜕𝑠 (𝑥, 𝑦, 𝑢, 𝑣, 𝑠) = 𝑢 𝜕𝐻 𝜕𝑥 (𝑥 + 𝑠𝑢, 𝑦 + 𝑠𝑣) + 𝑣 𝜕𝐻 𝜕𝑦 (𝑥 + 𝑠𝑢, 𝑦 + 𝑠𝑣)

,

  Vol. 1, No. 1, Article . Publication date: December 2023.

,

  yields in our case:|N 𝑥,𝑦,𝑢,𝑣 (𝑠)ŝ (𝑥, 𝑦, 𝑢, 𝑣)| ⩽ 𝑀 2 2𝑚 1 |𝑠ŝ (𝑥, 𝑦, 𝑢, 𝑣)| 2 .Therefore, by restricting the neighborhood 𝑈 of Γ(ℎ) ⊆ ŝ-1 (0) to:𝑈 ′ = (𝑥, 𝑦, 𝑢, 𝑣) ∈ 𝑈 , | ŝ (𝑥, 𝑦, 𝑢, 𝑣)| ⩽ 𝑚 1 𝑀 2we deduce that the Newton iterations N 𝑘 𝑥,𝑦,𝑢,𝑣 (0) converge quadratically fast to the desired limit ŝ (𝑥, 𝑦, 𝑢, 𝑣):|N 𝑘 𝑥,𝑦,𝑢,𝑣 (0)ŝ (𝑥, 𝑦, 𝑢, 𝑣)| ⩽ 𝑀 2 2𝑚 1 | ŝ (𝑥, 𝑦, 𝑢, 𝑣)| 2 𝑘 -1 | ŝ (𝑥, 𝑦, 𝑢, 𝑣)| ⩽ | ŝ (𝑥, 𝑦, 𝑢, 𝑣)| 2 2 𝑘 -1 . (B.1)• Composition with the initial guess. Let us fix 𝜂 > 0 such that any (𝑥 • , 𝑦 • ) and (𝑢 • , 𝑣 • ) 𝜂-close to (𝑥 * , 𝑦 * ) and ∇𝐻 (𝑥 * ,𝑦 * ) ∥ ∇𝐻 (𝑥 * ,𝑦 * ) ∥ 2 in (𝐶 0 2𝜋 ) 2 satisfy (𝑥 • (𝑡), 𝑦 • (𝑡), 𝑢 • (𝑡), 𝑣 • (𝑡)) ∈ 𝑈 ′ for all 𝑡 ∈ [0, 2𝜋]

Algorithm 8

 8 NewtonBall(𝑑, 𝜆) Input: defect 𝑑 ⩾ ∥T ( φ)φ ∥ and procedure 𝜆 s.t. T is 𝜆(𝑟 )-Lipschitz over 𝐵( φ, 𝑟 ) Output: (𝑟 + , 𝜆 + ) s.t. T is 𝜆 + -contracting over 𝐵( φ, 𝑟 + ) and 𝑑 + 𝜆 + 𝑟 + < 𝑟 + 1: compute 𝑟 max s.t. 𝜆(𝑟 max ) ≈ 0.5 2: push ((0, 𝜆(0)), (𝑟 max , 𝜆(𝑟 max ))) on a new empty stack 𝑆 ⊲ Bisection loop 3: while 𝑆 not empty do 4: ((𝑟 -, 𝜆 -), (𝑟 + , 𝜆 + )) ← pop(𝑆) ⊲ When interval subdivision is below user-defined 𝛿 𝜆 5: if 𝜆 + -𝜆 -⩽ 𝛿 𝜆 and 𝑑 + 𝜆 + 𝑟 + < 𝑟 + then 6: return (𝑟 + , 𝜆 + ) ⊲ Check whether interval [𝑟 -, 𝑟 + ] may contain a zero 7: else if 𝜆 + -𝜆 -> 𝛿 𝜆 and 𝑑 + 𝜆 -𝑟 + < 𝑟 + then 8: if 𝑑 + 𝜆 + 𝑟 + ⩽ 𝑟 + then clear(𝑆) end if ⊲ [𝑟 -, 𝑟 + ] necessarily contains a zero ⊲ Subdivide [𝑟 -, 𝑟 + ] and push resulting intervals on the stack 9: 𝑟 𝑚 ← (𝑟 -+ 𝑟 + )/2 10: 𝜆 𝑚 ← 𝜆(𝑟 𝑚 ) 11: push(((𝑟 𝑚 , 𝜆 𝑚 ), (𝑟 + , 𝜆 + )), 𝑆) 12:push(((𝑟 -, 𝜆 -), (𝑟 𝑚 , 𝜆 𝑚 )), 𝑆) No. 1, Article . Publication date: December 2023.

  Potential function 𝐻 , level ℎ, initial approximation (𝑥 • , 𝑦 • ) and degree 𝑁 ∈ N Output: Degree 𝑁 TPs 𝑢 • , 𝑣 • , š defining the approximate parameterization ( x, y) of Γ(ℎ) as in Eq. (3.6)

• , 𝑦 • , 𝑁 ) Input:

  1] 11: function 𝜆(𝑟 ) ⊲ Lipschitz constant of N over 𝐵( š, 𝑟 ) in 𝐶 0 𝑋 𝑟 ← 𝑋 + 𝑟𝑈 and 𝑌 𝑟 ← 𝑌 + 𝑟𝑉 𝜆 1 ← mag 𝑈 2 𝜕 2 𝐻 𝜕𝑥 2 (𝑋 𝑟 , 𝑌 𝑟 ) + 2𝑈𝑉 𝜕 2 𝐻 𝜕𝑥 𝜕𝑦 (𝑋 𝑟 , 𝑌 𝑟 ) + 𝑉 2 𝜕 2 𝐻 𝜕𝑦 2 (𝑋 𝑟 , 𝑌 𝑟 )

		2𝜋
	12:
	13:
	14:	return 𝜆 0 + 𝛼𝑟𝜆 1
	15: end function
	16: (𝑟

+ , 𝜆 + ) ← NewtonBall(𝑑, 𝜆)

then return 𝐹𝐴𝐼𝐿 end if

  𝜀 0 ) enclosing initial conditions such that it intersects 𝐻 -1 (ℎ) in the Γ(ℎ) component only, Algorithm OvalValid(𝐻, ℎ, 𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 , 𝜀 0 , 𝑥 𝑟 , 𝑦 𝑟 , 𝑥 • , 𝑦 • , 𝑢 • , 𝑣 • , š, 𝑁 ), if it does not fail, computes 𝜀 such that 𝑠 # of Definition 3.3 exists, satisfies Lemma 3.4 with winding number around (𝑥 𝑟 , 𝑦 𝑟 ) equal to 𝜍, and ∥

24: return 𝜀 Proposition 4.1. Given a point (𝑥 𝑟 , 𝑦 𝑟 ) inside the oval Γ(ℎ), and a ball 𝐵((𝑥 𝑖𝑛𝑖 , 𝑦 𝑖𝑛𝑖 ),

  1. Level curves of the potential function 𝐻

			our method			other methods
	𝑝	𝑁	(𝑎)	(𝑏)	CAPD [Johnson 2011] Arb/Sage
	4	55 0.32 1.6	0.44	7.9	0.012
	6	70 0.59 3.3	0.54	19.0	0.013
	8	100 0.60 4.8	0.67	49.4	0.014
	16	150 1.5	9.0	2.0	-	0.018
	24	225 1.5	9.1	2.2	-	0.021
	32	275 3.8 18.7	3.6	-	0.031
	48	425 3.8 18.9	6.4	-	0.040
	64	570 8.0 42.1	11.1	-	0.043
	96	870 9.4 46.1	23.8	-	0.078
	128	1165 18.8 96.2	41.7	-	0.084
	Table

Table 2 .

 2 Rigorous evaluation of ℑ(ℎ) with our algorithm and resulting sign alternations on small and big ovals. Computations are carried out with 256 bits of floating-point precision and the smallest degree 𝑁 for RTAs such that the rigorous enclosure 𝑰 𝑁 guarantees the sign of ℑ(ℎ). Tighter enclosures 𝑰 512 using a high degree 512 are also provided.

	-6]	-

Table

  Fig.6. Abelian integrals ℑ 2 , ℑ 3 , ℑ 7 in Li, Liu and Yang's system (7.4) with the parameters of Theorem 7.2, having 1, 5 and 1 zeros, respectively.

	oval	𝑟	𝑁	ℑ(ℎ)	𝛿	sign
	Γ 2	0.88	2000 [-9.54e-10, -9.53e-10] 5.45e-22	-
		0.92	2000 [5.55e-10, 5.56e-10] 3.92e-21	+
	Γ 3	0.0002 8000 [-5.21e-13, -5.20e-13] 4.16e-22	-
		0.001 8000 [8.45e-13, 8.46e-13] 3.24e-26	+
		0.0015 4000 [-2.16e-13, -2.15e-13] 5.11e-21	-
		0.0022 4000 [3.52e-13, 3.53e-13] 2.27e-22	+
		0.004 4000 [-6.13e-12, -6.12e-12] 6.76e-25	-
		0.0045 2000 [2.54e-12, 2.55e-12] 1.86e-20	+
	Γ 7	25	100	[2718.12, 2718.13]	3.63e-6	+
		30	100	[-3419.83, -3419.82]	6.00e-7	-

Florent Bréhard, Nicolas Brisebarre, Mioara Joldes, and Warwick Tucker

, Vol. 1, No. 1, Article . Publication date: December 2023.

√ X 0 , √ Y 0 Γ ++ Γ +-

ACKNOWLEDGMENTS

We wish to warmly thank J.-L. Figueras, G. Hanrot, S. Olver, R. M. Slevinsky and D. Wilczak for helpful discussions and advise. The first author wishes to thank Uppsala University and the CAPA group where part of this work was done. The fourth author wishes to thank Collegium de Lyon and the INRIA international chair programme for hosting and funding him during the time spent in Lyon while working on this paper.

This work was partially funded by the NuSCAP (ANR-20-CE48-0014) project of the French national research agency (ANR), by the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR), and by the Swedish Research Council (VR-2018-04265).