Ilya Loshchilov
email: ilya.loshchilov@epfl.ch

Marc Schoenauer
email: marc.schoenauer@inria.fr

Michèle Sebag
email: michele.sebag@inria.fr

Intensive Surrogate Model Exploitation in Self-adaptive Surrogate-assisted CMA-ES (saACM-ES)

Keywords: I.2.8 [Computing Methodologies]: Artificial Intelligence Problem Solving, Control Methods, and Search Evolution Strategies, CMA-ES, self-adaptation, surrogateassisted optimization, surrogate models, ranking support vector machine, black-box optimization

This paper presents a new mechanism for a better exploitation of surrogate models in the framework of Evolution Strategies (ESs). This mechanism is instantiated here on the self-adaptive surrogate-assisted Covariance Matrix Adaptation Evolution Strategy ( s * ACM-ES), a recently proposed surrogate-assisted variant of CMA-ES. As well as in the original s * ACM-ES, the expensive function is optimized by exploiting the surrogate model, whose hyper-parameters are also optimized online. The main novelty concerns a more intensive exploitation of the surrogate model by using much larger population sizes for its optimization.

The new variant of s * ACM-ES significantly improves the original s * ACM-ES and further increases the speed-up compared to the CMA-ES, especially on unimodal functions (e.g., on 20-dimensional Rotated Ellipsoid, s * ACM-ES is 6 times faster than aCMA-ES and almost by one order of magnitude faster than CMA-ES). The empirical validation on the BBOB-2013 noiseless testbed demonstrates the efficiency and the robustness of the proposed mechanism.

INTRODUCTION

Evolutionary Algorithms (EAs) have received a lot of attention regarding their potential to solve black-box optimiza-tion problems, where typically no additional information is available apart from the quality of evaluated solutions. By closely looking at the most successful EAs w.r.t. results on both artificial benchmark and real-world optimization problems, one may observe that invariance properties play a crucial role and represent a source of robustness. A good illustrative example is Covariance Matrix Adaptation Evolution Strategy (CMA-ES [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF]) which exhibits invariance w.r.t. rank-preserving transformations of the objective function and invariance w.r.t. orthogonal transformations of the search space (if the initial search point(s) are transformed accordingly). The algorithm has demonstrated the stateof-the-art performance on benchmark problems (see, e.g., CEC-2005 [START_REF] García | A study on the use of non-parametric tests for analyzing the evolutionary algorithms' behaviour: a case study on the CEC'2005 Special Session on Real Parameter Optimization[END_REF], BBOB-2010 [START_REF] Hansen | Comparing results of 31 algorithms from the black-box optimization benchmarking BBOB[END_REF]) and real-world optimization problems [START_REF] Hansen | References to CMA-ES Applications[END_REF].

One of the main limitations of EAs, the large number of function evaluations required for a reasonable accuracy of optimization, prevents EAs from being widely used on optimization problems which are expensive in terms of time (one evaluation takes several seconds or even hours) and/or in terms of money. This limitation is especially observable in the special, but quite common case of the unimodal noiseless continuous optimization, where gradient information is useful and quasi-Newton methods such as BFGS [START_REF] Shanno | Conditioning of Quasi-Newton Methods for Function Minimization[END_REF], proposed 40 years ago, usually outperform most of advanced EAs [START_REF] Hansen | Comparing results of 31 algorithms from the black-box optimization benchmarking BBOB[END_REF]. The latter is also the case for CMA-ES on some interesting (for benchmarking) optimization problems such as Rosenbrock function.

To address this limitation, numerous approaches to integrate surrogate/approximate models into EAs have been proposed (see, e.g., [START_REF] Jin | Surrogate-assisted evolutionary computation: Recent advances and future challenges[END_REF] for a recent overview of the stateof-the-art). In most of these approaches, a surrogate model f (x) (or simply f ) of the objective function f (x) (or simply f ) is built using the information about evaluated (training) solutions and corresponding objective values. Then, f is used to predict the quality of some promising test solutions in order to proceed towards the optimum of f faster in terms of number of function evaluations. The use of this approach usually leads to a loss of invariance w.r.t. rank-preserving transformations of f since the building of f is sensitive to these kinds of transformations. The property of invariance w.r.t. orthogonal transformations of the search space is also not preserved in most of the surrogate-assisted approaches, since the surrogate learning phase does not take into account the information about a coordinate system (if any) adapted during the search.

Until recently, surrogate-assisted CMA-ES algorithms preserved at most one of the above described invariance prop-erties (see, e.g., [START_REF] Runarsson | Ordinal Regression in Evolutionary Computation[END_REF][START_REF] Kern | Local Meta-Models for Optimization Using Evolution Strategies[END_REF]). The first surrogate-assisted CMA-ES which preserves both invariance properties, referred to as ACM-ES, was proposed in [START_REF] Loshchilov | Comparison-Based Optimizers Need Comparison-Based Surrogates[END_REF]. Its extension, referred to as self-adaptive surrogate-assisted CMA-ES ( s * ACM-ES [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF]), has demonstrated an ability to adapt during the search the hyper-parameters used to build the surrogate model, and this can be viewed as a sort of invariance. s * ACM-ES usually speed-ups the original CMA-ES by a factor between 2 and 4 on unimodal functions. However, it was surprisingly outperformed by ACM-ES algorithm on 10-dimensional Rosenbrock function [START_REF] Loshchilov | Comparison-Based Optimizers Need Comparison-Based Surrogates[END_REF]. An analysis of this observation led to an hypothesis that surrogate model exploitation is not sufficiently intensive in s * ACM-ES.

This paper investigates a simple mechanism of a more intensive surrogate model exploitation by using much larger population sizes, while keeping the default population size when optimizing the original expensive function. Such exploitation is useful when the surrogate model is sufficiently accurate and may lead to a certain divergence when the surrogate badly approximates the expensive function. Thus, a trade-off between between two cases should be respected.

The paper is organized as follows. First, Section 2 reviews surrogate-assisted Evolution Strategies and presents s * ACM-ES algorithm with a discussion of its main principles. Next, Section 3 introduces our new mechanism to surrogate model exploitation. The experimental validation of the proposed mechanism is reported and discussed in Section 4.1. Section 5 concludes the paper.

SURROGATE-ASSISTED EVOLUTION STRATEGIES

First, we discuss the various techniques used to learn and exploit surrogate models within Evolution Strategies. Then, we briefly describe s * ACM-ES algorithm whose mechanism of surrogate model exploitation will be further studied and modified in Section 3.

Background

One of the first surrogate-assisted (µ + λ)-ES and (µ, λ)-ES were proposed in [START_REF] Papadrakakis | Structural optimization using evolution strategies and neural networks[END_REF], where evaluations of an expensive structural optimization problem were replaced by a hiddenlayer Artificial Neural Network (ANN) trained by back-propagation. The authors suggested to re-learn the model at each iteration by adding new training points randomly drawn from a Gaussian distribution with the mean located in the center of the decision space.

Jin et al. [START_REF] Jin | Managing Approximate Models in Evolutionary Aerodynamic Design Optimization[END_REF] proposed individual-based and generationbased evolution control strategies for surrogate-assisted search with CMA-ES. Individual-based strategy corresponds to a pre-selection strategy, where λ ′ random (random strategy) or best (best strategy) controlled individuals out of λP re pre-selected individuals are evaluated with the expensive function f . In generation-based strategy, the whole population will be evaluated with the expensive function for η generations every κ generations, where η ≤ κ. The authors found that in both cases about 50% of individuals should be evaluated in order to have a good final speed-up, but they employed only the generation-based strategy, because they found it more suitable for parallel implementation. It is suggested to set the fraction of generations η κ that f is optimized to be proportional to E(k) Emax , where E(k) and Emax are current and maximum model errors, respectively.

Emmerich et al. [START_REF] Emmerich | Metamodel-Assisted Evolution Strategies[END_REF] studied Kriging model-based preselection strategy within (µ, κ, λ)-ES with µ = 15, λ = 100 and κ = 5, where the individuals that exceed the age of κ = 5 generations are eliminated from the selection procedure. The authors suggested to evaluate the best λ ′ = 10 individuals w.r.t. a criterion based on both the estimated value f (x) and the estimated local standard deviation σ(x) of the prediction Kriging model as in [START_REF] Trosset | Numerical optimization using computer experiments[END_REF]:

Sc(x) = f (x) -ρσ(x), (1)
where ρ ≥ 0 defines the selection trade-off between the most promising solutions with ρ = 0 and promising solutions in still unexplored search areas with ρ > 0. The experimental validation of the proposed surrogate-assisted algorithm showed that on unimodal functions the exploration (ρ = 1) does not harm, but leads to a better convergence on multimodal functions than f -based selection (ρ = 0).

Ulmer et al. [START_REF] Ulmer | Model-assisted steady-state evolution strategies[END_REF] studied RBF networks-based pre-selection strategy for (2,8)-ES with λP re = 30 and (1+1)-ES with λP re = 10 and concluded that the proposed meta-model assisted ES (MAES) performs better, especially on unimodal functions. Later, Ulmer et al. [START_REF] Ulmer | Evolution Strategies assisted by Gaussian Processes with Improved Pre-Selection Criterion[END_REF] analyzed GP-based preselection strategy similar to the one of [START_REF] Emmerich | Metamodel-Assisted Evolution Strategies[END_REF] for CMA-ES, where an exploration-based selection criterion was chosen. The results confirmed the observations of [START_REF] Emmerich | Metamodel-Assisted Evolution Strategies[END_REF] that on multimodal functions the exploration-based pre-selection should be preferred to "greedy" f -based pre-selection. Ulmer et al. [START_REF] Ulmer | Evolution strategies with controlled model assistance[END_REF] further studied Support Vector Regression (SVR)based pre-selection strategy for (µ, λ) Main Vector Adaptation (MVA [START_REF] Poland | Main Vector Adaptation: A CMA Variant with Linear Time and Space Complexity[END_REF], a CMA-ES variant with linear time and space complexity). They proposed to adjust λP re, depending on a model quality measure similar to the one proposed by [START_REF] Jin | Quality measures for approximate models in evolutionary computation[END_REF] and defined as summed rank of all correctly selected individuals. At each generation t, the actual model quality Q t is compared to a quality Q rand measured for the random model. The update procedure is controlled by a parameter δ λ P re as follows:

λ t+1 P re = λ t P re + Q max -Q t Q max -Q rand δ λ P re if Q t > Q rand λ t P re -Q rand -Q t Q rand δ λ P re otherwise, (2)
where Q max is the maximum possible quality. Experimental results of the proposed C-MAES on unimodal and multimodal functions confirmed that the adaptation works well even in dynamically changing noisy environment, where the speed-up up to a factor of 2 also can be achieved.

Runarsson [START_REF] Runarsson | Constrained Evolutionary Optimization by Approximate Ranking and Surrogate Models[END_REF] found that it is always desirable to evaluate with f at least one best point among λ evaluated with f , and proposed approximate ranking procedure which suggests to evaluate with f several additional points only if the model changes its ranking prediction of λ points.

Poland [START_REF] Poland | Explicit Local Models: Towards "Optimal" Optimization Algorithms[END_REF] proposed an ES-like algorithm with quadratic meta-models, which showed the best results (in that time) in the literature for an Evolutionary Algorithm on popular Rosenbrock benchmark problem, where it outperformed CMA-ES by a factor of 4. Unfortunately, the algorithm is not invariant to rank-preserving transformations of f and its performance probably will degrade significantly if f is scaled differently (also true for all regression to value-based surrogate-assisted optimizers).

Büche et al. [START_REF] Buche | Accelerating evolutionary algorithms with Gaussian process fitness function models[END_REF] proposed Gaussian Process Optimization Procedure (GPOP) which suggests to repeat the following procedure. First, build a GP-based model and directly optimize Eq. ( 1) by CMA-ES. When a local optimum for a given ρ is found, it can be evaluated with f and added to the training set. The training set consists of the union of the NC closest points to the current best solution and the NR most recently evaluated points. The search space for each local search is constrained by the hyper-rectangle around the current best solution. The GPOP outperforms CMA-ES on Sphere and Schwefel functions, but already for Rosenbrock the speed-up is less than 2.0 for n = 8, 16 and less than 0.4 for n = 32, where larger training sets are probably needed. On multimodal Rastrigin function GPOP in most cases is outperformed by CMA-ES.

Runarsson [START_REF] Runarsson | Ordinal Regression in Evolutionary Computation[END_REF] first proposed to use ordinal regression and Ranking Support Vector Machine (SVM) as surrogate models in Evolutionary Computation. The author studied the performance of Ranking SVM-based CMA-ES within the approximate ranking procedure and found that the surrogateassisted version outperforms the original CMA-ES for n ≤ 5, but shows no improvements for n > 5. The latter can be explained by the use of very small training set of only 60 training points. However, this paper has played an important role in the development of comparison-based surrogates (see, e.g., [START_REF] Loshchilov | Comparison-Based Optimizers Need Comparison-Based Surrogates[END_REF]).

Hoffmann and Holemann [START_REF] Hoffmann | Controlled Model Assisted Evolution Strategy with Adaptive Preselection[END_REF] suggested to keep λP re fixed, but adapt λ using the following formula, similar to (2):

λ t+1 = max(λ t -Q max -Q t Q rand δ λ , µ) if Q t > Q rand min(λ t + Q max -Q t Q max -Q rand δ λ , λP re) otherwise, (3)
Thus, the better the model -the smaller λ t+1 , and vice versa. The results of the proposed λ-controlled surrogateassisted λ-CMA-ES are similar, but not directly comparable to [START_REF] Ulmer | Evolution strategies with controlled model assistance[END_REF], where MVA was used as a baseline algorithm.

The local meta-model CMA-ES based on locally weighted regression was first proposed by [START_REF] Kern | Local Meta-Models for Optimization Using Evolution Strategies[END_REF] and later extended for large populations by [START_REF] Bouzarkouna | Investigating the Local-Meta-Model CMA-ES for Large Population Sizes[END_REF]. In lmm-CMA-ES, the surrogate model exploitation is controlled by the approximate ranking procedure [START_REF] Runarsson | Constrained Evolutionary Optimization by Approximate Ranking and Surrogate Models[END_REF], which suggests to evaluate with f only a fraction of λ new individuals if the ranking of some of λ individuals on surrogate models of the current and previous generations remains unchanged. The algorithm and its extensions have demonstrated a relatively good performance, but the time complexity usually scales between O(n 4 ) and O(n 6 ) depending on the version of the algorithm [START_REF] Kern | Fast Quadratic Local Meta-Models for Evolutionary Optimization of Anguilliform Swimmers[END_REF], limiting the range of application to problems with n 15 [START_REF] Kern | Fast Quadratic Local Meta-Models for Evolutionary Optimization of Anguilliform Swimmers[END_REF].

Self-adaptive Surrogate-assisted CMA-ES

In the following, we outline s * ACM-ES algorithm according to [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF], see Figure 1.

In s * ACM-ES, two optimization procedures are performed:

1. Optimization of the objective function f using CMA-ES (so-called "CMA-ES #1" in Figure 1) assisted by a surrogate model f α, built using a vector α of some surrogate model (learning) hyper-parameters.

2. Optimization of the surrogate model error Err(α) (using so-called "CMA-ES #2" in Figure 1) and corresponding surrogate model hyper-parameters α.

In both cases, the original generation/iteration procedure of CMA-ES is not modified, but called when necessary in order to change the state of CMA-ES in the search space, that Figure 1: Optimization loop of the s * ACM-ES (adopted from [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF]). also leads to a change of the internal parameters of CMA-ES (mean of mutation distribution m, covariance matrix C , step-size σ, evolution paths for σ and m).

First, the algorithm proceeds the optimization of f for a given number of generations gstart till a sufficient number of evaluated points is stored to build the surrogate model f α. The model is built using Ranking SVM and model hyperparameters α such that f α predicts only the ordering of test points, this allows f α-assisted search to be invariant w.r.t. rank-preserving transformations of f . The second invariance property, invariance w.r.t. orthogonal transformations of the search space, is preserved thanks to the following transformation applied to each of the training and test points x:

x ′ = C -1/2 (x -m), ( 4
)
where C is the covariance matrix adapted by CMA-ES. When the model is built, it is optimized by CMA-ES for a given number of generations n, where n is defined in [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF] by a linear function to be inversely proportional to a global model error Err(α). This global error is incrementally updated (with some relaxation factor) from a local error estimated with λ recently evaluated solutions stored in Λ as follows:

Err(α) = 2 |Λ|(|Λ| -1) |Λ| i=1 |Λ| j=i+1 wij .1 f α,i,j (5)
where 1 f α,i,j holds true if f α violates the f -based ordering constraint on pair (i, j) and wij defines the weight of this violation (wij = 1 by default). First (see Figure 1), s * ACM-ES i). builds the surrogate model f α, ii). optimizes it for n generations, then iii). optimizes the objective function f for 1 generation, and iv). estimates the surrogated model error by computing with Eq. ( 5) a fraction of comparison relations which were incorrectly predicted by f α with λ recently evaluated solutions, and finally v). adjusts n for the next iteration.

An important novelty of s * ACM-ES, apart from the preservation of the invariance properties of CMA-ES, is an ability to adapt hyper-parameters of surrogate learning procedure during the optimization of the objective function f . This procedure corresponds to the last step, illustrated in Figure 1, where an additional instance of CMA-ES (initialed before the loop of Figure 1) performs one optimization generation in the space of surrogate model hyper-parameters. In this step, CMA-ES generates λ hyp (λ hyp = 20 by default) different α-vectors of hyper-parameters and builds λ hyp corresponding surrogate models. These models are evaluated using Err(α) and µ hyp = ⌊λ hyp /2⌋ the most successful (with the smallest model error) out of λ hyp vectors of hyper-parameters are used to updates internal parameters of CMA-ES. The new mean of the mutation distribution of CMA-ES is then used as an estimate of the optimum in the space of surrogate model hyper-parameters. This estimate is finally used as α to build the surrogate model in the next iteration of s * ACM-ES, gradually adapting the surrogate model to the local topography of the objective function.

Finally, the algorithm optimizes the objective function f together with the hyper-parameters used to build its surrogate model f , that allows the user to define only the range of hyper-parameters and let s * ACM-ES to find their optimal values online.

NEW MECHANISM OF SURROGATE MODEL EXPLOITATION

As was mentioned before, in our experiments we observed that ACM-ES [START_REF] Loshchilov | Comparison-Based Optimizers Need Comparison-Based Surrogates[END_REF] outperforms s * ACM-ES on 10-dimensional Rosenbrock problem (but this is not the case in 20-D, where s * ACM-ES is better thanks to the adaptation of hyper-parameters). This result is quite surprising given that the surrogate model learning phases are very similar in both algorithms. These observations led to an hypothesis that the surrogate model exploitation is more intensive in ACM-ES on this particular problem and this is also beneficial because the surrogate is relatively accurate.

Intensive Exploitation

In the original ACM-ES algorithm, the surrogate model exploitation is independent on the model quality, that in certain cases may lead to a divergence of the algorithm when the surrogate provides almost random predictions. However, when the surrogate model is accurate enough, its exploitation by pre-selection substantially speed-ups CMA-ES [START_REF] Loshchilov | Comparison-Based Optimizers Need Comparison-Based Surrogates[END_REF]. In ACM-ES, λP re = 500 individuals are evaluated with f and only λ ′ = λ def ault 3

individuals among them (λ def ault is the default population size) are selected through a twostep selection procedure to be evaluated with the expensive function f . In s * ACM-ES, usually a smaller number of individuals than λP re is evaluated with f , but the change of the mean of the mutation distribution, and sometimes that is more important, of the step-size during n generations may lead to a larger divergence between the original and final multivariate normal distributions than with the pre-selection. Thus, s * ACM-ES is able to exploit f in the same way as CMA-ES exploits f , and if f ideally approximates f , then a speed-up of a factor of nmax (the maximum number of generations n when Err(α) = 0) is expected.

The first attempt to increase the intensity of exploitation by increasing nmax showed that larger values of nmax lead to moderate improvements on unimodal functions, but may increase the chance of premature divergence when the surrogate provides random predictions since a larger number of generations n is allowed. When the suggested n is much larger than some expected speed-up k (e.g., by a factor larger than 3), the search with the surrogate may start to deteriorate after ca. k generations by slowly becoming a random walk and, thus, destroying the evolution paths of CMA-ES.

Since longer runs may deteriorate the search, more attention should be paid to a single generation with f . In this generation, f can be more intensively exploited by estimating a local covariance matrix C loc constructed from a large number of points (e.g., in order of 10 5 ), sampled by CMA-ES and ranked according to f . This covariance matrix C loc usually is very similar to the one of CMA-ES, which was used for sampling, but in the same time, it stores "all" the information about covariances available from f for a given scale defined by the step-size. The use of C loc instead of the original estimate of CMA-ES allows to faster learn an appropriate covariance matrix, that in many cases leads to a faster convergence. However, the influence of C loc is somehow limited by the learning rate used to adapt the covariance matrix C of CMA-ES. In the same time, this influence should be controlled by the model error in order to avoid a possible degradation of C. The mechanism of this control might be very similar to the one already used in s * ACM-ES to control n.

Finally, we found a very simple exploitation mechanism which solves the above described problems and generalizes the second and third steps shown in Figure 1: 1. Optimize f for n generations by CMA-ES with population size λ = k λ λ def ault and number of parents µ = kµµ def ault , where k λ ≥ 1 and kµ ≥ 1.

2. Optimize f for 1 generation by CMA-ES with population size λ = λ def ault and number of parents µ = µ def ault .

This simple mechanism allows us to keep n in the order of the expected speed-up and, in the same time, to recover more information from f using a larger population size (when k λ > 1). In contrast to the estimation of C loc , the mean of the mutation distribution also changes (as well as other internal parameters) taking into account the information from λ = k λ λ def ault evaluated solutions. While, k λ and kµ can be adjusted depending on the model error, the control mechanism of n is already presented in s * ACM-ES.

Is should be noted that a more intensive exploitation occurs not when k λ is large and kµ = k λ , but rather when k λ is large and k λ ≫ kµ (e.g., kµ = 1).

Divergence Prevention

While the control mechanism of s * ACM-ES sets n to 0 when Err(α) reaches some surrogate error threshold value τerr (τerr= 0.45 by default [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF]), some oscillation between n = 0 and n = 1 is possible when Err(α) is very close to the threshold. In this case, a more intensive exploitation of f with a larger k λ might lead to a faster divergence of s * ACM-ES. The above described issue can be relevant for multimodal and noisy problems, where the surrogate model is typically inaccurate. To avoid the divergence, k λ can be adjusted according to Err(α) or n, e.g., k λ should increase together with n. Indeed, the optimal rule would be problemdependent and is difficult to choose a priori, but probably can be learned online. In this paper, in order to avoid a possible overfitting, we study a simple rule which states that

• k λ > 1 is used only if n ≥ nk λ ,
where nk λ is the number of generations which corresponds to a model error much smaller than τerr, such that the model is "accurate enough" to be intensively exploited.

EXPERIMENTAL VALIDATION

This Section presents an experimental validation of the proposed mechanism to surrogate model exploitation on a set of noiseless black-box optimization problems from the BBOB-2012 framework [START_REF] Hansen | Real-Parameter Black-Box Optimization Benchmarking[END_REF]. We compare the original CMA-ES [START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)[END_REF] and its so-called active version aCMA-ES [START_REF] Jastrebski | Improving Evolution Strategies through Active Covariance Matrix Adaptation[END_REF][START_REF] Hansen | Benchmarking a weighted negative covariance matrix update on the BBOB-2010 noiseless testbed[END_REF] in the IPOP and BIPOP scenarios of restarts [START_REF] Hansen | Benchmarking a weighted negative covariance matrix update on the BBOB-2010 noiseless testbed[END_REF][START_REF] Hansen | Benchmarking a BI-population CMA-ES on the BBOB-2009 function testbed[END_REF][START_REF] Loshchilov | Black-box Optimization Benchmarking of IPOP-saACM-ES and BIPOP-saACM-ES on the BBOB-2012 Noiseless Testbed[END_REF] to the results of s * ACM-ES in its original version [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF] and our modified version with the new mechanism of surrogate model exploitation. For all experiments, we use the Octave/Matlab source code provided1 by the authors of [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF] and, we make our modified version of this code also available online2 .

Experimental Setting

For all original algorithms the default parameters are used as given in [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF]. In our modified version of s * ACM-ES, apart from the k λ , kµ and nk λ , whose parameter settings will be discussed later in this Section, we also introduce an additional hyper-parameter to tune during the search: a stopping criterion niter of the quadratic programming solver of SVM. This parameter is set in [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF] to 1000Ntraining , where Ntraining is the number of training points. The online tuning of this parameter may potentially improve the quality of surrogate model learning and, in the same time, reduce the CPU complexity of the learning. A completed list of surrogate model hyper-parameters is summarized in Table 1, where offline tuned values and their ranges of variation for online tuning are given. A detailed description of these parameters can be found in [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF].

Parameter

Range for online tuning Offline tuned value Ntraining 4n, 2(40 + 4n

Results on Ellipsoid and Rosenbrock 20-D

First experiments were conducted on 20-dimensional Ellipsoid and Rosenbrock functions in order to investigate how the performance changes by changing k λ and kµ. We found that the greater the k λ , the better the performance is observed. Moreover, the results with kµ = 1 are usually better than with kµ = k λ , and this may be viewed as a more stronger pre-selection of individuals. This also can be interpreted as a trust-region method-based search, where the region is defined by the current covariance matrix and the step-size. In this context, n corresponds to the number of trust-region searches before the model is updated using new evaluated points.

Figure 2 shows the results for k λ = 1 (the original IPOPs * aACM-ES), k λ = 10, 100, 10000 (three modified versions with kµ = 1) as well as for the original and active versions of IPOP-CMA-ES. As can be clearly seen, the original IPOP-s * aACM-ES can be improved by a factor of about 2 on Ellipsoid function and by a factor of about 1.2 on Rosenbrock function for k λ ≥ 10. Thus, the new mechanism of surrogate model exploitation allows to be about 6 times faster than IPOP-aCMA-ES and almost by one order of magnitude faster than IPOP-CMA-ES. The algorithm demonstrates the speed-up of a factor of about 3.5 on Rosenbrock function, and, thus, performs as well as BFGS algorithm and by a factor between 1.1 and 1.4 slower than NEWUOA algorithm [START_REF] Hansen | Comparing results of 31 algorithms from the black-box optimization benchmarking BBOB[END_REF]. The performance of the latter algorithms, however, may degrade significantly for certain rank-preserving transformations of f (e.g., scaling of f ).

It should be noted that the computational complexity is growing linearly with k λ , and the evaluation of candidate solutions with f becomes relatively expensive (w.r.t. the cost of surrogate model learning) when k λ in 20-D is ca. 1000. The CPU time per function evaluation on F10 Rotated Ellipsoid 20-D (on 2.4 GHz processor and MatLab 2006) is 1.22, 0.313, 0.311, 0.367 and 1.52 sec. for k λ = 1, 10, 100, 1000 and 10000, respectively (note that with fixed hyper-parameters the cost would be λ hyp = 20 times smaller). Thus, k λ = 10 or k λ = 100 can be viewed as a reasonable trade-off for these problems given that the performance is (and if) growing only marginally for k λ > 10 (see Figure 2).

Figure 3 illustrates the dynamic of hyper-parameters adaptation on 20-dimensional Ellipsoid function for the original and modified (when k λ = 100) versions of IPOP-s * aACM-ES. As was expected, the number of training points also grows in the beginning, then becomes stable and finally decreases when the covariance matrix well approximates the optimal covariance matrix and the optimization of f becomes similar to the optimization of Sphere function. However, a smaller number of training points (and therefore CPU time) is usually needed for the new algorithm, that might be explained by the fact that the algorithm optimizes f faster and, therefore, outdated points are probably of a lesser importance.

Figure 3 also shows the adaptation of the number of iterations niter used to solve the quadratic programming problem of SVM learning. The results partially confirm our expectations that the offline tuned value of niter (see Table 1) was reasonably well chosen in [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF]. We also do not observe any moderate difference due to this adaptation on other BBOB benchmarks problems. This conjecture, however, should be more carefully verified by using a larger search range for ninit, and should not prevent from studying alternative strategies for SVM problem learning and corresponding stopping criteria.

The adaptation of other hyper-parameters is more difficult to interpret, however, in both cases, Cpow steadily increases as it was pointed out already in [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF].

Results on Noiseless BBOB Functions

The proposed modification of s * ACM-ES was implemented in BIPOP scenario of restarts and compared with BIPOP-CMA-ES [START_REF] Hansen | Benchmarking a BI-population CMA-ES on the BBOB-2009 function testbed[END_REF], BIPOP-aCMA-ES [START_REF] Loshchilov | Black-box Optimization Benchmarking of IPOP-saACM-ES and BIPOP-saACM-ES on the BBOB-2012 Noiseless Testbed[END_REF] and BIPOP-s * aACM-ES [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF] on 20-dimensional noiseless BBOB problems. In all experiments, the maximum number of function evaluations was set to 10 6 n. The number of function evaluations when surrogate models are used was set to 10 4 n to fit to the BBOB-2013 context of expensive optimization in contrast to 10 6 n used in BIPOP-s * aACM-ES [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF].

We found that the use of k λ = 10 for 10-D, k λ = 100 for 20-D and k λ = 1000 for 40-D (see [START_REF] Loshchilov | BI-population CMA-ES Algorithms with Surrogate Models and Line Searches[END_REF], here only 20-D case is shown) represents a reasonable trade-off between the speedup and the computational complexity (a simple formula to adjust k λ for other dimensions can be used). For problem dimensions smaller than 10, the use of large k λ may lead to a certain divergence even on unimodal functions. This can be interpreted as follows: when k λ is large, some of points sampled in 2-dimensional space might lie "far" from the current mean of the mutation distribution and "far" from the training set used to build the model. Thus, the model may incorrectly suggests these points, that may lead to a certain divergence. The use of some truncated sampling of points might resolve this issue, but this would lead to a change of mechanism of CMA-ES which expects that all points are sampled from multivariate normal distribution.

In order to avoid a possible divergence on multimodal functions, nk λ usually should be larger than 1. By some offline tuning we found that the use of nk λ = 3 and nk λ = 4 in most of the cases prevents the divergence (at least on tested BBOB problems), and we use nk λ = 4 in our experiments on BBOB. A simple alternative to this would be to decrease nmax or τerr.

Figure 4 shows the aggregated results for all 24 noiseless BBOB problems and different subgroups of these problems. The proposed version of BIPOP-s * aACM-ES is depicted as BIPOP-s * aACM-ES-k. The best improvements of results is observed on moderate and ill-conditioned functions, where the speed-up up to a factor of 2 is observed. On ill-conditioned problems, BIPOP-s * aACM-ES-k is about 3-6 times faster than BIPOP-aCMA-ES and by almost one order of magnitude faster than BIPOP-CMA-ES. The results of BIPOP-s * aACM-ES and BIPOP-s * aACM-ES-k are similar on multimodal functions, while BIPOP-s * aACM-ES slightly outperforms the latter for the budgets of function evaluations larger than 10 4 n, since the surrogate is still in use. The overall ranking of all algorithms is comparable and a difference is observed most likely due to the stochasticity of the process (at least for BIPOP-aCMA-ES and BIPOP-s * aACM-ES-k, which represent the same algorithm after 10 4 n function evaluations).

CONCLUSION AND PERSPECTIVES

We presented a new mechanism for surrogate model exploitation, where a larger population size is used to optimize the surrogate, while the default population size is used to optimize the expensive function. We found that while a more intensive exploitation of the surrogate is beneficial on unimodal functions, this also might lead to a certain divergence on multimodal functions. To avoid this potential issue, we suggested to more intensively exploit the surrogate only if it is sufficiently accurate. The proposed mechanism implemented in s * ACM-ES algorithm and its BIPOP extension demonstrates the speed-up of a factor up to 2 (compared with the original version of s * ACM-ES) on ill-conditioned BBOB problems. The resulting algorithm also preserves the properties of invariance of CMA-ES and might become a baseline version of s * ACM-ES in some near future.

The main perspective for further research is to replace the newly introduced offline-tuned parameters by some online procedure, which will be able to optimally control surrogate model exploitation. Such procedure might be based on some a posteriori analysis of what strategy was optimal in previous generations by analyzing (running experiments on) previous surrogate models and alternative exploitation strategies. Another perspective is to reduce the time complexity of the algorithm which limits its application to n 50 -100 and up to ca. 5000 training points (see, e.g., [START_REF] Loshchilov | Black-box Optimization Benchmarking of IPOP-saACM-ES and BIPOP-saACM-ES on the BBOB-2012 Noiseless Testbed[END_REF]).

Figure 2 :

 2 Figure 2: Comparison of the proposed surrogate-assisted versions of the original and active IPOP-CMA-ES algorithms on 20-dimensional Rotated Ellipsoid (Left) and Rosenbrock (Right) functions. The trajectories show the median of 15 runs.

Figure 3 :

 3 Figure 3: The median trajectories of normalized surrogate hyper-parameters estimated on 15 runs of the original (Left) and modified (Right) versions of IPOP-s * aACM-ES on 20-dimensional Rotated Ellipsoid.

Figure 4 :

 4 Figure 4: Bootstrapped empirical cumulative distribution of the number of objective function evaluations divided by dimension for 50 targets in 10 [-8..2] for all functions and subgroups in 20-D. The "best 2009" line corresponds to the best ERT observed during BBOB 2009 for each single target. The proposed algorithm is depicted as BIPOP-saACM-k. A detailed description how to read the figure is given in [8].

Table 1 :

 1 Surrogate hyper-parameters, default value and range of variation (extended from[START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF]).

	1.7 )	40 + 4n 1.7

https://sites.google.com/site/acmesgecco/

https://sites.google.com/site/newacmes/


Sylvain Chevillard
email: sylvain.chevillard@inria.fr

Marc Mezzarobba

Multiple-precision evaluation of the Airy Ai function with reduced cancellation

Keywords: Special functions, algorithm, numerical evaluation, arbitrary precision, Miller method, asymptotics, correct rounding, error bounds

The series expansion at the origin of the Airy function Ai(x) is alternating and hence problematic to evaluate for x > 0 due to cancellation. Based on a method recently proposed by Gawronski, Müller, and Reinhard, we exhibit two functions F and G, both with nonnegative Taylor expansions at the origin, such that Ai(x) = G(x)/F(x). The sums are now well-conditioned, but the Taylor coefficients of G turn out to obey an ill-conditioned three-term recurrence. We use the classical Miller algorithm to overcome this issue. We bound all errors and our implementation allows an arbitrary and certified accuracy, that can be used, e.g., for providing correct rounding in arbitrary precision.

Many mathematical functions (e.g., trigonometric functions, erf, Bessel functions) have a Taylor series of the form

y(x) = x s ∞ ∑ n=0 y n x dn , y n ∼ (-1) n λ α n n! κ (1)
with d, s ∈ Z and α, κ > 0. For large x > 0, the computation in finite precision arithmetic of such a sum is notoriously prone to catastrophic cancellation. Indeed, the terms |y n x n | are first growing before the series "starts to converge" when n κ ≥ αx.

In particular, when n κ ≈ αx, the terms y n x n usually get much larger than y(x). Eventually, their leading bits cancel out while lower-order bits that actually contribute to the first significant digits of the result get lost in roundoff errors. This cancellation phenomenon makes the direct computation by Taylor series impractical for large values of x. Often, the function y(x) admits an asymptotic expansion as x → +∞ that can be used very effectively to obtain numerical approximations when x is large, but might not provide enough accuracy (at least without resorting to sophisticated resummation methods) for intermediate values of x.

In the case of the error function erf(x), a classical trick going back at least to Stegun and Zucker [START_REF] Stegun | Automatic computing methods for special functions[END_REF] is to compute erf(x) as G(x)/F(x) where F(x) = e x 2 and [1, Eq. 7.6.2]

G(x) = e x 2 erf(x) = 2x √ π ∞ ∑ n=0 2 n 1 • 3 • • • (2n + 1) x 2n . (2)
The benefit of this transformation is that F and G are power series with nonnegative coefficients, and can thus be computed without cancellation. Algorithms based on (2) tend to behave well in some range a < x < b where x is large enough for cancellation to be problematic but small enough to make the use of asymptotic expansions at infinity inconvenient. Note that the obvious way to compute y(x) = e -x for x > 0 fits into the same framework, now with G(x) = 1 and F(x) = e x . Gawronski, Müller and Reinhard [START_REF] Gawronski | Reduced cancellation in the evaluation of entire functions and applications to the error function[END_REF], [START_REF] Reinhard | Reduced cancellation in the evaluation of entire functions and applications to certain special functions[END_REF] provide elements to understand where these rewritings "come from". They relate the amount of cancellation in the summation of a series [START_REF]Digital library of mathematical functions[END_REF] to the shape of the Phragmén-Lindelöf indicator of y, a classical tool from the theory of entire functions [START_REF] Levin | Lectures on entire functions[END_REF]. This description allows them to state criteria for choosing auxiliary series suitable for the evaluation of a given entire function in a given sector of the complex plane. They apply their method (called the "GMR method" in what follows) to obtain "reduced cancellation" evaluation algorithms for the error function and other related functions in various sectors.

In this article, we are interested in the evaluation for positive x of the Airy function Ai [START_REF]Digital library of mathematical functions[END_REF]Chap. 9]. The function Ai(x) satisfies the linear ordinary differential equation (LODE) Ai ′′ (x)x Ai(x) = 0 [START_REF] Brent | Modern computer arithmetic[END_REF] with initial values Ai(0) = A := 3 -2/3 Γ( 2 3 ) -1 , Ai ′ (0) = -B := -3 -1/3 Γ( 13 ) -1 . The classical existence theorem for LODE with complex analytic coefficients implies that Ai(x) is an entire function; and solving (3) by the method of power series yields the Taylor expansion Ai(x) = A f (x 3 ) -Bxg(x 3 ), where

f (x) = ∞ ∑ n=0 1•4• • •(3n -2) (3n)! x n , g(x) = ∞ ∑ n=0 2•5• • •(3n -1) (3n + 1)! x n .
Observe that while f and g are easy enough to evaluate individually, the difference A f (x 3 )-Bxg(x 3 ) causes catastrophic cancellation when computed in approximate arithmetic.

Using the GMR method, we derive a reduced cancellation algorithm for computing Ai(x). To our best knowledge, our algorithm for evaluating Ai(x) is new, and is the most efficient multiple-precision evaluation of Ai(x) when x is neither too small nor too large, while the precision is not large enough to make methods based on binary splitting [START_REF] Brent | Modern computer arithmetic[END_REF] competitive.

Besides the new application, the main difference between the present article and the work of Gawronski et al. is our setting of multiple-precision arithmetic "à la MPFR [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF]". Specifically, on the one hand, we are interested in arbitrary precision arithmetic rather than machine precision only. This makes it impossible, for instance, to tabulate the coefficients of auxiliary functions when these turn out to be hard to compute. Also, we are looking for rigorous error bounds instead of experimental error estimates. On the other hand, we restrict ourselves to numerical evaluation on a half line instead of a complex sector (though, in principle, the basic ideas generalize).

This article focuses on providing a complete algorithm in the specific case of Ai(x), x > 0. Yet, it should also be seen a case study, part of an effort to understand what the GMR method can bring in the context of multiple-precision computation, and, perhaps more importantly, how general and systematic it can be made. We discuss this last point further in Sec. VIII.

The rest of this text is organized as follows. In Sec. I, we use the GMR method to choose the functions F and G. Then, in Sec. II and III, we derive a few mathematical properties of these functions, including recurrences for their series expansions and various bounds. Sections IV to VI contain the details of our algorithm and its error analysis. Finally, in Sec. VII, we briefly describe our implementation of the algorithm.

I. The GMR method

We now review the GMR method and apply it to obtain candidate auxiliary series for the evaluation of the Ai function. Since the method itself is not crucial for our results, we summarize it in intuitive terms and refer the reader to the original works [START_REF] Gawronski | Reduced cancellation in the evaluation of entire functions and applications to the error function[END_REF], [START_REF] Reinhard | Reduced cancellation in the evaluation of entire functions and applications to certain special functions[END_REF] for more careful statements.

The starting point of the GMR method is the following observation. Let y(z) = ∑ n≥0 y n z n be an entire function. Assume that we have, in some intentionally vague sense,

y(re iθ ) ≈ exp(h(θ)r ρ ) (4)
for large r. (To make things precise, we would assume that y has finite order ρ, and that h is its indicator function with respect to ρ [START_REF] Levin | Lectures on entire functions[END_REF].)

We consider the computation of y(z) in floating-point arithmetic using its series expansion. It is well-known [START_REF] Chevillard | The functions erf and erfc computed with arbitrary precision and explicit error bounds[END_REF] that, if the sum is performed in floating-point arithmetic of precision t, the relative error between y(z) and the computed sum is roughly given by 2 Hence, the loss of precision by cancellation in the evaluation of y(re iθ ) is about

log 2 M(r) |y(re iθ )| ≈ [(max h) -h(θ)] r ρ .
For instance, when the y n all have the same complex argument, the maximum of h is reached for θ = 0, in accordance with the fact that the sum is optimally conditioned. In the case of the Airy Ai function, the following asymptotic equivalent holds as z tends to complex infinity in any open sector that avoids the negative real axis [1, Eq. 9.7.5]:

Ai(z) ∼ ̃ Ai(z) := exp(-2 3 z 3/2 ) 2 √ πz 1/4 . (5)
Additionally, Ai(x) is bounded for x < 0. Hence, we may take ρ = 3/2 and

h(θ) = -2 3 cos( 3 2 θ), -π < θ ≤ π (6)
(see Figure 1). The loss of precision is roughly proportional to 1 + cos( 3 2 θ). It is minimal in the directions of fastest growth θ = ± 2 3 π, and maximal for θ = 0. If now two entire functions y and F both satisfy conditions of the form (4) with the same ρ but different h (say h y and h F , respectively), we may expect that

G(z) = F(z)y(z) ≈ exp([h y (θ) + h F (θ)]r ρ ). ( 7
)
The GMR method consists in reducing the summation of the series y for z in some given sector to that of an auxiliary series F(z) and a modified series G(z) related by [START_REF] Gawronski | Reduced cancellation in the evaluation of entire functions and applications to the error function[END_REF]. The value of y(z) is then recovered as G(z)/F(z). The auxiliary series is chosen, based on the shape of h y , so that both h F and h G = h y + h F take values close to their maximum in the sector of interest.

There may be multiple choices, and it is not clear in general which one is better, except that the coefficients of F and G should be as easy to compute as possible. Gawronski et al. usually take F(z) = exp(az ρ ) and search for a value of a that makes (max h G )h G as small as possible on a whole complex sector. The choice of exponentials as auxiliary series is not appropriate in the case of Ai, since exp(z ρ ) is an entire function of z only for integer ρ.

However, as we are interested in one direction only, we can easily build a suitable auxiliary series from Ai itself. Indeed, we may "shift" the indicator function of Ai by 2π/3 to the left or to the right by changing z to j ±1 z, where j = e 2πi/3 . (Note that this is not the same as changing θ to θ ± 2 3 π in (6).) When we add such a shifted indicator to the original h Ai , one of the humps of the curve cancels out with the valley in the middle.

Using this idea, we set

F(x) = Ai( jx) Ai( j -1 x), G(x) = F(x) Ai(x). ( 8
)
The indicator functions of F and G are pictured on Figure 2. Based on their shapes, we expect that both series are optimally conditioned on the positive real axis. We shall prove that this is indeed the case in the next two sections. This second method of constructing auxiliary series seems to be new, and applies to many cases. For instance, applying it to the error function leads to

F(x) = -i erf ix = 2x √ π ∞ ∑ n=0 1 • 3 • • • (2n -1) (2n + 1)! 2 n x 2n , (9)
G(x) = F(x) erf x = 8x 2 π ∞ ∑ n=0 2 • 4 • • • (4n) (4n + 2)! 4 n x 4n , (10)
a slightly worse alternative to [START_REF] Bickley | Bessel functions[END_REF]. The advantage of (2) comes from the fact that e x 2 is faster to evaluate than [START_REF] Levin | Lectures on entire functions[END_REF].

II. The auxiliary series F

The functions F and G being chosen, we need to establish appropriate formulas to evaluate them, along with error bounds. Much of our analysis will be based on the following simple estimate [START_REF] Olver | Asymptotics and special functions[END_REF]Chap. 4,§4.1], where ̃ Ai was defined in Eq. [START_REF] Flajolet | Analytic combinatorics[END_REF].

Lemma 1: The Airy function Ai satisfies

|Ai(re iθ )/ ̃ Ai(re iθ ) -1| ≤ η 1 (θ)r -3/2 , |θ| < π,
where η 1 (θ) = 5 48 (cos θ 2 ) -7/2 . Now consider the power series expansion of F at the origin:

F(x) = Ai( jx) Ai( j -1 x) = ∞ ∑ n=0 F n x n . ( 11
)
Proposition 1: The coefficients F n are positive and satisfy the two-term recurrence relation

(n + 1)(n + 2)(n + 3)F n+3 -2(2n + 1)F n = 0 ( 12
)
with initial values

F 0 = 3 -4/3 Γ( 2 3 ) -2 , F 1 = (2 √ 3π) -1 , F 2 = 3 -2/3 Γ( 1 3 ) -2
. Proof: As a general fact, if two functions w and y each satisfy a homogeneous LODE with coefficients in Q(x), then their product wy satisfies an equation of the same class that can be explicitly computed [START_REF] Stanley | Enumerative combinatorics[END_REF]Sec. 6.4]. The functions Ai( j ±1 x) satisfy the same differential equation (3) as Ai(x) itself. Applying the procedure mentioned above to two copies of that equation yields F (3)

(x) -4xF ′ (x) -2F(x) = 0.
Similarly, when an analytic function y satisfies a homogeneous LODE over Q(x), we can compute a recurrence relation with coefficients in Q(n) on the coefficients y n of its power series expansion. In the case of F, we get [START_REF] Olver | Asymptotics and special functions[END_REF]. Finally, we compute the initial values F 1 , F 2 , F 3 from the first few terms of the Taylor expansion of Ai(x):

F(x) = (A -B jx + O(x 3 ))(A -B j -1 x + O(x 3 )) = A 2 -ABx + B 2 x 2 + O(x 3 ).
It is then apparent from ( 12) that F n > 0 for all n ∈ N. Thus, the coefficients of F(x) obey a two-term recurrence whose coefficients do not vanish for n ≥ 0. This allows one to compute them in a numerically stable way (see Sec. VI).

III. The modified series G

Recall that G(x) = Ai(x) Ai( jx) Ai( j -1 x), and set

G(x) = ̃ Ai(x) ̃ Ai( jx) ̃ Ai( j -1 x) = exp( 2 3 x 3/2 ) 8π 3/2 x 3/4 . ( 13
)
Proposition 2: The function G is an entire function with power series expansion of the form

G(z) = ∑ ∞ n=0 G n z 3n . The coefficient sequence (G n ) n∈N is determined from its first terms G 0 = A 3 = 1 9Γ(2/3) 3 , G 1 = A 3 2 -B 3 = 1 18Γ(2/3) 3 - 1 3Γ(1/3) 3 by the recurrence relation (n+1)(3n+4)(3n+5)(n+2)G n+2 -10(n+1) 2 G n+1 +G n = 0. ( 14
)
Proof: First, observe that G( jz) = G(z) and G(z) = G(z), so that the Taylor expansion of G at the origin is a power series in z 3 with real coefficients. The same routine reasoning as in the proof of Prop. 1 yields the LODE

G (4) (x) -10xG ′′ (x) -10G ′ (x) + 9x 2 G(x) = 0,
and from there the recurrence [START_REF] Reinhard | Reduced cancellation in the evaluation of entire functions and applications to certain special functions[END_REF]. The coefficients of ( 14) do not vanish for n ≥ 0, so that the sequence G n is indeed determined by G 0 , G 1 , and [START_REF] Reinhard | Reduced cancellation in the evaluation of entire functions and applications to certain special functions[END_REF].

Nonzero solutions of ( 14) decrease roughly as n! -2 for large n. Setting c n = n! 2 G n yields the "normalized" recurrence (3n + 4)(3n

+ 5) (n + 1)(n + 2) c n+2 -10c n+1 + c n = 0. ( 15
)
Letting n go to infinity in the coefficients of (15), we get a limit recurrence with constant coefficients whose characteristic polynomial 9α 2 -10α + 1 has two roots of distinct absolute value, namely 1 9 and 1. By the Perron-Kreuser theorem [21, Theorem B.10], it follows that any solution (v n ) n∈N of ( 14)

satisfies v n+1 /v n ∼ αn -2 with either α = 1 or α = 1 9 . Solutions (v n ) such that v n+1 /v n ∼ 1
9n 2 are called minimal and form a linear subspace of dimension 1 of the solutions of ( 14).

We shall prove that (G n ) actually is such a minimal solution of [START_REF] Reinhard | Reduced cancellation in the evaluation of entire functions and applications to certain special functions[END_REF]. But our analysis uses a bit more than the rough estimate G n ≈ n! -2 9 -n . Prop. 3 below provides a more precise estimate which implies the minimality. Before turning to it, we recall a standard bound on the tails of incomplete Gaussian integrals [1, Eq. 7.12.1] and state a second technical lemma.

Lemma 2: The complementary error function erfc

x = 1 - erf x = 2 √ π ∫ ∞ x e -t 2 dt satisfies 0 ≤ erfc x ≤ 1 √ πx e -x 2 for x > 0.
a(r) e -ϕ 2 dϕ ≤ ∫ ∞ a(r) e -ϕ 2 dϕ = √ π 2 erfc ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ √ 3r 1/8 √ 5 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ ,
and, by Lemma 2, I(r) ≤ 5 6 r 1/4 e -3 5 r 1/4 . When r ≥ 10, the last bound is decreasing and hence less than 0.51.

The following proposition is the main result of this section and the starting point of much of the error analysis that follows. Though somewhat technical, the proof is mostly routine.

Proposition 3: The sequence (G n ) satisfies

G n ∼ γ n = 1 4 √ 3π9 n n! 2 , n → ∞, with relative error |G n /γ n -1| ≤ 2.2n -1/4 (n ≥ 1). Notation 1: We write E(ε 1 , . . . , ε n ) = ∑ ∅ I⊂ 1,n ∏ i∈I ε i , so that, if |θ i | ≤ ε i for all i, then ∏ n i=1 (1 + θ i ) = 1 + θ with |θ| ≤ E(ε 1 , . . . , ε n ). Note that we obviously have E(α ε 1 , . . . , α ε n ) ≤ α E(ε 1 , . . . , ε n ) for 0 ≤ α ≤ 1.
Proof: The proof is a standard application of the saddle point method [5, §VIII.3], which we work out in some detail in order to get an explicit error bound.

We fix n > 10. We shall write z = re iθ in the following. The method prescribes to choose r = (3n + 3/4) 2/3 , so that G(z)

z 3n = G(r)
r 3n e O((z-r) 2 ) (

(i.e., d dz log( G(z)z -3n ) |z=r = 0). Now, guided by [START_REF] Smith | Efficient multiple-precision evaluation of elementary functions[END_REF] and G(z) ≈ G(z), we write

G n = 1 2πi ∮ G(z) z 3n+1 dz = 3 2π G(r)
r 3n I 1 , I 1 = ∫ π/3 -π/3 G(re iθ ) G(r) e -3niθ dθ.
Most of the weight of I 1 is concentrated around 0. We set

θ 0 = r -5/8 = (3n + 3/4) -5/12 . ( 17
)
Since n > 10, we remark that r ≥ 10 and 0 < θ 0 < 1/4 < π/3. We further let

I 2 = ∫ θ 0 -θ 0 G(re iθ ) G(r) e -3niθ dθ.
We first bound the error between I 1 and I 2 . When θ ∈ [0, π/3], using Lemma 1 and the connection formula [1, Eq. 9.2.12]

Ai(z) + j Ai( jz) + j -1 Ai( j -1 z) = 0, we see that |G(z)/ G(z)| is bounded by ⃒ ⃒ ⃒ ⃒ Ai(z)  Ai(z) ⃒ ⃒ ⃒ ⃒ ⃒ ⃒ ⃒ ⃒ Ai( j -1 z)  Ai( j -1 z) ⃒ ⃒ ⃒ ⃒ (⃒ ⃒ ⃒ ⃒ ̃ Ai(z)  Ai( jz) ⃒ ⃒ ⃒ ⃒ ⃒ ⃒ ⃒ ⃒ Ai(z)  Ai(z) ⃒ ⃒ ⃒ ⃒ + ⃒ ⃒ ⃒ ⃒ ̃ Ai( j -1 z)  Ai( jz) ⃒ ⃒ ⃒ ⃒ ⃒ ⃒ ⃒ ⃒ Ai( j -1 z)  Ai( j -1 z) ⃒ ⃒ ⃒ ⃒ ) ≤(1 + η 1 (θ)r -3/2 )(1 + η 1 (θ -2π 3 )r -3/2 ) ( e -4 3 cos 3θ 2 (1 + η 1 (θ)r -3 2 ) + 1 + η 1 (θ -2π 3 )r -3 2 ) .
It follows that |G(z)| ≤ 2.1| G(z)| when |z| ≥ 10. By symmetry, this inequality holds for -π/3 ≤ θ ≤ 0 too. Finally, we have

|I 1 -I 2 | ≤ 2 ∫ π/3 θ 0 |G(re iθ )| G(r) dθ ≤ 4.2 ∫ π/3 θ 0 | G(re iθ )| G(r) dθ = 4.2 ∫ π/3 θ 0 e 2 3 r 3/2 (-1+cos 3θ
2 ) dθ and hence, using Lemma 3,

|I 1 -I 2 | ≤ 2.15r -9/8 (18)
We now have to estimate I 2 . We write

I 2 = ∫ θ 0 -θ 0 G(re iθ ) G(re iθ ) • G(re iθ ) G(r) e -3niθ dθ.
On the one hand, Lemma 1 gives

G(z) = G(z)(1 + δ(θ)) with |δ(θ)| ≤ η 2 (θ)r -3/2 where η 2 (θ) = E(η 1 (θ), η 1 (θ + 2π 3 ), η 1 (θ -2π 3
)). Note that this bound increases with |θ|.

On the other hand, thanks to [START_REF] Smith | Efficient multiple-precision evaluation of elementary functions[END_REF], we can write G(re iθ )e -3niθ = G(r)e

2 3 r 3/2 (exp( 3 2 iθ)-1-3 2 iθ) = G(r)e (-3 4 θ 2 +u(θ))r 3/2 where |u(θ)| = 2 3 ⃒ ⃒ ⃒ ⃒ e 3 2 iθ -1 -3 2 iθ + 9 8 θ 2 ⃒ ⃒ ⃒ ⃒ ≤ 3 8 |θ| 3 . (19) Let v(θ) = G(re iθ )e -3inθ G(r) exp(-3 4 r 3/2 θ 2 ) -1.
Since θ 3 0 r 3/2 = r -3/8 by the choice (17), using [START_REF] Van Der Put | Galois theory of linear differential equations[END_REF] and the inequality |e z -1| ≤ |z|e |z| , we have, for |θ| ≤ θ 0 and r ≥ 10:

|v(θ)| ≤ 3 8 r -3/8 e 3 8 r -3/8 ≤ 0.44r -3/8 .
To sum up, we can now rewrite I 2 as

I 2 = ∫ θ 0 -θ 0 (1 + δ(θ))(1 + v(θ)) exp ( - 3 4 r 3/2 θ 2 ) dθ.
Now, define

I 3 = ∫ θ 0 -θ 0 exp ( - 3 4 r 3/2 θ 2 ) dθ, I 4 = 2 √ π √ 3 r -3/4 .
When θ ∈ [-θ 0 , θ 0 ], we have

|(1 + δ(θ))(1 + v(θ)) -1| ≤ E(η 2 (θ 0 )r -3/2 , 0.44r -3/8 ) ≤ 0.95r -3/8 , so we get |I 2 -I 3 | ≤ 0.95r -3/8 I 3 . (20
) Finally, I 3 is an incomplete Gaussian integral:

I 3 = 2 √ π √ 3r 3/4 erf( 1 2 √ 3r 3/4 θ 0 ) = I 4 erf( 1 2 √ 3r 1/8 ).
Lemma 2 yields

|I 3 /I 4 -1| ≤ 2r -1/8 e -3 4 r 1/4 √ 3π ≤ 0.66r -1/8 e -3 4 r 1/4 . (21)
Putting together [START_REF] Stegun | Automatic computing methods for special functions[END_REF][START_REF] Wasow | Asymptotic expansions for ordinary differential equations[END_REF][START_REF] Wimp | Computation with recurrence relations[END_REF], we obtain

I 1 = I 4 (1 + ) with | | ≤ |I 1 -I 2 | I 4 + ⃒ ⃒ ⃒ ⃒ ⃒ I 2 I 3 I 3 I 4 -1 ⃒ ⃒ ⃒ ⃒ ⃒ ≤ 1.06r -3/8 + E(0.95r -3/8 , 0.66r -1/8 e -3 4 r 1/4 ) ≤ 2.45r -3/8 ≤ 1.9n -1/4 .

Now, we can write

G n = 3 2π • G(r) r 3n I 4 (1 + ) = γ n H (1 + ) with H := √ 3 √ πr 3/4 γ n G(r) r 3n = n! 2 e 2n n 2n 2πn √ e (1 + 1 4n ) 2n+1
.

Stirling's formula in the form [START_REF] Robbins | A remark on Stirling's formula[END_REF]

⃒ ⃒ ⃒ ⃒ n!e n n n √ 2πn -1 ⃒ ⃒ ⃒ ⃒ ≤ s(n) = e 1 12n -1, combined with the bound ⃒ ⃒ ⃒ ⃒ √ e (1+1/(4n)) 2n -1 ⃒ ⃒ ⃒ ⃒ ≤ 1 16n , yields |H -1| ≤ E ( s(n), s(n), 1 16n , 1 4n ) ≤ 1 2n , n ≥ 11. ( 22
)
It follows that

⃒ ⃒ ⃒ ⃒ ⃒ G n γ n -1 ⃒ ⃒ ⃒ ⃒ ⃒ ≤ E(0.5n -1 , 1.9n -1/4 ) ≤ 2.2n -1/4 , n ≥ 11,
One easily checks that this bound is valid for 1 ≤ n ≤ 10. Note that the above bound is not the best we can get by this method. Indeed, by choosing the exponent of r in (17) closer to -3/4, we obtain |G n /γ n -1| = O(n -1/2+ ) for any > 0. This comes at the price of a larger constant factor and thus more terms to check separately.

A first consequence of Prop. 3 is that the series G(x) has nonnegative coefficients, as stated below. We also deduce several other technical results that will be used to bound various error terms. Recall that c n = n! 2 G n , and let τ = 3/20.

Corollary 1: The sequence (c n ) satisfies 0 ≤ c n+1 /c n ≤ τ for all n ≥ 0. Accordingly, we have 0 < G n+1 /G n ≤ τ(n + 1) -2 . In particular, the G n are positive.

Proof:

Prop. 3 implies that c n = (1 + θ n )/(4 √ 3π9 n ) with |θ n | ≤ 2.4n -1/4 . When n ≥ n 0 = 47610, we have |θ n | ≤ 0.148936. This implies c n ≥ 0 and c n+1 /c n ≤ (1/9)(1 + θ n+1 )/(1 + θ n ) ≤ τ.
The inequalities up to n = n 0 are checked by computing the corresponding terms with interval arithmetic.

Corollary 2: For any n ≥ 1, it holds that G n ≤ (e/(3n)) 2n .

Proof: Follows from Proposition 3 using n! ≥ (n/e) n .

Lemma 4: When x > 0 and N + 1 ≥ √ 2τx 3/2 , we have

N ∑ n=0 G n x 3n ≤ G(x) ≤ N-1 ∑ n=0 G n x 3n + 2G N x 3N .
Proof: The first inequality is obvious as G n ≥ 0. To show the second one, observe that a fortiori τx 3 /(n + 1) 2 ≤ 1 2 for all n ≥ N. Using Corollary 1, we deduce G n+1 x 3 /G n ≤ 1 2 , and hence

∑ n≥N G n x 3n ≤ G N x 3N (1 + 1 2 + 1 4 + • • • ) = 2G N x 3N . Lemma 5: The function G satisfies 0.01e 2 3 x 3/2 x -3/4 ≤ G(x) ≤ 0.04e 2 3 x 3/2 x -3/4 for all x ≥ 1/2. Proof: Let µ(x) = E(η 1 (0)x -3/2 , σx -3/2 , σx -3/2 ) where σ = η 1 ( 2π 3 ) = 5 6 √ 2. Lemma 1 implies |G(x)/ G(x) -1| ≤ µ(x) for x ≥ 0. We can check that µ(x) is a decreasing function and 0.01 ≤ 1 -µ(3) 8π 3/2 ,
1 + µ(3) 8π 3/2 ≤ 0.04, whence the desired inequality for x ≥ 3. For 0.5 < x ≤ 3, using the bounds 0 ≤ e x -(1 + x + 1 2 x 2 + 1 6 x 3 + 1 24 x 4 ) ≤ 5 48 x 4 (valid for 0 < x ≤ 2 3 3 3/2 ) and Lemma 4, we are reduced to checking explicit polynomial inequalities.

IV. Roundoff error analysis

We now turn to the floating-point implementation of the functions F(x) and G(x). To make the algorithm rigorous, we will use classical techniques of error analysis that we briefly recall here. We refer the reader, e.g., to Higham [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF], for proofs and complement of information.

We suppose that the precision of the floating-point format is t bits and that the exponent range is unbounded (in case it is bounded, it would probably be possible to rescale F(x) and G(x) by the same factor, to make them representable without changing the ratio F(x)/G(x)).

Notation 2: For x 0, we denote Exp(x) = ⌊log 2 |x|⌋ + 1, so that 2 Exp(x)-1 ≤ |x| < 2 Exp(x) .

Notation 3: If x ∈ R, ∘(x) denotes the floating-point number closest to x (ties can be decided either way). Circled operators such as ⊕ denote correctly rounded floating-point operations.

We always have ∘(x) = x (1 + δ) and x = ∘(x) (1 + δ ′ ) with |δ|, |δ ′ | ≤ 2 -t . We will also extensively use the relative error counter notation z ⟨k⟩.

Notation 4: We write ̂ z = z ⟨k⟩ when there exist δ

1 , . . . , δ k such that ̂ z = z ∏ k i=1 (1 + δ i ) ±1 with |δ i | ≤ 2 -t
for all i. Roughly speaking, each arithmetical operation adds one to the relative error counter of a variable. The overall error corresponding to an error counter can be bounded as follows.

Proposition 4: Suppose that we can write ̂ z = z ⟨k⟩ and that k

2 -t ≤ 1/2. Then ̂ z = z(1 + θ) with |θ| ≤ 2k • 2 -t .

V. Evaluation of the modified series

As we shall see in the next section, evaluating the auxiliary function F is fairly straightforward. The evaluation of G is more involved. Indeed, while ∑ G n x 3n is well-conditioned as a sum for x ≥ 0 (this is the whole point of the GMR method), the minimality of the sequence (G n ) among the solutions of [START_REF] Reinhard | Reduced cancellation in the evaluation of entire functions and applications to certain special functions[END_REF] implies that its direct recursive computation from the initial values G 0 and G 1 is numerically unstable (cf. [START_REF] Wimp | Computation with recurrence relations[END_REF]).

Algorithm 1: Evaluation of G

Input: a target precision p ≥ 1, a point x ≥ 0.5 Output: s such that |G(x) -s| ≤ 3 • 2 -p G(x) 1 Choose α, β, δ, γ s.t. α 3 e -1 x -3/2 , β (2/3) log 2 (e) x 3/2 , γ 1/ log 2 (20/3), δ (2/3) log 2 (e) ((20/3) 1/2 -1)

x 3/2 ; 2 N 0 ← max(1, ⌈(3/10) 1/2 x 3/2 -1⌉); 3 Choose N ≥ N 0 s.t. Exp((α N) 2N ) ≥ p + 9 + 3 4 Exp(x) -⌊β⌋; 4 Choose R ≥ max(N, (p + 2 + δ) γ); 5 Choose t s.t. 128 (N + 3) 2 -t ≤ 2 -p and (R + 2) 2 -t ≤ 2 -9 ; 6 for (a ← 1, b ← 0, i ← R -1; i ≥ 0; i ← i -1) do 7 c ← a; 8 a ← (10 ⊗ a) ⊖ (3i + 4)(3i + 5) ⊗ b ⊘ ((i + 1)(i + 2)); 9 b ← c; 10 if i = N -1 then s ′ ← a; 11 else if i < N -1 then s ′ ← a ⊕ (x 3 ⊗ s ′ ⊘ (i + 1) 2 ); 12 return s = (s ′ ⊘ a) ⊘ ∘(9 Γ(2/3) 3 );
There is a standard tool to handle this situation, namely Miller's backward recurrence method [START_REF] Bickley | Bessel functions[END_REF], [START_REF] Wimp | Computation with recurrence relations[END_REF]. Miller's method allows one to accurately evaluate the minimal solution (c n ) of a recurrence of the form

a 2 (n) u n+2 + a 1 (n) u n+1 + a 0 (n) u n = 0, a 0 (n)a 2 (n) 0. (23)
The idea is as follows: choose a starting index R and let (arbitrarily) u R+1 = 0 and u R = 1. Then compute u n as

-a 0 (n) -1 (a 1 (n) u n+1 + a 2 (n) u n+2 ), n = R -1, . . . , 1, 0.
It turns out that, for large R, the computed sequence (u n ) is close to a minimal solution of the forward recurrence. Since all minimal solutions are proportional to each other, we recover an approximation of c n as c n ≈ (c 0 /u 0 ) u n .

We use Miller's method to evaluate the minimal solution (c n ) of the normalized recurrence [START_REF] Robbins | A remark on Stirling's formula[END_REF], and we get an approximation of G(x) = ∑ n≥0 G n x 3n as S N = (c 0 /u 0 ) ∑ N-1 n=0 u n x 3n /(n! 2 ). The algorithm is summed up as Algorithm 1. The rest of this section is devoted to its proof of correctness, i.e., the proof that the value s it returns satisfies |G(x) -s| ≤ 3 • 2 -p G(x).

Proposition 5: With N as in Algorithm 1, the truncation error satisfies | ∑ ∞ n=N G n x 3n | ≤ 2 -p G(x) for all x ≥ 0.5. Proof: First, because of line 2 of the algorithm, we have

3 10 x 3 ≤ (N + 1) 2 , so that ⃒ ⃒ ⃒ ⃒ ∞ ∑ n=N G n x 3n ⃒ ⃒ ⃒ ⃒ ≤ 2G N x 3N ≤ 2 ( ex 3/2 3N ) 2N ≤ 2 (αN) -2N
by Lemma 4 and Corollary 2. Line 3 then ensures that

2 (αN) -2N ≤ 2 -p 128 e 2 3 x 3/2 x -3/4 ≤ 2 -p G(x),
the last inequality coming from Lemma 5.

There are two sources of error besides the truncation: first, (u n ) is not exactly proportional to (c n ), especially when n is close to R. Second, roundoff errors happen during the evaluation of (u n ). Rigorous bounds for both sources of error have been proposed by Mattheij and van der Sluis [START_REF] Mattheij | Error Estimates for Miller's Algorithm[END_REF]. We combine them with classical techniques (well-explained, e.g., in [START_REF] Chevillard | The functions erf and erfc computed with arbitrary precision and explicit error bounds[END_REF]) to choose the starting index R and the working precision t so as to guarantee the final accuracy.

We now recall Mattheij and van der Sluis' main result (adapted to our particular case, which simplifies the statement quite a bit). Consider a recurrence of the form (23). Denote by (c n ) a minimal solution that we wish to evaluate, and let (d n ) be the solution such that d 0 = d 1 = 1. Assume that (d n ) is a dominant solution and that the sequences (c n ), (d n ) and

(c n /d n ) are decreasing. Define H = d 0 c 0 ∑ R-1 i=0 c i d i and, for i ≤ R, U i = ( c i c i c i+1 c i d i+1 d i )
and

B i = ( -a 1 (i) a 0 (i) -a 2 (i) a 0 (i) 1 0
) .

Let v R ∈ R 2 be a column vector, and for i ≤ R -1, let v i be the result of the floating-point evaluation of B i v i+1 at precision t.

Write v i = (u i , u i+1 ) T . If all operations were exact, (u i ) would be the solution of the recurrence such that (u R , u R+1 ) T = v R . To take rounding errors into account, we write

v i = (B i +2 -t 𝒢 i )v i+1 for some matrix 𝒢 i instead. Define y R = (y R1 , y R2 ) T = U -1 R v R . Let ℱ i = ‖U -1
i 𝒢 i U i+1 ‖, the matrix norm being subordinate to the ‖ • ‖ ∞ norm for vectors, and let ℱ ≥ max i ℱ i .

Theorem 1: [10, Theorem 4.1] Provided that the quantities

ℱ R2 -t , c R d 0 d R c 0 ⃒ ⃒ ⃒ ⃒ ⃒ y R2 y R1 ⃒ ⃒ ⃒ ⃒ ⃒ , ‖y R ‖ ∞ |y R1 | (R + H)(1.3ℱ 2 -t )
are all bounded by 0.1, the approximate value u i computed by Miller's algorithm satisfies (c 0 /u 0 )

u i = c i (1 + θ i ) for some θ i such that |θ i | ≤ T i + R i ,
where

T i = 1.5 ( c R d i c i d R + c R d 0 c 0 d R ) ⃒ ⃒ ⃒ ⃒ ⃒ y R2 y R1 ⃒ ⃒ ⃒ ⃒ ⃒ , R i = 1.5 ‖y R ‖ ∞ y R1 ε (i + 2H), and ε = 1.3ℱ 2 -t .
Turning back to the special case c n = n! 2 G n , Theorem 1 applied to (15) yields the following. Recall that τ = 3/20.

Corollary 3: Set v R = (1, 0) T and

B i = ( 10 -r(i) 1 0
) where r(n) = (3n + 4)(3n + 5) (n + 1)(n + 2) .

Then, in the notation of Theorem 1, we have T i ≤ τ R-i and R i ≤ 76.5(i + 4)2 -t ≤ 76.5(N + 3)2 -t for all i < N.

Proving this corollary still requires some work. We postpone it for a bit to explore the consequences of this statement.

Observe that lines 8 and 9 of Algorithm 1 are equivalent to a floating-point evaluation of B i v i+1 where v i+1 = (a, b) T . Hence, at each loop turn, we have a = u i just after line 8. Lines 10 and 11 are a Horner-like evaluation scheme, so that s ′ ≈ ∑ N-1 i=0 u i x 3i /i! 2 at the end of the loop. More precisely, assuming x 3 is approximated by ∘(x) ⊗ ∘(x) ⊗ ∘(x) and division by (i + 1) 2 is performed as two successive divisions by (i + 1), an easy induction shows that one can write

s ′ = N-1 ∑ i=0 ( u i x 3i i! 2
) ⟨9(i + 1)⟩ .

Line 12 adds 3 to all error counters. The choice of t on line 5 ensures that (9(N + 1) + 3) • 2 -t ≤ 1/2. Using Prop. 4, we conclude that the sum s returned by Algorithm 1 satisfies

s = N-1 ∑ i=0 c 0 u i u 0 • x 3i i! 2 (1 + µ i ) = N-1 ∑ i=0 G i x 3i (1 + µ i )(1 + θ i ),
where

|µ i | ≤ 2 (9(i+1)+3)•2 -t ≤ 18(N +3) 2 -t and |θ i | ≤ T i +R i .
Since R ≥ N, we have T i ≤ τ for all i < N by Corollary 3. The choice of t also implies R i ≤ 76.5/256. Altogether, this ensures that

|θ i | ≤ 1. Writing (1 + µ i )(1 + θ i ) = 1 + δ i , we get |δ i | ≤ 2|µ i | + |θ i | ≤ 112.5(N + 3) 2 -t + τ R-i ≤ 2 -p + τ R-i , and therefore ⃒ ⃒ ⃒s - ∑ N-1 i=0 G i x 3i ⃒ ⃒ ⃒ ≤ 2 -p G(x) + τ R G(τ -1/3 x). Lemma 6: For x > 0.5, we have τ R G(τ -1/3 x) ≤ 2 -p G(x).
Proof: It follows from Lemma 5 (and τ < 1) that

G(x) G(τ -1/3 x) ≥ 1 4 exp ( 2 3 x 3/2 (1 -τ -1/2 )
) , and the algorithm ensures that

R log 2 (τ -1 ) ≥ p + 2 + 2 3 x 3/2 (τ -1/2 -1) log 2 (e), whence τ R ≤ 1 4 2 -p exp( 2 3 x 3/2 (1 -τ -1/2
)) and the result. We can now prove the correctness of Algorithm 1: Theorem 2: The value s returned by Algorithm 1 satisfies

|G(x) -s| ≤ 3 • 2 -p G(x).
Proof: It follows from Prop. 5 and the above discussion, since

|G(x) -s| ≤ |G(x) - ∑ N-1 i=0 G i x 3i | + |s - ∑ N-1 i=0 G i x 3i
|. The remainder of this section is devoted to the proof of Corollary 3. We begin with a crucial lemma. Let (d n ) be the solution of (15) defined by d 0 = d 1 = 1, let η(n) = 1/(3n 2 ), and let r(n) be as in Corollary 3.

Lemma 7: For all n ≥ 1, we have

d n+1 ≤ d n ≤ (1+η(n)) d n+1 .
Proof: We proceed by induction. Since d 2 = 9/10, the property is true for n = 1. Now, supposing it for an arbitrary n, we get (9

-η(n))d n+1 ≤ 10d n+1 -d n ≤ 9d n+1 , so 9 -η(n) r(n) d n+1 ≤ d n+2 ≤ 9 r(n) d n+1 .
We conclude by observing that 9/r(n) ≤ 1 and r(n)/(9-η(n)) ≤ 1 + η(n + 1) for n ≥ 1.

Corollary 4: For all n, we have 0.783 ≤ d n ≤ 1.

Proof: By Lemma 7, (d n ) is decreasing and d 0 = 1: this proves the right-hand side. For n ≥ 100, Lemma 7 implies

d 100 ≤ d n n-1 ∏ i=100 (1 + η(i)) ≤ d n p ∞ p 99 , p k = k ∏ i=1 (1 + η(i)).
Using the exact value

p ∞ = √ 3 π sinh( π √
3 ) [1, Eq. 4.36.1], we check that d n ≥ d 100 p 99 p -1 ∞ ≥ 0.783 for n ≥ 100. As (d n ) is decreasing, the inequality holds for n < 100 too.

This estimate, combined with Corollary 1 gives almost all we need to check the hypotheses of Theorem 1. We use the notation introduced for the statement of the theorem (specialized to the computation of c n = n! 2 G n using [START_REF] Robbins | A remark on Stirling's formula[END_REF], with (d n ) as above).

H = d 0 c 0 R-1 ∑ i=0 c i d i ≤ 1 0.783 R-1 ∑ i=0 τ i ≤ 2.
We have d 0 /d 1 = 1, d 1 /d 2 = 10/9, and d i /d i+1 ≤ 1+1/(3i 2 ) ≤ 10/9 for i ≥ 2, and, by definition of U R and v R ,

y R = U -1 R v R = c R det U R ( d R+1 /d R -1 -c R+1 /c R 1 ) ( 1 0 ) = c R det U R ( d R+1 /d R -c R+1 /c R ) .
It follows that |y R2 /y R1 | ≤ 10 9 τ = 1 6 , and hence ‖y R ‖ ∞ = |y R1 |. Finally, from the expression det

U i = c 2 i (d i+1 /d i -c i+1 /c i ), we obtain c -2 i det U i ≥ ( 9 10 -τ) =
|𝒢 i | ≤ ( 4 • 10 10 r(i) 0 0 ) ,
where, for matrices, |A| ≤ |B| means that the inequality holds entrywise. Since moreover

|U -1 i | ≤ |c i | | det U i | ( d i+1 /d i 1 c i+1 /c i 1 ) ≤ 4 3c i ( 1 1 3/20 1 ) , |U i+1 | ≤ c i+1 ( 1 1 c i+2 /c i+1 d i+2 /d i+1 ) ≤ c i+1 ( 1 1 3/20 1 ) ,
and

ℱ i ≤ |U -1 i | • |𝒢 i | • |U i+1 |, we get ℱ i ≤ 4 3 • 3 20
( 40 + 3 2 r(i) 40 + 10r(i) 3 20 (40 + 3 2 r(i)) 3 20 (40 + 10r(i))
) .

Hence ‖ℱ i ‖ ≤ 16 + 23 10 r(i). The result follows because r(i) ≤ 10 for all i ≥ 0.

We can finally prove Corollary 3.

Proof: To apply Theorem 1, we need to check that

(i) 39 R2 -t ≤ 0.1, (ii) c R d 0 d R c 0 ≤ 0.1 ⃒ ⃒ ⃒ ⃒ ⃒ y R1 y R2 ⃒ ⃒ ⃒ ⃒ ⃒ , (iii) (R + 2)(50.7 • 2 -t ) ≤ 0.1.
By definition of t, we have R + 2 ≤ 2 t-9 , so (iii) is satisfied, and (i) follows immediately. Corollary 5 combined with the inequalities d 0 /d R ≤ 1/0.783 and c R /c 0 ≤ τ R ≤ τ implies (ii).

In conclusion, we can apply Theorem 1, hence

T i = 1.5 ( c R d i c i d R + c R d 0 c 0 d R ) ⃒ ⃒ ⃒ ⃒ ⃒ y R2 y R1 ⃒ ⃒ ⃒ ⃒ ⃒ ≤ 1.5 0.783 ( τ R-i + τ R ) 1 6 . Therefore, T i ≤ 0.639 • τ R-i ≤ τ R-i as announced. Corollary 5 yields R i = 1.5 (50.7 • 2 -t )(i + 2H) ≤ 76.5 • (i + 4) 2 t .

VI. Evaluation of the auxiliary series

The implementation of the auxiliary series F is much easier than that of G. We limit ourselves to a sketch of the (fairly standard) algorithm.

A variable a 0 is used to successively evaluate F 0 , F 3 x 3 , F 6 x 6 , etc., using the recurrence [START_REF] Olver | Asymptotics and special functions[END_REF]. Accordingly, two variables a 1 and a 2 are used to evaluate the successive values of F 3i+1 x 3i+1 and F 3i+2 x 3i+2 . Each step adds at most 10 to the relative error counter of each variable. A variable s is used to accumulate the sum as the variables a k are updated. Therefore, after step K, we can write s = ∑ 3K-1 i=0 F i x i ⟨1 + 10⌊i/3⌋ + 3K -i⟩ (the term 3Ki representing the errors due to additions). Bounding all errors uniformly, we get

s = 3K-1 ∑ i=0 F i x i ⟨10K⟩ .
It is easy to see that q(i) = F i+3 /F i decreases for i ≥ 1, hence the loop can be stopped as soon as (i) q(3K)x 3 < 1/2, and (ii) a 0 , a 1 , a 2 < 2 Exp(s)-p-4 . These conditions ensure that the remainder ∑ i≥3K F i x i is bounded by 2(F 3K + F 3K+1 + F 3K+2 ) < 4(a 0 + a 1 + a 2 ) < 12 16 2 -p 2 Exp(s) . It is clear from Prop. 1 that F n+3 ≤ 4F n /n 2 for any n, hence (using n! ≈ (n/e) n ) we have F n ≈ (4e 2 /n 2 ) n/3 . This is most likely an overestimation of the true value. Moreover, we can approximate F(x) by 1 32 x -1/2 exp( 4 3 x 3/2 ) for x > 0.5. These estimates are used to get a rough overestimation of K. The working precision t is then chosen so that 20K • 2 -t ≤ 2 -3-p . Hence, we have | ∑ 3K-1 i=0 F i x i -s| ≤ 2 -3-p s ≤ 2 Exp(s)-3-p . The initial estimation of the truncation rank is very unlikely to be underestimated. In the case it would be smaller than the actual truncation rank decided on-the-fly by the above criterion, this is checked a posteriori and, if necessary, the evaluation is run again with an updated working precision.

We conclude by writing

|F(x) -s| ≤ |F(x) - ∑ 3K-1 i=0 F i x i | + |s - ∑ 3K-1 i=0 F i x i | ≤
| = |δ 3 + δ 4 (1 + δ 3 )| ≤ 7 • 2 -p ≤ 2 -(p-3) .
We developed a prototype implementation of this algorithm, based on the multiple-precision floating-point library MPFR [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF]. For simplicity, we supposed x ≥ 1/2 in the present paper, but our implementation is valid for any x ≥ 0.

We compared the results with those of the implementation of Ai(x) available in MPFR, run with a larger precision. Random tests using thousands of points x and target precisions p yield relative errors smaller than 2 -(p-3) between the result of our implementation and the result of MPFR, as predicted by theory.

We also ran our implementation and MPFR with the same accuracy to compare their performance (see Fig. 3). When x is large, our method is faster than MPFR, which uses the Taylor expansion of Ai at the origin. This is all the benefit of reducing the cancellation: for large x, the working precision of MPFR has a large overhead because of the bad condition number of the series. In contrast, when p is large compared to x, our implementation pays the cost of evaluating two series instead of one, with little benefit in terms of working precision.

To be completely fair, we should mention that MPFR does not implement the asymptotic expansion of Ai at infinity, which should be a better choice than our algorithm for large x and comparatively small p. On the other hand, our code currently uses only the naive series summation algorithm, while MPFR implements Smith's baby steps-giant steps technique [START_REF] Paterson | On the number of nonscalar multiplications necessary to evaluate polynomials[END_REF], [START_REF] Smith | Efficient multiple-precision evaluation of elementary functions[END_REF] that is more efficient. There is no theoretical obstacle to using Smith's method with our series: it only obfuscates a little bit the description of the algorithm and the roundoff error analysis. Once we will have implemented it, we will have a more complete picture with three areas, indicating what method between the Taylor series, the asymptotic series and our series is the most efficient, depending on (x, p).

VIII. Outlook: Towards a GMR algorithm?

As mentioned in the introduction, we see the present work as a case study. Indeed, in spite of the many technical details that occupy much of the space of this article, we really used few specific properties of the Airy function besides Equation (3). Looking back, our analysis essentially relies on the following ingredients.

(i) The ability to find auxiliary series. The indicator functions used in the GMR method depend only on the behaviour at (complex) infinity of the entire functions they are associated to. In the case of the solution of a LODE with analytic coefficients, this behaviour is entirely determined by the differential equation along with a finite number of "asymptotic initial values" [START_REF] Wasow | Asymptotic expansions for ordinary differential equations[END_REF], [START_REF] Van Der Put | Galois theory of linear differential equations[END_REF]. Once the indicator function is known, it remains to exhibit appropriate auxiliary series. Doing this in a truly general way remains an open problem. Yet, both the original GMR method and our variant apply to many cases, and it is likely that they can be combined and further generalized.

(ii) An efficient way to compute their coefficients. This seems to be considered a major limitation in the original GMR paper [START_REF] Gawronski | Reduced cancellation in the evaluation of entire functions and applications to the error function[END_REF]. But, as already mentioned, recurrences with polynomial coefficients automatically exist as soon as both the original function to evaluate and the auxiliary series satisfy differential equations with polynomial coefficients. Numerical stability is not much of an issue in the case of three-term recurrences, thanks to Miller's method, though many technical details must be settled. The situation is more complicated in general for recurrences of higher order. Observe, though, that we proved the minimality of (G n ) using essentially the same asymptotic properties that were exploited by the GMR method in the first place. The minimality may hence not be fortuitous and might generalize.

(iii) Upper and lower bounds on the coefficients and sums of the series F and G. All these bounds were derived, in a pretty systematic way, from the asymptotic expansion of Ai at infinity combined with Lemma 1. Bounds similar to that from Lemma 1 can themselves often be obtained from a LODE [START_REF] Olver | Asymptotics and special functions[END_REF].

(iv) Roundoff error analyses. Our error analyses follow a very regular pattern and could probably be abstracted to a more general case or automated.

In short, most steps of the present study could apparently be performed in a systematic way, starting from Eq. (3) plus a moderate amount of additional information. Systematizing the GMR method based on this observation seems a promising line of research. Solutions of LODE with polynomial coefficients are known in Computer Algebra as D-finite, or holonomic, functions. It would be interesting to isolate a subclass of D-finite functions to which the method applies in a truly systematic way, and attempt to automate it at least partially.

Figure 1 .

 1 Figure 1. The indicator function h of Ai.

  n≥0 |y n z n |)log 2 |y(z)|. Denote M(r) = sup |z|=r |y(z)| for all r > 0. Cauchy's formula implies max n (|y n |r n ) ≤ M(r), and under a suitable version of hypothesis (4), one can actually show that max n (|y n |r n ) ≈ M(r).

Figure 2 .

 2 Figure 2. Indicator functions for F (left) and G (right).

Lemma 3 : 8 ∫ 8 ∫ π/ 3 r 3 r 3 / 8

 3883338 The expression I(r) = r 9/satisfies 0 ≤ I(r) ≤ 0.51 for all r ≥ 10. Proof: We first use the inequality cos( 3 2 θ) ≤ 1 -9 10 θ 2 (valid for 0 ≤ θ ≤ π/3) followed by the change of variable ϕ = √ 3/5 r 3/4 θ to get I(r) ≤ r 9/-5/8 e -3 5 r 3/2 θ 2 dθ = √ 5 √ ∫ b(r) a(r) e -ϕ 2 dϕ, where a(r) = √ 3/5 r 1/8 and b(r) = π √ 15 r 3/4 . Now we have ∫ b(r)

Corollary 5 :

 5 The sequences (c n ), (d n ) and (c n /d n ) are decreasing. Moreover, the following inequalities hold:H ≤ 2, |y R2 /y R1 | ≤ 1 6 , ‖y R ‖ ∞ /|y R1 | = 1, and det U i ≥ 3 4 c 2 i . Proof: Corollary 1 shows that (c n ) is decreasing and Lemma 7 shows that (d n ) is decreasing. Together, they imply c n+1 d n+1 / c n d n ≤ τ(1 + η(n)) ≤ 1 5 for n ≥ 1. We check separately that c 1 /d 1 ≤ c 0 /d 0 . Hence (c n /d n ) is decreasing.Corollary 1 also implies c n ≤ τ n c 0 for any n, and Corollary 4 shows that d 0 = 1 ≤ d n /0.783 for any n. It follows that

3 4 .

 4 Lemma 8: A suitable value for ℱ is 39. Proof: The multiplication B i v i is performed on lines 8 and 9 of the algorithm. Denoting by a ′ and b ′ the new values of a and b, we can write b ′ = a anda ′ = 10a ⟨2⟩r(i)b ⟨5⟩ = 10a(1 + θ)r(i)b(1 + θ ′ )where (since t ≥ 5 by line 5) |θ| ≤ 4 • 2 -t and |θ ′ | ≤ 10 • 2 -t by Prop. 4. (This assumes that the multiplication by r(i) is performed through four successive multiplications and divisions). Therefore, we have

Figure 3 .

 3 Figure 3. Fastest method as a function of x and p (scales are logarithmic). MPFR 3.1.1, GMP 4.3.1, on an Intel Xeon at 2.67GHz.

  -t ( ∑ n≥0 |y n z n |)/|y(z)|. The sum ∑ n≥0 |y n z n | is larger than max n≥0 |y n z n |, and usually of the same order of magnitude.

	Therefore, the number of significant binary digits "lost by
	cancellation" is roughly
	log 2 (max

Acknowledgments: We thank Paul Zimmermann for pointing out the GMR article to us.