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

Guillaume Hanrot
email: guillaume.hanrot@ens-lyon.fr

Olivier Robert
email: olivier.robert@univ-st-etienne.fr

• O Robert

Exponential sums and correctly-rounded functions

The 2008 revision of the IEEE-754 standard, which governs floating-point arithmetic, recommends that a certain set of elementary functions should be correctly rounded. Successful attempts for solving the Table Maker's Dilemma in binary64 made it possible to design CRlibm, a library which offers correctly rounded evaluation in binary64 of some functions of the usual libm. It evaluates functions using a two step strategy, which relies on a folklore heuristic that is well spread in the community of mathematical functions designers. Under this heuristic, one can compute the distribution of the lengths of runs of zeros/ones after the rounding bit of the value of the function at a given floating-point number. The goal of this paper is to change, whenever possible, this heuristic into a rigorous statement. The underlying mathematical problem amounts to counting integer points in the neighborhood of a curve, which we tackle using so-called exponential sums techniques, a tool from analytic number theory.

INTRODUCTION

On most current computer systems, real numbers are approximated and represented by floating-point numbers [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]. For many years, floating-point arithmetic has been a mere set of cooking recipes, a situation described for instance in [START_REF] Kahan | Why do we need a floating-point standard?[END_REF], and numerical programs were not reliable nor portable.

The IEEE-754 [START_REF] Cody | A proposed radix-and-word-length-independent standard for floating-point arithmetic[END_REF], [START_REF]IEEE Standard for Binary Floating-Point Arithmetic[END_REF] standard, and its revision [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Standard[END_REF], [START_REF] Markstein | The New IEEE-754 Standard for Floating-Point Arithmetic[END_REF], for binary floating-point arithmetic (and the radix independent IEEE-854 [START_REF] Cody | A proposed radix and word length independent standard for floating-point arithmetic[END_REF], [START_REF]IEEE Standard for Radix Independent Floating-Point Arithmetic[END_REF] standard that followed) drastically improved the situation and put an end to this dangerous era. In particular, the IEEE-754 standard clearly specifies the formats of the floating-point representations of numbers, and the behavior of the four arithmetic operations and the square root.

And yet, as of today, the standard still does not rule the behavior of usual functions, such as the ones contained in the C mathematical library (libm), as precisely as it does for the four arithmetic operations and the square root. Part of a general effort to improve this situation [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF], this paper presents a novel rigorous approach for estimating the amount of so-called bad cases for rounding. As we shall see below, it is useful to design function evaluation routines with a higher quality of computation called correct rounding, and in particular, it plays a key role in order to estimate their performances.

Our arithmetic framework

Let 𝒟 𝑝 denote the set of exponent-unbounded, 𝑝-bit significand, radix-2 floating-point numbers (with 𝑝 1):

𝒟 𝑝 = { 𝑗 • 2 𝐸-𝑝+1 , 2 𝑝-1 |𝑗| 2 𝑝 -1, 𝑗, 𝐸 ∈ Z } ∪ {0}.
This set is not the set of available floating-point numbers on an existing system. It is an "ideal" system of floating-point numbers: a superset of an actual system, with no overflows, underflows, or subnormals [START_REF]IEEE Standard for Floating-Point Arithmetic. IEEE Standard[END_REF], [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF]. We prove results in 𝒟 𝑝 . These results remain true in an actual system, provided that no overflows, underflows or subnormals occur.

We call significand of a nonzero element 𝑥 = 𝑗 • 2 𝐸-𝑝+1 of 𝒟 𝑝 the number 𝑗/2 𝑝-1 and integral significand of 𝑥 the integer 𝑗. The exponent of an element 𝑗 •2 𝐸-𝑝+1 , with 2 𝑝-1 |𝑗| 2 𝑝 -1, is the integer 𝐸.

The result of an arithmetic operation whose input values belong to 𝒟 𝑝 may not belong to 𝒟 𝑝 (in general it does not). Hence that result must be rounded. The IEEE standard defines 4 different rounding modes; in the sequel, 𝑥 is any real number to be rounded:

• rounding towards +∞, or upwards: ∘ 𝑢 (𝑥) is the smallest element of 𝒟 𝑝 that is greater than or equal to 𝑥;

• rounding towards -∞, or downwards: ∘ 𝑑 (𝑥) is the largest element of 𝒟 𝑝 that is less than or equal to 𝑥;

• rounding towards 0: ∘ 𝑧 (𝑥) is equal to ∘ 𝑢 (𝑥) if 𝑥 < 0, and to ∘ 𝑑 (𝑥) otherwise;

• rounding to the nearest even: ∘ 𝑛 (𝑥) is the element of 𝒟 𝑝 that is closest to 𝑥. If 𝑥 is exactly halfway between two consecutive elements of 𝒟 𝑝 , ∘ 𝑛 (𝑥) is the one for which the integral significand 𝑗 is an even number.

The first three rounding modes are called directed rounding modes.

Correct rounding

The standard requires that the user should be able to choose one rounding mode among these ones, called the active rounding mode. An active rounding mode being chosen, when performing one of the 4 arithmetic operations, or when computing square roots, the obtained rounded result should be equal to the rounding of the exact result: this requirement on the quality of the computation is called correct rounding.

While the IEEE 754-1985 and 854-1987 standards required correctly rounded arithmetic operations, they did not do it for the most common mathematical functions, such as simple algebraic 1 functions like 1/ √ •, 3 √ •, . . . and also a few transcendental 2 functions like sine, cosine, exponentials, and logarithms of radices 𝑒, 2, and 10, etc. More generally, a natural target is the whole class of elementary functions 3 . A subset of these functions is usually available from the libms delivered with compilers.

This lack of requirement is mainly due to a difficult problem known as the Table Maker's Dilemma (TMD), a term coined by Kahan, which we present below.

Being able to provide correctly rounded functions is of utter interest:

• it greatly improves the portability of software;

• it allows one to design algorithms that use this requirement;

• this requirement can be used for designing formal proofs of pieces of software;

• one can easily implement interval arithmetic, or more generally one can get certain lower or upper bounds on the exact result of a sequence of arithmetic operations.

Implementation of correctly-rounded functions: Ziv's strategy

A. Ziv proposed a general methodology [START_REF] Ziv | Fast evaluation of elementary mathematical functions with correctly rounded last bit[END_REF], implemented in the libultim library 4 , that made possible the correctly rounded evaluation of the functions of the libm. It consists in iteratively improving the accuracy of the approximation until the correctly rounded value can be decided.

However, there was, up until recently, no practical bound on the termination time of Ziv's iteration: it may be proven to terminate for some transcendental functions, but the actual maximal accuracy required in the worst case was unknown. Note that in order to prove that Ziv's iteration terminates for a given real-valued function 𝜙 and a given floating-point number 𝑥, one has to guarantee a priori that 𝜙(𝑥) ̸ ∈ 𝒟 𝑝 (directed rounding modes) or 𝜙(𝑥) is not the middle of two consecutive elements of 𝒟 𝑝 (rounding to nearest mode). One has to deal separately with the exact cases such as exp(0), log(1), sin(0) or 3

√

𝑥 3 for 𝑥 ∈ 𝒟 ⌊𝑝/3⌋ etc. [START_REF] Lauter | Arrondi Correct de Fonctions Mathématiques[END_REF], [START_REF] Lauter | An efficient rounding boundary test for pow(𝑥, 𝑦) in double precision[END_REF], [START_REF] Muller | Elementary Functions, Algorithms and Implementation[END_REF].

Unfortunately, for many applications, this is not satisfactory, for several reasons:

• in libultim, the measured worst-case execution time is indeed three orders of magnitude larger than that of the usual libms;

1. We say that a function 𝜙 is algebraic if there exists 𝑃 ∈ Z[𝑥, 𝑦] ∖ {0} such that 𝑃 (𝑥, 𝜙(𝑥)) = 0.

2. A function is transcendental if it is not algebraic.

3. An elementary function is a function of one variable which is the composition of a finite number of arithmetic operations (+, -, ×, /), exponentials, logarithms, constants, and solutions of algebraic equations.

4. libultim was released by IBM. An updated version is now part of the GNU glibc and available under the GNU General Public License.

• A related problem is memory requirement, which is, for the same reason, unbounded in theory and much higher than the usual libms in practice;

• when the approximations are evaluated using a fast pipelined multiplier or multiplier-accumulator, the tests required by the previous strategy would require to wait many cycles before being able to restart the computation;

• for real-time applications, the delay of computation must be bounded.

Fast and cheap correctly-rounded function evaluation in binary64

For most of elementary functions, when evaluating them, one has to compute some approximation to the exact result, with an accuracy (which is then a rational integer power of the radix 2) somewhat higher than the "target" precision 𝑝. The Table Maker's Dilemma is the problem of determining what the accuracy of this approximation should be to make sure that rounding that approximation will always be equivalent to rounding the exact result. Ideally, we aim at getting the minimal possible accuracy.

On the one hand, there are theoretical results that yieldnot fully satisfactory -solutions for algebraic functions: the precision which the computations must be performed to is, in general, overestimated [START_REF] Iordache | On infinitely precise rounding for division, square root, reciprocal and square root reciprocal[END_REF], [START_REF] Lang | Bound on run of zeros and ones for algebraic functions[END_REF], [START_REF] Brisebarre | Correct rounding of algebraic functions[END_REF]. On the other hand, regarding transcendental functions, either no theoretical statement exists or they provide results that cannot be used in practical computations [START_REF] Nesterenko | On the approximation of the values of exponential function and logarithm by algebraic numbers[END_REF].

Therefore, algorithmic approaches to the TMD [START_REF] Lefèvre | Moyens Arithmétiques Pour un Calcul Fiable[END_REF], [START_REF] Lefèvre | Worst cases for correct rounding of the elementary functions in double precision[END_REF], [START_REF] Stehlé | Searching Worst Cases of a One-Variable Function Using Lattice Reduction[END_REF] had to be developed. They allowed for solving the TMD for the IEEE binary64 (also known as "double precision"). As a consequence, the revised IEEE-754 standard now recommends (yet does not require, due to the lack of results in the case of binary128) that the following functions should be correctly rounded: 𝑒 𝑥 , 𝑒 𝑥 -1, 2 𝑥 , 2 𝑥 -1, 10 𝑥 , 10 𝑥 -1, ln(𝑥), log 2 (𝑥), log 10 (𝑥), ln(1 + 𝑥), log 2 (1 + 𝑥), log 10 (1+𝑥), √ 𝑥 2 + 𝑦 2 , 1/ √ 𝑥, (1+𝑥) 𝑛 , 𝑥 𝑛 , 𝑥 1/𝑛 (𝑛 is an integer), sin(𝜋𝑥), cos(𝜋𝑥), arctan(𝑥)/𝜋, arctan(𝑦/𝑥)/𝜋, sin(𝑥), cos(𝑥), tan(𝑥), arcsin(𝑥), arccos(𝑥), arctan(𝑥), arctan(𝑦/𝑥), sinh(𝑥), cosh(𝑥), tanh(𝑥), sinh -1 (𝑥), cosh -1 (𝑥), tanh -1 (𝑥).

Thanks to these results, it is now possible to obtain correct rounding in binary64 in two Ziv steps only, which one may then optimize separately. This is the approach used in CRlibm 5 , a library which offers correctly rounded evaluation of the double-precision C99 standard elementary functions:

• the first quick step is as fast as a current libm, and provides an accuracy of 2 -52-𝑘 (𝑘 = 11 for the exponential function for instance), which is sufficient to round correctly to the 53 bits of binary64 in most cases;

• the second accurate step is dedicated to challenging cases. It is slower but has a reasonable bounded execution time, being tightly targeted at the hardestto-round cases computed by Lefèvre et al. [START_REF] Lefèvre | Towards correctly rounded transcendentals[END_REF], [START_REF] Lefèvre | Worst cases for correct rounding of the elementary functions in double precision[END_REF], [START_REF] Stehlé | Worst Cases and Lattice Reduction[END_REF], [START_REF] Stehlé | Searching Worst Cases of a One-Variable Function Using Lattice Reduction[END_REF], [START_REF] Stehlé | On the Randomness of Bits Generated by Sufficiently Smooth Functions[END_REF]. In particular, there is no need for arbitrary multiple precision anymore.

5. http://lipforge.ens-lyon.fr/www/crlibm/ This approach [START_REF] De Dinechin | Towards the postultimate libm[END_REF], [START_REF] De Dinechin | Fast and correctly rounded logarithms in double-precision[END_REF] leads to correctly-rounded function evaluation routines that are fast and have a reasonable memory consumption.

A heuristic probabilistic approach

This strategy relies on the following heuristic assumption about the first step: the proportion of "bad cases", that is to say the floating-point numbers whose evaluation cannot be correctly rounded if we use 2 -52-𝑘 as an intermediate accuracy, is around 2 1-𝑘 . One can find such an estimate used in [START_REF] Dunham | Feasibility of "perfect" function evaluation[END_REF] and a probabilistic study has been done in [START_REF] Gal | An accurate elementary mathematical library for the IEEE floating point standard[END_REF]. We now briefly recall the intuition behind this heuristic (see [START_REF] Muller | Elementary Functions, Algorithms and Implementation[END_REF], [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF] for a more detailed presentation).

Let 𝜙 be a real-valued function, assume that after the 𝑝 th bit, the bits of the significands of the values 𝜙(𝑥), where 𝑥 is a floating-point number, are sequences of independent random 0 or 1 with equal probability 1/2. The probability that after bit 𝑝, we have • in rounding to nearest mode, the bit sequence

100 • • • 0 ⏟ ⏞ 𝑘 bits or 011 • • • 1 ⏟ ⏞ 𝑘 bits •
or, in directed rounding modes, the bit sequence

00 • • • 0 ⏟ ⏞ 𝑘 bits or 11 • • • 1 ⏟ ⏞ 𝑘 bits is 2 -𝑘+1 .
Hence, if we have 𝑁 floating-point numbers in the domain being considered, the number of values 𝑥 for which we will have a bit sequence of the form indicated above is, under the probabilistic model stated above, around 𝑁 2 -𝑘+1 .

Goal and outline of the paper

The assumption on which this approach relies is a very strong one, and, actually, is wrong in full generality (one can think, for instance, of the exponential function around 0). And yet, it seems to be satisfied in most cases and as of today, proving it seems completely out of reach. The goal of this paper is to give solid theoretical foundations to this heuristic probability for certain ranges of values of 𝑘. We target in particular the values of 𝑘 that CRlibm uses in practice. For instance, our results prove in particular that:

• in binary64, for 𝑘 = 11, the number of bad cases for for the occurrence of the second accurate step.

-3 √ • over [1/2,
We shall start with a mathematical formalization of the problem we address. This will be done in Section 2. In Section 3, we shall present estimates, proved thanks to the so-called exponential sums techniques, a mathematical tool at the core of our study. These objects and the related theory are powerful tools mainly used in Analytic Number Theory [START_REF] Graham | Van der Corput's method of exponential sums[END_REF], [START_REF] Huxley | Area, lattice points, and exponential sums[END_REF], [START_REF] Montgomery | Ten lectures on the interface between analytic number theory and harmonic analysis[END_REF], [START_REF] Tenenbaum | Introduction to analytic and probabilistic number theory[END_REF], [START_REF] Krätzel | Lattice points[END_REF], [START_REF] Robert | On van der Corput's 𝑘-th derivative test for exponential sums[END_REF]. Then, in Section 4, we shall apply the results of Section 3 to elementary functions. This will enable us to estimate the number of occurrences of the second accurate step of the two Ziv step strategy, like the one used in the CRlibm design. We shall also analyze the constraints imposed by our approach (in particular, the ranges the parameter 𝑘 should belong to) and present in more detail the cases of the 3 √ • and exp functions. We shall report experiments with these two functions in binary64 and binary128 in order to evaluate the quality of our estimates.

Up to a renormalization (which associates a function 𝑓 to the function 𝜙, see ( 18) and ( 19)), our problem can be seen as a question about the number of integers 𝑚 in some interval of the form [1, 𝑀 ] such that the distance of 𝑓 (𝑚) to the integers is less than some 𝛿 > 0. This problem, which amounts to counting integer points (𝑚, 𝑛) in the 𝛿-neighborhood of the curve 𝑦 = 𝑓 (𝑥), will be treated in Section 3, which is the core of the paper. The main term of this number of points is expected to be 2𝑀 𝛿. We shall express the error term as a sum of periodic functions (namely fractional parts), and control this error term using an optimal truncation of the Fourier series introduced by Vaaler (Theorem 1). The problem then reduces to estimating exponential sums, for which we shall use van der Corput's inequality (Theorem 2).

Note that, actually, our results are valid for any 𝒞 2 function 6 , and not only for elementary functions.

Our work may be considered as complementary to the algorithmic determinations of the worst cases for correct rounding [START_REF] Lefèvre | Moyens Arithmétiques Pour un Calcul Fiable[END_REF], [START_REF] Lefèvre | Worst cases for correct rounding of the elementary functions in double precision[END_REF], [START_REF] Stehlé | Searching Worst Cases of a One-Variable Function Using Lattice Reduction[END_REF] in two ways:

• while their works allow for a worst case analysis of an elementary function implementation, our paper provides an attempt for a rigorous average case analysis of such an implementation. We shall come back to this in Remark 2.

• as the reader will notice in Subsection 4.1, our approach, roughly, works for values of 𝑘 up to 𝑝/3 bits. In this range, the aforementioned algorithmic approaches may prove pointless since that there are too many values to determine, whereas our approach makes it possible to estimate the amount of bad cases for rounding in a satisfactory way.

FORMALIZATION OF THE PROBLEM

Assume we wish to correctly round an elementary function 𝜙 (actually, our formalization is valid for any real-valued function). We consider that all input values are elements of 𝒟 𝑝 with the same exponent 𝑒 1 . A different analysis must be done for each possible value of 𝑒 1 .

6. Recall that a 𝒞 ℓ function over an interval 𝐼, ℓ 0, is a function that is ℓ-times differentiable and whose ℓ-th derivative is continuous.

Note that if 𝑥 is a bad case for 𝜙, then it is also a bad case for -𝜙 and -𝑥 is a bad case for 𝑡 ↦ → 𝜙(-𝑡) and 𝑡 ↦ → -𝜙(-𝑡). Hence we can assume that 𝑥 0 and 𝜙(𝑥) 0.

If the values of 𝜙(𝑥), for 𝑥 ∈ [2 𝑒1 , 2 𝑒1+1 ), are not all included in the binade 7 [2 𝑒2 , 2 𝑒2+1 ), we split the input interval into subintervals such that for each subinterval, there is a rational integer 𝑒 2 such that the values 𝜙(𝑥), for 𝑥 in the subinterval, are in [2 𝑒2 , 2 𝑒2+1 ).

We now consider the processing of one subinterval 𝐼 included in [2 𝑒1 , 2 𝑒1+1 ).

In rounding to nearest mode, the problem to be solved is the following: Problem 1 (Rounding to nearest mode). Given 𝑘 ∈ N, estimate the number of integers 𝑋, 2 𝑝-1 𝑋 2 𝑝 -1 (and, possibly, the restrictions implied by 𝑋/2 -𝑒1+𝑝-1 ∈ 𝐼) such that there exists

𝑌 ∈ N, 2 𝑝-1 𝑌 2 𝑝 -1, ⃒ ⃒ ⃒ ⃒ 𝜙 ( 𝑋 2 -𝑒1+𝑝-1 ) - 2𝑌 + 1 2 -𝑒2+𝑝 ⃒ ⃒ ⃒ ⃒ < 2 𝑒2-𝑝-𝑘 .
In directed rounding modes, the problem to be solved is the following: Problem 2 (Directed rounding modes). Given 𝑘 ∈ N, estimate the number of integers 𝑋, 2 𝑝-1 𝑋 2 𝑝 -1 (and, possibly, the restrictions implied by 𝑋/2 -𝑒1+𝑝-1 ∈ 𝐼) such that there exists

𝑌 ∈ N, 2 𝑝-1 𝑌 2 𝑝 -1, ⃒ ⃒ ⃒ ⃒ 𝜙 ( 𝑋 2 -𝑒1+𝑝-1 ) - 𝑌 2 -𝑒2+𝑝-1 ⃒ ⃒ ⃒ ⃒ < 2 𝑒2-𝑝-𝑘
. Now, we introduce the mathematical tools that will allow us to tackle Problems 1 and 2.

SOME EXPONENTIAL SUMS-BASED ESTIMATES

First, let us mention that the reader who is mainly interested in the applications to the design of mathematical functions should focus on Theorem 3, Lemma 4 and Corollary 1, which are the technical results that will be used in practice.

Several deep results (see [START_REF] Bombieri | On the order of 𝜁( 1 2 + 𝑖𝑡)[END_REF], [START_REF] Bombieri | Some mean-value theorems for exponential sums[END_REF], [START_REF] Bourgain | Decoupling, exponential sums and the Riemann zeta function[END_REF] for instance) in Analytic Number Theory rely on non-trivial bounding of the modulus of sums of the form ∑ 𝑎<𝑚 𝑎+𝑀 e(𝑓 (𝑚)), where

𝑎 ∈ Z, 𝑀 ∈ N, 𝑓 : [𝑎 + 1, 𝑎 + 𝑀 ] → R is a 𝒞 𝑘 function with 𝑘
1 and e(𝑡) := exp(2𝑖𝜋𝑡), 𝑡 ∈ R. The developed techniques prove very useful for estimating the number of integer points in the neighborhood of a curve, which actually is the usual context we face when we tackle the Table Maker's Dilemma. In-depth expositions of these techniques are presented in [START_REF] Graham | Van der Corput's method of exponential sums[END_REF], [START_REF] Huxley | Area, lattice points, and exponential sums[END_REF], [START_REF] Montgomery | Ten lectures on the interface between analytic number theory and harmonic analysis[END_REF].

In the present paper we shall use techniques based on the second derivative, and thus assume 𝑓 to be a 𝒞 2 function; first-order techniques (à la Kuzmin-Landau) do not apply to the present context, while higher order methods (when applicable) yield worse upper bounds, probably useless in the formats we are interested in (typically binary64 and binary128).

Given 0 < 𝜂 < 1/2, we want to establish estimates of the value

ℛ(𝑓 ; 𝜂) := # { 1 𝑚 𝑀, ‖𝑓 (𝑚)‖ < 𝜂 } , (1) 7
. A binade is an interval of the form [2 𝑘 , 2 𝑘+1 ) or (-2 𝑘+1 , -2 𝑘 ] for 𝑘 ∈ Z.
where ||𝑥|| denotes the distance of the real number 𝑥 to Z: ||𝑥|| = min({𝑥}, 1-{𝑥}), and {𝑥} = 𝑥-⌊𝑥⌋ is the fractional part of 𝑥.

When 𝜂 = 0, we set ℛ(𝑓 ; 0) := # {1 𝑚 𝑀, 𝑓 (𝑚) ∈ Z} .

In our work, we follow the approach developed by Vaaler in [START_REF] Vaaler | Some extremal functions in Fourier analysis[END_REF]. Let 𝜓(𝑡) denote the normalized fractional part

𝜓(𝑡) := 𝑡 -⌊𝑡⌋ -1 2 (𝑡 ∈ R).
The proof of our main result requires estimates of sums of fractional parts. Since our aim is to obtain upper bounds with suitable numerical constants, we shall use an appropriate approximation of 𝜓(𝑡), and a result due to Vaaler.

We consider the following function

 𝐽(𝑡) := ⎧ ⎨ ⎩ 1 if 𝑡 = 0, 𝜋𝑡(1 -|𝑡|)cot(𝜋𝑡) + |𝑡| if 0 < |𝑡| < 1, 0 if |𝑡| 1.
The (Fourier series type) approximation of 𝜓 we have in mind is the following:

𝜓 * 𝐻 (𝑡) := - ∑ 1 |ℎ|<𝐻 ̂ 𝐽 ( ℎ 𝐻 ) 2𝑖𝜋ℎ e(ℎ𝑡), 𝑡 ∈ R, 𝐻 ∈ N.
We start by recalling an upper bound for the approximation error:

Theorem 1. [START_REF] Vaaler | Some extremal functions in Fourier analysis[END_REF] Using previous notation, one has, for 𝑡 ∈ R, 𝐻 ∈ N,

|𝜓(𝑡) -𝜓 * 𝐻 (𝑡)| 1 2𝐻 ∑ |ℎ|<𝐻 ( 1 - |ℎ| 𝐻
) e(ℎ𝑡).

We shall now give an upper bound for ℛ(𝑓 ; 0) in terms of exponential sums.

Lemma 1. For any function 𝑓

: [1, 𝑀 ] → R, one has ℛ(𝑓 ; 0) 1 𝐻 ∑ |ℎ|<𝐻 ( 1 - |ℎ| 𝐻 ) 𝑀 ∑ 𝑚=1 e (ℎ𝑓 (𝑚)) (2)
for any 𝐻 ∈ N.

Proof: We note that for any 𝐻 ∈ N the term

1 𝐻 2 ⃒ ⃒ ⃒ ⃒ ⃒ 𝐻 ∑ ℎ=1 e(ℎ𝑡) ⃒ ⃒ ⃒ ⃒ ⃒ 2 = 1 𝐻 ∑ |ℎ|<𝐻 ( 1 - |ℎ| 𝐻 ) e(ℎ𝑡)
is nonnegative for any 𝑡 ∈ R, and is 1 for 𝑡 ∈ Z. It follows that

ℛ(𝑓 ; 0) 𝑀 ∑ 𝑚=1 1 𝐻 ∑ |ℎ|<𝐻 ( 1 - |ℎ| 𝐻
) e (ℎ𝑓 (𝑚)) , which gives the expected result.

We are now ready to state our main lemma.

Lemma 2. Let 𝛿 ∈ R >0 , and 𝑀, 𝐻 ∈ N. Then for any function

𝑓 : [1, 𝑀 ] → R one has |ℛ(𝑓 ; 𝛿) -2𝑀 𝛿| 𝐸 1 (𝐻) + 𝐸 2 (𝐻) + ℛ(𝑓 + 𝛿; 0) (3) and |ℛ(𝑓 ; 𝛿) -2𝑀 𝛿| 𝐸 1 (𝐻) + 2𝐸 2 (𝐻) (4)
where we have set

𝐸 1 (𝐻) := 2 𝐻-1 ∑ ℎ=1 | sin(2𝜋𝛿ℎ)| 𝜋ℎ ⃒ ⃒ ⃒ ⃒ ⃒ 𝑀 ∑ 𝑚=1 e (ℎ𝑓 (𝑚)) ⃒ ⃒ ⃒ ⃒ ⃒ , (5)
𝐸 2 (𝐻) := 1 𝐻 ∑ |ℎ|<𝐻 ( 1 - |ℎ| 𝐻 ) ⃒ ⃒ ⃒ ⃒ ⃒ 𝑀 ∑ 𝑚=1 e (ℎ𝑓 (𝑚)) ⃒ ⃒ ⃒ ⃒ ⃒ . (6)
Proof: We have

ℛ(𝑓 ; 𝛿) = 𝑀 ∑ 𝑚=1 #{𝜈 ∈ Z : 𝑓 (𝑚) -𝛿 < 𝜈 < 𝑓 (𝑚) + 𝛿} = 𝑀 ∑ 𝑚=1 (⌊𝑓 (𝑚) + 𝛿⌋ -⌊𝑓 (𝑚) -𝛿⌋) -ℛ(𝑓 + 𝛿; 0) = 2𝑀 𝛿 + 𝑀 ∑ 𝑚=1 (𝜓 (𝑓 (𝑚) -𝛿) -𝜓 (𝑓 (𝑚) + 𝛿))
-ℛ(𝑓 + 𝛿; 0).

Now let

𝐸 := 𝑀 ∑ 𝑚=1 (𝜓 (𝑓 (𝑚) -𝛿) -𝜓 (𝑓 (𝑚) + 𝛿))
and Now Theorem 1 implies that for 𝜂 = ±𝛿 one has

𝐸 * 𝐻 := 𝑀 ∑ 𝑚=1 (𝜓 * 𝐻 (𝑓 (𝑚) -𝛿) -𝜓 * 𝐻 (𝑓 (𝑚) + 𝛿)) .
𝑀 ∑ 𝑚=1 |𝜓 (𝑓 (𝑚) + 𝜂) -𝜓 * 𝐻 (𝑓 (𝑚) + 𝜂) | 1 2𝐻 ∑ |ℎ| 𝐻 ( 1 - |ℎ| 𝐻 ) 𝑀 ∑ 𝑚=1 e (ℎ (𝑓 (𝑚) + 𝜂)) 1 2𝐻 ∑ |ℎ| 𝐻 ( 1 - |ℎ| 𝐻 ) ⃒ ⃒ ⃒ ⃒ ⃒ 𝑀 ∑ 𝑚=1 e (ℎ𝑓 (𝑚)) ⃒ ⃒ ⃒ ⃒ ⃒ , so that |𝐸 -𝐸 * 𝐻 | 𝐸 2 (𝐻). Next, we prove that |𝐸 * 𝐻 | 𝐸 1 (𝐻). Indeed 𝐸 * 𝐻 = 𝑀 ∑ 𝑚=1 𝐻-1 ∑ |ℎ|=1 ̂ 𝐽 ( ℎ 𝐻 ) 2𝑖𝜋ℎ (e(ℎ(𝑓 (𝑚) -𝛿)) -e(ℎ(𝑓 (𝑚) + 𝛿))) = ∑ 1 |ℎ|<𝐻 ̂ 𝐽 ( ℎ 𝐻 ) 𝜋ℎ sin(-2𝜋𝛿ℎ) 𝑀 ∑ 𝑚=1 e (ℎ(𝑓 (𝑚)))
and the expected upper bound follows from the estimate

⃒ ⃒ ⃒ ̂ 𝐽(𝑡) ⃒ ⃒ ⃒ 1 for 𝑡 ∈ R.
Finally, writing

|ℛ(𝑓 ; 𝛿) -2𝑀 𝛿| |𝐸 * 𝐻 | + |𝐸 -𝐸 * 𝐻 | + ℛ(𝑓 + 𝛿; 0)
and using the previous estimates yields [START_REF] Bombieri | On the order of 𝜁( 1 2 + 𝑖𝑡)[END_REF]. Now, to prove (4), it is sufficient to check that ℛ(𝑓 + 𝛿; 0) 𝐸 2 (𝐻). This follows at once from (2).

Remark 1. The reader may have noticed that Lemma 2 provides two upper bounds for the same error |ℛ(𝑓 ; 𝛿) -2𝑀 𝛿|. The reason is the following: in practice, we shall set 𝛿 = 2 -𝑘 , 𝑘 ∈ N and there are functions, such as exp or cos for instance, for which we will be able to estimate precisely the quantity ℛ(𝑓 + 2 -𝑘 ; 0). In this case, inequality (3) will give us a better information than (4). However, the latter has the advantage of being general.

Theorem 3 and Corollary 1 are presented the same way, for the same reason.

We will need the following technical inequalities. Lemma 3. We have, for all 𝐻 ∈ N,

𝐻 ∑ ℎ=1 1 ℎ 1/2 2𝐻 1/2 ; (7) 𝐻-1 ∑ ℎ=1 ℎ 1/2 ( 1 - ℎ 𝐻 ) 4 15 𝐻 3/2 + 𝐻 1/2 ; (8) 𝐻-1 ∑ ℎ=1 1 ℎ 1/2 ( 1 - ℎ 𝐻 ) 4 3 𝐻 1/2 ; (9)
+∞ ∑ ℎ=1 | sin(2𝜋𝛿ℎ)| ℎ 3/2 5 2 (2𝜋𝛿) 1/2 + (2𝜋𝛿) 3/2 , 𝛿 > 0. ( 10
)
Proof: Note that the first three estimates are trivial for 𝐻 = 1. Hence in the sequel we suppose that 𝐻 2.

1) For ( 7), we have

𝐻 ∑ ℎ=1 1 ℎ 1/2 𝐻 ∑ ℎ=1 ∫ ℎ ℎ-1 d𝑡 𝑡 1/2 = ∫ 𝐻 0 d𝑡 𝑡 1/2 = 2𝐻 1/2 .
2) For (8), we have

𝐻-1 ∑ ℎ=1 ℎ 1/2 𝐻-1 ∑ ℎ=1 ∫ ℎ+1 ℎ 𝑡 1/2 d𝑡 = ∫ 𝐻 1 𝑡 1/2 d𝑡 2 3 𝐻 3/2 . Moreover 𝐻-1 ∑ ℎ=1 ℎ 3/2 𝐻 = -𝐻 1/2 + 𝐻 ∑ ℎ=1 ℎ 3/2 𝐻 -𝐻 1/2 + 𝐻 ∑ ℎ=1 ∫ ℎ ℎ-1 𝑡 3/2 𝐻 d𝑡 = -𝐻 1/2 + 1 𝐻 ∫ 𝐻 0 𝑡 3/2 d𝑡 = -𝐻 1/2 + 2 5 𝐻 3/2 .
Using these inequalities, we now have

𝐻-1 ∑ ℎ=1 ℎ 1/2 - 𝐻-1 ∑ ℎ=1 ℎ 3/2 𝐻 2 3 𝐻 3/2 - 2 5 𝐻 3/2 + 𝐻 1/2
which gives the expected result.

3) Similarly, for [START_REF] De Dinechin | Towards the postultimate libm[END_REF] we have

𝐻-1 ∑ ℎ=1 1 ℎ 1/2 ( 1 - ℎ 𝐻 ) = 𝐻 ∑ ℎ=1 1 ℎ 1/2 ( 1 - ℎ 𝐻 ) 𝐻 ∑ ℎ=1 ∫ ℎ ℎ-1 1 𝑡 1/2 ( 1 - 𝑡 𝐻 ) d𝑡 = ∫ 𝐻 0 1 𝑡 1/2 ( 1 - 𝑡 𝐻 ) d𝑡 = 4 3 𝐻 1/2 .

4)

We now estimate the last sum. We have

∑ ℎ 1/(2𝜋𝛿) | sin(2𝜋𝛿ℎ)| ℎ 3/2 2𝜋𝛿 ∑ ℎ 1/(2𝜋𝛿) 1 ℎ 1/2 4𝜋𝛿⌊1/(2𝜋𝛿)⌋ 1/2 2(2𝜋𝛿) 1/2
where we have used [START_REF] Cody | A proposed radix-and-word-length-independent standard for floating-point arithmetic[END_REF].

For the remaining terms, we write

∑ ℎ 1+⌊1/(2𝜋𝛿)⌋ | sin(2𝜋𝛿ℎ)| ℎ 3/2 ∑ ℎ 1+⌊1/(2𝜋𝛿)⌋ 1 ℎ 3/2 = 1 (1 + ⌊1/(2𝜋𝛿)⌋) 3/2 + ∑ ℎ 2+⌊1/(2𝜋𝛿)⌋ 1 ℎ 3/2 (2𝜋𝛿) 3/2 + ∑ ℎ 2+⌊1/(2𝜋𝛿)⌋ ∫ ℎ ℎ-1 d𝑡 𝑡 3/2 (2𝜋𝛿) 3/2 + ∫ +∞ 1/(2𝜋𝛿) d𝑡 𝑡 3/2 (2𝜋𝛿) 3/2 + 1 2 (2𝜋𝛿) 1/2 .
We now state an explicit version of van der Corput's inequality for exponential sums : this is Equation (6.14) of [35, p. 128].

Theorem 2. Let 𝑀 ∈ N, 𝜆 ∈ R >0 and 𝐶 ∈ [1, +∞). Let 𝑓 ∈ 𝒞 2 ( [1, 𝑀 ], R ) such that 𝜆 |𝑓 ′′ (𝑥)| 𝐶𝜆 for all 𝑥 ∈ [1, 𝑀 ]. Then ⃒ ⃒ ⃒ ⃒ ⃒ 𝑀 ∑ 𝑚=1 e(𝑓 (𝑚)) ⃒ ⃒ ⃒ ⃒ ⃒ 3𝐶𝑀 𝜆 1/2 + 6𝜆 -1/2 .
We should now explain why we focus on the case 𝑓 ∈ 𝒞 2 . The main term in the upper bound of Theorem 2 is the term 𝑀 𝜆 1/2 which provides a saving 𝜆 1/2 with respect to the trivial bound 𝑀 (at least when 𝜆 is small). More generally, van der Corput's method asserts that for 𝑘 2, when 𝑓 ∈ 𝒞 𝑘 on [1, 𝑀 ] and |𝑓 (𝑘) | has the order of magnitude 𝜆 𝑘 > 0 small, the saving is 𝜆 𝜃 𝑘 𝑘 with 𝜃 𝑘 = 1/(2 𝑘 -2) (see Theorem 3 of [START_REF] Robert | On van der Corput's 𝑘-th derivative test for exponential sums[END_REF]). In the applications we have in mind, |𝑓 (𝑘) | has a size close to 𝑀 1-𝑘 , which provides a saving 𝑀 (1-𝑘)𝜃 𝑘 . It is now plain that this saving is maximal for 𝑘 = 2.

Now we can prove the main result of this section: we give an explicit bound which controls the difference between ℛ(𝑓 ; 𝛿) and the expected main term 2𝑀 𝛿. In Section 4, we shall use it, jointly with Lemma 4, or Corollary 1 to legitimate the heuristic probabilistic estimate approach from Subsection 1.5: ℛ(𝑓 ; 𝛿) will count the number of bad cases and the main term 2𝑀 𝛿 will play the role of the probabilistic estimate. Note that we may obtain simpler expressions for the upper bounds in [START_REF] Dunham | Feasibility of "perfect" function evaluation[END_REF] and [START_REF] Fousse | MPFR: A Multiple-Precision Binary Floating-Point Library with Correct Rounding[END_REF]. However, this option has been ruled out in order to keep the best estimate possible. Corollary 1 actually gives more ready-to-use bounds.

Theorem 3. Let 𝑀 ∈ N ∖ {0}, 𝜆 ∈ R >0 and 𝐶 ∈ [1, +∞). Let 𝑓 ∈ 𝒞 2 ( [1, 𝑀 ], R ) such that 𝜆 |𝑓 ′′ (𝑥)| 𝐶𝜆 for all 𝑥 ∈ [1, 𝑀 ].
Then for any 𝛿 > 0, one has

|ℛ(𝑓 ; 𝛿) -2𝑀 𝛿| Δ 1 (𝑀, 𝜆) + ( 60 (2𝜋) 1/2 + 24(2𝜋) 1/2 𝛿 ) ( 𝛿 𝜆 ) 1/2 + ℛ(𝑓 + 𝛿; 0) (11)
and

|ℛ(𝑓 ; 𝛿) -2𝑀 𝛿| Δ 2 (𝑀, 𝜆) + ( 60 (2𝜋) 1/2 + 24(2𝜋) 1/2 𝛿 ) ( 𝛿 𝜆 ) 1/2 (12)
where we have set

𝛼 0 = 6 𝜋 + 4 5 , 𝛽 0 = 3 𝜋 + 4 5 , Δ 1 (𝑀, 𝜆, 𝐻) = 𝑀 𝐻 + 2𝛼 0 𝐶𝑀 𝜆 1/2 𝐻 1/2 + 6𝐶𝑀 𝜆 1/2 𝐻 1/2 + 16 𝜆 -1/2 𝐻 1/2 , Δ 2 (𝑀, 𝜆, 𝐻) = 2𝑀 𝐻 + 4𝛽 0 𝐶𝑀 𝜆 1/2 𝐻 1/2 + 12𝐶𝑀 𝜆 1/2 𝐻 1/2 + 32 𝜆 -1/2 𝐻 1/2 , Δ 1 (𝑀, 𝜆) = min 𝐻∈N Δ 1 (𝑀, 𝜆, 𝐻), Δ 2 (𝑀, 𝜆) = min 𝐻∈N Δ 2 (𝑀, 𝜆, 𝐻).
Proof: Using Lemma 2, we have, for any 𝐻 ∈ N

|ℛ(𝑓 ; 𝛿) -2𝑀 𝛿| 𝐸 1 (𝐻) + 𝐸 2 (𝐻) + ℛ(𝑓 + 𝛿; 0) and |ℛ(𝑓 ; 𝛿) -2𝑀 𝛿| 𝐸 1 (𝐻) + 2𝐸 2 (𝐻),
where 𝐸 1 (𝐻) and 𝐸 2 (𝐻) are defined by ( 5) and ( 6) respectively.

Our aim is now to bound each of the terms 𝐸 𝑗 using Theorem 2 for the inner sums, and Lemma 3 for the sums involved.

For 𝐸 1 , Theorem 2 yields

𝐸 1 (𝐻) 2 𝐻 ∑ ℎ=1 | sin(2𝜋𝛿ℎ)| 𝜋ℎ ( 3𝐶𝑀 𝜆 1/2 ℎ 1/2 + 6𝜆 -1/2 ℎ -1/2
) , and we have

2 𝐻 ∑ ℎ=1 | sin(2𝜋𝛿ℎ)| 𝜋ℎ 1/2 2 𝜋 𝐻 ∑ ℎ=1 1 ℎ 1/2 4 𝜋 𝐻 1/2 , 𝐻 ∑ ℎ=1 | sin(2𝜋𝛿ℎ)| 𝜋ℎ 3/2 ( 5 (2𝜋) 1/2 + 2(2𝜋) 1/2 𝛿 ) 𝛿 1/2
by using ( 7) and ( 10) respectively. This gives

𝐸 1 (𝐻) 𝐹 1 (𝐻) + ( 60 (2𝜋) 1/2 + 24(2𝜋) 1/2 𝛿 ) ( 𝛿 𝜆 ) 1/2
where we have set

𝐹 1 (𝐻) = 12 𝜋 𝐶𝑀 𝜆 1/2 𝐻 1/2 .
Similarly, for 𝐸 2 (𝐻), Theorem 2 gives

𝐸 2 (𝐻) 𝑀 𝐻 + 2 𝐻 𝐻-1 ∑ ℎ=1 ( 1 - |ℎ| 𝐻 ) ( 3𝐶𝑀 (ℎ𝜆) 1/2 + 6(ℎ𝜆) -1/2
) , and using ( 8) and ( 9) we have

𝐸 2 (𝐻) 𝐹 2 (𝐻)
where we have set

𝐹 2 (𝐻) = 𝑀 𝐻 + 8𝐶 5 𝑀 𝜆 1/2 𝐻 1/2 + 6𝐶𝑀 𝜆 1/2 𝐻 1/2 + 16 𝜆 -1/2 𝐻 1/2 .
Gathering the contribution of each 𝐸 𝑗 (𝐻), we have the bounds

|ℛ(𝑓 ; 𝛿) -2𝑀 𝛿| 𝐹 1 (𝐻) + 𝐹 2 (𝐻) + ℛ(𝑓 + 𝛿; 0) + ( 60 (2𝜋) 1/2 + 24(2𝜋) 1/2 𝛿 ) ( 𝛿 𝜆 ) 1/2 and |ℛ(𝑓 ; 𝛿) -2𝑀 𝛿| 𝐹 1 (𝐻) + 2𝐹 2 (𝐻) + ( 60 (2𝜋) 1/2 + 24(2𝜋) 1/2 𝛿 ) ( 𝛿 𝜆 ) 1/2
valid for any 𝐻 ∈ N. The expected result follows from taking the minimum in 𝐻 in both inequalities.

The form of the statement of Theorem 3 allows us to choose, for a given function, a value of 𝐻 that yields small values (and if possible, the smallest ones) of the upper bounds. But, it can prove convenient to give our reader some closed and ready-to-use expressions for the upper bounds for the error |ℛ(𝑓 ; 𝛿) -2𝑀 𝛿|. This is the purpose of Lemma 4 and Corollary 1.

Lemma 4. Let 𝑀 ∈ N ∖ {0}, 𝜆 ∈ R >0 and 𝐶 ∈ [1, +∞).
Then with the notation of Theorem 3, we have

Δ 1 (𝑀, 𝜆) 𝛼 2/3 0 𝐶 2/3 𝑀 𝜆 1/3 + 2𝛼 2/3 0 𝐶 2/3 max(𝑀 𝜆 1/3 , 𝑀 1/2 ) + 16𝛼 1/3 0 𝐶 1/3 𝑀 1/2 + (6 + 2𝛼 0 )𝛼 1/3 0 𝐶 4/3 𝑀 𝜆 2/3 and Δ 2 (𝑀, 𝜆) 2𝛽 2/3 0 𝐶 2/3 𝑀 𝜆 1/3 + 4𝛽 2/3 0 𝐶 2/3 max(𝑀 𝜆 1/3 , 𝑀 1/2 ) + 32𝛽 1/3 0 𝐶 1/3 𝑀 1/2 + (12 + 4𝛽 0 )𝛽 1/3 0 𝐶 4/3 𝑀 𝜆 2/3 . Proof:
We treat the cases 𝑀 𝜆 2/3 1 and 𝑀 𝜆 2/3 < 1 separately. In order to bound some of the terms below, we shall use the following inequality

√ 1 + ⌊𝑡⌋ √ 𝑡 + 1 √ 𝑡 (𝑡 > 0). ( 13
)
With the definition of 𝛼 0 , we recall that Theorem 3 gives

Δ 1 (𝑀, 𝜆) 𝑀 𝐻 + 2𝛼 0 𝐶𝑀 𝜆 1/2 𝐻 1/2 + 6𝐶 𝑀 𝜆 1/2 𝐻 1/2 + 16 𝜆 -1/2 𝐻 1/2
for any 𝐻 ∈ N.

If 𝑀 𝜆 2/3 1, we set 𝑥 := (𝛼 0 𝐶) -2/3 𝜆 -1/3 and we choose 𝐻 := 1 + ⌊𝑥⌋. In particular, 𝐻 𝑥, and (13) implies 𝐻 1/2 𝑥 1/2 + 𝑥 -1/2 , so that

Δ 1 (𝑀, 𝜆) 𝑀 𝑥 + 2𝛼 0 𝐶𝑀 𝜆 1/2 𝑥 1/2 + (6 + 2𝛼 0 )𝐶 𝑀 𝜆 1/2 𝑥 1/2 + 16 𝜆 -1/2 𝑥 1/2 3𝛼 2/3 0 𝐶 2/3 𝑀 𝜆 1/3 + (6 + 2𝛼 0 )𝛼 1/3 0 𝐶 4/3 𝑀 𝜆 2/3 + 16𝛼 1/3 0 𝐶 1/3 𝜆 -1/3
which implies the expected result since

3𝛼 2/3 0 𝐶 2/3 𝑀 𝜆 1/3 𝛼 2/3 0 𝐶 2/3 𝑀 𝜆 1/3 + 2𝛼 2/3 0 𝐶 2/3 max(𝑀 𝜆 1/3 , 𝑀 1/2 )
and 16𝛼

1/3 0 𝐶 1/3 𝜆 -1/3 16𝛼 1/3 0 𝐶 1/3 𝑀 1/2 .
If 𝑀 𝜆 2/3 1, we set 𝑥 := (𝛼 0 𝐶) -2/3 (𝑀 𝜆) -1 and we choose 𝐻 := 1 + ⌊𝑥⌋. In particular, 𝐻 𝑥, and (13) implies 𝐻 1/2 𝑥 1/2 + 𝑥 -1/2 , so that

Δ 1 (𝑀, 𝜆) 𝑀 𝑥 + 2𝛼 0 𝐶𝑀 𝜆 1/2 𝑥 1/2 + (6 + 2𝛼 0 )𝐶 𝑀 𝜆 1/2 𝑥 1/2 + 16 𝜆 -1/2 𝑥 1/2 𝛼 2/3 0 𝐶 2/3 𝑀 2 𝜆 + 2𝛼 2/3 0 𝐶 2/3 𝑀 1/2 + (6 + 2𝛼 0 )𝛼 1/3 0 𝐶 4/3 𝑀 3/2 𝜆 + 16𝛼 1/3 0 𝐶 1/3 𝑀 1/2 ,
which implies the expected result since

𝛼 2/3 0 𝐶 2/3 𝑀 2 𝜆 𝛼 2/3 0 𝐶 2/3 𝑀 𝜆 1/3 , 2𝛼 2/3 0 𝐶 2/3 𝑀 1/2 2𝛼 2/3 0 𝐶 2/3 max(𝑀 𝜆 1/3 , 𝑀 1/2 ) and (6 + 2𝛼 0 )𝛼 1/3 0 𝐶 4/3 𝑀 3/2 𝜆 (6 + 2𝛼 0 )𝛼 1/3 0 𝐶 4/3 𝑀 𝜆 2/3 .
For the estimates for Δ 2 (𝑀, 𝜆), we notice that Theorem 3 gives

Δ 2 (𝑀, 𝜆) 2 𝑀 𝐻 +2𝛽 0 𝐶𝑀 𝜆 1/2 𝐻 1/2 +6𝐶 𝑀 𝜆 1/2 𝐻 1/2 +16 𝜆 -1/2 𝐻 1/2
for any 𝐻 ∈ N. Reproducing the previous computations with 𝛽 0 instead of 𝛼 0 gives the expected result. This completes the proof. Combining Theorem 3 and Lemma 4, we obtain the following:

Corollary 1. We have, under the assumptions of Theorem 3, the two estimates

|ℛ(𝑓 ; 𝛿) -2𝑀 𝛿| 𝛼 2/3 0 𝐶 2/3 𝑀 𝜆 1/3 + 2𝛼 2/3 0 𝐶 2/3 max(𝑀 𝜆 1/3 , 𝑀 1/2 ) + 16𝛼 2/3 0 𝐶 2/3 𝑀 1/2 + (6 + 2𝛼 0 )𝛼 1/3 0 𝐶 4/3 𝑀 𝜆 2/3 + ( 60 (2𝜋) 1/2 + 24(2𝜋) 1/2 𝛿 ) ( 𝛿 𝜆 ) 1/2 + ℛ(𝑓 + 𝛿; 0) (14)
and

|ℛ(𝑓 ; 𝛿) -2𝑀 𝛿| 2𝛽 2/3 0 𝐶 2/3 𝑀 𝜆 1/3 + 4𝛽 2/3 0 𝐶 2/3 max(𝑀 𝜆 1/3 , 𝑀 1/2 ) + 32𝛽 2/3 0 𝐶 2/3 𝑀 1/2 + (12 + 4𝛽 0 )𝛽 1/3 0 𝐶 4/3 𝑀 𝜆 2/3 + ( 60 (2𝜋) 1/2 + 24(2𝜋) 1/2 𝛿 ) ( 𝛿 𝜆 ) 1/2 (15)
for any 𝛿 > 0.

As a final remark, let us point that the term (𝛿/𝜆) More generally, the same would hold for a function with 𝑓 (0) ∈ Z, 𝑓 ′ (0) = 𝑢/𝑣, and with derivatives of order 3 vanishing quickly enough so that 𝑓 is well approximated by a polynomial of degree 2 over (-(2𝛿/𝜆) 1/2 , (2𝛿/𝜆) 1/2 ). In terms of bad cases, this is, e.g., the well-known situation of the cosine function close to 0, or of the cube root close to 1. 8 . Let 𝐴 = 𝑎2 𝑝-1-𝑒1 and 𝐵 = 𝑏2 𝑝-1-𝑒1 ∈ N, and finally 𝑀 = 𝐵 -𝐴 + 1. In order to tackle Problems 1 and 2 of Section 2, we introduce ℐ(𝜙; 𝑝, 𝐴, 𝐵, 𝛿) :=

AN ESTIMATE OF THE NUMBER OF "BAD CASES" FOR THE FIRST QUICK STEP OF THE TWO STEP ZIV

STRATEGY Let 𝑝 1, 𝑒 1 , 𝑒 2 ∈ Z, 𝑎, 𝑏 ∈ 𝒟 𝑝 two floating-point numbers such that [𝑎, 𝑏] ⊂ [2 𝑒1 , 2 𝑒1+1 ), 𝛿 > 0 and 𝜙 : [𝑎, 𝑏] → R be a 𝒞 ℓ function, ℓ 2, such that 𝜙([𝑎, 𝑏]) ⊂ [2 𝑒2 , 2 𝑒2+1 )
• # { 𝐴 𝑗 𝐵, ⃦ ⃦ ⃦ ⃦ 2 𝑝-1-𝑒2 𝜙 ( 𝑗 2 𝑝-1-𝑒1 ) - 1 2 ⃦ ⃦ ⃦ ⃦ < 𝛿 } (16)
in rounding to nearest mode,

• # { 𝐴 𝑗 𝐵, ⃦ ⃦ ⃦ ⃦ 2 𝑝-1-𝑒2 𝜙 ( 𝑗 2 𝑝-1-𝑒1 )⃦ ⃦ ⃦ ⃦ < 𝛿 } ( 17
)
in directed rounding modes.

We call (𝐴, 𝐵, 𝛿)-bad cases (or simply bad cases when there is no ambiguity) the elements of these sets.

Our goal is to obtain fully effective estimates of the form ℐ(𝜙; 𝑝, 𝐴, 𝐵, 𝛿) = 2𝑀 𝛿 + 𝑜(𝑀 𝛿), or, at least, 8. As already said at the beginning of Section 2, considering this case is sufficient to address the general case

[𝑎, 𝑏] ⊂ ±[2 𝑒 1 , 2 𝑒 1 +1 ) and 𝜙([𝑎, 𝑏]) ⊂ ±[2 𝑒 2 , 2 𝑒 2 +1 ).
ℐ(𝜙; 𝑝, 𝐴, 𝐵, 𝛿) = 2𝑀 𝛿 + 𝑂(𝑀 𝛿). Though less accurate, the latter may still yield a fairly good upper bound for the number of bad cases. Recall that if 𝛿 = 2 -𝑘 , 𝑘 ∈ N, 2𝑀 𝛿 is the number of bad cases suggested by the heuristic approach of Subsection 1.5, so that ℐ(𝜙; 𝑝, 𝐴, 𝐵, 2 -𝑘 ) is exactly the number of bad cases for the first quick step of the two step Ziv strategy using intermediate accuracy 2 1-𝑝-𝑘 . In the sequel, we shall thus use 𝛿 = 2 -𝑘 for some 𝑘 ∈ N (the number of extra bits in the first step of Ziv strategy).

Remark 2. In the context of a two step Ziv strategy, a good estimate on ℐ(𝜙; 𝑝, 𝐴, 𝐵, 𝛿) allows one to give an average analysis of the cost of an evaluation of 𝜙(𝑥) for 𝑥 ∈ [𝑎, 𝑏].

Indeed, if the first step uses relative accuracy 2 1-𝑝-𝑘 , the second using relative accuracy 2 1-𝑝-𝑘 ′ (with 𝑘 ′ large enough so as to capture worst cases), the average cost of one evaluation of the function is 𝑐 𝜙,𝑝 (𝑘) + 𝑐 𝜙,𝑝 (𝑘 ′ )ℐ(𝜙; 𝑝, 𝐴, 𝐵, 2 -𝑘 )/𝑀 , where 𝑐 𝜙,𝑝 (𝑗) stands for the cost of one evaluation of 𝑓 with relative accuracy 2 1-𝑝-𝑗 . For choices of parameters for which our estimates yield ℐ(𝜙; 𝑝, 𝐴, 𝐵, 2 -𝑘 ) ≈ 2 1-𝑘 𝑀 , we obtain an average cost ≈ 𝑐 𝜙,𝑝 (𝑘)+2 1-𝑘 𝑐 𝜙,𝑝 (𝑘 ′ ) (the ≈ meaning here that we can determine a rigorous, small interval centered at 𝑐 𝜙,𝑝 (𝑘) + 2 1-𝑘 𝑐 𝜙,𝑝 (𝑘 ′ ) containing the average cost).

In order to connect these questions and the results obtained on (1),

•

In rounding to nearest mode, we introduce the function

𝑓 (𝑥) := 2 𝑝-1-𝑒2 𝜙 ( 𝐴 -1 + 𝑥 2 𝑝-1-𝑒1 ) - 1 2 , 1 𝑥 𝑀, (18)
and observe that ℛ(𝑓 ; 𝛿) = ℐ(𝜙; 𝑝, 𝐴, 𝐵, 𝛿).

•

Similarly, in directed rounding modes, we introduce the function

𝑓 (𝑥) := 2 𝑝-1-𝑒2 𝜙 ( 𝐴 -1 + 𝑥 2 𝑝-1-𝑒1 ) , 1 𝑥 𝑀, (19)
and observe again that ℛ(𝑓 ; 𝛿) = ℐ(𝜙; 𝑝, 𝐴, 𝐵, 𝛿).

We are thus in a position to apply Theorem 3 to derive rigorous and useful estimates for ℐ(𝜙; 𝑝, 𝐴, 𝐵, 𝛿). Note that, whatever the rounding mode we choose, only the values of the second derivative of 𝑓 matter in Theorem 3, hence the estimates that we obtain will be the same.

We now give the precise values of the constants involved in the statement of Theorem 3, namely 𝐶 and 𝜆.

𝜆 = min 𝑥∈[1,𝑀 ] |𝑓 ′′ (𝑥)| = 2 -𝑝+1+2𝑒1-𝑒2 min 𝑥∈[𝑎,𝑏] |𝜙 ′′ (𝑥)| .
Recall that 𝜆 has to be positive. Now,

𝐶𝜆 = max 𝑥∈[1,𝑀 ] |𝑓 ′′ (𝑥)| = 2 -𝑝+1+2𝑒1-𝑒2 max 𝑥∈[𝑎,𝑏] |𝜙 ′′ (𝑥)| , hence 𝐶 = max 𝑥∈[1,𝑀 ] |𝑓 ′′ (𝑥)| min 𝑥∈[1,𝑀 ] |𝑓 ′′ (𝑥)| = max 𝑥∈[𝑎,𝑏] |𝜙 ′′ (𝑥)| min 𝑥∈[𝑎,𝑏] |𝜙 ′′ (𝑥)| .
Note that this formula is legitimate since we assumed

min 𝑥∈[𝑎,𝑏] |𝜙 ′′ (𝑥)| = 2 𝑝-1-2𝑒1+𝑒2 𝜆 > 0.
Remark 3. In order to have 𝜆 > 0, one might have to split the interval into sub-intervals over which 𝜙 ′′ has no zero. More generally, if 𝜙 ′′ has huge variations, one might consider splitting the interval to avoid too large a value of 𝐶.

Before applying Theorem 3 in order to estimate the number of bad cases of the first quick step, we now wish to understand the constraints on the parameters 𝑀, 𝑘, 𝜆, 𝐶, and in particular, the number of bits 𝑘 for which we can derive a relevant result from Theorem 3. Therefore, we check the various constraints that the parameters must satisfy in order to make inequality (15) (from Corollary 1) relevant. This inequality is a consequence of Theorem 3 but, in our practical experiments in binary64 and binary128, the results it yields are very close to the one provided by Theorem (3).

Analysis of the results established in Section 3: constraints on the parameters

The main term 2𝑀 𝛿 should not be smaller than any of the terms of the right-hand side of [START_REF] Huxley | Area, lattice points, and exponential sums[END_REF]. So let us check the following conditions:

• we have 2𝑀 𝛿 4𝛽 2/3 0 𝐶 2/3 𝑀 𝜆 1/3 if and only if log 2 (1/𝜆) 3 - 2 log 2 (𝛽 0 𝐶) 3 -1 log 2 (1/𝛿) i.e. 𝑝 -2𝑒 1 + 𝑒 2 -log 2 ( min 𝑥∈[𝑎,𝑏] |𝜙 ′′ (𝑥)| ) -2 log 2 (𝛽 0 𝐶) -4 3𝑘. • We have 2𝑀 𝛿 32𝛽 2/3 0 𝐶 2/3 𝑀
2 log 2 (1/𝜆) 3 -log 2 (3 + 𝛽 0 ) - log 2 (𝛽 0 ) 3 - 4 log 2 (𝐶) 3 -1 log 2 (1/𝛿)
i.e.

2(𝑝-2𝑒

1 +𝑒 2 )-2 log 2 ( min 𝑥∈[𝑎,𝑏] |𝜙 ′′ (𝑥)| ) -3 log 2 (3+𝛽 0 ) -log 2 (𝛽 0 ) -4 log 2 (𝐶) -5 3𝑘.
• Now, we check the condition

2𝑀 𝛿 ( 60 (2𝜋) 1/2 + 24(2𝜋) 1/2 𝛿 ) ( 𝛿 𝜆 ) 1/2 .
Note that, as soon as 𝑘 2, 60/(2𝜋) 1/2 24(2𝜋) 1/2 𝛿. Therefore, we focus on the condition 2𝑀 𝛿 (60/(2𝜋)

1/2 ) ( 𝛿 𝜆 ) 1/2 , which is equivalent to 2 log 2 (𝑀 ) -log 2 (1/𝜆) -1 -2 log 2 (15) + log 2 (𝜋) 𝑘 i.e. 2 log 2 (𝑀 ) -𝑝 + 2𝑒 1 -𝑒 2 + log 2 ( min 𝑥∈[𝑎,𝑏] |𝜙 ′′ (𝑥)| ) -2 log 2 (15) + log 2 (𝜋) 𝑘.
Remark 4. It is probably worth commenting a little bit informally on the error term. Typically, for values of 𝑥 such that 𝑥 and 𝑓 (𝑥) have exponent close to 0, in our context, 𝜆 is going to be of the order of 2 -𝑝 ; assuming for the sake of simplicity that all constants including 𝐶 are 1, the upper bound is of the order of 𝑀 2 -𝑝/3 + 𝑀 1/2 + 2 (𝑝-𝑘)/2 (recall that 𝑀 is the number of floating-point numbers in the interval under consideration, so that the trivial upper bound would be 𝑀 ). As we can restrict to 𝑘 𝑝/3, the middle term is always absorbed by one of the other two. The error term is thus either roughly 𝑀 2 -𝑝/3 for 𝑀 2 2𝑝/3 ("large intervals") or 2 (𝑝-𝑘)/2 for 𝑀 < 2 2𝑝/3 ("small intervals"). Now, we illustrate our results on two examples: the cube root and the exponential functions. In both cases, we will first have to decompose the input binade in order for the output values to fit in one single binade, then find explicit expressions for 𝐶 and 𝜆 before being able to conduct actual experiments.

The 3 √ • function

In this case, we have 𝜙(𝑥) = 𝑥 1/3 . This function never overflows nor underflows in its exponent range. Note that it takes many exact values: for instance, if 𝑥 ∈ 𝒟 ⌊𝑝/3⌋ , 𝑦 = 𝑥 3 ∈ 𝒟 𝑝 has 𝑥 as an exact cube root in 𝒟 𝑝 .

Explicit values of the parameters

Let us apply our method on an interval [𝑎, 𝑏] enclosed in a single binade of exponent 𝑒 1 . Note that we can assume 𝑎, 𝑏 0 as the cube root function is odd; further, as

3 √ 2 3 • 𝑥 = 2 3
√ 𝑥, we can restrict our study to any 3 consecutive binades, for instance [1/2, 1), [START_REF]IEEE Standard for Binary Floating-Point Arithmetic[END_REF][START_REF]IEEE Standard for Radix Independent Floating-Point Arithmetic[END_REF], [START_REF]IEEE Standard for Radix Independent Floating-Point Arithmetic[END_REF][START_REF] Bombieri | Some mean-value theorems for exponential sums[END_REF].

We have

𝜙([𝑎, 𝑏]) ∩ 𝒟 𝑝 =[𝑎 1/3 , 𝑏 1/3 ] ∩ 𝒟 𝑝 ⊂[2 𝑒1/3 , 2 (𝑒1+1)/3 ) ∩ 𝒟 𝑝 ⊂[2 ⌊𝑒1/3⌋ , 2 ⌊𝑒1/3⌋+1
).

As such, we have 𝑒 2 = ⌊𝑒 1 /3⌋. Therefore, and as 𝑎, 𝑏 > 0, we have

𝜆 = 2 -𝑝+2+2𝑒1-⌊𝑒1/3⌋ 9𝑏 5/3 , 𝐶 = max 𝑥∈[𝑎,𝑏] |𝜙 ′′ (𝑥)| min 𝑥∈[𝑎,𝑏] |𝜙 ′′ (𝑥)| = (𝑏/𝑎) 5/3 .
Overall, we can always treat a full input binade at a time, so that 𝑀 = 2 𝑝-1 .

A simplified statement

We can now collect our results to illustrate the use of Lemma 4 to obtain a simple bound on the number of bad cases, depending only on the natural parameters 𝑝, 𝛿 of the problem. We actually give two bounds; the second one is more compact but loses a small constant factor as 𝑝 grows large.

Corollary 2. For the cube root function over the binade [1/2, 1), if 𝑝 11 the number of bad cases is in the interval [𝛿2 𝑝 -𝜌, 𝛿2 𝑝 +𝜌], where

𝜌 5.72 • 2 2𝑝/3 + 40.2 • 2 𝑝/2 + 19.7 • 2 𝑝/3 + 110 • 𝛿 1/2 2 𝑝/2 .
If, further, 𝑝 24, we have 𝜌 8.31

• 2 2𝑝/3 + 110 • 𝛿 1/2 2 𝑝/2 .
Proof: Follows from the second part of Lemma 4 by specializing the values of 𝐶, 𝜆 as seen before. Note that we have 𝜆 = 2 1-𝑝 /9(1 -2 -𝑝 ) 5/3 , which is 2 1-𝑝 /9 and 1.001 • 2 1-𝑝 /9 for 𝑝 11.

The second claim simply uses 2 𝑝/2 2 2𝑝/3 2 -4 and 2 𝑝/3 2 2𝑝/3 2 -8 for 𝑝 24.

Experimental results

We shall treat the binary64 (Table 1) and binary128 (Table 2) cases. We chose not to deal with the binary32 case as it is so easy, in that case, to compute bad cases exactly that our work seemed irrelevant. On the other hand, for binary64 and binary128, it is not realistic to do the computation of the exact values for a whole binade; we thus compare our bounds with the expected value under the heuristic probabilistic model. We shall also limit ourselves to the binade [1/2, 1), the results for the other two binades being slightly different, but very similar.

Our experiments are based on Theorem 3, but we checked that Lemma 4 gives the exact same tables (given the number of decimal figures displayed in the tables) in that setting.

The second column (expected) gives the expected number of bad cases under the probabilistic heuristic of Subsection 1.5, that is to say

2𝑀 𝛿 = 𝑀 2 -𝑘+1 .
The last column deserves further explanation: rather than giving the exact lower and upper bound, we prefer giving their "relative quality", which means the bound for |ℛ(𝑓 ; 𝛿)/2𝑀 𝛿 -1|. Thus, a bound close to 0% means that our bounds give a very precise order of magnitude, while a bound larger than 100% is not so good -though, while the corresponding lower bound is meaningless, the upper bound usually remains non trivial, even often useful as it still gives the correct order of magnitude. For instance, Table 1

The exp function

In this case, we have 𝜙(𝑥) = exp(𝑥). Note that in this section, we shall treat the exponential function both for positive and negative values of 𝑥. Equivalently, this can be seen as treating the two functions exp and 𝑥 ↦ → exp(-𝑥) for positive arguments, which is consistent with the formalization adopted in Section 2. The exp function overflows for large, positive 𝑥 and underflows for large, negative 𝑥; on the other hand, the only rational value of exp(𝑥), with 𝑥 a rational number, is exp(0) = 1, so that for the rational 𝛿 we use we shall ignore the term ℛ(𝑓 + 𝛿; 0) in ( 3), [START_REF] Dunham | Feasibility of "perfect" function evaluation[END_REF] and [START_REF] Graham | Van der Corput's method of exponential sums[END_REF]. We shall use that repeatedly in this section, together with the fact that log 2 is irrational so that 𝑥/ log 2 and 𝑥 log 2 can never be an integer when 𝑥 is a nonzero rational number.

𝑘

Explicit values of the parameters

Let [𝑎, 𝑏] be an interval such that ±[𝑎, 𝑏] ⊂ [2 𝑒1 , 2 𝑒1+1 ) and exp([𝑎, 𝑏]) ⊂ [2 𝑒2 , 2 𝑒2+1
). First observe that obviously, 𝜆 = 2 -𝑝+1+2𝑒1-𝑒2 exp(𝑎), 𝐶 = exp(𝑏 -𝑎).

We now discuss 𝑒 1 , 𝑒 2 . For exp([𝑎, 𝑏]) to be enclosed in a single binade of exponent 𝑒 2 , we need to have 2 𝑒2 exp(𝑎) exp(𝑏) < 2 𝑒2+1 , hence [𝑎, 𝑏] ⊂ [𝑒 2 log 2, (𝑒 2 + 1) log 2). For positive real numbers with 𝑒 1 -2, we have [2 𝑒1 , 2 𝑒1+1 ) ⊂ [0, log 2), hence we may deal with a whole binade at once and have 𝑒 2 = 0. The same is true for negative real numbers with 𝑒 1 -2, where we can take 𝑒 2 = -1. In both cases, we can take 𝑀 = 2 𝑝-1 .

For real numbers with 𝑒 1 -1, however, we have to split a given input binade into parts:

[2 𝑒1 , 2 𝑒1+1 ) = [2 𝑒1 , ⌈2 𝑒1 /log 2⌉ log 2) ∪ ⋃ ⌈ 2 𝑒 1 log 2 ⌉ ℓ< ⌊ 2 𝑒 1 +1 log 2 ⌋ [ℓ log 2, (ℓ + 1) log 2) ∪ [⌊ 2 𝑒1+1 /log 2 ⌋ log 2, 2 𝑒1+1 ) =:𝐼 0 (𝑒 1 ) ∪ ⋃ ℓ 𝐼 ℓ (𝑒 1 ) ∪ 𝐼 ∞ (𝑒 1 ).
Over 𝐼 0 (𝑒 1 ), we have 𝑒 2 = ⌊2 𝑒1 /log 2⌋, over 𝐼 ∞ (𝑒 1 ) we have 𝑒 2 = ⌊ 2 𝑒1+1 /log 2 ⌋ while over 𝐼 ℓ (𝑒 1 ) the corresponding exponent is ℓ. In this case, we are thus limited to values of 𝑀 of the order of 2 𝑝-1-𝑒1 log 2, or even smaller for 𝐼 0 and 𝐼 ∞ if 2 𝑒1 is close to a multiple of log 2.

A symmetric discussion holds similarly for negative 𝑥, where

𝐼 0 = [-⌈2 𝑒1 / log 2⌉ log 2, -2 𝑒1 ] gives 𝑒 2 = -⌈2 𝑒1 / log 2⌉, 𝐼 ℓ = [-(ℓ + 1) log 2, -ℓ log 2) gives 𝑒 2 = -1 -ℓ, 𝐼 ∞ = (-2 𝑒1+1 , - ⌊ 2 𝑒1+1 /
(𝑥) 𝑒 1 𝐼 𝑒 2 -2 -1 < 0 𝐼 0 -⌈2 𝑒1 /log 2⌉ -1 𝐼 ℓ -1 -ℓ 𝐼 ∞ -⌈2 𝑒1+1 /log 2⌉ -2 0 > 0 𝐼 0 ⌊2 𝑒1 /log 2⌋ -1 𝐼 ℓ ℓ 𝐼 ∞ ⌊2 𝑒1+1 /
𝜌 1 2.05 • 2 2𝑝/3 + 11.2 • 2 𝑝/2 + 10.1 • 2 𝑝/3 + 31.5 • 𝛿 1/2 2 𝑝/2 . Over [2 log 2, 2) we have 𝑎 ∈ [2 log 2, 2 log 2 + 2 1-𝑝 ] and 𝑏 = 2 -2 1-𝑝 , 𝑒 1 = 0, 𝑒 2 = 2; this gives 𝜆 ∈ [2 -𝑝+1 , 1.001 • 2 -𝑝+1 ], 𝐶 𝑒 2 /4, 𝑀 2 𝑝 (1 -log 2)
, and then

𝜌 2 3.4 • 2 2𝑝/3 + 26 • 2 𝑝/2 + 17.6 • 2 𝑝/3 + 36.7 • 𝛿 1/2 2 𝑝/2 .
Adding the two estimates yields the claimed result.

Experimental results

We now include similar tables as in Subsection 4.2. Those tables are obtained by using the decomposition above, and applying Theorem 3 for each of those intervals -mimicking the strategy of Corollary 3 (except that the slightly sharper Theorem 3 is used instead of the ready-to-use Lemma 4).

In order to demonstrate our bounds in various settings, we treat three different binades with a rather different behavior: [1, 2), [START_REF] Stehlé | On the Randomness of Bits Generated by Sufficiently Smooth Functions[END_REF]64) and [1/64, 1/32) in order to illustrate the behavior of our bounds for large/small operands of the exponential function. We do not illustrate our bounds for large, negative 𝑥 as in that case we get large negative 𝑒 2 , which yield a value of 𝜆 1 and thus a bound 𝑀 , hence empty.

Before going into this direction, we first mention that we stress-tested our bounds with the binade [2 -26 , 2 -25 ) where, depending on 𝑘 and the rounding mode, either none or all the floating-point numbers give rise to a hard-to-round case; as expected in this setting, our bounds give the trivial result that the number of bad cases is in [0, 2 𝑝-1 ].

Our results are collected in Table 3 and Table 4. We obtain, for instance, in binary128, for 𝑘 = 29, that the number of bad cases for exp over [START_REF]IEEE Standard for Binary Floating-Point Arithmetic[END_REF][START_REF]IEEE Standard for Radix Independent Floating-Point Arithmetic[END_REF] We conclude this discussion by a binade which is of a particular interest for the CRlibm library, as the latter starts by performing an argument reduction to restrict to the interval [0, log(2) • 2 -12 ) ; we thus give a few values of our bounds for the binade [2 -14 , 2 -13 ) in Table 5

Discussion

We now try to summarize the teachings that our experiments bring regarding our bounds.

The good situations

In all situations presented above, our experiments show that, as expected, our bounds provide orders of magnitude which go from very precise to reasonable when 𝑘 approaches 𝑝/3. As expected, the relative quality of our bounds improve a lot with 𝑝 for fixed 𝑘 -Remark 4 states that we expect the relative quality to be of the order 2 𝑘-𝑝/3 when input and output exponents are close to 0.

Until this paper only the probabilistic model (and no rigorous result -except for binary32 where complete enumeration is easy) was available. Note in particular that Tables 1 and 3 cover ranges for 𝑘 which correspond to the values of 𝑘 used commonly for the implementation of the quick phase of CRLibm functions -as such, demonstrate that results are sufficiently sharp to allow for a rigorous analysis of Ziv's method in that setting (see also Table 5). One might wonder also how our bounds work for MPFR [START_REF] Fousse | MPFR: A Multiple-Precision Binary Floating-Point Library with Correct Rounding[END_REF] implementation of Ziv's strategy; MPFR code suggests that the number of guard bits, depending on the function, is usually of the order of 𝑎 • log 2 (𝑝) + 𝑏 for 𝑎 ∈ {1, 2}, and not-too-large 𝑏; this corresponds to 𝛿 ≈ 2 -𝑏 𝑝 -𝑎 , a range where our bounds will be even better than for CRlibm except for very small intervals or very small/large output exponents, see below.

The bad situations

In this subsection we specifically target a few examples where we expect our bounds to behave poorly. In view of Remark 4 and Section 4.1, it should come as no surprise that the bounds behave poorly in the situations where • Phenomenon 1: 𝜆 is either larger than the expected 2 -𝑝 (which means either large 𝑒 1 or large negative 𝑒 2 )

• Phenomenon 2: or, conversely, 𝜆 much smaller than expected -in which case the term (𝛿/𝜆) 1/2 may grow wild (which means large negative 𝑒 1 or large 𝑒 2 ).

• Phenomenon 3: Also, if 𝑀 is much smaller than 2 𝑝-1 , the term, independent of 𝑀 , (𝛿/𝜆) 1/2 may become larger than the main term 2𝑀 𝛿 which has a linear dependency in 𝑀 ; to fix ideas, when exponents are close to zero, this is bound to happen whenever 𝑀 becomes smaller than 2 (𝑝+𝑘)/2 . Note that Phenomena 2 and 3 are related with the somewhat unavoidable term (𝛿/𝜆) 1/2 , see the end of Section 3.

Phenomena 1 & 2 do not occur for the cube root function, as we can always restrict to exponents in the range [0, 2].

Phenomenon 1 is obvious on the case of the binade (-64, -32] of the exponential function where the results are terribly bad: in binary64 with 𝛿 = 2 -15 , we get a (ridiculous) bound for the number of bad cases of the order of 2 128 , whereas the right order of magnitude is expected to be 2 38 .

Phenomenon 2 depends on the relative order of magnitudes of 𝑝 vs. 𝑒 1 , 𝑒 2 and would be somewhat visible on the binades [1/64, 1/32) and [START_REF] Stehlé | On the Randomness of Bits Generated by Sufficiently Smooth Functions[END_REF]64) in binary32. It becomes visible, for instance, in binary64 on the binade [2 -25 , 2 -24 ) where for 𝛿 = 2 -5 the bound is already trivial -and, in order to compare with previous tables, the relative radius is already ≈ 3660%.

What happens on the binade [START_REF] Stehlé | On the Randomness of Bits Generated by Sufficiently Smooth Functions[END_REF]64) in binary64 and binary128 is of a different nature and comes mostly from the fact that we have to split the interval (i.e., essentially Phenomenon 3) in many sub-intervals over which our bound is much worse. We come back to Phenomenon 2 in the final remark below.

Note that Phenomenon 1 corresponds to large (in absolute value) 𝑥 or zeros of 𝜙, and case Phenomenon 2 to small 𝑥 (again in absolute value) or 𝑥 close to a singularity of 𝜙. Note also that if 𝑒 1 and 𝑒 2 are simultaneously large or small with 𝑒 2 ≈ 2𝑒 1 , neither Phenomenon 1 nor 2 occurs; this may happen if 𝑥 is a zero of order 2 of 𝜙, or if 𝜙(𝑥) tends to infinity as 𝑥 2 .

It should probably be pointed that in the bad cases problem and Ziv's strategy, for most elementary functions the situations of very large/small input/output exponents are handled in a specific way, and we feel that those bad situations have thus little importance on the relevance of our results.

One further comment

Often in our tables, the relative radius seems to roughly get multiplied by 2 when 𝑘 increases by 1. It comes from the fact that in practice the 𝑀 𝜆 1/3 term of the error bound often seems to dominate over the term (𝛿/𝜆) 1/2 even for moderately small values of 𝛿. As the former term is independent of 𝛿, this means that the error bound is close to constant when 𝛿 varies, explaining why the percentages displayed double when 𝑘 increases by one -the error term does not change a lot while the main term decreases by a factor of 2.

Note that when Phenomenon 2 occurs (i.e. (𝛿/𝜆) 1/2 dominates) this no longer happens and a similar analysis shows that the relative radius then increases by a factor of √ 2 when 𝑘 increases by 1.

CONCLUSION

In this paper, we gave the first rigorous version of a heuristic commonly used by the designers of mathematical libraries.

To do so, we established some theoretical statements, based on so-called exponential sums techniques. Our results apply to any 𝒞 2 function and, though we only dealt with the radix 2 case, our results would remain valid for any other radix.

Our experiments demonstrate that, even though our bounds have their weaknesses, they provide sharp results for the analysis of the first phase of Ziv's strategy. They legitimate the two Ziv step strategy used for instance by the designers of the CRlibm library: it makes it possible to evaluate some elementary functions with correct rounding while keeping an average performance and a memory overhead due to correct rounding negligible [START_REF] De Dinechin | Towards the postultimate libm[END_REF], [START_REF] De Dinechin | Fast and correctly rounded logarithms in double-precision[END_REF].

The current main limit in our work is the fact that, if 𝑝 denotes the precision of the floating-point environment that we use, we can roughly estimate the numbers of bad cases up to 𝑘 ≈ 𝑝/3 bits but not beyond. Note that the nature of our techniques makes doubtful the fact that we shall ever be able to tackle the "worst cases" (which would correspond to 𝑘 ≈ 𝑝).

There are two reasons for the current 𝑘 ≈ 𝑝/3 limit: we lose the parameter 2 -𝑘 in some terms of the result stated in Lemma 2; and the use of Theorem 2. However, the latter has been shown to be optimal for 𝒞 2 functions [START_REF] Robert | On van der Corput's 𝑘-th derivative test for exponential sums[END_REF].

In the future, we shall focus on tackling some improvements of these two results, for instance in the more restrictive, and yet quite relevant and useful, framework of elementary and special functions. We shall also try to overcome to some extent some of the phenomena pointed in Subsection 4.4.2, in order for instance to be able to cope with large arguments of trigonometric functions.

Then

  |𝐸 -𝐸 * 𝐻 | 𝑀 ∑ 𝑚=1 |𝜓 (𝑓 (𝑚) -𝛿) -𝜓 * 𝐻 (𝑓 (𝑚) -𝛿) | + 𝑀 ∑ 𝑚=1 |𝜓 (𝑓 (𝑚) + 𝛿) -𝜓 * 𝐻 (𝑓 (𝑚) + 𝛿) |.

  sign

  1/2 , which may look somewhat singular in our bounds, has a very natural interpretation. Indeed, if 𝑓 is a polynomial of the form 𝑓 (𝑥) = 𝑎 + 𝑢 𝑣 𝑥 + 𝜆 2 𝑥 2 , with 𝑎, 𝑢, 𝑣 ∈ Z, 𝑣 small, then for any 𝑥 ∈ (-(2𝛿/𝜆) 1/2 , (2𝛿/𝜆) 1/2 ) divisible by 𝑣, we shall have ‖𝑓 (𝑥)‖ < 𝛿, which gives 2

	𝑣	( 2𝛿 𝜆	) 1/2 bad cases in
	this interval.

  tells us that, in binary64, if 𝑘 = 14, the number of bad cases for 3 √ • over [1/2, 1) lies between (1 -0.46) • 2 42 and (1 + 0.46) • 2 42 , which leads to a probability between (1 -0.46) • 2 -10 and (1 + 0.46) • 2 -10 for the occurrence of the second accurate step.

	𝑘	expected	rel. radius
	5	281474976710656	0.1%
	6	140737488355328	0.18%
	7	70368744177664	0.36%
	8	35184372088832	0.72%
	9	17592186044416	1.5%
	10	8796093022208	2.9%
	11	4398046511104	5.8%
	12	2199023255552	11.5%
	13	1099511627776	23%
	14	549755813888	46%
	15	274877906944	92%
		TABLE 1
	Our bounds for the cube root function in the binary64 case, 𝑥 ∈ [1/2, 1)

  ⌉ . We sum up the previous discussion concerning 𝑒 1 and 𝑒 2 in the following table.

			log 2 ⌋	log 2) gives 𝑒 2 =
	-	⌈	2 𝑒1+1 / log 2

  For the exponential function over the binade[START_REF]IEEE Standard for Binary Floating-Point Arithmetic[END_REF][START_REF]IEEE Standard for Radix Independent Floating-Point Arithmetic[END_REF], if 𝑝 11 the number of bad cases is in the interval [𝛿2 𝑝 -𝜌, 𝛿2 𝑝 +𝜌], where𝜌 5.45 • 2 2𝑝/3 + 37.2 • 2 𝑝/2 + 27.7 • 2 𝑝/3 + 68.2 • 𝛿 1/2 2 𝑝/2 .If, further, 𝑝 24, we have 𝜌 7.89• 2 2𝑝/3 + 68.2 • 𝛿 1/2 2 𝑝/2 . Proof: We split [1, 2) as [1, 2 log 2) ∪ [2 log 2, 2).Over [1, 2 log 2), we have 𝑎 = 1 and 𝑏 2 log 2, 𝑒 1 = 0, 𝑒 2 = 1; this gives 𝜆 = 2 -𝑝 𝑒, 𝐶 4/𝑒, 𝑀 2 𝑝 (log 2 -1/2) + 1 0.194 • 2 𝑝 for 𝑝 11. The first part of Lemma 4 then yields

	log 2⌋
	4.3.2 A simplified statement
	Again, we illustrate Lemma 4 to get a global bound over the
	binade [1, 2) of a compact form.
	Corollary 3.

  lies between (1 -0.014) • 2 84 and (1 + 0.014) • 2 84 , which leads to a probability between (1 -0.014) • 2 -28 and (1 + 0.014) • 2 -28 for the occurrence of the second accurate step.

	𝑘	expected	rel. rad. rel. rad.	rel. rad.
			[1, 2)	[32, 64)	[1/64, 1/32)
	5	281474976710656	0.09%	1%	0.02%
	6	140737488355328	0.17%	1.9%	0.02%
	7	70368744177664	0.34%	3.7%	0.03%
	8	35184372088832	0.67%	7.3%	0.05%
	9	17592186044416	1.3%	15%	0.09%
	10	8796093022208	2.7%	29%	0.2%
	11	4398046511104	5.4%	58%	0.3%
	12	2199023255552	11%	116%	0.6%
	13	1099511627776	22%	231%	1.1%
	14	549755813888	43%	462%	2.0%
	15	274877906944	86%	924%	3.9%
			TABLE 3
	Our bounds for the exponential function in the binary64 case
	𝑘	expected	rel. rad.	rel. rad.	rel. rad.
			[1, 2)	[32, 64)	[1/64, 1/32)
	8	4.0564819e31	6.4e-7%	6.9e-6 %	2.8e-8
	11 5.0706024e30	5.1e-6%	5.6e-5%	2.2e-7
	14 6.3382530e29	4.1e-5%	4.5e-4%	1.8e-6
	17 7.9228162e28	3.3e-4%	3.6e-3%	1.5e-5
	20 9.9035203e27	2.7e-3%	2.9e-2%	1.2e-4
	23 1.2379400e27	2.1e-2%	2.3e-1%	9.0e-4
	26 1.5474250e26	1.7e-1%	1.9%	7.2e-3
	29 1.9342813e25	1.4%	15%	5.8e-2
	32 2.4178516e24	11%	116%	0.46%
	35 3.0223145e23	86%	924%	3.7%
			TABLE 4
	Our bounds for the exponential function in the binary128 case

  .

	𝑘	rel. radius
		[2 -14 , 2 -13 )
	11	14%
	12	19%
	13	26%
	14	38%
	15	54%
		TABLE 5
	Some more CRlibm-specific bounds for the exponential function,
	binary64 format

ACKNOWLEDGEMENTS

The authors wish to warmly thank Christoph Lauter and Jean-Michel Muller for their very helpful answers to our various questions. Also, the careful rereading of the anonymous reviewers was a real help in improving the paper.

• This work was partly supported by the TaMaDi and FastRelax projects of the French Agence Nationale de la Recherche.


PhD Julien Zwang

Abstract

A needs assessment study was conducted including the participation of over 400 persons in three districts of Rakhine State, Myanmar. The study was composed of a representative quantitative survey of 334 persons with disabilities, and a qualitative study including eight focus group discussions and 20 key informant indepth interviews. Persons with disabilities were identified through snowball sampling stratified by age and sex using a reference population with disabilities. No study had ever been conducted on people with disabilities in this State. A majority (55%) of the respondents had a severe disability (could 'not do at all') and 99% 'had a lot of difficulties' to do at least one of the six categories of the Washington Group short set of questions to assess the type and the severity of the disability. The most frequent severe functional impairment was mobility (36% of the respondents). Higher psychological distress scores were detected in respondents under 50 years old, from Sittwe and having no access to bathroom facilities. Despite a primary health service coverage of 69%, the proportion of persons with disabilities who reported unmet needs in medical services was 75% (95%CI 69%-81%) and related to a lack of financial resources, expensive medical services, availability and distance. The persons with the most recent onset of disabilities and severe hearing impairment were the most at risks for not receiving appropriate medical care. In persons with some mobility impairment, 78% were not using any assistive device. In respondents aged of 15 years old or more, the unemployment rate was 78% and higher in females, persons with severe disability, and from Myauk-U district. Only 3% received a vocational training and they were 4 fold times more likely to currently work compared to other persons with disabilities. Only 3% of the persons with disabilities were receiving a pension, principally males over 50 years old. Overall, 81% (95%CI 77%-85%) of the respondents did not receive any assistance or social service within the last two years in particular in Myauk-U and Thandwe districts. The most frequent assistive devices or services needed to improve the function of the respondent were related to mobility (46%), physiotherapy, glasses, and hearing aids. The highest priorities reported for improving the living conditions were engaging in business, micro-credit, employment, physical rehabilitation and medication. In a context of extreme poverty and high unmet social needs and where females were the most vulnerable, forthcoming interventions should focus on persons with the most recent onset of disabilities as they have the most acute health and rehabilitation needs, access to bathrooms and clean water, food aid, inclusive vocational training programs and livelihood opportunities.

Study context

Rakhine State is believed to be one of the poorest regions of the Republic of the Union of Myanmar (Myanmar) and is marked by historical and ethno-political tensions. The development of Rakhine State has long lagged behind much of the rest of the country, leading to a situation in which both the local Rakhine and other ethnic minority communities have considerable unaddressed health and socio-economic needs. The 2012 intercommunal violence, resulting in the segregation between the Rakhine and Muslim communities, has an ongoing negative impact on the local political and economic situation and has aggravated the condition of the most vulnerable populations particularly persons with disabilities.

During the census conducted in Myanmar in 2014, the estimated prevalence rate of people with disabilities in Rakhine State was 5.3%, greater than the national average of 4.6%. However, a considerable segment of the population was left out of the exercise and was not enumerated because of ethno-political tensions and the demand of many local people to self-identify as Rohingya not conceded by the local authorities (UN, 2015).

So far, no study, assessment, or situational analysis had ever been conducted on persons with disabilities in Rakhine State. The limited knowledge about the situation and the needs of people with disabilities has prevented the implementation of inclusive interventions and/or targeted services.

As a result, there is a lack of advocacy initiatives addressed to policy makers, international donors, and humanitarian and development agencies to develop inclusive policies and programs for persons with disabilities in humanitarian and development actions.

Moreover, there is a very limited understanding of the protection risks that are faced by girls, boys, women, and men with disabilities in Rakhine State, especially due to its unique socio-political situation. Therefore, a needs assessment study on persons with disabilities in Rakhine State was required.

Objectives

The overall aim of the study was to assess the situation of persons with disabilities, their needs in relationship to basic services and identify barriers to access these services in selected districts of the Rakhine State.

Methods

Sampling methodology and case finding

The needs assessment of persons with disabilities in Rakhine State was achieved by conducting a quantitative and qualitative population based assessment study. Given the large size of the State, 5% of the country (Figure 1), the anthropological approach and technical non-probabilistic snowball sampling was used for participant identification [START_REF] Zwang | Social Context of Premarital Fertility in Rural South-Africa[END_REF]. This technique enabled us to explore a network of relationships from one local community contact. Once the person with disabilities was identified in the field, the quantitative questionnaire was administered. As in the Myanmar Census 2014 (Department of Population Ministry of Immigration and Population, May 2015), the Washington Group (WG) short set of questions was used to identify and to assess the type and severity of the functional limitations of the respondent (CDC, 2010).

Sampling for the eight focus discussion groups and the 20 key informant in-depth interviews, using four different tools, were made from a network of community based organizations (CBO). Opting for CBOs was appropriate insofar as these organizations play a mediation and information supply role for person with disabilities. The interviews were carried out in Burmese, Rakhine, and local languages used by the Muslim population, recorded and transcribed verbatim. In order to ensure the representativeness of the quantitative survey, a randomized sample stratified by age and sex using a reference population was employed. An age profile analysis on 4 million persons with disabilities with censuses data from five ASEAN countries (Cambodia, Indonesia, Malaysia, Philippines and Vietnam), with disability data, (IPUMS international) was conducted to define the age and sex profile of the population to be surveyed. The Myanmar census 2014 disability prevalence rates was used to select the districts to be surveyed.

Statistical analysis

Continuous data were reported as mean (standard deviation) or median with range (minimum, maximum) for non-normally distributed data. Categorical data were compared using the chisquare or the Fisher exact test as appropriate. The risks of the persons with disabilities were analyzed using multivariate logistic regression presented as adjusted risks (AOR). Multivariate linear regression was used for analyzing continuous variables. Confidence intervals (CI) were calculated at 95% (95%CI) and comparisons considered significant when p <0.05. Data were analyzed using Stata v13 (Stata corp.).

Ethics

Ethical approval for the study was sought and obtained from Handicap International's Technical Resource Department. Informed oral consent was obtained from all survey participants. Informed oral consent was obtained after explanation of the interview content. For people under 16 years old or those with intellectual impairment, consent was sought from a family member, who was present during all interviews. All interviews were anonymous and full names not recorded.

Results

Location

Three out of the five districts of the Rakhine State were surveyed (Figure 1). Overall, 31 Rakhine villages and 3 cities in each district were surveyed (70% of the respondents). Muslim communities customarily stay in the north of the Rakhine State and were interviewed (30% of the respondents) in 2 villages in Sittwe and 2 villages in Myauk-U districts and in one internally displaced persons (IDP) camp in Sittwe.

Respondents characteristics

The study, including training and data collection, was conducted between March and May 2016. The number of persons interviewed and included in the quantitative survey was 334 of whom 144 were from Sittwe district (43%), 124 from Myauk-U (37%), and 66 from Thandwe (20%). The total number of male was 169 and of female 165, resulting in a sex ratio (male/female) of 1.02 (Figure 2). The mean age of the population interviewed was 37 years ranging from 3 to 98 years old; without differentiation between sexes (p=0.849). Overall (excluding the IDP camp), the age distribution of the sample was representative of the population with disabilities in Rakhine State (p=0.210). The overall sex ratio (p=0.717) was not different from the target nor was the place of residence (86% surveyed from rural area compared to 83% targeted, p=0.350).

The eight focus group discussions (FGD), included groups with representation of persons with all types of impairment as well groups with specific criteria. One FGD was conducted with persons with visual impairment, one with people with mobility impairment, one in an IPD camp, two FGD were held in Myauk-U with one group representing the Rakhine community and one group representing the Muslim community and two FGD in Rakhine villages in Sittwe district. All FGD were gender balanced except one specific FGD held with women only. The 20 in-depth interviews were conducted in key informants from Nongovernmental Organizations (NGO), schools, and social action and health services.

Type and onset of disabilities

The mean age at onset of disability was 16 years old and the mean duration of disabilities was 21 years (Figure 3). The breakdown by cause of disability was disease (46%), congenital (27%), accident (21%), other (4%) and don't know (1%). No differentiation between sexes was detected (p=0.569). The proportion of respondents with chronic disease was 22% with no differentiation between sexes (p=0.803).

Figure 3: Lifespan with disabilities by type of impairment Overall, 55% of the respondents could 'not do at all' and 99% 'had a lot of difficulties' to do at least one of 6 categories of the Washington Group set of questions to assess the type and severity of the disability. No differentiation between sexes was detected (p=0.846, p=0.119, respectively). Respondents from Thandwe district were at lower risk for any severe disability (AOR 0.14) compared to other districts. Severe functional impairment ('cannot do at all') was 36% for walking, 12% for bathing, 12% for visual, 9% for communicating, 6% for hearing and 5% for remembering (Figure 4). Persons with poliomyelitis were at higher risks (AOR 2.69) for walking impairment. There was no difference observed in children under five years old compared to older respondents (p=0.395).

Female (AOR 2.31) compared to male and respondent over 35 years old (AOR 4.92) were significantly more at risk for severe visual impairment. No difference was detected for the other types of functional impairment.

Individuals' characteristics

In the 15 years old or more population, 52% of the respondents were married, no differentiation between sexes detected (p=0.801) and 47% were not able to read or write. Female (AOR 2.04) compared to male, persons with congenital disabilities (AOR 2.47) and persons with severe hearing impairment (AOR 2.05) were at higher risks for illiteracy.

Female, 27: I never attended school. My parents did not want me to go to school because they thought I could not walk properly and the school was located far away.

Livelihood

The average household size was 5 persons and the median monthly income of a household with a person with a disability was 45,000 MMK (around 38 USD) ranging from 0 to 200,000

MMK. As a result, the average income per person was 8 USD a month.

A

Psychological impact

Asking questions about five different types of psychological distress using a score ranging from never (0) to very often (4), the mean score was 5.5 ranging from 0 to 19. Respondents under 50 years old (p=0.005), from Sittwe (p=0.001) and persons having no access to a bathroom (p=0.001) had a significant higher psychological distress score compared to other groups.

94% of the respondents have never received any psychological support to deal with their current life situation. 2% received support from an NGO, 1% from a social worker, 1% from a peer support group and the other 2% from family, neighbors, psychologist or local organization for people with disabilities.

Persons with severe disability were more likely to receive any kind of support (AOR 3.26).

The main system to cope with disability was the family (51%), 'had no other option' (18%), 'overcame the situation alone' (16%) and 'my spirituality' (9%).

Female, 22: I feel depressed sometimes when I cannot do what I want. But when it happens my family usually encourages me.

Access to services

The proportion of person with disabilities reporting there was a school in their community was 81% (95%CI 77%-85%) at an estimated average 30 minutes distance; the proportion was 69% (95%CI 65%-74%) for primary health services at an estimated average of 40 minutes distance; 13% for social services at an estimated average 40 minutes distance; 8% for traditional healers at an estimated average one hour distance; 3% for psychological support at an estimated average more than one hour distance; 2% for vocational services at an estimated average 30 minutes distances and 0% for rehabilitation service.

Male, 15: I think that the medical service is not good for us and our Rakhine state is very poor.

Health services and rehabilitation

Overall, 75% (95%CI 69%-81%) of the persons with disabilities who needed any health facilities could not go there as much as needed. These unmet needs were composed of persons who could not attend at all any health facility and others that should have gone more often. 43% of the respondents went to any health facilities while 36% needed to go more often. Of the 57% of the respondents who did not go, 68% needed to attend a service.

The respondents most at risks for unmet needs in health facilities or to 'see the doctor' within the last 2 years were those with recent onset of disabilities (under 3 years, AOR 2.20; from 4 to 6 years, AOR 1.83) compared to those with longer onset, and those with severe hearing impairment (AOR 2.87), no other significant risks factors were detected.

Female, 65: Once I met a man at the bank and I received financial charity aid from him. With this money I bought a hearing aid from a shop. I never received anything from the hospital for my ears.

The most frequent assistive devices or services needed to improve the functional abilities of the respondent were related to mobility. 46% of respondents needed a cane, wheelchair, crutches or prostheses.

In persons with mobility impairment, 78% were not using any assistive device. Among those using a device (usually homemade), 65% used a cane, 22% used crutches, 7% used prostheses and one person used a wheel chair. The most at risk for not using any assistive device were female (AOR 2.14), aged below 50 years old (AOR 3.24) and persons living in Myauk-U (AOR 5.07) or Thandwe (AOR 3.48).

Female, 33: I think that the quality of medical care is not good. I have never heard that there is a rehabilitation center for people with disabilities.

The needs for a physiotherapist (35%), prostheses (24%), and social services (29%) were all mainly restricted by the availability (86%, 62%, and 78%, respectively).

Male, 52:

We have difficulty to access the hospital as it is far away. Before, we could not go to the hospital because of the conflict but after some time, we were allowed to go with permission of the authorities but it is still a burden.

Community participation

64% (95%CI 59%-69%) of the persons with disabilities felt abandoned or hated by their by their community and that life is harder compared to others because they cannot do things like before. The most at risk for these opinions were persons between 25 to 34 years of age (AOR 3.13), living in Sittwe (AOR 1.80), with severe disability (AOR 2.91), with a higher psychological distress score (AOR 1.20) and those who have unmet needs in health services (AOR 2.05).

The proportion of respondents who feel like other people look at them differently due to their impairment was 59%. The most at risk for this perception were under 35 years old (AOR 2.08) and with a higher psychological distress score (AOR 1.12).

Female, 33: Our community does not have superstitious beliefs anymore about persons with disabilities but my neighbors gossip about me and insult me because I am a woman with disabilities and I want to get married.

31% (95%CI 26%-36%) of the respondents reported that people use pejorative terms to refer to them because of their impairment. The most frequently discriminated persons were under 50 years old (AOR 2.60), from Sittwe (AOR 4.21) or the IDP camp (AOR 5.74) or with a higher psychological distress score (AOR 1.12).

Male, 48: I feel depressed when people make fun of me when I go out of my house. Sometimes, I want to give up my life.

In a list of four different types of community participation (family visits, community festivities, shopping at the market, and worshipping), the proportion of persons with disabilities never participating in any of all the four community events was 28% (95%CI 24%-33%). Persons with severe disabilities (AOR 7.13), aged of 25 years old or more (AOR 1.77) and females (AOR 1.95) were more at risk of never attending any of these community activities.

Female, 33: I never go to the festival with my husband because I am afraid that he will lose face when we walk together.

The proportion of respondents who ever had a leader role in the community was 16%. Females (AOR 7.81) and persons with congenital disabilities (AOR 2.45), compared to other type of disabilities, were less likely to have a role in the community. For those who never had a leader role, a majority (51%) thought they would never become a community leader.

Awareness about rights

In Rakhine State, 94% of the persons with disabilities had never heard about laws for persons with disabilities and 94% had never heard about the Convention on the Rights of Persons with Disabilities (CRPD).

Female, 25: It would be very good for us if we had a meeting place to talk together once a year.

The proportion of respondents who were member of an organization for people with disabilities for one year or more was 2%. However, the proportion increased to 38% for persons who were members less than one year, with the average membership duration being two weeks. This finding can be attributed to study participants deciding to register with the local disabled persons' organization and study partner, after being identified but prior to being interviewed for the study. The most frequent answers of the persons with disabilities for improving their living conditions were business or micro-credit (53%), employment (49%), physical rehabilitation (49%), medication (39%), others (mainly financial support, 17%) and education (16%).

Male, 43: I wish I had a permanent job, to stand on my own feet without relying on others. I feel having little value. I am depressed because I only do casual work.

Conclusions and recommendations

The

Gender issues

Overall, females with disabilities were more vulnerable than male with disabilities. They were at higher risks for illiteracy, unemployment, access to clean water, access to a bathroom, not receiving any assistance, lower community participation, not using assistive device, and having more financial difficulties to access medical services compared to male with disabilities. Focus on enabling women and girls to access services is recommended for all service providers, including Government, NGOs and CBOs.

Awareness

The study showed very low awareness about rights of people with disabilities (6%) and low affiliation to a CBO or DPO.

The creation and strengthening of self-help groups (SHG) should be applied to increase social cohesion, dissemination of information, and awareness-raising on the rights of persons with disabilities.

Mapping of available services within the State and regionally is recommended. Information campaigns on access to services would be an asset for persons with disabilities to promote access to health, rehabilitation, economic and education services.

Persons with disabilities must be included in the provision of services, especially services essential to basic needs (primary healthcare, food distribution, WASH). This includes the consideration of persons with disabilities in a humanitarian crisis. Physical and communication accessibility must also be considered.

The recent increase in the proportion of members from an organization for persons with disabilities was a direct result of the study design using snowball sampling technique. This indicates a high desire to engage with one's peers and participate in a social organization. Emphasis on creating greater social opportunities for persons with disabilities is recommended.

Community differences

The mean age of persons with disabilities in the IDP camp was significantly lower compared to other people with disabilities. Further studies should be conducted to explain this difference.

While people with disabilities from the IDP camp had a lower livelihood condition (lower item score) compared to others, they had better access to clean water, a bathroom, and food aid compared to other persons with disabilities in the Rakhine State.

Muslim communities were at higher risks for illiteracy, and Muslims from the villages had a higher need for mobility assistive devices.

Study limitations

Despite a randomization under age and sex constraint, the snowball sampling technique might have detection limits for people with sensorial and intellectual impairments. Without any robust population reference distribution by type of disabilities, it is difficult to distinguish whether the technique used here underestimated some types of disabilities or not. Therefore, the snowball sampling technique might be influenced by the community perception to identify persons with disabilities, and physical impairment related to mobility might have been overrated compared to sensorial and intellectual impairments.

Further analysis on the distribution by type of impairment should be conducted to further define the population with disabilities according to age, sex, and the local environment.

Figure 1 :

 1 Figure 1: Map of the districts surveyed, Rakhine State, Myanmar

Figure 2 :

 2 Figure 2: Respondents' age group distribution by sex

Figure 4 :

 4 Figure 4: Proportion of respondents with functional impairment by type and severity

Figure 5 :

 5 Figure 5: Priorities to improve living conditions, persons with disabilities, Rakhine State, Myanmar

  large majority of the respondents did not know what the monthly income was (85%) probably because most of the persons with disabilities do not work and do not generate income. Moreover, most of the household workers have casual work (farming) subject to seasonal variability or small scale activities for basic self-sustainability.In a list of 22 basic items in a household (air conditioner, bank account, bicycle, car, cattle, cellphone, chicken, couch, electricity, fan, fridge, internet connection, land, motorbike, landline phone, radio, rice cooker, solar power system, sewing machine, stove, table and chairs, TV), the mean number of items owned by the household was 3.4, ranging from 0.8 in the IDP camp, 1.8 in Myauk-U and, 3.8 in Sittwe to 6.6 in Thandwe.

	15% (95%CI 11%-19%) of the respondents were sometimes (3 81% (95%CI 77%-85%) of the respondents have not received any
	times a month or more) not able to eat for 24 hours within the assistance or social service within the last two years. The persons
	last month, the risks being higher in children under 15 years old the most in need of assistance were female (AOR 1.82), from
	(AOR 2.02) and persons with severe disabilities (AOR 4.21). Myauk-U (AOR 4.87) or Thandwe (AOR 5.71) as well as persons
	with less severe (AOR 2.49) compared to those with severe 78% of the persons with disabilities aged of 15 years or more did disability. 100% of the persons from the IDP camp are receiving not currently work. The most at risk for unemployment were the persons with a severe disability (AOR 4.27), female (AOR 2.04) a monthly ration of rice.
	and those from Myauk-U district (AOR 2.09).
	3% of the respondent aged of 15 years or more received a
	vocational training and they were almost 4 times more likely to
	currently work compared to those who never received
	vocational training (71% vs. 20%, respectively).
	Female, 22: I am interested in vocational training but I don't
	know what kind of vocational training is suitable for me.
	While 100% in the IDP camp had access to clean water, 34%
	(95%CI 29%-39%) of the other respondents did not have access 50% of the respondents aged of 15 years or more could not to clean water; the most at risk were female (AOR 2.41), as well engage in any household activities by themselves, in particular as respondents from Thandwe (AOR, 8.81) and Myauk-U (AOR persons with severe disability (AOR 1.94) and from Myauk-U 27.26) compared to Sittwe district. (AOR 2.01).
	51% (95%CI 45%-56%) of the respondents did not have access to 3% of the respondents receive a pension; in particular male (AOR a bathroom with the most at risk were female (AOR 1.58) and 10.05) and respondents over 50 years old (AOR 21.40) compared residents from Myauk-U (AOR 2.14). to other groups.

  needs assessment study of persons with disabilities in Rakhine State, Myanmar, used a representative sample stratified on sex and age with snowball sampling to identify persons with disabilities for a sample of over 400 persons (including both quantitative and qualitative surveys). The study design was also designed to identify and include persons with severe disabilities. The unmet needs in health and social services and high unemployment rate indicated that the persons with disabilities are prone to extreme poverty in Rakhine State.

	Unmet service needs Water, sanitation, hygiene (WASH)
	In Rakhine State, the unmet needs rate in health services for Half of the persons with disabilities did not have access to a
	persons with disabilities was 75%. The persons most in need or bathroom and the study showed that persons without access
	not receiving appropriate medical care were those with recent were at higher risk for psychological distress. In order to
	disabilities and with severe hearing impairment. decrease psychological distress, and increase access to
	bathrooms, further sanitation strategies should be implemented The training of DSW workers on disability identification and the with a focus on females and respondents living in Myauk-U needs of the most recently impaired persons, basic psychosocial district. Additionally, toilet, water points, and bathing areas support, cares and referral mechanisms should be implemented. must be physically accessible for persons with disabilities, Moreover, taking into consideration that persons with severe especially those persons with mobility issues, assistive devices, hearing impairment were at higher risk for illiteracy, DSW workers should be trained for sign language and access to and visual impairments.
	education and accessible communication formats/languages is
	recommended. Vocational training
	Despite a primary health services coverage of 69%, access to a
	Overall, awareness raising on the study results and subsequent training on disability inclusive development, with State level general practitioner, specialist, medication, physiotherapist, or While the persons with disabilities were almost 4 times more vocational training were restricted by a lack of financial likely to currently work if they had vocational training, the need resources, availability of the service and distance (most of the for vocational training is high but restricted by availability. The services are only available in Yangon for the whole country). access to vocational training was also frequently refused authorities (district, township and village) is recommended, 81% of the respondents did not receive any assistance or social because of the individual's disability. Forthcoming interventions especially as it relates to livelihoods and the access and inclusion service within the last two years. Upcoming assistance strategies should focus on awareness raising and inclusive policies for of persons with disabilities in educational and health programs. should focus on women and persons living in Myauk-U or persons with disabilities in vocational training programs in
	Thandwe districts where the aid was more scarce. particular those with a severe disability, female, and residing in
	Myauk-U district where the employment rate was the lowest.

The proportion of respondents in need of general health practitioners was 68% but access was limited by a lack of financial resource (32%), the availability of the service (24%) or the distance (19%). The need for medication was 64%, limited by a lack of financial resource (29%), availability (28%) or distance (13%); and the need for specialists was 55%, limited by the availability (52%), a lack of financial resource (23%) or the distance (13%).

Male, 17: According to doctors, my leg would be in a better situation if I was under regular medical care but because of financial difficulties my family cannot afford to follow the physician's recommendations.

The need for vocational training was 46%, restricted by the availability (52%) and because access was refused due to the impairment of the person (31%).

The need for school was 38% but the access was mainly refused because of the impairment of the respondents (73%).

Extreme poverty

The livelihood score of the persons with disabilities (in a list of basic items owned by the person with disabilities' household) and an average income per person of 8 USD a month revealed very low living conditions associated with no access to clean water (31%) or a bathroom (51%).

Regarding