# Efficient and Validated Numerical Evaluation of Abelian Integrals

Florent Bréhard, Nicolas Brisebarre, Mioara Joldeş, Warwick Tucker

► To cite this version:

Florent Bréhard, Nicolas Brisebarre, Mioara Joldeş, Warwick Tucker. Efficient and Validated Numerical Evaluation of Abelian Integrals. ACM Transactions on Mathematical Software, In press, 10.1145/3637550. [hal-03561096v2](https://hal.science/hal-03561096v2)

HAL Id: hal-03561096

https://hal.science/hal-03561096v2

Submitted on 12 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Efficient and Validated Numerical Evaluation of Abelian Integrals

FLORENT BRÉHARD, Univ. Lille, CNRS, Centrale Lille, UMR 9189 CRISTAL, France

NICOLAS BRISEBARRE, CNRS, LIP, Inria AriC, Université de Lyon, France

MIOARA JOLDES, LAAS-CNRS, Université de Toulouse, CNRS, France

WARWICK TUCKER, School of Mathematics, Monash University, Australia

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.

CCS Concepts: • Mathematics of computing → Quadrature; Ordinary differential equations; Approximation; Interval arithmetic; Interpolation; Arbitrary-precision arithmetic.

Additional Key Words and Phrases: Hilbert's 16th problem, Abelian integral, limit cycles, rigorous numerics, trigonometric polynomial interpolation, Newton-like operator

**ACM Reference Format:**

Florent Bréhard, Nicolas Brisebarre, Mioara Joldes, and Warwick Tucker. 2023. Efficient and Validated Numerical Evaluation of Abelian Integrals. 1, 1 (December 2023), 38 pages. https://doi.org/10.1145/nnnnnnn.nnnnnnn

# 1 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 [Lanford 1982], the existence of chaos in the Kuramoto–Sivashinsky equations [Wilczak 2003], and the (almost) finiteness of relative equilibria for the 5-body problem of celestial mechanics [Albouy and Kaloshin 2012].

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

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).

Authors' addresses: Florent Bréhard, Univ. Lille, CNRS, Centrale Lille, UMR 9189 CRISTAL, F-59000, Lille, France, florent.brehard@univ-lille.fr; Nicolas Brisebarre, CNRS, LIP, Inria AriC, Université de Lyon, Lyon, France, nicolas.brisebarre@ens-lyon.fr; Mioara Joldes, LAAS-CNRS, Université de Toulouse, CNRS, Toulouse, France, joldes@laas.fr; Warwick Tucker, School of Mathematics, Monash University, VIC 3800, Australia, warwick.tucker@monash.edu.
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 $H : \mathbb{R}^2 \to \mathbb{R}$, a rational rescaling factor $\mu : \mathbb{R}^2 \to \mathbb{R}$, a polynomial perturbation $(P, Q) : \mathbb{R}^2 \to \mathbb{R}^2$ and a regular 1 level value $h$, compute the Abelian integral along the oval $\Gamma(h)$ in the plane:

$$ \mathfrak{I}(h) = \int_{\Gamma(h)} \frac{P(x, y)dy - Q(x, y)dx}{\mu(x, y)}, \qquad (1.1) $$

where an oval $\Gamma(h)$ denotes a bounded connected component of the level set $H(x, y) = h$ for a regular value $h^2$.

We propose a validated numerical algorithm, which, given polynomial or rational functions $H, \mu, P, Q$, together with the level $h$, computes an interval that contains the exact value of $\mathfrak{I}(h)$. 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 $\mathfrak{I}(h)$ – 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 implementation<fnref n="3" /> in Julia<fnref n="4" /> [Bezanson et al. 2017]. 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 [Bertot and Castéran 2013] in the near future. The end goal is to provide an efficient and fully certified calculator for the infinitesimal Hilbert's 16th problem.

## 1.1 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 [Hilbert 1900] 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 $\mathcal{H}(n)$ – 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 $n$ can display. Note that the bound $\mathcal{H}(n)$ should be uniform, that is, it should not depend on the particular polynomial vector field, only on its degree $n$. As of today, this question is not resolved even in the simplest case $n = 2$. Even finding non-trivial lower bounds for $\mathcal{H}(n)$ appears to be very hard.

<footnote n="1">¹It means that $\nabla H(x, y) \neq 0$ for all $(x, y)$ s.t. $H(x, y) = h$.</footnote>

<footnote n="2">²This makes $\Gamma(h)$ a Jordan curve, $C^\infty$-diffeomorphic to the unit circle.</footnote>

<footnote n="3">³The repositories for our code are available from https://gitlab.inria.fr/abintvalid.</footnote>

<footnote n="4">⁴https://julialang.org/</footnote>
In light of the lack of progress regarding bounds for $H(n)$, in the mid-seventies, V.I. Arnold [Arnold 1977, 1990] 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 $\mathbb{R}^2$:

$$ egin{cases} \dot{x} = -\frac{\partial H}{\partial y}(x, y) + \varepsilon P(x, y), \\ \dot{y} = \frac{\partial H}{\partial x}(x, y) + \varepsilon Q(x, y), \end{cases} \qquad (1.2) $$

where $H(x, y)$ is a polynomial of degree at most $n+1$. The polynomials $P$ and $Q$ of degree at most $n$, and the positive number $\varepsilon > 0$ define the small perturbation of the Hamiltonian system.

Is there a bound $Z(n)$ on the number of limit cycles the system (1.2) can have (for small $\varepsilon$), that only depends on the degree $n$?

For the infinitesimal problem, significant progress has been made, notably the proof of the finiteness of $Z(n)$ for all $n$ [Binyamini et al. 2010], and the uniform bound for the quadratic case $Z(2) = 2$ [Chen et al. 2006]

## 1.2 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 $\Sigma$, that is an open interval of a line crossing all orbits it encounters non-tangentially. This interval may be parameterized by the parameter $h$ (the energy level) of the unperturbed system for some domain $h^{-} < h < h^{+}$. By continuity, for every $h \in (h^{-}, h^{+})$ and sufficiently small $\varepsilon$, the trajectory originating from the point $\Sigma(h)$ on the transversal in the perturbed Hamiltonian system will cross $\Sigma$ again. The Poincaré return map $\Pi$ associates to $h$ and $\varepsilon$ the parameter $\Pi(h, \varepsilon)$ corresponding to the first point of return to $\Sigma$. We call $d(h, \varepsilon) = \Pi(h, \varepsilon) - h$ the displacement function. Clearly, the point $\Sigma(h)$ belongs to a periodic orbit of the perturbed system if and only if $d(h, \varepsilon) = 0$, and this is a limit cycle if and only if the zero of $d(\cdot, \varepsilon)$ is isolated.

The Poincaré-Pontryagin theorem roughly states that, for a Hamiltonian (or more generally, integrable) system perturbed by $\varepsilon$-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 $\varepsilon$ of the Poincaré return map. We now state the following adaptation of this result (see [Christopher and Li 2007] for a proof) that will prove useful in the sequel.

**THEOREM 1.3 (GENERALIZED POINCARÉ-PONTRYAGIN THEOREM).** Let $H: \mathbb{R}^2 \to \mathbb{R}$ be a real analytic potential function, $P, Q: \mathbb{R}^2 \to \mathbb{R}$ real analytic functions, $\varepsilon > 0$, and $\mu: \mathbb{R}^2 \to \mathbb{R}$ an analytic rescaling factor. Consider the perturbed integrable system:

$$ egin{cases} \dot{x} = -\mu(x, y) \frac{\partial H}{\partial y}(x, y) + \varepsilon P(x, y), \\ \dot{y} = \mu(x, y) \frac{\partial H}{\partial x}(x, y) + \varepsilon Q(x, y). \end{cases} \qquad (1.3) $$
Let $\Gamma(h)$ be an oval of H of level h over which $\mu$ does not vanish. Then the displacement function $d(h, \varepsilon)$ is approximated as

$$d(h, \varepsilon) = \varepsilon\mathfrak{I}(h) + O(\varepsilon^2), \quad \text{as } \varepsilon \to 0,$$

where the Abelian integral $\mathfrak{I}(h)$ was defined in (1.1).

In particular, $d(h, \varepsilon)$ and $\mathfrak{I}(h)$ have the same sign for small $\varepsilon > 0$. Hence, the number of isolated zeros of $\mathfrak{I}(h)$ where a change of sign occurs (in particular, simple zeros of $\mathfrak{I}(h)$) provides a lower bound for the number of limit cycles of (1.3) that exist for small $\varepsilon > 0$. Considering zeros of $\mathfrak{I}(h)$ 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) $\Gamma(h)$, see [Christopher and Li 2007].

A pessimistic, yet constructive, upper bound for $\mathcal{Z}(n)$ was obtained by the authors of [Binyamini et al. 2010] by bounding the number of zeros of $\mathfrak{I}(h)$ in terms of the degree $n$ only, using the Picard-Fuchs differential equations satisfied by $h \mapsto \mathfrak{I}(h)$. 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 $\mathcal{H}(n)$ [Li et al. 2009] or $\mathcal{Z}(n)$ [Johnson and Tucker 2010].

## 1.3 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., [Li et al. 2009] for $\mathcal{H}(3) \ge 13$), the numerical evaluation of $\mathfrak{I}(h)$ 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 [Moore 1966; Tucker 2011], 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 $\mathfrak{I}(h)$ 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 [Kapela et al. 2021]<fnref n="5" /> using the celebrated CAPD<fnref n="6" /> [Kapela et al. 2020] 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 [Johnson and Tucker 2011], by one of the authors of this article, computes the Abelian integral by subdividing the interior of the oval $\Gamma(h)$ and enclosing its border using parallelotopes. A similar approach can be found in [Martin et al. 2013] 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 [Beltrán and Leykin 2013; Xu et al. 2018]. 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 $\Gamma(h)$. 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

<footnote n="5">⁵See also the tutorial slides: https://kam.mff.cuni.cz/conferences/swim2015/slides/wilczak.pdf</footnote>

<footnote n="6">⁶http://capd.ii.uj.edu.pl/</footnote>
mostly restricted to curves of genus 0 [Abhyankar and Bajaj 1988; Hong and Schicho 1998], while approximate methods (piecewise rational, splines, etc.) developed for more general curves most of the time come without guaranteed error bounds (see e.g., [Bajaj and Xu 1997; Gao and Li 2004; Yang et al. 2010]). These methods usually target the standard binary 64 precision (i.e. double precision), which is not sufficient for our problem. A notable exception is the computation of certain elliptic integrals [Bost and Mestre 1988], with quadratic convergence due to the arithmetic-geometric mean.

## 1.4 Our approach and contributions

We propose a higher-order method where the oval $\Gamma(h)$ 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 I for $\mathfrak{I}(h)$ (cf. (1.1)) in $O(N \log N)$ arithmetic operations, where N is the degree used for trigonometric polynomials to approximate the curve $\Gamma(h)$. As $N \to \infty$, the diameter of I tends exponentially fast to 0.*

The theorem says that in order to compute $\mathfrak{I}(h)$ rigorously to $k$ digits of accuracy, one has to use a degree $N = O(k)$, where the constant hidden in this notation depends on the Hamiltonian system and $h$, but not on $k$.

Our approach is based on an *a posteriori* validation scheme, where we first approximate the oval $\Gamma(h)$ 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 $\Gamma(h)$ is approximated by the two routines OVALAPPROXINIT and OVALAPPROXREFINE (both described in Section 3). The result is a pair of degree $N$ trigonometric polynomials $(\hat{x}, \hat{y})$, such that $t \in [0, 2\pi] \mapsto (\hat{x}(t), \hat{y}(t))$ is a smooth approximate parameterization of $\Gamma(h)$.

(2) Next, Algorithm OVALVALID (described in Section 4) computes an *a posteriori* error bound for the distance between $(\hat{x}, \hat{y})$ and $\Gamma(h)$, 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 $\mathfrak{I}(h)$ 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.

# 2 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.

## 2.1 Interpolation of periodic functions

Throughout this paper, we shall need to approximate parameterizations of the ovals $\Gamma(h)$. Seeing that these parameterizations are continuous and $2\pi$-periodic, a natural choice for approximating them is via trigonometric polynomials [Zygmund 2002].
**Definition 2.1 (Trigonometric Polynomials).** A degree $N$ Trigonometric Polynomial (TP) $\tilde{f}$ is a real-valued trigonometric polynomial over $[0, 2\pi]$:

$$ \tilde{f}(t) = a_0 + \sum_{k=1}^{N} (a_k \cos(kt) + b_k \sin(kt)), $$

where $a_N$ or $b_N$ is nonzero. We represent $\tilde{f}$ in terms of its $2N + 1$ coefficients.

Let $\text{TP}_N(\mathbb{R})$ denote the vector space of trigonometric polynomials with real coefficients and degree at most $N$, and let $(C_{2\pi}^0, \|\cdot\|_{\infty})$ denote the vector space of real continuous $2\pi$-periodic functions, equipped with the supremum norm: $\|f\|_{\infty} = \max_{0 \le x \le 2\pi} |f(x)|$.

A usual way to obtain quasi-optimal trigonometric polynomial approximants is to consider the Fourier series expansion of the function $f \in C_{2\pi}^0$ truncated to order $N$. And yet, from a computational point of view, it is often more interesting to consider the degree (at most) $N$ trigonometric polynomial that interpolates the function under consideration at the equispaced points $(\frac{2k\pi}{2N+1})^{2N}$ [Wright et al. 2015]. We then define the interpolation operator

$$ \begin{gathered} \tilde{I}_N : C_{2\pi}^0 \to \text{TP}_N(\mathbb{R}) \\ f \mapsto \tilde{f} \quad \text{s.t.} \quad \tilde{f}\left(\frac{2k\pi}{2N+1}\right) = f\left(\frac{2k\pi}{2N+1}\right) \quad \text{for } k=0, \dots, 2N. \end{gathered} \tag{2.1} $$

The Lebesgue constant $\Lambda_N$ [Powell 1981] associated to the trigonometric interpolation operator $\tilde{I}_N$ is the operator norm of $\tilde{I}_N$, that is to say $\Lambda_N = \sup_{f \in C_{2\pi}^0, \|f\|_{\infty}=1} \|\tilde{I}_N(f)\|_{\infty}$. A first nice feature of these interpolation trigonometric polynomials is that $\Lambda_N = O(\log N)$ [Ehlich and Zeller 1966], which makes them quasi-optimal approximations with respect to $\|\cdot\|_{\infty}$.

A second nice feature is that we can compute the coefficients of the interpolation polynomial by the Fast Fourier Transform, hence in $O(N \log(N))$ arithmetic operations. Indeed, if $\tilde{f}(t) = a_0 + \sum_{k=1}^N (a_k \cos(kt) + b_k \sin(kt))$ and $f_j = \tilde{f}\left(\frac{2jk\pi}{2N+1}\right)$ for $j = 0, \dots, 2N$, we have, for $k = 1, \dots, N$,

$$ a_0 = \frac{1}{2N+1} \sum_{j=0}^{2N} f_j, \quad a_k = \frac{1}{2N+1} \sum_{j=0}^{2N} f_j \cos\left(\frac{2jk\pi}{2N+1}\right), \quad b_k = \frac{1}{2N+1} \sum_{j=0}^{2N} f_j \sin\left(\frac{2jk\pi}{2N+1}\right). $$

Having this in mind, we define

$$ \begin{gathered} \text{FFT}_N : \mathbb{R}^{2N+1} \to \text{TP}_N(\mathbb{R}) \\ (f_j)_{j=0}^{2N} \mapsto \tilde{f} \quad \text{s.t.} \quad \tilde{f}(t_j) = f_j \quad \text{for } 0 \le j \le 2N, \end{gathered} $$

and its inverse transform

$$ \begin{gathered} \text{IFFT}_N : \text{TP}_N(\mathbb{R}) \to \mathbb{R}^{2N+1} \\ \tilde{f} \mapsto (\tilde{f}(t_j))_{j=0}^{2N} \quad \text{where} \quad t_j = \frac{2\pi j}{2N+1} \quad \text{for } 0 \le j \le 2N, \end{gathered} $$

which allows us to interpolate a function and to evaluate a TP on the equispaced grid in only $O(N \log(N))$ arithmetic operations. By combining both routines, we can compute the degree $2N$ product of two degree $N$ TPs with the same asymptotic complexity $O(N \log(N))$. In practice, this outperforms the “naive” multiplication algorithm in $O(N^2)$ arithmetic operations for degrees $N$ of the order of magnitude of 100 or more. We assume this $O(N \log(N))$ asymptotic complexity for the product of degree $N$ TPs in what follows.
When working with real trigonometric approximations, a convenient way to bound the supremum norm of functions is to use the $|| \cdot ||_{\ell^1}$ norm directly defined from the Fourier coefficients $a_k, b_k$ of $f$:

$$ \|f\|_{\ell^1} = |a_0| + \sum_{k=1}^{\infty} (|a_k| + |b_k|). $$

This norm is well-defined as long as the sum converges, and then satisfies $||f||_{\ell^1} \ge ||f||_{\infty}$. Note that the sum does not necessarily converge for $f \in C_{2\pi}^0$. However, this becomes true under very mild regularity assumptions, and a fortiori<fnref n="7" /> for analytic $f$ as assumed throughout this article. Finally, we have that $||fg||_{\ell^1} \le ||f||_{\ell^1}||g||_{\ell^1}$.

## 2.2 Rigorous trigonometric approximations

By analogy with intervals used to represent real numbers rigorously on computers [Moore 1966], (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-arithmetic* 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 [Brisebarre and Joldeş 2010] 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., [Figueras and de la Llave 2017; Hungria et al. 2016]). 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 $N$ *Rigorous Trigonometric Approximation* (RTA) is a pair $f = (\check{f}, \varepsilon)$, with $\check{f}$ a degree $N$ TP and $\varepsilon \ge 0$ representing the closed $\varepsilon$-ball around $\check{f}$ in the Banach space $C_{2\pi}^0$:

$$ \bar{B}(\check{f}, \varepsilon) = \{ f \in C_{2\pi}^0 \mid \|f - \check{f}\|_{\infty} \le \varepsilon \}. $$

Therefore, we say that a degree $N$ RTA $f = (\check{f}, \varepsilon)$ represents a function $f \in C_{2\pi}^0$ (or $f$ contains $f$, or simply $f \in f$) if there exists a function $s \in C_{2\pi}^0$ with $\|s\|_{\infty} \le \varepsilon$ such that $f(t) = \check{f}(t) + s(t)$ for all $t \in \mathbb{R}$.

Linear operations are trivially defined on RTAs to match the (Banach) linear space structure of $(C_{2\pi}^0, || \cdot ||_{\infty})$: for all degree $N$ RTAs $(\check{f}, \varepsilon)$, $(\check{g}, \eta)$ and $\lambda \in \mathbb{R}$, we have

$$ (\check{f}, \varepsilon) + (\check{g}, \eta) = (\check{f} + \check{g}, \varepsilon + \eta), $$

$$ (\check{f}, \varepsilon) - (\check{g}, \eta) = (\check{f} - \check{g}, \varepsilon + \eta), $$

$$ \lambda(\check{f}, \varepsilon) = (\lambda\check{f}, |\lambda|\varepsilon). $$

Also, $\text{BOUND}(f)$, defined as:

$$ \text{BOUND}((\check{f}, \varepsilon)) = \sum_{k=0}^{N} |a_k| + \sum_{k=1}^{N} |b_k| + \varepsilon. $$

<footnote n="7">⁷This is a straightforward corollary of the exponential decay of Fourier coefficients for analytic functions, see e.g. [Wright et al. 2015, Thm. 4.1].</footnote>
satisfies BOUND(f) ≥ ||f||ℓ¹ + ε ≥ ||f||∞ + ε, so that ||f||∞ ≤ BOUND(f) for all f ∈ f.

The product of two degree N RTAs is the degree 2N RTA defined by

$$ (\tilde{f}, \varepsilon) \times (\check{g}, \eta) = (\tilde{f}\check{g}, \text{BOUND}(\tilde{f})\eta + \text{BOUND}(\check{g})\varepsilon + \varepsilon\eta), $$

which can be computed in O(N log(N)) arithmetic operations using the fast multiplication algorithm for TPs mentioned previously.

We shall also define the integral of f = (f̃, ε) with f̃ = a₀ + ΣNk=1 (aₖ cos(kt) + bₖ sin(kt)) over one period:

$$ \int_{0}^{2\pi} f(t)dt = 2\pi (a_{0} + [-\varepsilon, \varepsilon]). $$

LEMMA 2.3. If the RTAs f = (f̃, ε) and g = (g̃, η) represent functions f and g, then the RTA λf represents λf for any λ ∈ ℝ, and the RTA f★g represents f★g for ★ ∈ {+, -, ×}. Also, ∫2π₀ f(t)dt ∈ ∫2π₀ f(t)dt for all f ∈ f.

PROOF. We prove the least trivial case only, namely the multiplication. By hypothesis, there exist functions u and v in C02π with ||u||∞ ≤ ε and ||v||∞ ≤ η, such that f = f̃ + u and g = g̃ + v. Therefore,

$$ \begin{aligned} ||fg - \tilde{f}\tilde{g}||_{\infty} &\le ||\tilde{f}v + \check{g}u + uv||_{\infty} \le ||\tilde{f}||_{\infty}||v||_{\infty} + ||\check{g}||_{\infty}||u||_{\infty} + ||u||_{\infty}||v||_{\infty} \\ &\le \text{BOUND}(\tilde{f})\eta + \text{BOUND}(\check{g})\varepsilon + \varepsilon\eta, \end{aligned} $$

which by definition implies fg ∈ f × g.

□

REMARK 2.4. By analogy with a globally set precision for floating-point and interval arithmetic, we can choose to fix a global degree N 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:

$$ T_{\text{RUNC}}(f, N) = \left( a_0 + \sum_{k=1}^{N} (a_k \cos(kt) + b_k \sin(kt)), \sum_{k=N+1}^{m} (|a_k| + |b_k|) + \varepsilon \right), $$

for f = (a₀ + Σm k=1 (aₖ cos(kt) + bₖ sin(kt)), ε) when m > N. Clearly, if f ∈ f, then f ∈ T_RUNC(f, N).

These elementary operations on RTAs are the building blocks of more complex ones. For instance, if p(x) is a polynomial and f an RTA for f, then substituting x by f and applying an evaluation scheme of choice with arithmetic operations on RTAs gives an RTA p(f) for p(f).

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

## 2.3 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 Γ(h). 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 $\varepsilon \ge \|\dot{\varphi} - \varphi^*\|_{\infty}$ is rigorously computed by expressing $\varphi^*$ as a fixed point of a contracting operator $\mathcal{T}$, and applying the Banach fixed-point theorem (see Theorem 2.5 below). The result is an RTA $\varphi = (\dot{\varphi}, \varepsilon)$ such that $\varphi^* \in \varphi$.

**THEOREM 2.5.** (see, e.g., [Yamamoto 1998, Thm. 2.1]) Let $(E, \|\cdot\|)$ be a Banach space, $\dot{\varphi} \in E$, and $\mathcal{T}: E \to E$. If we can find $r, \lambda \ge 0$ such that:

* $\mathcal{T}$ is $\lambda$-contracting over the closed ball $\bar{B}(\dot{\varphi}, r)$, i.e., $\lambda < 1$ and:

$$ \|\mathcal{T}(\varphi_1) - \mathcal{T}(\varphi_2)\| \le \lambda \|\varphi_1 - \varphi_2\| \quad \text{for all } \varphi_1, \varphi_2 \in \bar{B}(\dot{\varphi}, r); $$

* We have $d + \lambda r \le r$, where $d := \|\mathcal{T}(\dot{\varphi}) - \dot{\varphi}\|$ (called the defect);

Then $\mathcal{T}$ admits a unique fixed point $\varphi^*$ in the ball $\bar{B}(\dot{\varphi}, r)$, and we have the following enclosure for the approximation error:

$$ \frac{d}{1 + \lambda} \le \|\dot{\varphi} - \varphi^*\| \le \frac{d}{1 - \lambda}. $$

Assume the functional equation is of the form $\mathcal{F}(\varphi) = 0$ with $\mathcal{F}: C_{2\pi}^0 \to C_{2\pi}^0$ being of class $C^1$, and let $\dot{\varphi}$ be an approximate zero of $\mathcal{F}$. A classical method in rigorous numerics (see e.g., [Hungria et al. 2016; Yamamoto 1998]) to obtain an equivalent fixed-point equation $\mathcal{T}(\varphi) = \varphi$ is to construct $\mathcal{T}: C_{2\pi}^0 \to C_{2\pi}^0$ as a Newton-like operator:

$$ \mathcal{T}(\varphi) = \varphi - \mathcal{A}(\mathcal{F}(\varphi)). $$

Here, $\mathcal{A}: C_{2\pi}^0 \to C_{2\pi}^0$ is an invertible bounded linear operator approximating $(\mathcal{D}\mathcal{F}(\dot{\varphi}))^{-1}$, the inverse of the differential of $\mathcal{F}$ at $\dot{\varphi}$, making $\mathcal{T}$ a local contraction in a neighborhood of $\dot{\varphi}$.

**2.3.1 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

$$ c_r(t)\varphi(t)^r + c_{r-1}(t)\varphi(t)^{r-1} + \dots + c_1(t)\varphi(t) + c_0(t) = 0, \quad t \in [0, 2\pi], $$

with $c_0, \dots, c_r \in C_{2\pi}^0$, that is to say, $\mathcal{F} \in C_{2\pi}^0[X]$ is a polynomial with coefficients in $C_{2\pi}^0$.

In this case, the differential $\mathcal{D}\mathcal{F}(\dot{\varphi})$ coincides with the multiplication by $\mathcal{F}'(\dot{\varphi}) \in C_{2\pi}^0$, where the prime symbol denotes the usual differentiation in a ring of polynomials, here $C_{2\pi}^0[X]$. Therefore, the linear operator $\mathcal{A}$ is set to be the multiplication by a TP $\acute{a}$ constructed by interpolation to approximate $1/\mathcal{F}'(\dot{\varphi})$:

$$ \mathcal{A}(\varphi) = \acute{a}\varphi, \quad \text{with} \quad \acute{a}(t) \approx \frac{1}{\mathcal{F}'(\acute{a}(t))}. $$

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

**LEMMA 2.6.** If $\mathcal{F} \in C_{2\pi}^0[X]$, $\mathcal{A}(\varphi)(t) = \acute{a}(t)\varphi(t)$, and $\mathcal{T}(\varphi) = \varphi - \mathcal{A}(\mathcal{F}(\varphi))$, then $\mathcal{T}$ is $\lambda(r)$-Lipschitz over $\bar{B}(\dot{\varphi}, r)$ with $\lambda(r) = \lambda_0 + \alpha r \lambda_1(r)$, obtained from the (rigorously computed) bounds:

$$ \lambda_0 \ge \|1 - \acute{a} \mathcal{F}'(\dot{\varphi})\|_{\infty}, \quad \alpha \ge \|\acute{a}\|_{\infty}, \quad \lambda_1(r) \ge \sup_{\varphi \in \bar{B}(\dot{\varphi}, r)} \|\mathcal{F}''(\varphi)\|_{\infty}. $$

**PROOF.** To get an upper bound for the Lipschitz constant for $\mathcal{T} \in C_{2\pi}^0[X]$ over $\bar{B}(\dot{\varphi}, r)$, we simply bound $\|\mathcal{T}'(\varphi)\|_{\infty}$ for $\|\varphi - \dot{\varphi}\|_{\infty} \le r$ by the triangle inequality. First, $\|\mathcal{T}'(\dot{\varphi})\|_{\infty} = \|1 - \acute{a} \mathcal{F}'(\dot{\varphi})\|_{\infty} \le \lambda_0$ by definition of $\lambda_0$. Then, $\|\mathcal{T}'(\varphi) - \mathcal{T}'(\dot{\varphi})\|_{\infty} = \|\acute{a} (\mathcal{F}'(\varphi) - \mathcal{F}'(\dot{\varphi}))\|_{\infty} \le \alpha r \lambda_1(r)$ using the mean value theorem, by definition of $\alpha$ and $\lambda_1$, and since $\varphi \in \bar{B}(\dot{\varphi}, r)$. Adding both bounds, we obtain the expected upper bound for the Lipschitz constant $\lambda(r)$. $\square$
Note that $\mathcal{A}$ being invertible is a byproduct of the second assumption of Theorem 2.5: except in the trivial case $d = 0$, the inequality implies $\lambda_0 < 1$, so $\check{a}$ never vanishes.

**EXAMPLE 2.7 (DIVISION OF RTAs).** For two functions $g, h \in C_{2\pi}^0$ with $h(t) \neq 0$ for all $t$, represented by RTAs $g$ and $h$, Algorithm RTAD<sub>IV</sub>($g, h, N$) computes a degree $N$ TP $\check{\varphi}$ approximating the quotient $\varphi^* = g/h$ by trigonometric interpolation, and validates it using Lemma 2.6, since $\varphi^*$ is the unique zero of $\mathcal{F} \in C_{2\pi}^0[X]$ defined by $\mathcal{F}(\varphi) = h\varphi - g$. When a global degree $N$ for the RTA computations is fixed, we may use the notation $g/h$ for RTAD<sub>IV</sub>($g, h, N$).

<caption>**Algorithm 1 RTAD<sub>IV</sub>(g, h, N)**</caption>

**Input:** RTAs $g = (\check{g}, \delta)$ and $h = (\check{h}, \eta)$, and approximation degree $N \in \mathbb{N}$
**Output:** Degree $N$ RTA $\varphi$ representing the quotient $g/h$

▷ Build degree $N$ candidate TP $\check{\varphi}$ by interpolation

1: $(g_j)_{j=0}^{2N} \leftarrow \text{IFFT}_N(\check{g})$ and $(h_j)_{j=0}^{2N} \leftarrow \text{IFFT}_N(\check{h})$

2: $\check{\varphi} \leftarrow \text{FFT}_N\left(\left(\frac{g_j}{h_j}\right)_{j=0}^{2N}\right)$

▷ Build TP $\check{a}$ defining the Newton-like operator and bound the Lipschitz constant

3: $\check{a} \leftarrow \text{FFT}_N\left(\left(\frac{1}{h_j}\right)_{j=0}^{2N}\right)$

4: $\lambda \leftarrow \text{BOUND}(1 - \check{a}h)$

▷ Compute a posteriori error bound, if possible

5: **if** $\lambda < 1$ **then**

6:     $d \leftarrow \text{BOUND}(\check{a}(h\check{\varphi} - g))$

7:     $\varepsilon \leftarrow \frac{d}{1-\varepsilon}$

8:     **return** $\varphi = (\check{\varphi}, \varepsilon)$

9:     **else**

10: **return** FAIL

11: **end if**

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

**LEMMA 2.8.** Let $g = (\check{g}, \delta)$ and $h = (\check{h}, \eta)$ be RTAs. Then,

(i) For all $N \in \mathbb{N}$, if RTAD<sub>IV</sub>(g, h, N) returns an RTA $\varphi$, then for all $g \in g$ and $h \in h$, $h(t) \neq 0$ for all $t$ and $g/h \in \varphi$.

(ii) If $\check{h}(t) \neq 0$ for all $t$ and $\eta < \|\check{h}^{-1}\|_{\ell^1}^{-1}$, then there exists an $N_0 \in \mathbb{N}$ such that for all $N \ge N_0$, Algorithm RTAD<sub>IV</sub>(g, h, N) does not fail, and the remainder of $\varphi = (\check{\varphi}, \varepsilon)$ satisfies

$$ \varepsilon \le \frac{\|\check{h}^{-1}\|_{\ell^1} \delta + \|\check{g}/\check{h}^2\|_{\ell^1} \eta}{1 - \|\check{h}^{-1}\|_{\ell^1} \eta} + O(\kappa^{-N}), \quad \text{as } N \to \infty, $$

for some $\kappa > 1$ depending on $\check{g}$ and $\check{h}$.

(iii) RTAD<sub>IV</sub>(g, h, N) runs in $O(N' \log(N'))$ arithmetic operations, where $N' = \max(\deg g, \deg h, N)$.

**PROOF.** (i) Soundness. The Newton-like operator associated to $\mathcal{F}$ is defined by $\mathcal{T}(\varphi) = \varphi - \check{a}(h\varphi - g)$ with $\check{a}$ a degree $N$ TP approximating $1/\mathcal{F}'(\check{\varphi}) = 1/h(t)$ (line 3). Since $\mathcal{F}$ (and hence $\mathcal{T}$) has degree 1,
$\mathcal{F}'' = 0$. Hence the bound $\lambda_1(r)$ in Lemma 2.6 is identically 0, and we have $\lambda(r) = \lambda_0 = \|1 - \hat{a}h\|_{\infty}$, rigorously upper-bounded by $\lambda$ in line 4. Hence, as soon as $\lambda < 1$, we can set $r$ arbitrary large in Theorem 2.5 and return the bound implemented in lines 6 and 7:

$$ \|\check{\varphi} - \varphi^*\|_{\infty} \le \frac{\|\mathcal{T}(\check{\varphi}) - \varphi^*\|_{\infty}}{1 - \lambda} \le \frac{\|\check{a}(h\check{\varphi} - g)\|_{\infty}}{1 - \lambda} \le \frac{\text{BOUND}(\check{a}(h\check{\varphi} - g))}{1 - \lambda}. $$

Therefore, the degree $N$ RTA $\varphi = (\check{\varphi}, \varepsilon)$ returned in line 8 contains $g/h$ for all $g \in g$ and $h \in h$. Note moreover that the nonvanishing of all $h \in h$ is a byproduct of the inequality $\|1 - \hat{a}h\|_{\infty} \le \lambda < 1$.

(ii) *Convergence*. Suppose that $\check{h}(t) \ne 0$ for all $t$, and that $\eta < \|\check{h}^{-1}\|_{\ell^1}^{-1}$. In particular, $\eta < \|\check{h}^{-1}\|_{\infty}^{-1} = \min_{t \in \mathbb{R}} |\check{h}(t)|$. Since $\check{h}^{-1}$ (resp. $\check{g}/\hat{h}$) is analytic over $\mathbb{R}$, then $\check{a} = \mathcal{I}_N(\check{h}^{-1})$ (resp. $\check{\varphi} = \mathcal{I}_N(\check{g}/\hat{h})$) converges to $\check{h}^{-1}$ (resp. $\check{g}/\hat{h}$) w.r.t. $\|\cdot\|_{\ell^1}$ exponentially fast as $N \to \infty$ (see e.g. [Wright et al. 2015, Thm. 4.2]$^8$).

Now for the bound $\lambda$ computed in line 4, we have for all $h \in h$,

$$ \|1 - \hat{a}h\|_{\ell^1} \le \|1 - \check{a}\|_{\ell^1} \|h - \check{h}\|_{\ell^1} \le \|\check{a} - \check{h}^{-1}\|_{\ell^1} \|\check{h}\|_{\ell^1} + \|\check{a}\|_{\ell^1} \eta \le \|\check{h}^{-1}\|_{\ell^1} \eta + O(\kappa^{-N}), $$

for some $\kappa > 1$, as $N \to \infty$. In particular, Algorithm RTADiv always terminates for sufficiently large $N$. Similarly, for the bound $d$ computed in line 6, for all $g \in g$ and $h \in h$,

$$ \|\check{a}(h\check{\varphi} - g)\|_{\ell^1} \le \underbrace{\|\check{a}(\check{h}\check{\varphi} - \check{g})\|_{\ell^1}}_{\to 0} + \underbrace{\|\check{a}\|_{\ell^1}}_{\to \|\check{h}^{-1}\|_{\ell^1}} \delta + \underbrace{\|\check{a}\check{\varphi}\|_{\ell^1}}_{\to \|\check{g}/\check{h}^2\|_{\ell^1}} \eta \le \|\check{h}^{-1}\|_{\ell^1} \delta + \|\check{g}/\check{h}^2\|_{\ell^1} \eta + O(\kappa^{-N}), $$

for some $\kappa > 1$, as $N \to \infty$. This finally yields the expected estimate for $\varepsilon$ computed in line 7.

(iii) *Complexity*. The asymptotic complexity is determined by the FFT$_N/IFFT$_N routines used for evaluation on the equispaced grid, interpolation and TP/RTA multiplications. $\square$

**REMARK 2.9.** In the case where the polynomial operator $\mathcal{F}$ has degree greater than 1 (which will be the case for Algorithm OVALVALID in Section 4 as long as the potential function $H(x, y)$ has degree at least 2), an explicit radius $r$ satisfying the hypothesis of Theorem 2.5 must be rigorously computed. Given the function $r \mapsto \lambda(r)$, this amounts to locating the smallest positive zero of $(1 - \lambda(r))r - d$. A simple bisection method NEWTONBALL is proposed in Appendix C for this purpose.

### 3 APPROXIMATION OF THE OVAL $\Gamma(h)$ AND NUMERICAL INTEGRATION

This section is devoted to the computation of a smooth approximate parameterization $t \in [0, 2\pi] \to (\check{x}(t), \check{y}(t))$ of the oval $\Gamma(h)$ with $\check{x}, \check{y}$ TPs of degree $N$.

First, an implicit parameterization for $\Gamma(h)$ is given by the unscaled, unperturbed Hamiltonian system itself (i.e., System (1.2) with $\varepsilon = 0$):

$$ \left\{ \begin{aligned} \dot{x} &= -\frac{\partial H}{\partial y}(x, y), \\ \dot{y} &= \frac{\partial H}{\partial x}(x, y). \end{aligned} \right. \qquad (3.1) $$

For a level $h$ on the portion of $\Sigma$ of interest, the trajectory of (3.1) starting at $(x_{ini}, y_{ini}) = \Sigma \cap H^{-1}(h)$ is periodic since it exactly follows the oval $\Gamma(h)$. Call $t_{end}$ the time of first return onto $\Sigma$. By a linear change of the independent variable $t$, we may assume that $t_{end} = 2\pi$, giving a periodic analytic parameterization

$$ t \in [0, 2\pi] \mapsto (x^*(t), y^*(t)) \text{ of } \Gamma(h), \qquad (3.2) $$

with a winding number equal to $\pm 1$ with respect to any point inside $\Gamma(h)$.

<footnote>$^8$This theorem provides convergence bounds w.r.t. the $\|\cdot\|_{\infty}$ norm, but the same bounds actually hold for the $\|\cdot\|_{\ell^1}$ norm since the proof consists in bounding the Fourier coefficients.</footnote>
The approximation of $\Gamma(h)$ consists of two steps:

* First, Algorithm OVALAPPROXINIT presented in Section 3.1 computes an initial approximation $(x^\circ, y^\circ) \in (C_{2\pi}^0)^2$ of $(x^*, y^*)$, with TPs $x^\circ, y^\circ$ of a low degree $N_0$. The accuracy of this initial guess (and therefore $N_0$) should ensure the convergence of the Newton iterations in the next routine OVALAPPROXREFINE. It only depends on the system under consideration (i.e., $H, \Sigma$ and $h$) 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 $\tilde{x}, \tilde{y}$ of higher degree $N$. As explained in Section 3.2, this is done by applying Newton iterations to the initial guess so as to "project" it onto $\Gamma(h)$. We provide a complexity analysis and an asymptotic convergence estimate of $(\tilde{x}, \tilde{y})$ to $\Gamma(h)$ in terms of the degree $N$.

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.

## 3.1 Algorithm OVALAPPROXINIT: initial approximation of the curve

Algorithm OVALAPPROXINIT is given an initial point $(x_{ini}, y_{ini}) \in \Gamma(h) \cap \Sigma$ on the transversal. It first calls a numerical ODE solver to approximately solve the differential system (3.1) until the first return onto $\Sigma$. It provides us with a continuous function $t \mapsto (\tilde{x}(t), \tilde{y}(t))$, which we assume to be rescaled over the time interval $[0, 2\pi]$. Finally, in order to obtain smooth approximations to $(x^*, y^*)$, cf. (3.2), we interpolate $(\tilde{x}, \tilde{y})$ at evenly spaced points and obtain low degree TPs $x^\circ, y^\circ$ (lines 2–5) representing our analytic initial guess for $\Gamma(h)$. The proof of the following Proposition is postponed to the end of this section.

**PROPOSITION 3.1.** *Under Assumption 3.2, Algorithm OVALAPPROXINIT($H, \Sigma, x_{ini}, y_{ini}, N_0$) computes degree $N_0$ TPs $x^\circ$ and $y^\circ$ that converge uniformly to the rescaled parameterization $(x^*, y^*) \in (C_{2\pi}^0)^2$ of $\Gamma(h)$:*

$$ \|x^\circ - x^*\|_{\infty}, \|y^\circ - y^*\|_{\infty} \to 0 \quad \text{as} \quad N_0 \to \infty. $$

### Algorithm 2 OVALAPPROXINIT($H, \Sigma, x_{ini}, y_{ini}, N_0$)

**Input:** Potential function $H$, transversal $\Sigma$, initial point $(x_{ini}, y_{ini})$, and approximation degree $N_0 \in \mathbb{N}$
**Output:** Degree $N_0$ TPs $x^\circ, y^\circ$ approximating $\Gamma(h)$

* **Approximation of $\Gamma(h)$, computed with precision parameter $N_0^{-1}$**

1: $(\tilde{x}, \tilde{y}) \leftarrow DSOLVE(H, \Sigma, x_{ini}, y_{ini}, N_0^{-1})$

2: **for** $j = 0$ **to** $2N_0$ **do**

3: &nbsp;&nbsp;&nbsp;&nbsp; $x_j^\circ \leftarrow \tilde{x} \left( \frac{2j\pi}{2N_0 + 1} \right)$, $y_j^\circ \leftarrow \tilde{y} \left( \frac{2j\pi}{2N_0 + 1} \right)$

4: **end for**

* **Interpolate the points**

5: $x^\circ \leftarrow FFT_{N_0}((x_j^\circ)^{2N_0})$, $y^\circ \leftarrow FFT_{N_0}((y_j^\circ)^{2N_0})$

6: **return** $(x^\circ, y^\circ)$
**3.1.1 The numerical solver.** We need a numerical routine `DSOLVE(H, Σ, x<sub>ini</sub>, y<sub>ini</sub>, δ)` that integrates the ODE system (3.1) with some precision parameter `δ`, starting from the initial point `(x<sub>ini</sub>, y<sub>ini</sub>) ∈ Γ(h) ∩ Σ`. This routine has to detect numerically the first return onto `Σ` and stop at the corresponding time `t<sub>end</sub>`. Moreover, we assume that this routine rescales the time interval `[0, t<sub>end</sub>]` 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 $\delta$ 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
$t \mapsto (\tilde{x}(t), \tilde{y}(t))$.

**3.1.2 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 $\Sigma$ is trivially recast as a linear equation in the Cartesian coordinates. Note, however, that the event is already satisfied at $t=0$ (up to numerical errors), since $(x_{ini}, y_{ini}) \in \Sigma$, 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 $(x_r, y_r)$ of choice inside $\Gamma(h)$, that is the ODE system:

$$
\begin{equation}
\begin{cases}
x(t) = x_r + \rho(t) \cos \theta(t), \\
y(t) = y_r + \rho(t) \sin \theta(t),
\end{cases}
\quad \text{with} \quad
\begin{cases}
\dot{\rho} = \sin \theta \frac{\partial H}{\partial x}(x, y) - \cos \theta \frac{\partial H}{\partial y}(x, y), \\
\dot{\theta} = \frac{\cos \theta \frac{\partial H}{\partial x}(x, y) + \sin \theta \frac{\partial H}{\partial y}(x, y)}{\rho}.
\end{cases}
\end{equation}
$$

The stopping condition is when $\theta(t)$ reaches $\theta(0) + \zeta 2\pi$, with $\zeta = \pm 1$ the winding number of $(x^*, y^*)$ with respect to $(x_r, y_r)$.

**3.1.3 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 $(\tilde{x}, \tilde{y})$, computed with a precision parameter $\delta = 1/N_0$, converges fast enough to $(x^*, y^*)$ to balance the overestimation factor due to the degree $N_0$ trigonometric interpolation. In practice, such a mild convergence assumption in $O(\delta)$ is guaranteed by most iterative schemes, even a simple Euler scheme with time step $\delta$.

ASSUMPTION 3.2. The numerical routine DSOLVE(H, Σ, x<sub>ini</sub>, y<sub>ini</sub>, δ) computes a continuous approximation t ∈ [0, 2π] ↦→ (x̃(t), ỹ(t)) of the rescaled parameterization (x<sup>*</sup>, y<sup>*</sup>), cf. (3.2), of Γ(h), such that:

$$
\|\tilde{x} - x^*\|_{\infty}, \|\tilde{y} - y^*\|_{\infty} = O(\delta) \quad \text{as} \quad \delta \to 0.
$$

We can now prove the convergence of Algorithm OVALAPPROXINIT.

*Proof of Proposition 3.1.* Let the symbol *z* denote either *x* or *y*. Let η = *z*̃ − *z*\* be the error between *z*\* and the approximate solution *z̃* computed by DSOLVE in line 1. Thanks to the linearity of the operator *I*<sub>*N*<sub>0</sub></sub> (see (2.1)), the error between the interpolant *z*° computed in lines 2–5 and *z*\* can be bounded from above by the error of interpolation for *z*\*, plus the interpolation of the error η :

$$
\|z^{\circ} - z^{*}\|_{\infty} = \|\mathcal{I}_{N_0}(\tilde{z}) - z^{*}\|_{\infty} \le \|\mathcal{I}_{N_0}(z^{*}) - z^{*}\|_{\infty} + \|\mathcal{I}_{N_0}(\eta)\|_{\infty}.
$$

• *Error of interpolation.* Since *z*\* is 2π-periodic and continuously differentiable, the trigonometric interpolants *I*<sub>*N*<sub>0</sub></sub>(*z*\*) converge uniformly to *z*\* [Zygmund 2002, Chap. X].
• Interpolation of error. By Assumption 3.2, $η$ is a continuous function over $[0, 2π]$ and $||η||_∞$ converges uniformly to 0 in $O(N_0^{-1})$ as the precision parameter $δ = N_0^{-1}$ tends to 0. From Section 2.1, we have

$$||\mathcal{I}_{N_0}(η)||_∞ \leq Λ_{N_0}||η||_∞ = O(\log(N_0)/N_0) \to 0 \text{ as } N_0 \to \infty. \quad \square$$

## 3.2 Algorithm OVALAPPROXREFINE: refining the curve by projection

As its name suggests, the role of Algorithm OVALAPPROXREFINE is to refine the initial guess $(x^°, y^°)$ with TPs $(\tilde{x}, \tilde{y})$ of higher degree $N \geq N_0$, yielding an analytic approximate parameterization of $Γ(h)$ of very high accuracy.

By solving algebraic equations only, it is in general not possible to recover the parameterization $(x^*, y^*) \in (C_{2π}^0)^2$ satisfying the differential system (3.1). However, the value of $\mathfrak{H}(h)$ is independent of the parameterization used for $Γ(h)$. Therefore, we define another one, $(x^#, y^#) \in (C_{2π}^0)^2$, called the projected parameterization, which Algorithm OVALAPPROXREFINE approximates using Newton's method on points sampled on the initial curve $(x^°, y^°)$.

**Definition 3.3.** Let $(x^°, y^°) \in (C_{2π}^0)^2$ be an approximate parameterization of $Γ(h)$, $(u^°, v^°) \in (C_{2π}^0)^2$, and $ε > 0$. If the equation

$$H(x^°(t) + s(t) u^°(t), y^°(t) + s(t) v^°(t)) = h, \quad t \in [0, 2\pi], \quad (3.3)$$

has a unique solution $s^# \in C_{2π}^0$ with $||s^#||_∞ \leq ε$, called the shift, then the curve $(x^#, y^#) \in (C_{2π}^0)^2$,

$$\begin{cases} x^#(t) = x^°(t) + s^#(t) u^°(t), \\ y^#(t) = y^°(t) + s^#(t) v^°(t), \end{cases} \quad (3.4)$$

is called the projection of $(x^°, y^°)$ onto $Γ(h)$ w.r.t. the direction $(u^°, v^°)$.

The following Lemma, whose proof is postponed to Appendix B, establishes the good properties of the projected parameterization when the initial guess $(x^°, y^°)$ returned by OVALAPPROXINIT is sufficiently close to $Γ(h)$.

**LEMMA 3.4.** There exists $η > 0$ (depending on $H$ and $h$) such that if $(x^°, y^°)$ and $(u^°, v^°)$ are analytic and $η$-close<fnref n="9" /> to $(x^*, y^*)$ and $\frac{\nabla H(x^*, y^*)}{||\nabla H(x^*, y^*)||_2}$, respectively, in $(C_{2π}^0)^2$, then:

(i) $(x^#, y^#) \in (C_{2π}^0)^2$ is well-defined and it is an analytic parameterization of $Γ(h)$ with same orientation as $(x^*, y^*)$.

(ii) For all $t \in [0, 2\pi]$, letting $(x, y) = (x^°(t), y^°(t))$ and $(u, v) = (u^°(t), v^°(t))$, the Newton operator associated to Equation (3.3) at $t$ is well-defined:

$$N_{x,y,u,v}(s) = s - \frac{H(x + su, y + sv) - h}{u \frac{\partial H}{\partial x}(x + su, y + sv) + v \frac{\partial H}{\partial y}(x + su, y + sv)}, \quad (3.5)$$

and the Newton iterations $s^{(k)} = N_{x,y,u,v}^k(0)$ converge quadratically fast to $s^#(t)$.

Based on those properties, Algorithm OVALAPPROXREFINE discretizes the initial guess $(x^°, y^°)$ over the equispaced grid $(t_j)$ of size $2N + 1$ (line 1), then applies several Newton iterations (3.5), implemented in lines 6–8, and finally interpolates the obtained values $s_j$ in order to construct a degree $N$ TP $š$ approximating the shift $s^#$. This defines the approximate parameterization $(\tilde{x}, \tilde{y})$ of $Γ(h)$:

$$\left\{ \begin{array}{l} \tilde{x}(t) = x^°(t) + š(t) u^°(t), \\ \tilde{y}(t) = y^°(t) + š(t) v^°(t), \end{array} \right. \quad t \in [0, 2\pi]. \quad (3.6)$$

<footnote n="9">⁹It means that $\max(||x^° - x^*||_∞, ||y^° - y^*||_∞) \leq η$.</footnote>
**Algorithm 3** OVALAPPROXREFINE(H, h, x°, y°, N)

**Input:** Potential function H, level h, initial approximation (x°, y°) and degree N ∈ $\mathbb{N}$

**Output:** Degree N TPs u°, v°, š defining the approximate parameterization (x̃, ỹ) of Γ(h) as in Eq. (3.6)

► Sample initial points

1: $(x_j)_{j=0}^{2N} \leftarrow \text{IFFT}_N(x^\circ), (y_j)_{j=0}^{2N} \leftarrow \text{IFFT}_N(y^\circ)$

► Apply $[\log N]$ Newton iterations

2: **for** j = 0 to 2N **do**

3:   $(u_j, v_j) \leftarrow \frac{\nabla H(x_j, y_j)}{\|\nabla H(x_j, y_j)\|_2}$

4:   $s_j \leftarrow 0$

5:   **for** k = 1 to $[\log N]$ **do**

6:   $\delta \leftarrow \frac{h - H(x_j, y_j)}{u_j \frac{\partial H}{\partial x}(x_j, y_j) + v_j \frac{\partial H}{\partial y}(x_j, y_j)}$

7:   $s_j \leftarrow s_j + \delta$

8:   $(x_j, y_j) \leftarrow (x_j, y_j) + \delta (u_j, v_j)$

9: **end for**

10: **end for**

► Interpolate

11: $u^\circ \leftarrow \text{FFT}_N \left( (u_j)_{j=0}^{2N} \right), \quad v^\circ \leftarrow \text{FFT}_N \left( (v_j)_{j=0}^{2N} \right)$

12: $\check{s} \leftarrow \text{FFT}_N \left( (s_j)_{j=0}^{2N} \right)$

13: **return** $(u^\circ, v^\circ)$, $\check{s}$

The following proposition establishes the exponential convergence of this method w.r.t. the approximation degree N and its quasi-linear arithmetic complexity.

**PROPOSITION 3.5.** There exists $\eta > 0$ (depending on H and h) such that if the analytic initial guess $(x^\circ, y^\circ) \in (C_{2\pi}^0)^2$ is $\eta$-close to $(x^*, y^*)$, then there is $\kappa > 1$ such that for any $N \in \mathbb{N}$, Algorithm OVALAPPROXREFINE(H, h, x°, y°, N) computes a degree N TP š in O(N log N) arithmetic operations, and the corresponding curve (x̃, ỹ), cf. (3.6), converges exponentially fast to the projected parameterization $(x^#, y^#)$ of $\Gamma(h)$:

$$ \|\check{s} - s^\#\|_\infty, \|\tilde{x} - x^\#\|_\infty, \|\tilde{y} - y^\#\|_\infty = O(\kappa^{-N}) \quad \text{as} \quad N \to \infty. $$

**PROOF.** Let $\eta' > 0$ denote the $\eta$ of Lemma 3.4. Clearly, one can chose an $\eta \in (0, \eta']$ such that if $(x^\circ, y^\circ)$ is $\eta$-close to $(x^*, y^*)$, then $(u^\circ, v^\circ) = \mathcal{I}_N \left( \frac{\nabla H(x^\circ, y^\circ)}{\|\nabla H(x^\circ, y^\circ)\|_2} \right)$ (lines 3 and 11) is $\eta'$-close to $\frac{\nabla H(x^*, y^*)}{\|\nabla H(x^*, y^*)\|_2}$, since by assumption $\nabla H$ does not vanish over $\Gamma(h)$. Then the analytic shift $s^\#$ and the projected parameterization $(x^#, y^#)$ 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:

$$ \|\check{s} - s^\#\|_\infty \le \|\mathcal{I}_N(s^\#) - s^\#\|_\infty + \Lambda_N \sup_{0 \le j \le 2N} |\check{s}(t_j) - s^\#(t_j)|, $$

where $t_j = \frac{2j\pi}{2N+1}$ are the equispaced interpolation points.

*Error of interpolation.* Since $s^\# \in C_{2\pi}^0$ is $2\pi$-periodic and analytic over $\mathbb{R}$ by Lemma 3.4, there is a $\rho > 0$ such that the open horizontal strip $\mathbb{R} + (-\rho, \rho)i$ avoids all possible singularities of $s^\#$ in the
complex plane. Then the interpolants $\mathcal{I}_N(s^#)$ converge uniformly to $s^#$ over $\mathbb{R}$ in $O(\kappa^{-N})$ for any $\kappa \in (1, e^\rho)$ (see [Wright et al. 2015, Thm. 4.2]).

*Interpolation of the error.* Since the interpolation in line 12 uses approximations $s_j$ of the exact points $s^#(t_j)$, the resulting “interpolation of the error” depends on the accuracy reached after $k = \lceil \log N \rceil$ Newton iterations (lines 2–10). By Eq. (B.1) in the proof of Lemma 3.4, we have:

$$|s_j - s^#(t_j)| \le |\mathcal{N}_{x_j, y_j, u_j, v_j}^k(0) - s^#(t_j)| \le \frac{|s^#(t_j)|}{2^{2k-1}} \le 2^{-N} \frac{\|s^#\|_\infty}{2}.$$

*Combining both errors, we conclude that* $\|s - s^#\|_\infty = O(\kappa^{-N})$ for any $\kappa \in (1, \min(e^\rho, 2))$. It is also clear from Eq. (3.4) and (3.6) that the same convergence happens for $(\tilde{x}, \tilde{y})$ to $(x^#, y^#)$.

*Complexity.* The sampling of initial points and reinterpolation through $\text{FFT}_N/\text{IFFT}_N$ requires $O(N \log N)$ arithmetic operations. For the projection step (lines 2–10), each of the $2N + 1$ points is refined through $\lceil \log N \rceil$ Newton iterations, where the cost of a single iteration is $O(1)$. This justifies the overall complexity in $O(N \log N)$ arithmetic operations. □

# 4 A POSTERIORI VALIDATION OF THE OVAL $\Gamma(h)$

The purpose of Algorithm OVALVALID presented in this section is to enclose the oval $\Gamma(h)$ in a validated tube $T(\tilde{x}, \tilde{y}, u^\circ, v^\circ, \varepsilon)$ around the approximate parameterization $(\tilde{x}, \tilde{y})$ (3.6) containing $2\pi$-periodic curves:

$$T(\tilde{x}, \tilde{y}, u^\circ, v^\circ, \varepsilon) = \left\{ (x, y) \in (C_{2\pi}^0)^2 \mid \exists s \in C_{2\pi}^0, \|s\|_\infty \le \varepsilon \text{ and } \begin{cases} x(t) = \tilde{x}(t) + s(t)u^\circ(t) \\ y(t) = \tilde{y}(t) + s(t)v^\circ(t) \end{cases} \right\}.$$

The width $\varepsilon$ 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 $\tilde{s}$ computed by OVALAPPROXREFINE and the exact shift $s^#$ defining the parameterization $(x^#, y^#)$ of $\Gamma(h)$ (see Definition 3.3). Under a few additional conditions also checked by OVALVALID (winding number, etc. $^{10}$), this ensures that $(x^#, y^#)$ exists, satisfies Lemma 3.4, and lies in the tube $T(\tilde{x}, \tilde{y}, u^\circ, v^\circ, \varepsilon)$.

Using Equation (3.3), we construct a polynomial Newton-like operator $\mathcal{N} \in C_{2\pi}^0[X]$ acting on a function $s \in C_{2\pi}^0$:

$$\mathcal{N}(s) = s - \tilde{a} \left( H(x^\circ + s u^\circ, y^\circ + s v^\circ) - h \right), \quad (4.1)$$

for which $s^#$ is the unique fixed point in some small neighborhood of $\tilde{s}$. Here,

$$\tilde{a}(t) \approx \frac{1}{u^\circ(t) \frac{\partial H}{\partial x}(\tilde{x}(t), \tilde{y}(t)) + v^\circ(t) \frac{\partial H}{\partial y}(\tilde{x}(t), \tilde{y}(t))},$$

is a TP built by interpolation to make the derivative of $\mathcal{N}$,

$$\mathcal{N}'(s) = 1 - \tilde{a} \left( u^\circ \frac{\partial H}{\partial x}(x^\circ + s u^\circ, y^\circ + s v^\circ) + v^\circ \frac{\partial H}{\partial y}(x^\circ + s u^\circ, y^\circ + s v^\circ) \right),$$

small in the neighborhood of $s = \tilde{s}$ (that is, geometrically, in the neighborhood of the approximate curve $(\tilde{x}, \tilde{y})$).

Contrary to the a posteriori validation for the division presented in Section 2.3.1, the Newton-like operator $\mathcal{N}$ in (4.1) is not of degree 1 in $s$ unless $H$ itself has degree 1. As a consequence, the value of $\|\mathcal{N}'(s)\|_\infty$ depends on $s$, so that the conditions of Theorem 2.5 must be satisfied with a finite radius $r > 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:

<footnote>$^{10}$Although 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.</footnote>
* a procedure computing an upper bound $\lambda(r)$ for $\sup_{\|s-\tilde{s}\|_{\infty}\le r} \|N'(s)\|_{\infty}$ for any $r \ge 0$,

* an upper bound $d$ for the defect $\|N(\tilde{s}) - \tilde{s}\|_{\infty}$,

NEWTONBALL($d, \lambda$) either computes a radius $r$ and a Lipschitz constant $\lambda = \lambda(r)$ satisfying the hypotheses of Theorem 2.5, or returns an error.

**Algorithm 4 OVALVALID($H, h, x_{ini}, y_{ini}, \varepsilon_0, x_r, y_r, \zeta, x^\circ, y^\circ, u^\circ, v^\circ, \tilde{s}, N$)**

**Input:** $H, h, \zeta = \pm 1$, $(x_{ini}, y_{ini}) \in \mathbb{R}^2, \varepsilon_0 > 0$, $(x_r, y_r) \in \mathbb{R}^2$, and degree $N$ TPs $x^\circ, y^\circ, u^\circ, v^\circ, \tilde{s}$

**Output:** a bound $\varepsilon > 0$ defining a tube $T(\tilde{x}, \tilde{y}, u^\circ, v^\circ, \varepsilon)$ around $(\tilde{x}, \tilde{y})$ containing $(x^\#, y^\#)$

▷ Build TP ă defining the Newton-like operator N

1: $(x_j)_{j=0}^{2N} \leftarrow \text{IFFT}_N(\tilde{x})$ and $(y_j)_{j=0}^{2N} \leftarrow \text{IFFT}_N(\tilde{y})$

2: $(u_j)_{j=0}^{2N} \leftarrow \text{IFFT}_N(u^\circ)$ and $(v_j)_{j=0}^{2N} \leftarrow \text{IFFT}_N(v^\circ)$

3: ă $\leftarrow \text{FFT}_N\left(\left(\frac{1}{u_j \frac{\partial H}{\partial x}(x_j, y_j) + v_j \frac{\partial H}{\partial y}(x_j, y_j)}\right)_{j=0}^{2N}\right)$

▷ Bound operator norm of N, using rigorous operations with intervals and RTAs

4: $x \leftarrow x^\circ + \tilde{s}u^\circ$ and $y \leftarrow y^\circ + \tilde{s}v^\circ$ ▷ $(x, y)$ RTAs for $(\tilde{x}, \tilde{y})$

5: $x' \leftarrow x^{\circ'} + \tilde{s}'u^\circ + \tilde{s}'u^{\circ'}$ and $y' \leftarrow y^{\circ'} + \tilde{s}'v^\circ + \tilde{s}'v^{\circ'}$

6: $d \leftarrow \text{BOUND}(s)$ where $s \leftarrow \tilde{a}(h - H(\mathbf{x}, \mathbf{y}))$ ▷ defect

7: $\lambda_0 \leftarrow \text{BOUND}(\ell_0)$ where $\ell_0 \leftarrow 1 - \tilde{a}\left(u^\circ \frac{\partial H}{\partial x}(\mathbf{x}, \mathbf{y}) + v^\circ \frac{\partial H}{\partial y}(\mathbf{x}, \mathbf{y})\right)$

8: $\alpha \leftarrow \text{BOUND}(\tilde{a})$

9: $X \leftarrow \text{BOUND}(\mathbf{x}) \cdot [-1, 1]$ and $Y \leftarrow \text{BOUND}(\mathbf{y}) \cdot [-1, 1]$

10: $U \leftarrow \text{BOUND}(u^\circ) \cdot [-1, 1]$ and $V \leftarrow \text{BOUND}(v^\circ) \cdot [-1, 1]$

11: function $\lambda(r)$ ▷ Lipschitz constant of N over $\bar{B}(\tilde{s}, r)$ in $C_{2\pi}^0$

12: $X_r \leftarrow X + rU$ and $Y_r \leftarrow Y + rV$

13: $\lambda_1 \leftarrow \text{mag} \left( U^2 \frac{\partial^2 H}{\partial x^2} (X_r, Y_r) + 2UV \frac{\partial^2 H}{\partial x \partial y} (X_r, Y_r) + V^2 \frac{\partial^2 H}{\partial y^2} (X_r, Y_r) \right)$

14: return $\lambda_0 + \alpha r \lambda_1$

15: end function

16: $(r^+, \lambda^+) \leftarrow \text{NEWTONBALL}(d, \lambda)$ ▷ Find stable ball

▷ Validated error bound $\varepsilon$

17: $\varepsilon \leftarrow d/(1 - \lambda^+)$

▷ Check initial conditions

18: **if** $(\mathbf{x}(0) - \varepsilon u^\circ(0) - x_{ini})^2 + (\mathbf{y}(0) - \varepsilon v^\circ(0) - y_{ini})^2 > \varepsilon_0^2$ **or**

19: $(\mathbf{x}(0) + \varepsilon u^\circ(0) - x_{ini})^2 + (\mathbf{y}(0) + \varepsilon v^\circ(0) - y_{ini})^2 > \varepsilon_0^2$, **then return** FAIL **end if**

▷ Check the winding number of the curve

20: $\mathbf{r}2 \leftarrow (\mathbf{x} - x_r)^2 + (\mathbf{y} - y_r)^2$

21: **if** $\neg \text{ISPOSITIVE} (\mathbf{r}2 - \varepsilon^2 \text{BOUND}(u^\circ + v^\circ))$ **then return** FAIL **end if**

22: $\mathcal{S} \leftarrow \frac{1}{2\pi} \int_0^{2\pi} \frac{(\mathbf{x}(t) - x_r)\mathbf{y}' - (\mathbf{y}(t) - y_r)\mathbf{x}'}{\mathbf{r}2(t)} dt$

23: **if** $\mathcal{S} \cap \mathbb{Z} \neq \{\zeta\}$ **then return** FAIL **end if**

24: return $\varepsilon$

**PROPOSITION 4.1.** Given a point $(x_r, y_r)$ inside the oval $\Gamma(h)$, and a ball $\bar{B}((x_{ini}, y_{ini}), \varepsilon_0)$ enclosing initial conditions such that it intersects $H^{-1}(h)$ in the $\Gamma(h)$ component only, Algorithm OVALVALID($H, h$,
$x_{ini}, y_{ini}, \varepsilon_0, x_r, y_r, x^\circ, y^\circ, u^\circ, v^\circ, \check{s}, N$), if it does not fail, computes $\varepsilon$ such that $s^\#$ of Definition 3.3 exists, satisfies Lemma 3.4 with winding number around $(x_r, y_r)$ equal to $\zeta$, and $\|\check{s} - s^\#\|_\infty \le \varepsilon$.

Moreover, for a sufficiently large degree N, the process never fails and returns a bound $\varepsilon$ that converges exponentially fast to 0 w.r.t. N, in $O(N \log(N))$ arithmetic operations.

**PROOF.** We apply Theorem 2.5 over the space $C_{2\pi}^0$ to compute an upper bound $\varepsilon$ for $\|s^\#\|_\infty$, where $s^\# \in C_{2\pi}^0$ 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 $\check{a}$ interpolating the function which represents the inverse differential. Next, the algorithm computes a rigorous bound $d$ on the defect $\|N(\check{s}) - \check{s}\|_\infty$ (line 6), and produces a procedure $\lambda$ s.t. $\lambda(r)$ is a Lipschitz constant for $N$ over $\bar{B}(\check{s}, r)$ in $C_{2\pi}^0$ w.r.t. to the $\|\cdot\|_\infty$ norm (lines 7–15). The formula for $\lambda(r)$ in line 14 comes from Lemma 2.6, with $U, V$ enclosures for the range of $u^\circ, v^\circ$, and $X_r, Y_r$ enclosures for the range of $x^\circ + su^\circ = \check{x} + (s - \check{s})u^\circ$ and $y^\circ + sv^\circ = \check{y} + (s - \check{s})v^\circ$, for all $s \in \bar{B}(\check{s}, r)$. After that, the auxiliary routine NEWTONBALL (line 16) described in Appendix C rigorously computes a suitable radius $r^+$ satisfying the hypotheses of Theorem 2.5, if possible. We thus deduce the existence and uniqueness of $s^\# \in C_{2\pi}^0$ (and hence $(x^\#, y^\#)$) together with a validated bound $\varepsilon$ (line 17).

*   **Initial conditions.** By hypothesis, the ball $\bar{B}((x_{ini}, y_{ini}), \varepsilon_0)$ in $\mathbb{R}^2$ intersects $H^{-1}(h)$ in its component $\Gamma(h)$ only. Hence, by continuity, checking that $(x^\#(0), y^\#(0))$ lies in that ball is sufficient to ensure that $(x^\#(t), y^\#(t))$ belongs to the component $\Gamma(h)$ for all $t \in [0, 2\pi]$. 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 $T(\check{x}, \check{y}, u^\circ, v^\circ, \varepsilon)$ all avoid the point $(x_r, y_r)$ inside $\Gamma(h)$. Then, the winding number of $(\check{x}, \check{y})$ is rigorously computed in line 22 and compared to the reference value $\zeta$ at line 23. Since the tube avoids $(x_r, y_r)$, the winding number remains the same for all curves in it (so in particular $(x^\#, y^\#)$).

*   **Analyticity.** The Banach fixed-point based argument used previously only guarantees the continuity of $(x^\#, y^\#)$. To gain more regularity, we need to consider a larger function space, namely the Banach space of continuous $2\pi$-periodic functions defined over a horizontal strip $S_\rho := \mathbb{R} + [-\rho, \rho]i$ of $\mathbb{C}$ for some small $\rho > 0$, equipped with the norm: $\|\varphi\|_{S_\rho} = \sup_{z \in S_\rho} |\varphi(z)|^{11}$. The a posteriori validation argument based on Theorem 2.5 applies the same way, except that we need to replace $\|\cdot\|_\infty$ with $\|\cdot\|_{S_\rho}$ in all the bounds defining $d$ and $r \mapsto \lambda(r)$. If we choose $\rho$ small enough, the inequalities $\lambda(r^+) < 1$ and $d + r^+\lambda(r^+) < r^+$ guaranteed by NEWTONBALL still hold for the same $r^+$, which means that Theorem 2.5 now proves the existence of a unique continuous periodic solution to (3.3) extending $s^\#$ over $S_\rho$. Moreover, the iterations $N^k(0)$ are all trigonometric polynomials (hence analytic) and they converge uniformly to $s^\#$ over $S_\rho$ as $k \to \infty$. As a consequence, the limit $s^\#$ is analytic over $\mathbb{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 $\check{s}$ computed by OVALAPPROXREFINE tends exponentially fast to $s^\#$. This holds for the $\|\cdot\|_\infty$ norm, but also for the $\|\cdot\|_{\ell^1}$ norm since they are TPs of degree $O(N)$. As a consequence, the computed bounds $d$ and $\lambda_0$ tend to 0 exponentially fast, while the other bounds $\alpha, X, Y, U, V$ and $\lambda_1$ (for given $r$) tend to the mathematical values obtained by replacing the approximations $\check{s}, \check{x}, \check{y}$ by $s^\#, x^\#, y^\#$. Therefore, in the routine NEWTONBALL, for $N$ sufficiently large, the initial radius $r_{max}$ computed in line 1 is bounded away from 0 and satisfies $d + \lambda(r_{max})r_{max} < r_{max}$. Hence, NEWTONBALL necessarily returns an $r$ satisfying $d + \lambda(r)r < r$, and then OVALVALID computes $\varepsilon = d/(1 - \lambda^+) \le 2d$, which tends exponentially fast to 0. To conclude, similar arguments show that the remaining tests (initial conditions and winding number) never fail for $N$ sufficiently large.

<footnote n="11">¹¹Note that this argument is purely theoretical: there is no need to reimplement Algorithm OVALVALID to actually compute with this norm $\|\cdot\|_{S_\rho}$.</footnote>
• *Complexity.* The asymptotic complexity w.r.t. the degree *N* is determined by the *FFT*<sub>*N*</sub>/<i>IFFT</i><sub>*N*</sub> 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 *N*, since the function *r* ↦→ λ(*r*) (lines 11–15) involves intervals rather than RTAs. □

5 RIGOROUS ENCLOSURE OF THE ABELIAN INTEGRAL ℑ(h)

According to Proposition 4.1, the tube T(\x, \y, u, v, \epsilon) computed by OVALVALID contains a unique analytic parameterization (x^#, y^#) of Γ(h) with correct winding number, so that we have:

$$
\mathfrak{I}(h) = \int_{0}^{2\pi} \frac{P(x^{\#}(t), y^{\#}(t))y^{\prime\#}(t) - Q(x^{\#}(t), y^{\#}(t))x^{\prime\#}(t)}{\mu(x^{\#}(t), y^{\#}(t))} \mathrm{d}t. \quad (5.1)
$$

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

$$
I_{\tilde{s}} = \int_{0}^{2\pi} \frac{P(\tilde{x}(t), \tilde{y}(t))\tilde{y}'(t) - Q(\tilde{x}(t), \tilde{y}(t))\tilde{x}'(t)}{\mu(\tilde{x}(t), \tilde{y}(t))} \mathrm{d}t, \quad (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:

$$
c(x, y) \mathrm{d}x \wedge \mathrm{d}y = \left[ \frac{\partial(P/\mu)}{\partial x}(x, y) + \frac{\partial(Q/\mu)}{\partial y}(x, y) \right] \mathrm{d}x \wedge \mathrm{d}y, \quad (5.3)
$$

over the region located between the two curves (x̃, ỹ) and (x^#, y^#), which is safely overapproximated
by the tube T(x̃, ỹ, u^o, v^o, ε).

Algorithm 5 ABINTVALIDQUAD(P, Q, μ, x°, y°, u°, v°, š, ε, N)

Input: P, Q, μ as in Problem 1.1, degree N TPs x°, y°, u°, v°, š, and ε > 0 computed by OVALVALID

Output: interval enclosure I for ℑ(h) (1.1)

► Compute rigorous integral along approximate curve using RTAs

1: x ← x° + š u° and y ← y° + š v° ► (x, y) RTAs for (x̃, ỹ)

2: x' ← x°' + š' u° + š' u°' and y' ← y°' + š' v° + š' v°'

3: P ← P(x, y), Q ← Q(x, y), and μ ← μ(x, y)

4: E ← (Py - Qx) / μ

5: I<sub>š</sub> ← ∫<sub>0</sub><sup>2π</sup> E(t) dt

► Compute integration error in the tube using Stokes' theorem

6: x<sub>ε</sub> ← (x° + š u°, ε BOUND(u°)) and y<sub>ε</sub> ← (y° + š v°, ε BOUND(v°))

7: c ← $\frac{\partial(P/\mu)}{\partial x}$ + $\frac{\partial(Q/\mu)}{\partial y}$ ► symbolically as a rational fraction

8: C ← BOUND(c(x<sub>ε</sub>, y<sub>ε</sub>))

9: B<sub>1</sub> ← BOUND(x'u° - y'u°)

10: B<sub>2</sub> ← BOUND(u°v° - v°u°)

11: δ ← 2πε(2B₁ + εB₂)C

► Return interval enclosure

12: return I ← I<sub>š</sub> + [-δ, δ]
**PROPOSITION 5.1.** Let P, Q, μ be as in Problem 1.1. If ABINTVALIDQUAD(P, Q, μ, x°, y°, u°, v°, š, ε, N) with ε > 0 computed by OVALVALID does not fail, then it returns an interval enclosure I of ℑ(h).
Moreover, if N is sufficiently large, then it does never fail and computes an interval I whose diameter tends exponentially fast to 0 w.r.t. N, in O(N log(N)) arithmetic operations.

**PROOF.** First, an interval enclosure $I_{\tilde{s}}$ of the integral $I_{\tilde{s}}$ (5.2) along the approximate curve $(\tilde{x}, \tilde{y})$ is computed rigorously by using RTA representations and operations for the integrand (lines 3 and 4), and then integrating it over $[0, 2\pi]$ (line 5). Next, lines 6–11 compute a bound $\delta$ for the error between $I_{\tilde{s}}$ and the exact value $\mathfrak{I}(h)$ taken along $(x^{\#}, y^{\#})$ of $\Gamma(h)$, so that ABINTVALIDQUAD returns a validated enclosure $I = I_{\tilde{s}} + [-\delta, \delta]$ for $\mathfrak{I}(h)$. The strategy to obtain $\delta$ is detailed below.

• *Error bound $\delta$ with Stokes’ theorem.* Let $(x, y)$ denote the standard coordinates of the plane $\mathbb{R}^2$, and add a third one, $z$, living in $\mathbb{R}/2\pi\mathbb{Z}$, corresponding to the time $t$. The approximate parameterization $(\tilde{x}, \tilde{y})$ traces a smooth simple closed curve $t \in \mathbb{R}/2\pi\mathbb{Z} \mapsto (x^{\circ}(t) + \tilde{s}(t)u^{\circ}(t), y^{\circ}(t) + \tilde{s}(t)v^{\circ}(t), t)$ in $\mathbb{R}^2 \times \mathbb{R}/2\pi\mathbb{Z}$. Similarly, $(x^{\#}, y^{\#})$ traces $t \in \mathbb{R}/2\pi\mathbb{Z} \mapsto (x^{\circ}(t) + s^{\#}(t)u^{\circ}(t), y^{\circ}(t) + s^{\#}(t)v^{\circ}(t), t)$. Since $s^{\#}$ is analytic, it only has a finite number of zeros in $[0, 2\pi]$. Hence, since $(u^{\circ}, v^{\circ})$ never vanishes, the two curves only intersect each other a finite number of times. The “strip” $S$ between them,

$$S = \{(x^{\circ}(t) + su^{\circ}(t), y^{\circ}(t) + sv^{\circ}(t), t), \quad t \in \mathbb{R}/2\pi\mathbb{Z}, s \in [\tilde{s}(t), s^{\#}(t)] \text{ or } [s^{\#}(t), \tilde{s}(t)]\},$$

is a finite collection of two-dimensional cells (with two singularities each). Choosing the appropriate orientation for each cell, the difference between $\mathfrak{I}(h)$ and $I_{\tilde{s}}$ is equal to the integral of the 1-form $(P/\mu)dy - (Q/\mu)dx$ along the oriented boundary $\partial S$ of $S$:

$$\mathfrak{I}(h) - I_{\tilde{s}} = \int_{\partial S} \frac{P(x, y)dy - Q(x, y)dx}{\mu(x, y)}.$$

Using Stokes' theorem [Lee 2013, Thm. 14.9], and since P and Q do not depend on z, we obtain:

$$|\mathfrak{I}(h) - I_{\tilde{s}}| \le \int_S |c(x, y)| dx \wedge dy,$$

with $c(x, y)$ defined as in Equation (5.3).

A bound $C$ for the integrand $|c(x, y)|$ is computed in lines 6–8. Moreover, the integration domain $S$ can be safely overapproximated by the geometric realization $T_{\varepsilon}$ of the tube $T(\tilde{x}, \tilde{y}, u^{\circ}, v^{\circ}, \varepsilon)$ in $\mathbb{R}^2 \times \mathbb{R}/2\pi\mathbb{Z}$, which is the two-dimensional manifold parameterized by:

$$\sigma: \mathbb{R}/2\pi\mathbb{Z} \times [-\varepsilon, \varepsilon] \longrightarrow T_{\varepsilon} \\ (t, s) \longmapsto \begin{pmatrix} \tilde{x}(t) + s u^{\circ}(t) \\ \tilde{y}(t) + s v^{\circ}(t) \\ t \end{pmatrix}.$$

Hence, the difference $\mathfrak{I}(h) - I_{\tilde{s}}$ is bounded by:

$$|\mathfrak{I}(h) - I_{\tilde{s}}| \le C \int_{T_{\varepsilon}} |dx \wedge dy| \\ = C \int_{0}^{2\pi} \int_{-\varepsilon}^{\varepsilon} \left| \begin{pmatrix} \tilde{x}'(t) & u^{\circ}(t) \\ \tilde{y}'(t) & v^{\circ}(t) \end{pmatrix} \right| + s \left| \begin{pmatrix} u^{\circ'}(t) & u^{\circ}(t) \\ v^{\circ'}(t) & v^{\circ}(t) \end{pmatrix} \right| ds dt \\ \le 2\pi\varepsilon (2B_1 + \varepsilon B_2) C = \delta,$$

$$\text{with } B_1 \ge \sup_{0 \le t \le 2\pi} \left| \begin{pmatrix} \tilde{x}'(t) & u^{\circ}(t) \\ \tilde{y}'(t) & v^{\circ}(t) \end{pmatrix} \right| \text{ and } B_2 \ge \sup_{0 \le t \le 2\pi} \left| \begin{pmatrix} u^{\circ'}(t) & u^{\circ}(t) \\ v^{\circ'}(t) & v^{\circ}(t) \end{pmatrix} \right|,$$

computed in lines 9 and 10. Hence, $\mathfrak{I}(h) = I_{\tilde{s}} + (\mathfrak{I}(h) - I_{\tilde{s}}) \in I_{\tilde{s}} + [-\delta, \delta] = I$.
• *Convergence*. When $N$ is large enough, $(\tilde{x}, \tilde{y})$ is sufficiently close to $\Gamma(h)$, and since $\mu$ does not vanish in a neighborhood of this oval, all the functions involved in the computation of $I_{\tilde{s}}$ (lines 3–5) are analytic; hence the diameter of this interval tends exponentially fast to 0 w.r.t. $N$.

The same arguments hold for the computation of Stokes' error bound $\delta$ in lines 6–11. Note that $B_1, B_2$ and $C$ are bounded since they converge to the corresponding $\| \cdot \|_{\ell^1}$ norms with $(\tilde{x}, \tilde{y})$ and $(u^\circ, v^\circ)$ replaced by $(x^\#, y^\#)$ and $\frac{\nabla H(x^\circ, y^\circ)}{\|\nabla H(x^\circ, y^\circ)\|_2}$, respectively. Since $\varepsilon$ tends to 0 exponentially fast w.r.t. $N$ according to Proposition 4.1, so does $\delta$ (line 11), and consequently the total diameter of $I$ in line 12.

• *Complexity*. The asymptotic complexity w.r.t. the globally set RTA degree $N$ is determined by the FFT$_N$/$IFFT_N$ routines used for TP/RTA multiplications and divisions. $\square$

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.

<caption>**Algorithm 6** ABINTVALID($H, h, \Sigma, P, Q, \mu, x_{ini}, y_{ini}, \varepsilon_0, x_r, y_r, \zeta, N_0, N$)</caption>

**Input:** Potential $H$, level $h$, transversal $\Sigma$, $P, Q, \mu$ as in Problem 1.1, $(x_{ini}, y_{ini}) \in \mathbb{R}^2, \varepsilon_0 > 0$, $(x_r, y_r) \in \mathbb{R}^2, \zeta = \pm 1, N_0 \in \mathbb{N}$ and $N \ge N_0$

**Output:** interval enclosure $I$ of $\mathfrak{I}(h)$

1: $(x^\circ, y^\circ) \leftarrow \text{OVALAPPROXINIT}(H, \Sigma, x_{ini}, y_{ini}, N_0)$

2: $(u^\circ, v^\circ), \check{s} \leftarrow \text{OVALAPPROXREFINE}(H, h, x^\circ, y^\circ, N)$

3: $\varepsilon \leftarrow \text{OVALVALID}(H, h, x_{ini}, y_{ini}, \varepsilon_0, x_r, y_r, \zeta, x^\circ, y^\circ, u^\circ, v^\circ, \check{s}, N)$

4: **return** $I \leftarrow \text{ABINTVALIDQUAD}(P, Q, \mu, x^\circ, y^\circ, u^\circ, v^\circ, \check{s}, \varepsilon, N)$

**THEOREM 5.2.** Let $H, P, Q, \mu$ be as in Problem 1.1, $h$ be a regular level value on a transversal $\Sigma$, the ball of center $(x_{ini}, y_{ini})$ and radius $\varepsilon_0$ contain a unique point in $H^{-1}(h) \cap \Sigma$, the point $(x_r, y_r)$ be in the interior of the oval $\Gamma(h)$, and $\zeta \in \{-1, 1\}$ be the winding number of $\Gamma(h)$ w.r.t. $(x_r, y_r)$ in the orientation prescribed by the Hamiltonian. Then for a sufficiently large degree $N_0$ (depending on the preceding input data but not on the target accuracy), as the approximation degree $N \to \infty$, Algorithm ABINTVALID($H, h, \Sigma, P, Q, \mu, x_{ini}, y_{ini}, \varepsilon_0, x_r, y_r, \zeta, N_0, N$) never fails and returns an interval enclosure $I$ of $\mathfrak{I}(h)$ in $O(N \log N)$ arithmetic operations. The diameter of $I$ tends to 0 exponentially fast, i.e., there exists $\kappa > 1$ depending on $H, h, \Sigma, P, Q, \mu, x_{ini}, y_{ini}, \varepsilon_0, x_r, y_r, \zeta, N_0$ such that $\text{diam}(I) = O(\kappa^{-N})$.

# 6 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.

## 6.1 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 package<fnref n="12" /> which provides a very modular framework to approximate and manipulate functions in

<footnote n="12">¹²https://juliaapproximation.github.io/ApproxFun.jl/latest/</footnote>
various bases, including Fourier: A degree $N$ TP is stored as the length $2N + 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 library<sup>13</sup>). Since this data structure is parameterized by the type of the coefficients, TPs can also be defined with interval coefficients<sup>14</sup>. Therefore, in *RigorousFourier.jl*, RTAs are implemented as pairs made up of a TP with interval coefficients together with a floating-point 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 $N$ is set by the user so that TPs and RTAs are automatically truncated to degree $N$ after each operation. Note that the degree $N' = 2^{\lceil \log_2(N) \rceil} \ge N$ is used internally for FFT based operations (such as multiplication) for the sake of efficiency.

**REMARK 6.1.** Using floating-point coefficients for TPs $\tilde{f}$ adds rounding errors to the approximation error w.r.t. to the represented functions $f$. Similarly, using floating-point interval coefficients for RTAs $f$ enlarges the set of functions contained in $f$, due to the growth of the interval widths. However, it can be easily proved that with $p \approx N$ 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 $O(N \cdot 2^{-N})$ for some $c > 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 [Cyranka 2014; Lessard 2018], overestimation of tail coefficients [Lessard 2018], wrapping effect when input coefficients are given as large intervals [Liu and Kreinovich 2010]). 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:

(a) 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., [Brisebarre et al. 2020]). Using $p \approx N$ bits of precision, the claimed bounds are in $O(N^2 2^{-N})$, 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.

(b) 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 (a). We however also provide the timings for this method (b) since it is the only fully rigorous one in our prototype code.

<footnote n="13"><sup>13</sup>https://www.mpfr.org/</footnote>

<footnote n="14"><sup>14</sup>Intervals are provided by *IntervalArithmetic.jl*, https://juliaintervals.github.io/IntervalArithmetic.jl/latest/intro/.</footnote>
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` package<fnref n="15" /> for multivariate polynomials used in the data structures representing integrable systems, oval families and perturbations;

* the `DifferentialEquations.jl` package<fnref n="16" /> that provides a numerical ODE solver to instantiate the `DSOLVE` routine used by Algorithm `OVALAPPROXINIT` to compute the initial guess $(x^°, y^°)$. As of today, it uses by default a Tsitouras 5/4 Runge-Kutta method with 4th order interpolant [Tsitouras 2011]. 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 (a)), the convergence results of the preceding sections in $O(κ^{-N})$ for some $κ ∈ (1, 2)$ still hold for the floating-point implementation of `ABIINTVALID`. Combining this with the bit complexity of $O(N \log N)$ for arithmetic operations on N-bit floating-point numbers [Harvey and van der Hoeven 2021], a refined version of Theorem 5.2 says that Algorithm `ABIINTVALID`, using degree N TPs/RTAs and floating-point arithmetic with N bits of precision (neglecting considerations about the exponent range), computes a rigorous enclosure of $\mathfrak{I}(h)$ with $O(N)$ bits of accuracy in $O(N^2 \log^2 N)$ bit operations.

## 6.2 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 [Johnson and Tucker 2011], 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 [Neumaier 2003]. Our strategy to evaluate an Abelian integral using CAPD is to augment the plane $(x, y)$ with a third dimension $I$ (the integral), and then add a third component to the vector field, representing $\dot{I}$, which is equal to the integrand of (1.1), where $dx$ and $dy$ are replaced by the values of $\dot{x}$ and $\dot{y}$ 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 library<fnref n="17" /> [Johansson 2017], 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 [Johansson 2018]. Although it cannot rigorously approximate general algebraic functions yet (in order to approximate the oval),

<footnote n="15">¹⁵https://github.com/JuliaAlgebra/TypedPolynomials.jl.</footnote>

<footnote n="16">¹⁶https://diffeq.sciml.ai/v2.0/</footnote>

<footnote n="17">¹⁷https://arblib.org/</footnote>
we used it, as a SageMath<fnref n="18" /> 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 [Laírez et al. 2019]. 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.

# 7 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 [Johnson 2011] to claim $H(4) \ge 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 $H$ makes it possible to divide the ovals into arcs $x \mapsto y(x)$ and $y \mapsto x(y)$ with simple formulas, so that $\mathfrak{I}(h)$ can easily be computed as a classical integral over a segment of $\mathbb{R}$ with an explicit integrand.

* The second example (Section 7.2) is a cubic Hamiltonian system used in [Li et al. 2009] to demonstrate $\mathcal{Z}(3) \ge 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.

## 7.1 Johnson's symmetric quartic system revisited

In [Johnson 2011], T. Johnson constructed a perturbed quartic pseudo-Hamiltonian vector field of the form:

$$ \begin{cases} \dot{x} = -4y^2(y^2 - Y_0), \\ \dot{y} = 4xy(x^2 - X_0) + \epsilon g(x, y), \end{cases} \qquad (7.1) $$

where $X_0 = \frac{9}{10}$, $Y_0 = \frac{11}{10}$, and $g$ is a degree 4 perturbation,

$$ g(x, y) = \alpha_{00} + \alpha_{20}x^2 + \alpha_{22}x^2y^2 + \alpha_{40}x^4 + \alpha_{04}y^4, \qquad (7.2) $$

with well-chosen coefficients $\alpha_{ij}$. Using a rigorous validation integration routine, he claimed to prove the existence of 26 limit cycles, thus surpassing the previously known record $H(4) \ge 22$ [Christopher and Li 2007]. 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 $\alpha_{ij}$, one can prove the following Theorem establishing $H(4) \ge 24^{19}$.

**THEOREM 7.1.** The quartic system (7.1) with $X_0 = \frac{9}{10}$, $Y_0 = \frac{11}{10}$, and the degree 4 perturbation $g$ defined as in (7.2) with coefficients:

$$ \begin{aligned} \alpha_{00} &= -0.78622148667854837664, & \alpha_{20} &= 0.87723523612653436051, & \alpha_{22} &= 1, \\ \alpha_{40} &= 0.23742713894293038223, & \alpha_{04} &= -0.21823846173078863753, \end{aligned} $$

<footnote n="18">¹⁸https://www.sagemath.org/</footnote>

<footnote n="19">¹⁹In the meantime, Prohens and Torregrosa [Prohens and Torregrosa 2019] showed $H(4) \ge 28$ by a different kind of computer-aided proof (and for a different system).</footnote>
has at least 24 limit cycles.

A particularity of the potential function H associated to this system is that the ovals can be parameterized explicitly using arcs x(y) and y(x) (see Section 7.1.1 below). This was the strategy adopted in [Johnson 2011] 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.

**7.1.1 Potential function, ovals and symmetries.** The quartic system (7.1) admits the following potential function:

$$H(x, y) = (x^2 - X_0)^2 + (y^2 - Y_0)^2, \qquad (7.3)$$

together with the rescaling factor $\mu(x, y) = y$.

The level set associated to the parameter $h \ge 0$, represented in the $(x^2, y^2)$ plane, is the portion of the circle of center $(X_0, Y_0)$ and radius $\sqrt{h}$ located in the positive quadrant (see Figure 1a). In the $(x, y)$ plane, this corresponds to the ovals depicted in Figure 1b. Only those not intersecting the x-axis (over which $\mu$ vanishes) must be considered. We call them *small* and *big* ovals, respectively:

* When $h \in (0, X_0^2)$, the circle in the $(x^2, y^2)$-plane entirely lies in the positive quadrant. In the $(x, y)$ plane, this results into four symmetric *small* ovals:

$$\Gamma^{++}(h) = H^{-1}(h) \cap \mathbb{R}_{\ge 0} \times \mathbb{R}_{\ge 0}, \qquad \Gamma^{+-}(h) = H^{-1}(h) \cap \mathbb{R}_{\ge 0} \times \mathbb{R}_{\le 0}, \\ \Gamma^{--}(h) = H^{-1}(h) \cap \mathbb{R}_{\le 0} \times \mathbb{R}_{\ge 0}, \qquad \Gamma^{--}(h) = H^{-1}(h) \cap \mathbb{R}_{\le 0} \times \mathbb{R}_{\le 0}.$$

* When $h \in (X_0^2, Y_0^2)$, a portion of this circle in the $(x^2, y^2)$-plane crosses the y-axis, yielding two symmetric *big* ovals:

$$\Gamma^{+}(h) = H^{-1}(h) \cap \mathbb{R} \times \mathbb{R}_{\ge 0}, \qquad \Gamma^{-}(h) = H^{-1}(h) \cap \mathbb{R} \times \mathbb{R}_{\le 0}.$$

In what follows, the notation $\Gamma(h)$ stands for the small oval $\Gamma^{++}(h)$ when $h \in (0, X_0^2)$, and for the big oval $\Gamma^{+}(h)$ when $h \in (X_0^2, Y_0^2)$. The corresponding Abelian integral is:

$$\mathfrak{I}(h) = - \int_{\Gamma(h)} \frac{g(x, y)}{y} dx.$$

By symmetry and the fact that $g(x, y)$ contains only even powers of x and y, it is clear that the Abelian integrals taken along the other ovals of $H^{-1}(h)$ are equal, up to the sign.

Due to the very specific form of the potential function H in (7.3), we can divide the ovals into four arcs $x \mapsto y_{up}(x)$, $x \mapsto y_{down}(x)$, $y \mapsto x_{left}(y)$ and $y \mapsto x_{right}(y)$, 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.

**7.1.2 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 [Johnson 2011] 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 [Johnson 2011] is limited by the binary64 (double precision) format and the finite order quadrature used to compute $\mathfrak{I}(h)$.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Level curves of the potential function H</figcaption></figure>

<table><caption>Table 1. Timings in seconds to rigorously evaluate $\mathfrak{I}(0.25)$ associated to (7.1) with relative error at most $10^{-p}$, using our method with degree $N$ RTAs and floating-point FFT (a) or interval FFT (b) for RTA multiplication (see Remark 6.2), and other software: the CAPD library, the original code of [Johnson 2011] and the Arb library.</caption><thead><tr><th rowspan="2">p</th><th colspan="3">our method</th><th colspan="3">other methods</th></tr><tr><th>N</th><th>(a)</th><th>(b)</th><th>CAPD</th><th>[Johnson 2011]</th><th>Arb/Sage</th></tr></thead><tbody><tr><td>4</td><td>55</td><td>0.32</td><td>1.6</td><td>0.44</td><td>7.9</td><td>0.012</td></tr><tr><td>6</td><td>70</td><td>0.59</td><td>3.3</td><td>0.54</td><td>19.0</td><td>0.013</td></tr><tr><td>8</td><td>100</td><td>0.60</td><td>4.8</td><td>0.67</td><td>49.4</td><td>0.014</td></tr><tr><td>16</td><td>150</td><td>1.5</td><td>9.0</td><td>2.0</td><td>—</td><td>0.018</td></tr><tr><td>24</td><td>225</td><td>1.5</td><td>9.1</td><td>2.2</td><td>—</td><td>0.021</td></tr><tr><td>32</td><td>275</td><td>3.8</td><td>18.7</td><td>3.6</td><td>—</td><td>0.031</td></tr><tr><td>48</td><td>425</td><td>3.8</td><td>18.9</td><td>6.4</td><td>—</td><td>0.040</td></tr><tr><td>64</td><td>570</td><td>8.0</td><td>42.1</td><td>11.1</td><td>—</td><td>0.043</td></tr><tr><td>96</td><td>870</td><td>9.4</td><td>46.1</td><td>23.8</td><td>—</td><td>0.078</td></tr><tr><td>128</td><td>1165</td><td>18.8</td><td>96.2</td><td>41.7</td><td>—</td><td>0.084</td></tr></tbody></table>

### 7.1.3 Computer-assisted proof of $\mathcal{H}(4) \ge 24$
In Table 2, we compute rigorous interval enclosures $I_N$ for $\mathfrak{I}(h)$, for specific values of $h$, using our algorithm with RTAs of degree $N$ as small as possible, as long as the resulting interval guarantees the sign of $\mathfrak{I}(h)$. The existence of sufficiently many simple zeros needed to prove Theorem 7.1 results from the obtained sign alternation. Figures 2 and 3 provide a graphical representation of the sign alternations on small and big ovals.
<table><caption>Table 2. Rigorous evaluation of $Im(h)$ 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 $N$ for RTAs such that the rigorous enclosure $I_N$ guarantees the sign of $Im(h)$. Tighter enclosures $I_{512}$ using a high degree 512 are also provided.</caption><thead><tr><th>ovals</th><th>r</th><th>h</th><th>N</th><th>I<sub>N</sub></th><th>I<sub>512</sub></th><th>sign(Im(h))</th></tr></thead><tbody><tr><td rowspan="6">small</td><td>0.5</td><td>0.25</td><td>17</td><td>[3.5963e-5, 9.6953e-5]</td><td>[6.6457e-5, 6.6458e-5]</td><td>+</td></tr><tr><td>0.78</td><td>0.6084</td><td>37</td><td>[-1.3730e-4, -1.6580e-5]</td><td>[-7.6939e-5, -7.6938e-5]</td><td>-</td></tr><tr><td>0.88</td><td>0.7744</td><td>93</td><td>[1.6407e-9, 3.1821e-8]</td><td>[1.6730e-8, 1.6731e-8]</td><td>+</td></tr><tr><td>0.89</td><td>0.7921</td><td>100</td><td>[-3.8327e-8, -1.1066e-9]</td><td>[-1.9717e-8, -1.9716e-8]</td><td>-</td></tr><tr><td>0.895</td><td>0.801025</td><td>102</td><td>[9.5936e-10, 1.1267e-07]</td><td>[5.6812e-8, 5.6813e-8]</td><td>+</td></tr><tr><td>0.8987</td><td>0.80766169</td><td>111</td><td>[-5.1962e-7, -5.2103e-8]</td><td>[-2.8586e-7, -2.8585e-7]</td><td>-</td></tr><tr><td rowspan="3">big</td><td>0.901</td><td>0.811801</td><td>197</td><td>[-1.2620e-5, -3.4042e-7]</td><td>[-6.4798e-6, -6.4797e-6]</td><td>-</td></tr><tr><td>0.93</td><td>0.8649</td><td>128</td><td>[2.0334e-5, 6.2143e-4]</td><td>[3.2088e-4, 3.2089e-4]</td><td>+</td></tr><tr><td>0.95</td><td>0.9025</td><td>140</td><td>[-1.9667e-4, -7.0461e-6]</td><td>[-1.0186e-4, -1.0185e-4]</td><td>-</td></tr></tbody></table>

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2. Rigorously computed interval enclosures of $Im(h)$ on small ovals, for $h \in \{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.</figcaption></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. Rigorously computed interval enclosures of $Im(h)$ on big ovals, for $h \in \{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.</figcaption></figure>
Proof of Theorem 7.1. Let $h_i$ ($1 \le i \le 9$) denote the 9 values of $h$ in Table 2, taken in increasing order. The rigorous intervals enclosures given in this table are sufficient to guarantee the sign of $\Im(h_i)$. According to Theorem 1.3, there exists for each $h_i$, an $\varepsilon_i > 0$ such that $d(h_i, \varepsilon)$ and $\Im(h_i)$ share the same (strict) sign whenever $0 < \varepsilon \le \varepsilon_i$. Hence, with $\varepsilon^* = \min_{0 \le i \le 9} \varepsilon_i > 0$, we have that the displacement function $h \mapsto d(h, \varepsilon)$ alternates sign at least 5 times on $(0, X_0^2)$ and at least 2 times on $(X_0^2, Y_0^2)$, for each fixed $0 < \varepsilon \le \varepsilon^*$, 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 \times 4 + 2 \times 2 = 24$ limit cycles in the quartic system (7.1) whenever $0 < \varepsilon \le \varepsilon^*$. $\square$

## 7.2 Li, Liu and Yang's cubic system with 13 limit cycles

<figure><img src="image_7.png" /><figcaption>Fig. 4. Ovals of potential function $H$ (7.5), with $\lambda = 0.7$ and $\kappa = 1.1$.</figcaption></figure>

In [Li et al. 2009], C. Li, C. Liu and J. Yang showed that $\mathcal{Z}(3) \ge 13$ using the following perturbed cubic Hamiltonian system:

$$ \begin{cases} \dot{x} = -y^3 + \kappa^2 y, \\ \dot{y} = x^3 + (1 - \lambda)x^2 - \lambda x + \varepsilon y(\alpha_1 + \alpha_2 x + \alpha_3 x^2 + \alpha_4 y^2), \end{cases} \quad (7.4) $$
with well-chosen<fnref n="20" /> 0 < $\lambda$ < 1, $\kappa$ > 1, and coefficients $\alpha_1, \alpha_2, \alpha_3, \alpha_4 \in \mathbb{R}$ for the $\varepsilon$-small perturbation.
The critical points of the corresponding first integral,

$$H(x, y) = \underbrace{\frac{x^4}{4} + \frac{1-\lambda}{3}x^3 - \frac{\lambda}{2}x^2 + \frac{y^4}{4} - \frac{\kappa^2 y^2}{2}}_{F(x)}, \qquad (7.5)$$

depicted in Figure 4, are ranked according to their level value:

$$\begin{array}{l@{\quad}l@{\quad}l@{\quad}l@{\quad}l@{\quad}l}
H(P_3) & < & H(P_1) & < & H(P_2) & < & H(P_6) & < & H(P_4) & < & H(P_5). \\
\parallel & & \parallel & & \parallel & & & & & \\
H(P_9) & & H(P_7) & & H(P_8) & & & & & 0
\end{array}$$

The resulting families of ovals can be considered as a (horizontally) dissymmeterized version of the ovals of [Johnson 2011] discussed in the previous section: $\Gamma_2$ (and symmetrically $\Gamma_5$) when $h \in (H(P_3), H(P_2))$, $\Gamma_1$ (and sym. $\Gamma_4$) when $h \in (H(P_1), H(P_2))$, $\Gamma_3$ (and sym. $\Gamma_6$) when $h \in (H(P_2), H(P_6))$, $\Gamma_8$ when $h \in (H(P_6), H(P_4))$, $\Gamma_9$ when $h \in (H(P_4), H(P_5))$, and $\Gamma_7$ when $h > H(P_4)$.

The Abelian integrals under consideration are:

$$\mathfrak{I}_j(h) = \int_{\Gamma_j(h)} y(\alpha_1 + \alpha_2 x + \alpha_3 x^2 + \alpha_4 y^2) \, dx, \quad \text{for } 1 \le j \le 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 $\mathfrak{I}(3) \ge 13$ with effective values for the coefficients $\lambda, \kappa, \alpha_1, \alpha_2, \alpha_3, \alpha_4$ in System (7.4).

**THEOREM 7.2.** Consider System (7.4) with parameters:

$$\begin{array}{l@{\quad}l@{\quad}l}
\lambda = 0.17, & \alpha_1 = 8.808855593098, & \alpha_3 = -2.599597169555 \cdot 10^{-5}, \\
\kappa = 20, & \alpha_2 = -2.078279433211 \cdot 10^{-5}, & \alpha_4 = -7.340712733831 \cdot 10^{-3}.
\end{array}$$

Then, over their respective interval of definition:

* $\mathfrak{I}_2$ (and symmetrically $\mathfrak{I}_5$) has at least one simple zero;

* $\mathfrak{I}_3$ (and symmetrically $\mathfrak{I}_6$) has at least 5 simple zeros;

* $\mathfrak{I}_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 [Li et al. 2009], based on the careful analysis of limit cycles for Liénard systems by F. Dumortier and C. Li [Dumortier and Li 2003].

### 7.2.1 Timings.
We benchmark our implementation on ovals $\Gamma_3$ (non symmetric, non convex) and $\Gamma_8$ (since it is not even star-shaped). We set $\lambda = \frac{7}{10}$ and $\kappa = \frac{11}{10}$ (the same values were used to plot Figure 4), and $\alpha_1 = \alpha_2 = \alpha_3 = \alpha_4 = 1$. Let $Q_1$ (resp. $Q_2$) denote the topmost intersection of the homoclinic connection of $P_6$ (resp. $P_4$) with the y-axis. Then, for $\Gamma_3$, we consider the periodic orbit starting from $\frac{1}{10}P_2 + \frac{9}{10}Q_1$, and for $\Gamma_8$, we take as initial point $\frac{1}{2}Q_1 + \frac{1}{2}Q_2$ (this corresponds exactly to the ovals $\Gamma_3$ and $\Gamma_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

<footnote n="20">²⁰Note that it is assumed in [Li et al. 2009] that $\kappa > 10$.</footnote>
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 [Johnson and Tucker 2011]. 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><caption>Table 3. Timings in seconds to rigorously evaluate an Abelian integral along ovals from the $\Gamma_3$ and $\Gamma_8$ family with relative error at most $10^{-p}$, using our implementation with degree $N$ RTAs and floating-point FFT (a) or interval FFT (b) for RTA multiplication (see Remark 6.2), and CAPD.</caption><thead><tr><th colspan="5">oval Γ<sub>3</sub></th><th colspan="5">oval Γ<sub>8</sub></th></tr><tr><th>p</th><th>N</th><th>(a)</th><th>(b)</th><th>CAPD</th><th>p</th><th>N</th><th>(a)</th><th>(b)</th><th>CAPD</th></tr></thead><tbody><tr><td>4</td><td>100</td><td>0.45</td><td>2.2</td><td>0.25</td><td>4</td><td>220</td><td>0.90</td><td>4.6</td><td>0.51</td></tr><tr><td>6</td><td>125</td><td>0.46</td><td>2.2</td><td>0.38</td><td>6</td><td>280</td><td>1.7</td><td>11.6</td><td>0.83</td></tr><tr><td>8</td><td>170</td><td>1.0</td><td>5.8</td><td>0.54</td><td>8</td><td>385</td><td>1.8</td><td>11.7</td><td>1.2</td></tr><tr><td>16</td><td>400</td><td>2.1</td><td>11.6</td><td>1.7</td><td>16</td><td>860</td><td>4.4</td><td>23.7</td><td>3.8</td></tr><tr><td>24</td><td>625</td><td>4.3</td><td>23.5</td><td>4.1</td><td>24</td><td>1330</td><td>8.9</td><td>47.5</td><td>9.8</td></tr><tr><td>32</td><td>855</td><td>5.0</td><td>25.5</td><td>8.3</td><td>32</td><td>1805</td><td>10.8</td><td>52.7</td><td>25.3</td></tr><tr><td>48</td><td>1310</td><td>11.0</td><td>52.0</td><td>22.1</td><td>48</td><td>2750</td><td>22.3</td><td>122.6</td><td>64.5</td></tr><tr><td>64</td><td>1765</td><td>11.1</td><td>55.4</td><td>46.2</td><td>64</td><td>3695</td><td>22.8</td><td>136.2</td><td>139.5</td></tr><tr><td>96</td><td>2675</td><td>27.4</td><td>136.3</td><td>158.5</td><td>96</td><td>5585</td><td>59.7</td><td>295.8</td><td>485.3</td></tr><tr><td>128</td><td>3580</td><td>30.0</td><td>137.1</td><td>361.6</td><td>128</td><td>7475</td><td>60.7</td><td>298.8</td><td>1084.0</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Obtaining numerical $\beta_i$ realizing 5 zeros for $\tilde{\mathfrak{F}}_3$ in Dumortier and Li's system [Dumortier and Li 2003].</figcaption></figure>

### 7.2.2 Computer-assisted proof of $\mathfrak{Z}(3) \ge 13$. Before proving Theorem 7.2, we present a heuristic process, directly inspired by [Li et al. 2009], that yields potential relevant parameters. To do so, we focus on $\mathfrak{F}_2$ and $\mathfrak{F}_3$ (once the parameters fixed, the result on $\mathfrak{F}_7$ will be straightforward to establish
<figure><img src="image_9.png" /><figcaption>Fig. 6. Abelian integrals $\mathfrak{I}_2$, $\mathfrak{I}_3$, $\mathfrak{I}_7$ in Li, Liu and Yang's system (7.4) with the parameters of Theorem 7.2, having 1, 5 and 1 zeros, respectively.</figcaption></figure>

<table><caption>Table 4. Rigorous evaluations of $\mathfrak{I}_2(h)$, $\mathfrak{I}_3(h)$ and $\mathfrak{I}_7(h)$ with the parameters of Theorem 7.2 using ABINTVALID. The ovals $\Gamma_2(h)$, $\Gamma_3(h)$ and $\Gamma_7(h)$ are parameterized with a variable $r$ denoting a position on a transversal, with $h = H(-1 + r, \kappa)$, $h = H(0, \kappa + r)$, and $h = H(r, 0)$, respectively. $N$ is the degree used for RTAs, and $\delta$ is the diameter of the computed interval enclosure for $\mathfrak{I}(h)$.</caption><thead><tr><th>oval</th><th>r</th><th>N</th><th>Im(h)</th><th>δ</th><th>sign</th></tr></thead><tbody><tr><td rowspan="2">Γ<sub>2</sub></td><td>0.88</td><td>2000</td><td>[-9.54e-10, -9.53e-10]</td><td>5.45e-22</td><td>-</td></tr><tr><td>0.92</td><td>2000</td><td>[5.55e-10, 5.56e-10]</td><td>3.92e-21</td><td>+</td></tr><tr><td rowspan="6">Γ<sub>3</sub></td><td>0.0002</td><td>8000</td><td>[-5.21e-13, -5.20e-13]</td><td>4.16e-22</td><td>-</td></tr><tr><td>0.001</td><td>8000</td><td>[8.45e-13, 8.46e-13]</td><td>3.24e-26</td><td>+</td></tr><tr><td>0.0015</td><td>4000</td><td>[-2.16e-13, -2.15e-13]</td><td>5.11e-21</td><td>-</td></tr><tr><td>0.0022</td><td>4000</td><td>[3.52e-13, 3.53e-13]</td><td>2.27e-22</td><td>+</td></tr><tr><td>0.004</td><td>4000</td><td>[-6.13e-12, -6.12e-12]</td><td>6.76e-25</td><td>-</td></tr><tr><td>0.0045</td><td>2000</td><td>[2.54e-12, 2.55e-12]</td><td>1.86e-20</td><td>+</td></tr><tr><td rowspan="2">Γ<sub>7</sub></td><td>25</td><td>100</td><td>[2718.12, 2718.13]</td><td>3.63e-6</td><td>+</td></tr><tr><td>30</td><td>100</td><td>[-3419.83, -3419.82]</td><td>6.00e-7</td><td>-</td></tr></tbody></table>

with our approach). In this case, [Li et al. 2009] uses the change of variable $\tilde{y} = \frac{y^2 - \kappa^2}{2}$ to obtain a simpler potential function,

$$ \tilde{H}(x, \tilde{y}) = \tilde{y}^2 + F(x), \qquad (7.6) $$

whose families of ovals $\tilde{\Gamma}_1$, $\tilde{\Gamma}_2$ and $\tilde{\Gamma}_3$ are depicted in Figure 5a. The corresponding Abelian integrals $\tilde{\mathfrak{I}}_j$ for $1 \le j \le 3$ are:

$$ \tilde{\mathfrak{I}}_j(h) = \int_{\tilde{\Gamma}_j(h)} \tilde{y}(\beta_1 + \beta_2x + \beta_3x^2 + \beta_4\tilde{y}^2) dx = \beta_1 \tilde{\mathfrak{I}}_{j1}(h) + \beta_2 \tilde{\mathfrak{I}}_{j2}(h) + \beta_3 \tilde{\mathfrak{I}}_{j3}(h) + \beta_4 \tilde{\mathfrak{I}}_{j4}(h). \quad (7.7) $$

The authors of [Li et al. 2009] establish, as soon as $\kappa$ is large enough, an equivalence between the numbers of zeros of $\mathfrak{I}_2$ and $\tilde{\mathfrak{I}}_2$, resp. $\mathfrak{I}_3$ and $\tilde{\mathfrak{I}}_3$.

Zeros of $\tilde{\mathfrak{I}}_2$ and $\tilde{\mathfrak{I}}_3$. By combining the main existence result of [Dumortier and Li 2003] with some asymptotic analysis, C. Li, C. Liu and J. Yang prove [Li et al. 2009, Lemma 1] that if $\lambda$ is close to a certain critical value $\lambda^*$, then there exist coefficients $\beta_1, \beta_2, \beta_3, \beta_4$ such that $\tilde{\mathfrak{I}}_2$ has at least one simple zero and $\tilde{\mathfrak{I}}_3$ has at least 5 simple zeros. The reason why 5 zeros can be obtained from the 4-term linear combination $\tilde{\mathfrak{I}}_3$ is that for this critical value $\lambda^*$, the Wronskian $\tilde{W}_3(h) = \det(\tilde{\mathfrak{I}}_{3j}^{(i-1)}(h))_{1 \le i, j \le 4}$ has a
double zero $\tilde{h}_0 = \tilde{r}_0^2$ (see Figure 5b). Then, an appropriate perturbation of $\lambda$ around $\lambda^*$ creates two simple zeros for $\tilde{W}(h)$, from which one hopes to deduce the existence of values for the $\beta_i$ realizing 5 zeros$^{21}$.

In contrast with the non-constructive existence proof in [Li et al. 2009], we confirm this heuristic in a constructive way, by determining candidate values for $\lambda$ and the $\beta_i$ purely numerically. To do so, we compute Chebyshev approximations of the Wronskian $r \mapsto \tilde{W}_3(r^2)^{22}$ for several values of $\lambda$ 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 value$^{23}$. The value $\lambda = 0.17$ is found to fulfill those requirements, that is 2 zeros $\tilde{h}_1 = \tilde{r}_1^2 < \tilde{r}_2^2 = \tilde{h}_2$ for $\tilde{W}_3$ (at least numerically, see Figure 5b), with $\tilde{r}_1$ and $\tilde{r}_2$ not too close for better numerical stability.

Now, to obtain values for the $\beta_i$, we fix $\beta_1 = 1$, and numerically solve for $\beta_2, \beta_3$ and $\beta_4$ the equations $\tilde{\mathfrak{I}}_3(\tilde{r}_1^2) = 0$ for $\tilde{r}_1, \tilde{r}_2$ and the midpoint $\tilde{r}_3 = \frac{\tilde{r}_1+\tilde{r}_2}{2}$. As desired, the numerical plot of the resulting combination $r \mapsto \tilde{\mathfrak{I}}_3(r)$ in Figure 5b exhibits 2 additional zeros $\tilde{r}_4$ and $\tilde{r}_5$, making a total of 5 zeros for $\tilde{\mathfrak{I}}_3$. Finally, a numerical plot suggests that $\tilde{\mathfrak{I}}_2$ also has one simple zero.

### Zeros of $\tilde{\mathfrak{I}}_2$, $\tilde{\mathfrak{I}}_3$ and $\tilde{\mathfrak{I}}_7$.
As claimed by [Li et al. 2009, Lemma 2], by selecting a sufficiently large parameter $\kappa > 1$ in (7.4) and defining the coefficients $\alpha_i$ from the $\beta_i$ using the following formulas:

$$ \alpha_1 = \frac{\beta_1}{\kappa^5} + \frac{3\beta_4}{2\kappa}, \quad \alpha_2 = \frac{\beta_2}{\kappa^5}, \quad \alpha_3 = \frac{\beta_3}{\kappa^5}, \quad \alpha_4 = -\frac{\beta_4}{2\kappa^3}, $$

the 5 zeros of $\tilde{\mathfrak{I}}_3$ and the single zero of $\tilde{\mathfrak{I}}_2$ are recovered in $\tilde{\mathfrak{I}}_3$ and $\tilde{\mathfrak{I}}_2$. Moreover, with these coefficients, an additional zero can be found for $\tilde{\mathfrak{I}}_7$, according to [Li et al. 2009, Lemma 3].

Fixing $\kappa = 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 $\alpha_i$, 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 $\mathfrak{Z}(3) \ge 13$.

### Proof of Theorem 7.2.
We perform rigorous evaluations of Abelian integrals $\tilde{\mathfrak{I}}_2$, $\tilde{\mathfrak{I}}_3$ and $\tilde{\mathfrak{I}}_7$ on the points given in Table 4 using Algorithm ABINTVALID, with the values of $\lambda, \kappa$ and the $\alpha_i$ given in the statement of the theorem. By continuity of these three functions, it follows that $\tilde{\mathfrak{I}}_2$, $\tilde{\mathfrak{I}}_3$ and $\tilde{\mathfrak{I}}_7$ have at least 1, 5 and 1 simple zeros, respectively. Moreover, by the symmetry $(x, y) \mapsto (x, -y)$, $\tilde{\mathfrak{I}}_5$ and $\tilde{\mathfrak{I}}_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 $\epsilon > 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 $\tilde{\mathfrak{I}}_2$ and $\tilde{\mathfrak{I}}_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 $(x_r, y_r)$. Alternative methods for this should be investigated in the future to avoid this phenomenon.

<footnote>$^{21}$Indeed, the case of a single (simple) zero for the Wronskian (together with additional non-vanishing conditions on lower-dimensional Wronskians of the system) would imply that $\tilde{\mathfrak{I}}_3$ cannot exhibit more than 4 zeros (see [Novaes and Torregrosa 2017, Cor. 1.4]).</footnote>

<footnote>$^{22}$For convenience, we use $r = \sqrt{h}$ instead of $h$ as the independent variable.</footnote>

<footnote>$^{23}$For 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.</footnote>
## 7.3 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 [Brisebarre et al. 2020], 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(n)$ and $Z(n)$ bounds.

## 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.

## REFERENCES

S. S. Abhyankar and C. L. Bajaj. 1988. Computations with Algebraic Curves. In *Symbolic and Algebraic Computation, International Symposium ISSAC'88, Rome, Italy, July 4-8, 1988, Proceedings (Lecture Notes in Computer Science, Vol. 358)*, Patrizia M. Gianni (Ed.). Springer, 274-284. https://doi.org/10.1007/3-540-51084-2_26

A. Albouy and V. Kaloshin. 2012. Finiteness of central configurations of five bodies in the plane. *Ann. of Math.* (2) 176, 1 (2012), 535-588. https://doi.org/10.4007/annals.2012.176.1.10

V. I. Arnold. 1977. Loss of stability of self-induced oscillations near resonance, and versal deformations of equivariant vector fields. *Funkcional. Anal. i Priložen.* 11, 2 (1977), 1-10, 95. Translated in English in *Functional Anal. Appl.* **11** (1977), no. 2, 85-92.

V. I. Arnold. 1990. *Ten problems*. Adv. Soviet Math., Vol. 1. Amer. Math. Soc., Providence, RI. 1-8 pages.

C. L. Bajaj and G. L. Xu. 1997. Piecewise rational approximations of real algebraic curves. *J. Comput. Math.* 15, 1 (1997), 55-71.

C. Beltrán and A. Leykin. 2013. Robust certified numerical homotopy tracking. *Found. Comput. Math.* 13, 2 (2013), 253-295. https://doi.org/10.1007/s10208-013-9143-2

Y. Bertot and P. Castéran. 2013. *Interactive Theorem Proving and Program Development: Coq'Art: The Calculus of Inductive Constructions*. Springer Science & Business Media.

J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah. 2017. Julia: a fresh approach to numerical computing. *SIAM Rev.* 59, 1 (2017), 65-98. https://doi.org/10.1137/141000671

G. Binyamini, D. Novikov, and S. Yakovenko. 2010. On the number of zeros of Abelian integrals. *Invent. Math.* 181, 2 (2010), 227-289. https://doi.org/10.1007/s00222-010-0244-0

J.-B. Bost and J.-F. Mestre. 1988. Moyenne arithmético-géométrique et périodes des courbes de genre 1 et 2. *Gaz. Math.* 38 (1988), 36-64.

N. Brisebarre and M. Joldeş. 2010. Chebyshev interpolation polynomial-based tools for rigorous computing. In *ISSAC 2010—Proceedings of the 2010 International Symposium on Symbolic and Algebraic Computation*. ACM, New York, 147-154. https://doi.org/10.1145/1837934.1837966

N. Brisebarre, M. Joldeş, J.-M. Muller, A.-M. Naneş, and J. Picot. 2020. Error Analysis of Some Operations Involved in the Cooley-Tukey Fast Fourier Transform. *ACM Trans. Math. Softw.* 46, 2, Article 11 (may 2020), 27 pages. https://doi.org/10.1145/3368619

F. Bréhard. 2011. *Certified numerics in function spaces: Polynomial approximations meet computer algebra and formal proof*. Ph. D. Dissertation. École normale supérieure de Lyon – Université de Lyon, Lyon, France. https://tel.archives-ouvertes.fr/tel-02337901

H. Cartan. 1995. *Elementary theory of analytic functions of one or several complex variables*. Dover Publications, Inc., New York. 228 pages. Translated from the French, Reprint of the 1973 edition.

F. Chen, C. Li, J. Llibre, and Z. Zhang. 2006. A unified proof on the weak Hilbert 16th problem for $n = 2$. *Journal of Differential Equations* 221, 2 (2006), 309-342. https://doi.org/10.1016/j.jde.2005.01.009

C. Christopher and C. Li. 2007. *Limit cycles of differential equations*. Birkhäuser Verlag, Basel. viii+171 pages.

J. Cyranka. 2014. Efficient and generic algorithm for rigorous integration forward in time of dPDEs: Part I. *J. Sci. Comput.* 59, 1 (2014), 28-52. https://doi.org/10.1007/s10915-013-9749-1

F. Dumortier and C. Li. 2003. Perturbation from an elliptic Hamiltonian of degree four. IV. Figure eight-loop. *J. Differential Equations* 188, 2 (2003), 512-554. https://doi.org/10.1016/S0022-0396(02)00111-0
H. Ehlich and K. Zeller. 1966. Auswertung der Normen von Interpolationsoperatoren. *Math. Ann.* 164 (1966), 105–112.
https://doi.org/10.1007/BF01429047

C. Epstein, W. L. Miranker, and T. J. Rivlin. 1982a. Ultra-arithmetic. I. Function data types. *Math. Comput. Simulation* 24, 1 (1982), 1–18. https://doi.org/10.1016/0378-4754(82)90045-3

C. Epstein, W. L. Miranker, and T. J. Rivlin. 1982b. Ultra-arithmetic. II. Intervals of polynomials. *Math. Comput. Simulation* 24, 1 (1982), 19–29. https://doi.org/10.1016/0378-4754(82)90046-5

J.-L. Figueras and R. de la Llave. 2017. Numerical computations and computer assisted proofs of periodic orbits of the Kuramoto-Sivashinsky equation. *SIAM J. Appl. Dyn. Syst.* 16, 2 (2017), 834–852.

X.-S. Gao and M. Li. 2004. Rational quadratic approximation to real algebraic curves. *Comput. Aided Geom. Design* 21, 8 (2004), 805–828. https://doi.org/10.1016/j.cagd.2004.07.009

D. Harvey and J. van der Hoeven. 2021. Integer multiplication in time $O(n\log n)$. *Annals of Mathematics* 193, 2 (2021), 563–617. https://doi.org/10.4007/annals.2021.193.2.4

D. Hilbert. 1900. Mathematische Probleme. Vortrag, gehalten auf dem internationalen Mathematiker-Congress zu Paris 1900. *Nachr. Ges. Wiss. Göttingen, Math.-Phys. Kl.* 1900 (1900), 253–297.

H. Hong and J. Schicho. 1998. Algorithms for trigonometric curves (simplification, implicitization, parameterization). *J. Symbolic Comput.* 26, 3 (1998), 279–300. https://doi.org/10.1006/jsco.1998.0212

A. Hungria, J.-P. Lessard, and J. D. Mireles James. 2016. Rigorous numerics for analytic solutions of differential equations: the radii polynomial approach. *Math. Comp.* 85, 299 (2016), 1427–1459.

F. Johansson. 2017. Arb: efficient arbitrary-precision midpoint-radius interval arithmetic. *IEEE Trans. Comput.* 66, 8 (2017), 1281–1292. https://doi.org/10.1109/TC.2017.2690633

F. Johansson. 2018. Numerical integration in arbitrary-precision ball arithmetic. In *Mathematical Software – ICMS 2018. Springer Lecture Notes in Computer Science*, 255–263. https://doi.org/10.1007/978-3-319-96418-8_30

T. Johnson. 2011. A quartic system with twenty-six limit cycles. *Exp. Math.* 20, 3 (2011), 323–328.

T. Johnson and W. Tucker. 2010. An improved lower bound on the number of limit cycles bifurcating from a quintic Hamiltonian planar vector field under quintic perturbation. *Internat. J. Bifur. Chaos Appl. Sci. Engrg.* 20, 1 (2010), 63–70.
https://doi.org/10.1142/S0218127410025405

T. Johnson and W. Tucker. 2011. On a computer-aided approach to the computation of Abelian integrals. *BIT* 51, 3 (2011), 653–667. https://doi.org/10.1007/s10543-011-0318-4

T. Kapela, M. Mrozek, D. Wilczak, and P. Zgliczyński. 2020. CAPD::DynSys: a flexible C++ toolbox for rigorous numerical analysis of dynamical systems. *Commun. Nonlinear Sci. Numer. Simul.* (2020), 105578.

T. Kapela, D. Wilczak, and P. Zgliczyński. 2021. Recent advances in rigorous computation of Poincaré maps. *arXiv:2104.08046 [math.NA]*

P. Lairez, M. Mezzarobba, and M. Safey El Din. 2019. Computing the Volume of Compact Semi-Algebraic Sets. In *Proceedings of the 2019 on International Symposium on Symbolic and Algebraic Computation*. ACM, New York, NY, USA, 259–266. https://doi.org/10.1145/3326229.3326262

O. E. Lanford, III. 1982. A computer-assisted proof of the Feigenbaum conjectures. *Bull. Amer. Math. Soc.* (N.S.) 6, 3 (1982), 427–434. https://doi.org/10.1090/S0273-0979-1982-15008-X

J. M. Lee. 2013. *Introduction to smooth manifolds* (second ed.). Graduate Texts in Mathematics, Vol. 218. Springer, New York. xvi+708 pages.

J.-P. Lessard. 2018. Computing discrete convolutions with verified accuracy via Banach algebras and the FFT. *Appl. Math.* 63, 3 (2018), 219–235. https://doi.org/10.21136/AM.2018.0082-18

C. Li, C. Liu, and J. Yang. 2009. A cubic system with thirteen limit cycles. *J. Differential Equations* 246, 9 (2009), 3609–3619. https://doi.org/10.1016/j.jde.2009.01.038

G. Liu and V. Kreinovich. 2010. Fast convolution and fast Fourier transform under interval and fuzzy uncertainty. *J. Comput. System Sci.* 76, 1 (2010), 63–76. https://doi.org/10.1016/j.jcss.2009.05.006

K. Makino and M. Berz. 2003a. Taylor models and other validated functional inclusion methods. *Int. J. Pure Appl. Math.* 4, 4 (2003), 379–456.

K. Makino and M. Berz. 2003b. Taylor models and other validated functional inclusion methods. *Int. J. Pure Appl. Math.* 6, 3 (2003), 239–316.

B. Martin, A. Goldsztejn, L. Granvilliers, and C. Jermann. 2013. Certified parallelotope continuation for one-manifolds. *SIAM J. Numer. Anal.* 51, 6 (2013), 3373–3401. https://doi.org/10.1137/130906544

R. E. Moore. 1966. *Interval analysis*. Vol. 4. Prentice-Hall Englewood Cliffs.

A. Neumaier. 2003. Taylor forms—use and limits. *Reliab. Comput.* 9, 1 (2003), 43–79.

D. D. Novaes and J. Torregrosa. 2017. On extended Chebyshev systems with positive accuracy. *J. Math. Anal. Appl.* 448, 1 (2017), 171–186. https://doi.org/10.1016/j.jmaa.2016.10.076

M. J. D. Powell. 1981. *Approximation theory and methods*. Cambridge University Press.
R. Prohens and J. Torregrosa. 2019. New lower bounds for the Hilbert numbers using reversible centers. *Nonlinearity* 32, 1 (2019), 331–355. https://doi.org/10.1088/1361-6544/aae94d

C. Tsitouras. 2011. Runge-Kutta pairs of order 5(4) satisfying only the first column simplifying assumption. *Comput. Math. Appl.* 62, 2 (2011), 770–775. https://doi.org/10.1016/j.camwa.2011.06.002

W. Tucker. 2011. *Validated numerics: a short introduction to rigorous computations*. Princeton University Press.

D. Wilczak. 2003. Chaos in the Kuramoto-Sivashinsky equations—a computer-assisted proof. *J. Differential Equations* 194, 2 (2003), 433–459. https://doi.org/10.1016/S0022-0396(03)00104-9

G. B. Wright, M. Javed, H. Montanelli, and L. N. Trefethen. 2015. Extension of Chebfun to periodic functions. *SIAM J. Sci. Comput.* 37, 5 (2015), C554–C573. https://doi.org/10.1137/141001007

J. Xu, M. Burr, and C. Yap. 2018. An Approach for Certifying Homotopy Continuation Paths: Univariate Case. In *Proceedings of the 2018 ACM International Symposium on Symbolic and Algebraic Computation* (New York, NY, USA) (ISSAC '18). Association for Computing Machinery, New York, NY, USA, 399–406.

N. Yamamoto. 1998. A numerical verification method for solutions of boundary value problems with local uniqueness by Banach's fixed-point theorem. *SIAM J. Numer. Anal.* 35, 5 (1998), 2004–2013.

H. Yang, B. Jüttler, and L. Gonzalez-Vega. 2010. An evolution-based approach for approximate parameterization of implicitly defined curves by polynomial parametric spline curves. *Math. Comput. Sci.* 4, 4 (2010), 463–479.

A. Zygmund. 2002. *Trigonometric series. Vol. I, II* (third ed.). Cambridge University Press, Cambridge. xii; Vol. I: xiv+383 pp.; Vol. II: viii+364 pages.

## A ALGORITHM ISPOSITIVE: CHECKING POSITIVITY OF RTAS

In rigorous numerics, it is common to determine if a function $f$ 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 = (\tilde{f}, \epsilon)$ is preconditioned by an approximate inverse $\tilde{g} \approx \frac{1}{f}$. If $f$ does not vanish, then the product $\tilde{g}f$ is likely to be close to 1, which is measured by the quantity $\|1 - \tilde{g}f\|_\ell^\epsilon$.

**Algorithm 7 ISPOSITIVE(f, N)**

**Input:** RTA $f = (\tilde{f}, \epsilon)$ and approximation degree $N$

**Output:** Boolean true or false, with true guaranteeing that $f > 0$ for all $f \in f$

1: $a_0 \leftarrow$ constant coefficient of $\tilde{f}$

2: **if** $a_0 \le 0$ **then**

3: **return** false

4: **else if** BOUND($f - a_0$) < $a_0$ **then**

5: **return** true

6: **else**

7: $\tilde{g} \leftarrow$ FFT$_N((\frac{1}{f_j})_{j=0}^{2N})$, where $(f_j)_{j=0}^{2N} \leftarrow$ IFFT$_N(\tilde{f})$

8: **if** $\tilde{g}$ is well defined and BOUND$(1 - \tilde{g}f) < 1$ **then**

9: **return** true

10: **else**

11: **return** false

12: **end if**

13: **end if**

**LEMMA A.1.** Let $f = (\tilde{f}, \epsilon)$ be an RTA. Then,

(i) For any approximation degree $N$, if `ISPOSITIVE(f, N)` returns true, then $f > 0$ for all $f \in C_{2\pi}^0$ represented by $f$.
(ii) If moreover $\tilde{f} > 0$ and $\varepsilon < \|\tilde{f}^{-1}\|_{\ell^1}^{-1}$ (which implies $f > 0$ for all $f \in f$), then there is an $N_0$ such that $\text{IsPOSITIVE}(f, N)$ returns true for all $N \ge N_0$.

(iii) $\text{IsPOSITIVE}(f, N)$ runs in $O(N' \log(N'))$ arithmetic operations, where $N' = \max(\deg f, N)$.

PROOF. (i) Soundness. Assume $\text{IsPOSITIVE}(f, N)$ returns true. A first possibility is that the condition in line 4 is true, implying that $a_0 > 0$ and that any function $f \in f$ can be written as $a_0 + h$ with $\|h\|_{\infty} < a_0$, and is therefore positive. The other possibility (line 8) is that we obtained a $\tilde{g}$ such that for any $f \in f$, $\|1 - \tilde{g}f\|_{\infty} < 1$, implying that $f(t) \ne 0$ for all $t$. By continuity, $f$ has constant sign over $\mathbb{R}$. In particular, this is true for $\tilde{f}$. Since the condition in line 2 was false, i.e., $a_0 = \int_0^{2\pi} \tilde{f} dt > 0$, we get $\tilde{f} > 0$. Finally, since the ball in $C_{2\pi}^0$ denoted by $f$ is convex, the same statement $f > 0$ is true for all $f \in f$.

(ii) Completeness. Suppose conversely that $\tilde{f} > 0$ and that $\varepsilon < \|\tilde{f}^{-1}\|_{\ell^1}^{-1}$. In particular, $\varepsilon < \|\tilde{f}^{-1}\|_{\infty}^{-1} = \min_{t \in \mathbb{R}} \tilde{f}(t)$, implying that $f > 0$ for all $f \in f = (\tilde{f}, \varepsilon)$. We prove that $\text{IsPOSITIVE}(f, N)$ returns true for $N$ sufficiently large.

First, the condition in line 2 is false since $a_0 = \int_0^{2\pi} \tilde{f} dt > 0$. Then, the algorithm either returns true if the condition in line 4 is true, or it computes $\tilde{g} \leftarrow \mathcal{I}_N(\tilde{f}^{-1})$ in line 7. Since $1/\tilde{f}$ is analytic, $\|\tilde{g} - \tilde{f}^{-1}\|_{\ell^1}$ converges to 0 as the approximation degree $N$ tends to infinity. Hence,

$$ \text{BOUND}(1 - \tilde{g}f) = \|1 - \tilde{g}\tilde{f}\|_{\ell^1} + \varepsilon \text{BOUND}(\tilde{g}) \to 0 + \varepsilon \|\tilde{f}^{-1}\|_{\ell^1} < 1 \quad \text{as } N \to \infty. $$

Therefore, there is an $N_0$ such that for all $N \ge N_0$, the degree $N$ TP $\tilde{g}$ computed by trigonometric interpolation makes the condition in line 8 true, so that $\text{IsPOSITIVE}(f, N)$ returns true.

(iii) Complexity. The asymptotic complexity is determined by the $FFT_N/IFFT_N$ routines used for evaluation on the equispaced grid, interpolation and TP/RTA multiplication. $\square$

## 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 $\hat{s} : U \to \mathbb{R}$ over a neighborhood $U \subseteq \mathbb{R}^4$ of

$$ \hat{\Gamma}(h) = \left\{ (x, y, u, v) \in \mathbb{R}^4 \text{ s.t. } (x, y) \in \Gamma(h) \text{ and } (u, v) = \frac{\nabla H(x, y)}{\|\nabla H(x, y)\|_2} \right\}, $$

that serves as a shift to project any point $(x, y)$ onto $\Gamma(h)$ in the direction given by $(u, v)$. Then, we prove that over a possibly smaller neighborhood $U' \subseteq U$ of $\hat{\Gamma}(h)$, the Newton iterations $\mathcal{N}_{x,y,u,v}^k(0)$ (see Eq. (3.5)) converge quadratically fast to $\hat{s}(x, y, u, v)$ for all $(x, y, u, v) \in U'$. Finally, we compose by $t \in [0, 2\pi] \to (x^\circ(t), y^\circ(t))$ to complete the proof.

* Analytic projection. Consider the equation:

$$ F(x, y, u, v, s) = H(x + su, y + sv) - h = 0. $$

By definition, $F(x, y, u, v, 0) = 0$ for all $(x, y, u, v) \in \hat{\Gamma}(h)$. Moreover, the derivative with respect to $s$,

$$ \frac{\partial F}{\partial s}(x, y, u, v, s) = u \frac{\partial H}{\partial x}(x + su, y + sv) + v \frac{\partial H}{\partial y}(x + su, y + sv) $$

satisfies $\frac{\partial F}{\partial s}(x, y, u, v, 0) = \|\nabla H(x, y)\|_2 > 0$ for all $(x, y, u, v) \in \hat{\Gamma}(h)$, since $h$ is a regular value of $H$. By the analytic implicit function theorem [Cartan 1995, Chap. IV, Prop. 6.1], there exists a neighborhood $U \subseteq \mathbb{R}^4$ of $\hat{\Gamma}(h)$, a neighborhood $V \subseteq \mathbb{R}^5$ of $\hat{\Gamma}(h) \times \{0\}$, and a unique analytic function $\hat{s} : U \to \mathbb{R}$ such that:

$$ ((x, y, u, v, s) \in V \wedge F(x, y, u, v, s) = 0) \Leftrightarrow ((x, y, u, v) \in U \wedge s = \hat{s}(x, y, u, v)). $$

In particular, $\hat{s}(x, y, u, v) = 0$ if and only if $(x, y) \in \Gamma(h)$.
• *Convergence of Newton iterations.* Since $\hat{\Gamma}(h)$ is compact, we may assume (up to restricting the neighborhoods $U$ and $V$) that $\frac{\partial F}{\partial s}(x, y, u, v, s) \ge m_1$ for some $m_1 > 0$, and $|\frac{\partial^2 F}{\partial s^2}(x, y, u, v, s)| \le M_2$ for some $M_2 > 0$, for all $(x, y, u, v, s) \in V$. Now for any $(x, y, u, v, s) \in V$, the classical error analysis of Newton’s method [Tucker 2011, §5.1.2] yields in our case:

$$|N_{x,y,u,v}(s) - \hat{s}(x, y, u, v)| \le \frac{M_2}{2m_1} |s - \hat{s}(x, y, u, v)|^2.$$

Therefore, by restricting the neighborhood $U$ of $\hat{\Gamma}(h) \subseteq \hat{s}^{-1}(0)$ to:

$$U' = \{(x, y, u, v) \in U, |\hat{s}(x, y, u, v)| \le \frac{m_1}{M_2}\},$$

we deduce that the Newton iterations $N_{x,y,u,v}^k(0)$ converge quadratically fast to the desired limit $\hat{s}(x, y, u, v)$:

$$|N_{x,y,u,v}^k(0) - \hat{s}(x, y, u, v)| \le \left(\frac{M_2}{2m_1}\right)^{2k-1} |\hat{s}(x, y, u, v)| \le \frac{|\hat{s}(x, y, u, v)|}{2^{2k-1}}. \quad (B.1)$$

• *Composition with the initial guess.* Let us fix $\eta > 0$ such that any $(x^\circ, y^\circ)$ and $(u^\circ, v^\circ)$ $\eta$-close to $(x^*, y^*)$ and $\frac{\nabla H(x^*, y^*)}{\|\nabla H(x^*, y^*)\|_2}$ in $(C_{2\pi}^0)^2$ satisfy $(x^\circ(t), y^\circ(t), u^\circ(t), v^\circ(t)) \in U'$ for all $t \in [0, 2\pi]$. By the property of composition of analytic functions, the map $t \mapsto s^\#(t) = \hat{s}(x^\circ(t), y^\circ(t), u^\circ(t), v^\circ(t))$ is analytic, and so is $(x^\#(t), y^\#(t))$ by Eq. (3.4). Also, by choosing $\eta > 0$ sufficiently small, we can make $(x^\#(t), y^\#(t))$ sufficiently close to $(x^*, y^*)$, which ensures that the winding number with respect to a given point inside $\Gamma(h)$ remains the same ($\pm 1$ depending on the orientation). This concludes the proof that $(x^\#(t), y^\#(t))$ is an analytic parameterization of $\Gamma(h)$ with same orientation as $(x^*, y^*)$, although not equal to $(x^*, y^*)$ in general. Finally, the quadratic convergence of the Newton iteration scheme for any point $(x^\circ(t), y^\circ(t))$ w.r.t. direction $(u^\circ(t), v^\circ(t))$ follows from the preceding analysis. $\square$

### 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 $r$ and a Lipschitz constant $\lambda$ satisfying the hypotheses of Theorem 2.5. To this aim, given a function $r \mapsto \lambda(r)$ providing a guaranteed upper bound for the Lipschitz constant for the fixed-point operator $\mathcal{T}$ over $\bar{B}(\tilde{\phi}, r)$, and an upper bound for the defect $d \ge \|\mathcal{T}(\tilde{\phi}) - \tilde{\phi}\|$., Algorithm NewtonBall implements a simple bisection method to identify the smallest zero of $r \mapsto f(r) := (1 - \lambda(r))r - d$. Specifically, it maintains a stack of subintervals $[r^-, r^+]$ to be investigated, each of them satisfying $f(r^-) \le 0$. Under the reasonable assumption that $\lambda(r)$ is an increasing function of $r$, the range of $f$ over $[r^-, r^+]$ can be overapproximated by $[(1 - \lambda^+)r^- - d, (1 - \lambda^-)r^+ - d]$ where $\lambda^\pm := \lambda(r^\pm)$, thus giving a simple test (line 7) to detect the possible existence of a zero of $f$ in $[r^-, r^+]$. If moreover $f(r^+) \ge 0$, then $[r^-, r^+]$ necessarily contains a zero, and the stack of remaining intervals is cleared (line 8). When the process successfully terminates, it returns a pair $(r^+, \lambda^+)$ 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 $d + \lambda^+ r^+ < r^+$. By doing so, we make sure that this inequality remains valid under small perturbations of $\lambda(r)$, 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.
**Algorithm 8 NEWTONBALL(d, λ)**

**Input:** defect $d \ge \|\mathcal{T}(\dot{\varphi}) - \ddot{\varphi}\|$ and procedure $\lambda$ s.t. $\mathcal{T}$ is $\lambda(r)$-Lipschitz over $\bar{B}(\ddot{\varphi}, r)$

**Output:** $(r^+, \lambda^+)$ s.t. $\mathcal{T}$ is $\lambda^+$-contracting over $\bar{B}(\ddot{\varphi}, r^+)$ and $d + \lambda^+r^+ < r^+$

1: compute $r_{max}$ s.t. $\lambda(r_{max}) \approx 0.5$

2: push $((0, \lambda(0)), (r_{max}, \lambda(r_{max})))$ on a new empty stack $S$

▷ Bisection loop

3: **while** $S$ not empty **do**

4:   $((r^{-}, \lambda^{-}), (r^{+}, \lambda^{+})) \leftarrow \text{pop}(S)$

▷ When interval subdivision is below user-defined $\delta_{\lambda}$

5:   **if** $\lambda^{+} - \lambda^{-} \le \delta_{\lambda}$ **and** $d + \lambda^{+}r^{+} < r^{+}$ **then**

6:   **return** $(r^{+}, \lambda^{+})$

▷ Check whether interval $[r^{-}, r^{+}]$ may contain a zero

7:   **else if** $\lambda^{+} - \lambda^{-} > \delta_{\lambda}$ **and** $d + \lambda^{-}r^{+} < r^{+}$ **then**

8:   **if** $d + \lambda^{+}r^{+} \le r^{+}$ **then** clear(S) **end if**

▷ Subdivide $[r^{-}, r^{+}]$ and push resulting intervals on the stack

9:   $r_{m} \leftarrow (r^{-} + r^{+})/2$

10:   $\lambda_{m} \leftarrow \lambda(r_{m})$

11:   push(((r_m, \lambda_m), (r^+, \lambda^+)), S)

12:   push(((r^{-}, \lambda^{-}), (r_m, \lambda_m)), S)

13:   **end if**

14:   **end while**

▷ Empty stack

15:   **return** FAIL