# Formally Verified Roundoff Error Bounds on LogSumExp-based Computations

Paul Bonnot, Benoît Boyer, Florian Faissole, Claude Marché, Raphaël Rieu-Helft

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

## To cite this version:

Paul Bonnot, Benoît Boyer, Florian Faissole, Claude Marché, Raphaël Rieu-Helft. Formally Verified Roundoff Error Bounds on LogSumExp-based Computations. 2025. hal-04991131

HAL Id: hal-04991131

https://inria.hal.science/hal-04991131v1

Preprint submitted on 14 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Formally Verified Roundoff Error Bounds on LogSumExp-based Computations

Paul Bonnot¹, Benoît Boyer², Florian Faissole²,
Claude Marché¹, Raphaël Rieu-Helft³

¹Université Paris-Saclay, CNRS, ENS Paris-Saclay, Inria, LMF,
Gif-sur-Yvette, 91190, France.

²Mitsubishi Electric R&D Centre Europe, Rennes, France.

³TrustInSoft, Paris, 75014, France.

## Abstract

We study the numerical accuracy of some specific computer programs performing numerical computations. Such a numerical accuracy is expressed in terms of a bound on the difference between the floating-point computations and the corresponding rounding-free computations using mathematical real numbers. We do not only seek to discover such bounds “on paper” but we aim at obtaining computer-assisted formal proofs that these bounds are correct for any possible inputs. The functions we study come from the domain of machine learning: a function computing the logarithm of the sum of exponentials of a sequence, and a similar function related to the computation of *mutual information*. The bounds obtained are original results, because they are parameterized by assumed error bounds of the underlying implementations of the logarithm and exponential functions. The methodology we follow to conduct our formal proofs is also original, using a combination of the Why3 environment for deductive verification, an original modeling of floating-point computations using *unbounded* numbers, and the J³ environment for proving properties on C source code.

**Keywords:** Floating-Point, Mathematical Functions, LogSumExp, Roundoff Error, Formal specification and proof, Why3 environment, C programs

## 1 Introduction

Software is involved in many industrial systems nowadays. In cyber-physical systems in a broad sense, the software controlling a system must perform numerical computations, that are typically based on floating-point arithmetic. The floating-point representation of numbers, and the operations on them, are standardized by the IEEE-754 standard [45]. Despite of these rules, guessing the accuracy of a program that
<figure><img src="image_3.png" /><figcaption>Fig. 1 Graphical visualization of LSE (on the left) compared to the max function (on the right).</figcaption></figure>

compounds thousands of elementary operations, without software assistance, becomes almost impossible. However, recent history has shown that underestimating these errors could have catastrophic consequences [61]. That explains the recent interest for the formal verification of floating-point properties of numerical programs in proof assistants [12, 59] or deductive verification environments [18]. Formally proving the accuracy of floating-point computations is a complex topic addressed by different approaches in the scientific literature. Recent overviews of this topic can be found in the Handbook of Floating-Point Arithmetic [58], or surveys by Melquiond [54] and Boldo et al. [17].

In contrast with elementary operations, the accuracy of mathematical functions, say as the exponential and the logarithm as provided by the libm library for C code, is not enforced by the standard and could depend on the target architecture and on a specific library. The same goes for more complex mathematical algorithms, for which the accuracy strongly depends on the accuracies of underlying libm functions. For these reasons, there was very little effort to formally verify programs that call mathematical functions implementations.

An example of algorithm making use of libm functions is the Log-Sum-Exp algorithm, abbreviated as LSE. It applies to an n-dimensional vector $a = (a_0, ..., a_{n-1})$ and computes the logarithm of the sum of exponentials over its components:

$$ \text{LSE}(a) = \log \left( \sum_{0 \le i < n} \exp(a_i) \right) $$

It is typically used in machine-learning applications [24, 37, 56]. This algorithm is a smooth approximation of the max function, as illustrated on Figure 1 for a 2d vector, the smoothness property of inner functions being a pre-requisite for the efficiency of machine-learning algorithms. This function is frequently used in implementations of statistical classifiers [44, 50] and satisfies the property: $\max_{0 \le i < n}(a_i) \le \text{LSE}(a) \le \max_{0 \le i < n}(a_i) + \log(n)$. Blanchard et al. [10] showed that the roundoff error of a floating-point implementation of the LSE function can be bounded by relatively tight values as long as no overflow or underflow occurs. They present a pen-and-paper
proof taking advantage of both floating-point arithmetic specificities and mathemati-cal properties of the log and exp functions. They assume that the implementations of these two functions are correctly rounded, that is, their relative errors are bounded by the ε unit roundoff (see Section 2.1). Therefore, they provide the best accuracy that can be achieved in the considered floating-point format. Compared to Blanchard et al. [10], one of the contribution of this paper is to consider the possibility of using less accurate but more efficient implementations of exp and log. Such implementations may be useful in contexts involving energy-saving small devices like IoT [40]. In this paper we also consider the other function, abbreviated as SLSE, defined as

$$ \mathrm{SLSE}(a, \rho) = \sum_{0 \le i < n} \log_2 \left( \sum_{0 \le j < n} \exp \left( - \frac{(a_i + \rho - a_j)^2}{2} \right) \right) $$

where ρ is some extra parameter, and log₂ is the logarithm in base 2. This is one of the functions involved in computing mutual information [44, 50]. As for LSE, we provide bounds on the roundoff error of implementations of this function.

Another contribution of this paper is to give formal proofs of our results. Generally speaking, obtaining formal proofs on the accuracy of floating-point programs is not a simple task. One can start from a software environment for proving functional properties of programs, and augment it with a formalization of floating-point arithmetic, typically via a library built on top of a formalization of real numbers that provides a rounding function and its logical properties. It was done for example in 2007 by Boldo and Filliâtre [16] using the Rocq proof assistant (formerly Coq), augmented with the Flocq library [20] for floating-point arithmetic, allowing to prove programs using the Rocq general-purpose environment. These proofs typically require a large amount of manually written proof steps. To obtain a higher degree of automation, Ayad and Marché [5] proposed a setting making use of the Frama-C [48] environment for static analysis on C source code, with a dedicated library of ACSL [9] specifications that allows to discharge the proofs to various theorem provers, in particular SMT solvers but also the Gappa [30] solver specialized in floating-arithmetic. Boldo and Marché [18] presented an overview of what could be achieved on increasingly complex codes, using combination of automated solvers and the Rocq proof assistant for the most complex proof obligations. With the addition, later on, of some built-in support for floating-point operations in SMT solvers, Fumex et al. [36] showed that this methodology can reach a fairly high amount of automation.

Initially we planned to follow the methodology above to directly prove a C code for LSE. Yet, to achieve the proofs in a reasonably simple manner, we identified two important aspects that should be emphasized: first, the use of the Why3 verification environment and its intermediate language WhyML, and second the use of unbounded floating-point numbers. This article is an extended version of our FMCAD paper [23], which already contained the accuracy result for LSE, which we now extend also to the SLSE function. We also give the pen-and-paper proofs that were not given in the FMCAD paper, and present more details on the WhyML formalization. The rest of this paper is organized as follows. Section 2 is a preliminary section where we introduce the notion of unbounded floating-point numbers (in Section 2.2), introduce the
Why3 setting for formalization and illustrate it on unbounded floats (Section 2.3), expose some first accuracy results about compound summations in Section 2.4 and formalize them in Section 2.5. In Section 3 we expose our results bounding the accuracy of LSE, parameterized by assumed accuracies of implementations of exp and log. Their formalization is shown in Section 3.4. We then derive a formal proof of a C code computing LSE in Section 3.5. Section 4 is devoted to the SLSE function, whose accuracy result requires to state and prove a few preliminary lemmas. We conclude in Section 5 with a discussion on related work and future work. The code formalizing our results is publicly available on the Toccata gallery of verified programs [51].

# 2 Floating-Point Numbers and Accuracy of Summations

In this preliminary section we first give a quick overview of floating-point arithmetic (Section 2.1) before introducing our notion of unbounded floating-point numbers in Section 2.2. Towards the formalization, we introduce WhyML briefly in Section 2.3 and present how floating-point arithmetic is formalized, including the unbounded version. Section 2.4 presents accuracy results about compound summations, and Section 2.5 shows how these results are formalized and proved.

## 2.1 Floating-Point Arithmetic in a Nutshell

The IEEE-754 standard [45] defines several formats of representation of floating-point numbers. A format is characterized by a precision *p* as well as upper and lower bounds $e_{max}$ and $e_{min}$ for the exponent. A floating-point number is either a value among $+\infty$, $-\infty$, and NaN, or a value $\pm m \times 2^{e-p+1}$ where $m, e \in \mathbb{Z}$, $0 \le m \le 2^p - 1$ and $e_{min} \le e \le e_{max}$. The largest representable number in this format is $\text{maxf} = (2 - 2^{-p-1}) \times 2^{e_{max}}$, and the smallest positive representable number is $2^{e_{min} - p+1}$.

In this paper we are not interested in a particular format since our proof methodology is independent of the format used, however only two formats are currently supported in our formal proofs: single format (32 bits) where $p = 24$, $e_{max} = 127$ and $e_{min} = -126$, and double format (64 bits) where $p = 53$, $e_{max} = 1023$ and $e_{min} = -1022$. We use the symbol *rnd* to denote the *rounding* of a real number to a floating-point number. The IEEE-754 standard defines several rounding modes. Here we consider only the mode *nearest-ties-to-even*: when a real number *x* lies within an interval $[x_1; x_2]$ of two consecutive floating-point numbers, then *rnd*(*x*) is either $x_1$ or $x_2$: the one of these which is closest to *x*, or in case *x* is exactly in the middle, the one among $x_1$ and $x_2$ whose mantissa is even. Also, when *x* is too large (larger than or equal to the middle of $\text{maxf}$ and $2 \times 2^{e_{max}}$), *rnd*(*x*) is $+\infty$. We use the symbols $\oplus$, $\ominus$, $\otimes$, $\oslash$ to denote the basic operations of addition, subtraction, multiplication and division of floating-point numbers. As specified by IEEE-754, all these operations must use the best possible rounding, that is $x \oplus y = \text{rnd}(x + y)$ and similarly for the three other operations.
The main property of the rounding function that we use in this paper is the following: for any real number $x$ such that $|x| \le \max f$, $\text{rnd}(x)$ is finite and

$$|\text{rnd}(x) - x| \le \varepsilon|x| + \eta \quad (1)$$

where $\varepsilon = \frac{2^{-p}}{1+2^{-p}}$ and $\eta = 2^{e_{\min}-p}$. This property can be considered as well-known and folklore in the literature, see for example Jeannerod and Rump [46]. In seminal publications, such as the Handbook of Computer Arithmetic [58] or Higham's survey [43], the simpler term $\varepsilon = 2^{-p}$ is used instead of $\frac{2^{-p}}{1+2^{-p}}$, inducing a slightly larger bound. Jeannerod and Rump [46, Theorem 2.1] showed that the refined bound is actually optimal in the sense that there exist some inputs values and floating-point formats (with certain conditions) for which it is attained. In most cases, the precision gain obtained using this optimal bound instead of $2^{-p}$ is small. Anyway, the latter results and proofs simply use the symbol $\varepsilon$ to denote either of the bounds.

As remarked by Jeannerod and Rump [46], Property (1) can be refined in the special case of addition because underflowing additions are exact:

$$|(x \oplus y) - (x + y)| \le \varepsilon|x + y| \quad (2)$$

that is, the term $\eta$ can be removed from Formula (1). Moreover, it should be noted that (see for example the Handbook [58])

$$|(x \oplus y) - (x + y)| \le |x| \quad (3)$$

and symmetrically

$$|(x \oplus y) - (x + y)| \le |y| \quad (4)$$

The combination of the properties (2), (3) and (4) is used later on to obtain bounds on compound summations.

## 2.2 Unbounded Floating-Point Numbers

The notion of unbounded floating-point number is somewhat simple, and is in fact not original: it is commonly used in the literature on numerical programs [58] and also in advanced formalization such as Flocq [20]. Roughly speaking, unbounded floating-point number are very much like standard IEEE floating-point numbers, except that their exponent can be arbitrarily large: it is a value $\pm m \times 2^{e-p+1}$ where $0 \le m \le 2^p-1$, $e_{\min} \le e$ without upper bound on $e$. As a consequence, unlike standard floating-point numbers, the four basic operations on unbounded floating-point numbers *never overflow*. There is no need for special values for infinities to represent the result of unbounded floating-point operations. On the other hand, notice that unbounded floats include sub-normal numbers. There is an injection from finite IEEE float numbers to unbounded float numbers. The properties (2), (3) and (4) indeed hold for unbounded floating-point numbers. This fact allows us to separate the proofs concerning functional behavior of numerical programs from the proof of absence of overflow: to prove a property on floating-numbers it suffices to prove the same on unbounded floats, and
separately prove that each floating-point operation involved does not overflow nor
produces NaN values.

## 2.3 WhyML Formalization of Unbounded Floating-Point Arithmetic

Our formal proof methodology makes use of the unbounded floating-point numbers introduced in Section 2.2 and heavily relies on WhyML. WhyML is the language of Why3 [11], a general-purpose environment for deductive verification. Why3 is used as an intermediate tool by several front-ends including Frama-C [48] for C code, Spark for Ada code [53] and Creusot [35] for Rust code. The Why3 environment allows the user to access a large set of different provers, including the SMT solvers which support the SMT-LIB floating-point theory, such as CVC4 [7], cvc5 [6], Z3 [34] and Alt-Ergo [26]. Why3 also provides to discharge proofs using Rocq and solvers dedicated to real arithmetic and floating-point arithmetic, such as Gappa but also dReal [38] and Metitarski [1]. For the formal proofs of this paper it appears we do not rely upon the latter solvers dedicated to arithmetic, but only use the SMT solvers. Yet, relying only on SMT solvers to achieve the proofs is not possible for the complex steps: Why3 provides the user with a large set of techniques and tools to achieve complex proofs, for example via the use of lemma functions, which are, roughly speaking, a way to construct a proof by writing a program. This approach is heavily used below.

A starting point of our formalization was to design a new WhyML theory for unbounded floats. An excerpt of that theory is given in Figure 2. In this theory, the type `udouble` is abstract, and only assumed to be given a function `to_real` which returns the real number represented by any `udouble`. The rounding function `around` that rounds any real number to a `udouble` is also declared abstractly, axiomatized by property (1). The basic operations are defined as the rounding of the real operations. The properties (2), (3) and (4) are stated as axioms in the theory. This library is part of the standard library of Why3 since its 1.8 release. To provide guarantees that this theory is consistent with IEEE floats, it is realized using Rocq and its Flocq library.

## 2.4 Accuracy of Compound Summations

The compound sum of a vector $(a_0, \dots, a_{n-1})$ of floating-point numbers is the sum of all $a_i$. Defining it properly is more complex than the compound sum of real numbers because $\oplus$ is not associative. It is thus necessary to choose the order in which the additions are done. We make the choice to associate to the left, meaning that we define the compound sum from $a_0$ (included) to $a_k$ (excluded), denoted by $\bigoplus_{0 \le i < k} a_i$, by the following recursive equations:
$$ \bigoplus_{0 \le i < k} a_i = 0 \text{ if } k \le 0, \text{ and } \bigoplus_{0 \le i < k} a_i = \left( \bigoplus_{0 \le i < k-1} a_i \right) \oplus a_{k-1} $$
when $0 < k$. Associating to the left is important because it impacts the final result of a sum. Yet the bounds we prove in the following are invariant by permutation of the element of the input vector. In other words, the same bounds could be proved when the sum is performed in any other order.

The following theorem states a bound on compound summations. It is not a new result of us but only a slight reformulation of a theorem by Jeannerod and Rump [46].
module UFloat

type u
  (** The type of unbounded floats *)
function to_real u : real (** injection of ufloats to real numbers *)

function uround mode real : u (** The rounding function *)
constant eps:real          (** axiomatized by property (1) *)
constant eta:real
axiom uround_rne:
forall x:real. abs (uround RNE x - to_real x) ≤ eps × abs x + eta

function uadd (x y:udouble) : udouble (** addition *)
= uround RNE (to_real x + to_real y)

(** properties (2), (3) and (4) *)
axiom add_rounding :forall x y:udouble.
abs (to_real (uadd x y) - (to_real x + to_real y))
≤ abs (to_real x + to_real y) × eps
axiom add_bound_left:forall x y:udouble.
abs (to_real (uadd x y) - (to_real x + to_real y)) ≤ abs (to_real x)
axiom add_bound_right:forall x y:udouble.
abs (to_real (uadd x y) - (to_real x + to_real y)) ≤ abs (to_real y)
[...]
end (* module UFloat *)

module UDouble
constant eps_double : real = 0x1p-53 / (1.0 + 0x1p-53)
constant eta_double : real = 0x1p-1075
clone UFloat with eps = eps_double, eta = eta_double
end

<caption>Fig. 2 Generic theory of unbounded doubles in WhyML (excerpt): conversion to real numbers, rounding, addition and its properties. It is instantiated to unbounded doubles with particular values for ε and η.</caption>

**Theorem 1** (Accuracy of compound summation) For any vector a of unbounded doubles, and any n ≥ 0:

$$
\left| \bigoplus_{0 \le i < n} a_i - \sum_{0 \le i < n} a_i \right| \le (n-1)\varepsilon \sum_{0 \le i < n} |a_i|
$$

*Proof* The proof is done by induction on *n*. The base cases *n* = 0 or 1 are trivial since the summations are both equal to 0 and *a*₀ respectively, without any rounding involved. For the inductive case, we assume *n* > 1 and we assume the property holds for *n* − 1. We first prove the additional property

$$
\left| \left( \bigoplus_{0 \le i < n} a_i \right) - \left( \bigoplus_{0 \le i < n-1} a_i + a_{n-1} \right) \right| \le \varepsilon \left( \sum_{0 \le i < n-1} |a_i| + (n-1)|a_{n-1}| \right) \quad (5)
$$
by distinguishing two cases: If $|a_{n-1}| \le \varepsilon \sum_{0 \le i < n-1} |a_i|$ : by Property (4) we have

$$ \left| \left( \bigoplus_{0 \le i < n} a_i \right) - \left( \bigoplus_{0 \le i < n-1} a_i + a_{n-1} \right) \right| \le |a_{n-1}| \le \varepsilon \sum_{0 \le i < n-1} |a_i| \quad \text{which proves Formula (5). Otherwise, that is when } \varepsilon \sum_{0 \le i < n-1} |a_i| < |a_{n-1}|, \text{ we use Property (2) to state that} $$

$$ \left| \left( \bigoplus_{0 \le i < n} a_i \right) - \left( \bigoplus_{0 \le i < n-1} a_i + a_{n-1} \right) \right| \le \varepsilon \left| \bigoplus_{0 \le i < n-1} a_i + a_{n-1} \right| \quad \text{and we have} $$

$$ \begin{aligned} \left| \bigoplus_{0 \le i < n-1} a_i + a_{n-1} \right| &= \left| \bigoplus_{0 \le i < n-1} a_i - \sum_{0 \le i < n-1} a_i + \sum_{0 \le i < n-1} a_i + a_{n-1} \right| \\ &\le \left| \bigoplus_{0 \le i < n-1} a_i - \sum_{0 \le i < n-1} a_i + \sum_{0 \le i < n-1} |a_i| + |a_{n-1}| \right| \quad \text{(triangular inequality)} \\ &\le (n-2)\varepsilon \sum_{0 \le i < n-1} |a_i| + \sum_{0 \le i < n-1} |a_i| + |a_{n-1}| \quad \text{(induction hypothesis)} \\ &\le (n-2)|a_{n-1}| + \sum_{0 \le i < n-1} |a_i| + |a_{n-1}| \quad \text{(case assumption and } n-2 \ge 0) \\ &= \sum_{0 \le i < n-1} |a_i| + (n-1)|a_{n-1}| \end{aligned} $$

which proves Formula (5). We obtain the final result as follows, starting by applying triangular inequality

$$ \left| \bigoplus_{0 \le i < n} a_i - \sum_{0 \le i < n} a_i \right| \le \left| \left( \bigoplus_{0 \le i < n} a_i \right) - \left( \bigoplus_{0 \le i < n-1} a_i + a_{n-1} \right) \right| + \left| \left( \bigoplus_{0 \le i < n-1} a_i + a_{n-1} \right) - \sum_{0 \le i < n} a_i \right| $$

and then on the left part we use Formula (5) and on the right part, the term $a_{n-1}$ cancels so we can apply the induction hypothesis, hence

$$ \begin{aligned} \left| \bigoplus_{0 \le i < n} a_i - \sum_{0 \le i < n} a_i \right| &\le \varepsilon \left( \sum_{0 \le i < n-1} |a_i| + (n-1)|a_{n-1}| \right) + (n-2)\varepsilon \sum_{0 \le i < n-1} |a_i| \\ &= \varepsilon \left( (n-1) \sum_{0 \le i < n-1} |a_i| + (n-1)|a_{n-1}| \right) = \varepsilon(n-1) \sum_{0 \le i < n} |a_i| \end{aligned} $$

□

From the previous theorem, we deduce the following corollary, which will be useful in the proofs we perform afterwards, in particular for proving absence of overflows.

**Corollary 2 (Bound on summations)** For any constant $M_a$, any vector a such that $|a_i| \le M_a$ for all $0 \le i < n$ we have

$$ \left| \bigoplus_{0 \le i < n} a_i \right| \le n(1 + \varepsilon(n-1))M_a $$

*Proof* From Theorem 1 we get $\left| \bigoplus_{0 \le i < n} a_i \right| \le \left| \sum_{0 \le i < n} a_i \right| + (n-1)\varepsilon \sum_{0 \le i < n} |a_i| \le nM_a + (n-1)\varepsilon nM_a$ and conclude easily. □

Theorem 4 below is our generalization of Theorem 1 to the case where the vector a itself is already an approximation of a real vector, with some assumed roundoff error. We start with an auxiliary lemma on summations of real numbers.
let rec function u_sum (a: int → udouble) (n: int) : udouble
variant { n }
= if n <= 0 then uzero else uadd (u_sum a (n-1)) (a (n-1))

<caption>Fig. 3 WhyML definition of compound summations of unbounded doubles</caption>

**Lemma 3** (Composition of errors in a sum) *Given any vectors a and $\hat{a}$ and any parameters A and B such that for all i, $|\hat{a}_i - a_i| \le A|a_i| + B$ we have:*

$$ \left| \sum_{0 \le i < n} \hat{a}_i - \sum_{0 \le i < n} a_i \right| \le A \sum_{0 \le i < n} |a_i| + nB $$

*Proof* Straightforward by induction on n. □

**Theorem 4** (Accuracy of summations, generalized) *Given any vectors a and $\hat{a}$ such that for all i, $|\hat{a}_i - a_i| \le A|a_i| + B$ we have:*

$$ \left| \bigoplus_{0 \le i < n} \hat{a}_i - \sum_{0 \le i < n} a_i \right| \le (A + (n-1)\varepsilon(1+A)) \sum_{0 \le i < n} |a_i| + Bn(1+(n-1)\varepsilon) $$

*Proof* By Theorem 1, we have $\bigoplus_{0 \le i < n} \hat{a}_i \le \sum_{0 \le i < n} \hat{a}_i + (n-1)\varepsilon \sum_{0 \le i < n} |\hat{a}_i|$. By Lemma 3 we have $\sum_{0 \le i < n} \hat{a}_i \le \sum_{0 \le i < n} a_i + A \sum_{0 \le i < n} |a_i| + nB$ and $\sum_{0 \le i < n} |\hat{a}_i| \le \sum_{0 \le i < n} |a_i| + A \sum_{0 \le i < n} |a_i| + nB$. Therefore $\bigoplus_{0 \le i < n} \hat{a}_i \le \sum_{0 \le i < n} a_i + A \sum_{0 \le i < n} |a_i| + nB + (n-1)\varepsilon \left( \sum_{0 \le i < n} |a_i| + A \sum_{0 \le i < n} |a_i| + nB \right) \le \sum_{0 \le i < n} a_i + (A + (n-1)\varepsilon(1+A)) \sum_{0 \le i < n} |a_i| + Bn(1+(n-1)\varepsilon)$ which proves the upper bound, a similar reasoning works for the lower bound. □

## 2.5 Formalization of Accuracy Results on Summations

We start by formally defining the compound operator $\oplus$ in WhyML. Given that its mathematical definition above is recursive, the idiomatic way to define it in WhyML is defining this function as a recursive program that we prove terminating, as shown in Figure 3. Note that the function here is defined on our new type `udouble` for unbounded doubles. The vector is provided as a function `a` from integers to floating-point numbers. We make use here of the higher-order features provided by Why3.

The WhyML statement corresponding to Theorem 1 is given by the function `u_sumookv` of Figure 4, with attribute `ghost` telling it is a lemma function, not intended for program execution but for writing proofs. Note how the pre-condition and the post-condition correspond respectively to the hypothesis and the statement of the theorem. It means that if someone provides an implementation of that function and proves that this implementation satisfies the contract, then for all values of the parameters, that is for all `a` and `n`, satisfying the pre-condition, the post-condition holds: this is exactly the statement of the expected accuracy property.
```h
function f_to_real (a:int → udouble) : int → real
    = fun i → to_real (a i)
function f_abs_to_real (a:int → udouble) : int → real
    = fun i → abs (to_real (a i))

let rec ghost u_sumendor (a:int → udouble) (n:int)
    requires { 0 <= n }
    variant { n }
    ensures { abs (to_real (u_sum a n) - sum (f_to_real a) 0 n)
               ≤ from_int (n-1) × eps × sum (f_abs_to_real a) 0 n }
```

<caption>Fig. 4 WhyML statement corresponding to Theorem 1.</caption>

We provide a body for this function as shown in Figure 5. We use recursive calls so as to simulate a proof by induction. The structure of the code follows quite closely the structure of the pen-and-paper proof of Theorem 1 of Section 2.4.

The WhyML source code can then be fed to Why3, and provers executed on the generated VCs. We run in parallel the provers Alt-Ergo 2.6.0, CVC4 1.8, cvc5 1.2.0 and Z3 4.12.2, with a time limit of 5 seconds. On each VC we retain the prover which answer the fastest. For Theorem 1, whose statement is shown in Figure 4 and proof is shown on Figure 5, we have 11 VCs, 4 are proved by Alt-Ergo, 3 by CVC4 and 4 by cvc5. For Corollary 2, we have a total of 9 VCs, 4 proved by Alt-Ergo, 1 by CVC4 and 4 by cvc5. For Theorem 4, we have a total of 9 VCs, 6 proved by Alt-Ergo, 2 by cvc5 and 1 by Z3. The details on which prover is proving which VC can be found on the dedicated page of Toccata gallery [51].

# 3 Accuracy of Log-Sum-Exp

## 3.1 Parametric Approximations of exp and log

In implementations using floating-point numbers, not only the sum is subject to rounding, but also the computations of functions exp and log. Here, we do not discuss any particular implementations of these two functions. We do not want to rely, as Blanchard et al. [10] do, on perfectly rounded implementations of exponential and logarithm. Instead we assume we have implementations that are possibly less precise, the precision of them being specified as parameters.

Concerning exponential first, we assume given an implementation $\widehat{exp}$ which satisfies the following property: for any real $x$ such that $|x| \le M_{exp}$, $|\widehat{exp}(x) - \exp(x)| \le E_{exp} \exp(x)$ where $M_{exp}$ and $E_{exp}$ are two positive parameters. In the following we need to assume $E_{exp} \le 0.5$, a reasonable assumption, which in particular implies that $\widehat{exp}(x)$ is always non-negative. Concerning logarithm we assume similarly an implementation satisfying the following property: for any real $x$ such that $0 < x \le M_{log}$ $|\widehat{log}(x) - \log(x)| \le E_{log} |\log(x)|$ where $M_{log}$ and $E_{log}$ are positive parameters.

Notice that we do not claim that there exist implementations of approximations of exponential and logarithm, satisfying the properties above, for any value of the parameters $E_{exp}$, $M_{exp}$, $E_{log}$ and $M_{log}$. We just assume we are given some. Indeed
```c
if n <= 1 then (* base cases, trivial *) return ();
(* induction *)

u_sumendor a (n-1); (** recursive call to obtain induction hypothesis *)
let exact_sum = sum (f_to_real a) 0 (n-1) in (* sum in real numbers *)
let abs_sum = sum (f_abs_to_real a) 0 (n-1) in (* sum of absolute values *)
begin ensures { abs exact_sum ≤ abs_sum }
    sum_abs (f_to_real a) 0 (n-1);
end;

let psum : udouble = u_sum a (n-1) in (* the previous sum, in udouble *)
let res : udouble = uadd psum (a (n-1)) in (* the final sum, in udouble *)
let s : real = to_real psum in (* the previous sum, in real *)
let anm1 = to_real (a (n-1)) in (* the last element as a real *)
let delta = to_real res - (s + anm1) in (* the difference to bound *)
begin (* additional property (5) *)
    ensures { abs delta ≤ eps × (abs_sum + from_int (n-1) × abs anm1) }
    if abs anm1 ≤ eps × abs_sum then (** case 1 of the proof, trivial *) () else
        begin
            (* case 2 of the proof *)
            (* we know eps * abs_sum < abs anm1 *)
            assert { abs (s + anm1) ≤ abs (s - exact_sum) + abs abs_sum + abs anm1 };
            assert { from_int (n-2) × (eps × abs_sum) ≤ from_int (n-2) × abs anm1 }
        end
end;

assert { abs (to_real (u_sum a n) - sum (f_to_real a) 0 n)
           ≤ abs delta + abs (s + anm1 - sum (f_to_real a) 0 n) }
```

<caption>Fig. 5 WhyML body of `u_sumendor`, corresponding to the proof of Theorem 1.</caption>

it is known in the literature that such implementations exist for double<fnref n="1" /> precision, with a correct rounding, that is with $E_{\text{exp}} = E_{\text{log}} = \varepsilon$, $M_{\text{log}} = \text{maxf}$, and $M_{\text{exp}}$ at most 708 (for larger values the exponential overflows): see for example Daramy et al. [29] and the implementations provided by the CORE-MATH project [60]. Notice also that we assume only some relative error ($E_{\text{exp}}$ and $E_{\text{log}}$) but no absolute error. For exponential, this is not needed because for an argument at least $-708$ the result is never a sub-normal. For a completely different reason, the logarithm do not need to return any sub-normal either, because the logarithm of the floating-point successor of 1 is around $2^{-52}$, larger than a sub-normal too (and similar for the predecessor).

## 3.2 Analysis of Accuracy of LSE

Our main result concerning the accuracy of the computation of $\tilde{\text{LSE}}$ is given by Theorem 6 below. See Section 5 for a comparison with the previous result of Jean-nerod and Rump. To prove this theorem we need to establish first an auxiliary lemma telling how errors are propagated through the logarithm.

<footnote n="1">¹for simple precision, $M_{\text{exp}}$ should be at most 88.</footnote>
**Lemma 5 (Error propagation through logarithm)** For any positive real numbers $x$ and $\hat{x}$ such that $|\hat{x}-x| \le Ax$, for some parameter $A$ with $0 \le A < 1$ we have $|\log \hat{x} - \log x| \le -\log(1-A)$.

*Proof* We have $x(1-A) \le \hat{x} \le x(1+A)$, so $\log(x(1-A)) = \log(x) + \log(1-A) \le \log(\hat{x}) \le \log(x(1+A)) = \log(x) + \log(1+A)$. Since $0 \le A < 1$, we have$^2 \log(1+A) \le -\log(1-A)$, therefore $\log(x) + \log(1-A) \le \log(\hat{x}) \le \log(x) - \log(1-A)$. $\square$

**Theorem 6 (Accuracy of LSE)** For any $n \ge 1$, such that $(n-1)\varepsilon \le \frac{1}{4}$, for any vector $a$ of size $n$ such that for all $i$, $|a_i| \le M_a$ for some $M_a \le M_{\text{exp}}$, and assuming that

$$n(1 + \varepsilon(n - 1)) \exp(M_a)(1 + E_{\text{exp}}) \le M_{\log} \quad (6)$$

we have

$$|\widehat{\text{LSE}}(a) - \text{LSE}(a)| \le E_{\log}|\text{LSE}(a)| - \log(1 - (E_{\text{exp}} + (n-1)\varepsilon(1 + E_{\text{exp}}))(1 + E_{\log}))$$

As shown in the proof below, the hypothesis (6) above is required to call the $\widehat{\log}$ function on the proper interval of definition. A bound on the size of the input, relatively to $\varepsilon$, is needed to apply Lemma 5 with $A = E_{\text{exp}} + (n-1)\varepsilon(1 + E_{\text{exp}})$: to show that $A$ is smaller than 1, together with the hypothesis $E_{\text{exp}} \le 0.5$, the condition $(n-1)\varepsilon \le \frac{1}{4}$ suffices. In double format it means that $n$ should be no larger than $\frac{1}{4\varepsilon} = 2^{51}$, while in single format, $n$ should be no larger than $\frac{1}{4\varepsilon} = 2^{20}$.

*Proof* We apply Theorem 4 with $\widehat{\exp}(a)$ for $\hat{a}$, $\exp(a)$ for $a$, $E_{\text{exp}}$ for $A$ and $0$ for $B$. The hypothesis of Theorem 4 is satisfied due to the hypothesis $|a_i| \le M_{\text{exp}}$ that is needed to call $\widehat{\exp}$ on each $a_i$. Therefore

$$\left| \bigoplus_{0 \le i < n} \widehat{\exp}(a_i) - \sum_{0 \le i < n} \exp(a_i) \right| \le (E_{\text{exp}} + (n-1)\varepsilon(1 + E_{\text{exp}})) \sum_{0 \le i < n} \exp(a_i)$$

We then use Lemma 5 with $\bigoplus_{0 \le i < n} \widehat{\exp}(a_i)$ for $\hat{x}$ and $\sum_{0 \le i < n} \exp(a_i)$ for $x$ and $A = E_{\text{exp}} + (n-1)\varepsilon(1+E_{\text{exp}})$. Since $E_{\text{exp}}$ is assumed not larger than $0.5$, hence $A = E_{\text{exp}}+(n-1)\varepsilon(1+E_{\text{exp}}) \le 0.5 + \frac{1}{4} \times 1.5 < 1$. By applying Lemma 5, we get

$$\left| \log \left( \bigoplus_{0 \le i < n} \widehat{\exp}(a_i) \right) - \log \left( \sum_{0 \le i < n} \exp(a_i) \right) \right| \le - \log(1 - (E_{\text{exp}} + (n-1)\varepsilon(1 + E_{\text{exp}})))$$

that is

$$\left| \log \left( \bigoplus_{0 \le i < n} \widehat{\exp}(a_i) \right) - \text{LSE}(a) \right| \le - \log(1 - (E_{\text{exp}} + (n-1)\varepsilon(1 + E_{\text{exp}})))$$

Therefore

$$\begin{aligned} & \left| \widehat{\text{LSE}}(a) - \text{LSE}(a) \right| \\ & \le \left| \widehat{\text{LSE}}(a) - \log \left( \bigoplus_{0 \le i < n} \widehat{\exp}(a_i) \right) \right| + \left| \log \left( \bigoplus_{0 \le i < n} \widehat{\exp}(a_i) \right) - \text{LSE}(a) \right| \\ & = \left| \log \left( \bigoplus_{0 \le i < n} \widehat{\exp}(a_i) \right) - \log \left( \bigoplus_{0 \le i < n} \widehat{\exp}(a_i) \right) \right| + \left| \log \left( \bigoplus_{0 \le i < n} \widehat{\exp}(a_i) \right) - \text{LSE}(a) \right| \end{aligned}$$

<footnote>$^2$because $\log(1+A) + \log(1-A) = \log((1+A)(1-A)) = \log(1-A^2) \le 0$.</footnote>
At this point we want to apply the assumption of the accuracy of $\widehat{\log}$. For that, we need to show that the argument $\bigoplus_{0 \le i < n} \widehat{\exp}(a_i)$ is in the assumed domain, that is, is positive and no larger than $M_{\log}$. It is positive as a sum of positive numbers, and according to Corollary 2 it is bounded by $n(1 + \varepsilon(n-1))M_a$ where $M_a$ is a bound for all $\widehat{\exp}(a_i)$. Such a bound can be obtained from the assumption on the accuracy of $\widehat{\exp}$ giving

$$ \widehat{\exp}(a_i) \le \exp(a_i)(1 + E_{\exp}) \le \exp(M_{\exp})(1 + E_{\exp}) $$

Our sum is thus bounded by $n(1 + \varepsilon(n-1))\exp(M_{\exp})(1 + E_{\exp})$ which is smaller or equal to $M_{\log}$ using the hypothesis (6). We can continue our proof with

$$ \begin{aligned} |\widehat{\text{LSE}}(a) - \text{LSE}(a)| & \le E_{\log} \times \left| \log \left( \bigoplus_{0 \le i < n} \widehat{\exp}(a_i) \right) \right| + \left| \log \left( \bigoplus_{0 \le i < n} \widehat{\exp}(a_i) \right) - \text{LSE}(a) \right| \\ & \le (E_{\log} + 1) \times \left| \log \left( \bigoplus_{0 \le i < n} \widehat{\exp}(a_i) \right) - \text{LSE}(a) \right| + E_{\log} \times |\text{LSE}(a)| \\ & \le E_{\log} |\text{LSE}(a)| - \log \left(1 - (E_{\exp} + (n-1)\varepsilon(1 + E_{\exp}))\right) (1 + E_{\log}) \end{aligned} $$

which is the expected bound. □

## 3.3 Discussion on the Variations of the Bound on Accuracy

The error bound of $\widehat{\text{LSE}}$ has two parts: a relative part $E_{\log}$, and a constant part $-\log(1 - (E_{\exp} + (n-1)\varepsilon(1 + E_{\exp}))(1 + E_{\log}))$. Note that $-\log(1-x) < 2x$ for $x \le \frac{1}{2}$. We can therefore bound the constant error by $2 \times (E_{\exp} + (n-1)\varepsilon(1 + E_{\exp}))(1 + E_{\log})$. The factors that dominate the constant bound are $E_{\exp}$ and $\varepsilon \times (n-1)$. In Section 5.1, we compare this bound with the one proposed by Blanchard et al. [10]. The error bound grows linearly with $E_{\log}$, $E_{\exp}$ and $n$. Since it is possible to choose an implementation of $\widehat{\log}$ and $\widehat{\exp}$ with specific bounds, having the error bound of $\widehat{\text{LSE}}$ depending on $E_{\exp}$ and $E_{\log}$ is useful in order to control the error. To give some instances of the obtained bound, we can choose specific values for the parameters $E_{\log}$, $E_{\exp}$, $M_{\exp}$, $M_{\log}$, $M_a$ and $n$. Let us assume reasonable bounds in practice on $n$ and $M_a$ that is $2^{10} = 1024$ and $M_a = 25$.

Let us assume first we have some correctly rounded implementations of exponential and logarithm, that is $E_{\exp} = 2^{-53}$ and $E_{\log} = 2^{-53}$. Then, to ensure that hypothesis (6) of Theorem 6 holds, it suffices to have $M_{\log}$ larger than $n(1 + \varepsilon(n-1))\exp(M_a)(1 + E_{\exp}) \le 2^{10}(1 + 2^{-53} \times 1023)\exp(25)(1 + 2^{-53}) \le 7.38 \times 10^{13}$. Assuming thus that the implementation of $\widehat{\log}$ is correctly rounded on the domain given by the bound $M_{\log}$ above, the relative error on $\widehat{\text{LSE}}$ is $E_{\log} = 2^{-53}$ and the absolute error is bounded by $2 \times (E_{\exp} + (n-1)\varepsilon(1 + E_{\exp}))(1 + E_{\log}) \le 2 \times (2^{-53} + 1023 \times 2^{-53}(1 + 2^{53}))(1 + 2^{-53}) \le 2.28 \times 10^{-13}$.

Let us assume less precise implementations of exponential and logarithm with $E_{\exp} = 2^{-40}$ and $E_{\log} = 2^{-36}$. These are some bounds for efficient implementations of exponential and logarithm that empirically seemed sufficiently accurate and energy-saving for industrial applications like IoT systems or deep-learning frameworks [40]. Then, to ensure that hypothesis (6) of Theorem 6 holds, it suffices to have $M_{\log}$ larger than $n(1 + \varepsilon(n-1))\exp(M_a)(1 + E_{\exp}) \le 2^{10}(1 + 2^{-53} \times 1023)\exp(25)(1 + 2^{-40}) \le 7.38 \times 10^{13}$ that is roughly the same bound as above with correct rounding on $\widehat{\exp}$ and $\widehat{\log}$. In other words, the required bound on the input domain of logarithm depends mostly on the bound on inputs and the number of elements in the input sequence. The relative
constant exp_max_value : real (* constant M<sub>exp</sub> *)
axiom exp_max_value_spec: 0.0 < exp_max_value ≤ 708.0

constant exp_error: real (* constant E<sub>exp</sub> *)
axiom exp_error_bound : 0.0 < exp_error ≤ 0.5

function u_exp (x:udouble) : udouble (* function $\widehat{exp}$ *)
axiom u_exp_spec :forall x:udouble. abs (to_real x) ≤ exp_max_value →
abs (to_real (u_exp x) − exp (to_real x)) ≤ exp (to_real x) × exp_error

[...]

function u_lse (a: int → udouble) (n:int) : udouble
= u_log (u_sum (fun i → u_exp (a i)) n)

<caption>Fig. 6 Declaration of $\widehat{exp}$ in WhyML, with assumed accuracy.</caption>

error on the computation of LSE is now $E_{\log} = 2^{-36}$ and the absolute error is bounded by $2 \times (E_{\exp} + (n-1)\varepsilon(1+E_{\exp})) (1+E_{\log}) \le 2 \times (2^{-40} + 1023 \times 2^{-53}(1+2^{53})) (1+2^{-36}) \le 2.05 \times 10^{-12}$. This absolute error is roughly twice the absolute error of the case with correct rounded implementations of exp and log.

Notice that if the size of the sequence is significantly larger than the assumed bound 1024, then the required value for $M_{\log}$ gets significantly larger, indeed it increases roughly linearly with this size.

## 3.4 Accuracy of LSE proved in WhyML

To start with, we need to declare the approximations of exp and log that we consider. The declaration of $\widehat{exp}$ is shown in Figure 6 and the one of $\widehat{log}$ is similar. These declarations are axiomatic so as to make them parametric in the values of $M_{\exp}$, $E_{\exp}$ and such. The definition of $u_lse$ for $\widehat{LSE}$ follows naturally.

The WhyML statement corresponding to Theorem 6 is given in Figure 7. Preconditions on lines 2–3 express assumptions on the size of the array. The precondition on line 4 expresses the bound on the array elements, and the precondition on lines 7–8 expresses Hypothesis (6). The post-condition on lines 10–12 expresses the bound on accuracy given by Theorem 6. The body of the WhyML function `lse State` constitutes a formal proof of Theorem 6. It is 38 lines long and available at on the dedicated page of Toccata gallery [51]. For this function, Why3 generates 27 VCs, 13 of them are proved by Alt-Ergo, 1 by CVC4, 11 by cvc5, and 2 by Z3. Of course, that proof requires Lemma 5 which needs to be proved first. Its proof generates 8 VCs, 4 of them are proved by Alt-Ergo, 2 by CVC4, 1 by cvc5 and 1 by Z3.

## 3.5 Proving a C code implementing LSE

We aim to achieve proofs on concrete C code. For that, we use the environment TIS-kernel, a fork of Frama-C, and its J³ plug-in for deductive verification, which is a
<figure><img src="image_4.png" /><figcaption>Fig. 7 Statement of Theorem 6 in WhyML.</figcaption></figure>

prototype under development. Alternatively, there should be no technical difficulty to
achieve the proofs of our C code using the regular Frama-C environment and its Wp
plug-in for deductive verification.

Our C code computing the LSE function is given on Figure 8. In a first step, let's ignore the potential floating-point overflow, and focus on proving the accuracy property. To specify the intended behavior and its properties, we build a bridge to WhyML definitions (see our report [22] for technical details), so that for example we can use the WhyML definitions of `u_exp` and `u_log` in the ACSL annotations. It provides in particular a function `to_udouble` that promotes a regular C double to an unbounded double. We declare and specify the auxiliary C functions for computing approximations of exp and log, as shown on Figure 9.

The first post-condition, on lines 10–11 of Figure 8, thus expresses that the result of the C function is equal to the LSE function defined in WhyML. The second post-condition, on lines 12–16 of Figure 8, expresses the expected bounding property, as stated by Theorem 6. The second post-condition is going to be proved easily from the first one, and the accuracy result on `u_lse` already proved in WhyML. The precondition on line 3 is required to ensure that all calls to the exponential are within its correct domain. The precondition on lines 4–6 expresses the required hypothesis 6 of Theorem 6.

The first post-condition on lines 10–11 is an easy consequence of the definition of the LSE function, and the loop invariant given on lines 21–22. The post-condition on lines 12–16 is proved by invoking the proof of the same statement already done in WhyML. Together with the simple loop invariants on lines 20 and 24, all the VCs are proved, in particular the expected post-conditions, except two of them. The first unproved VC is related to line 47 where it is required to show that `s` fits in the expected domain of the approximated logarithm. The second unproved VC is related to line 42 where it is required to show the absence of floating-point overflow when performing addition. To prove the VC on line 47 and thus prove that the sum `s` on line 47 fits in the expected range of the log, we need to state the additional loop invariants on lines 24 and 25–26 to bound the sum. To prove that these invariants hold, we need
1 /*@ requires 0 < size <= 1024 && max_a <= exp_max_value;
2 @ requires \initialized (&a[0..size-1]);
3 @ requires \forallall integer i; 0 <= i < size ==> \abs(a[i]) <= max_a;
4 @ // the hypothesis (6) of Theorem 6
5 @ requires \exp(max_a) * (1.0 + exp_error) * size *
6 @ (1.0 + (eps * (size-1))) <= log_max_value;
7 @ // additional requirements to prevent overflow on addition
8 @ requires max_a <= 701.0;
9 @ requires log_max_value <= 0x1p1023;
10 @ // result is equal to the WhyML def of LSE on udouble
11 @ ensures to_udouble(\result) == u_lse(a, size);
12 @ // the accuracy property
13 @ ensures
14 @ \abs(\result - lse_exact(a, size)) <= log_error * \abs(lse_exact(a,size))
15 @ - \log(1 - (exp_error + eps * (size-1) * (1+exp_error)))
16 @ * (1+log_error);
17 @*/
18 double log_sum_exp(size_t size) {
19 int i; double s = 0.0;
20 /*@ loop invariant 0 <= i <= size;
21 @ loop invariant // to prove the first post-condition
22 @ to_udouble(s) == u_sum_of_u_exp(a, 0, i);
23 @ // for proving s is the domain of the log
24 @ loop invariant (i == 0 ? s == 0.0 : 0.0 < s);
25 @ loop invariant \forallall integer j; 0 <= j < i ==>
26 @ \abs(to_real(u_exp(to_udouble(a[j])))) <= \exp(max_a) * (1.0+exp_error) ;
27 @ loop assigns i, s;
28 @ loop variant size-i;
29 @*/
30 for (i = 0; i < size; i++) {
31 /*@ assert 0.0 <= to_real(u_exp(to_udouble(a[i])));
32 @ assert to_real(to_udouble(a[i])) <= max_a ;
33 @ assert \exp(to_real(to_udouble(a[i])))) <= \exp(max_a) ;
34 @ assert
35 @ \abs(to_real(u_exp(to_udouble(a[i])))) - \exp(to_real(to_udouble(a[i]))))
36 @ <= \exp(max_a) * exp_error ;
37 @ assert to_real(u_exp(to_udouble(a[i])))) <= \exp(max_a) * (1.0+exp_error) ;
38 @ assert // invocation of Corollary (2)
39 @ usum_double_bound(u_sum_of_u_exp(a, 0, i),
40 @ \exp(max_a)*(1.0+exp_error), size);
41 @*/
42 s += expApprox(a[i]);
43 }
44 /*@ assert // another invocation of Corollary (2)
45 @ usum_double_bound(to_udouble(s), \exp(max_a)*(1.0+exp_error), size);
46 @*/
47 return logApprox(s);
48 }

Fig. 8 C code for computing LSE, annotated with ACSL specifications.

<caption>Fig. 8 C code for computing LSE, annotated with ACSL specifications.</caption>
/*@ requires \abs(x) <= exp_max_value;
@ ensures to_udouble(\result) == u_exp(to_udouble(x));
@ assigns \nothing;
@*/
extern double expApprox(double x);

/*@ requires 0 < x <= log_max_value;
@ ensures to_udouble(\result) == u_log(to_udouble(x));
@ assigns \nothing;
@*/
extern double logApprox(double x);

<caption>Fig. 9 External C functions for $\widehat{exp}$ and $\widehat{log}$, specified in ACSL.</caption>

to invoke Corollary 2 again. Achieving this proof is a bit involved, requiring all the extra intermediate assertions on lines 31–39. The last VC remaining to prove is the absence of numerical overflow when computing the addition on line 42. It is indeed expected since the C code operates on true IEEE floating-point numbers and not the unbounded ones. To achieve this, the given pre-conditions are not enough, we need to assume extra bounds on the inputs. So far we assumed the input numbers smaller than $M_{exp}$, for which no upper bound is assumed so far. Yet, we add the exponentials of these numbers, and summing up to say 1024 of these numbers, we can indeed have an overflow. A tighter bound must be assumed. We assume here, on lines 8–9 of Figure 8, that $M_a$ is smaller than 701 and $M_{log}$ is smaller than $2^{1023}$. With these extra assumptions, and thanks to the already stated and proved loop invariants on lines 24 and 25–26, the VC is proved. In all, 56 VCs are generated. 50 of them are proved by Alt-Ergo, within a 5 seconds time limit, and CVC4 and cvc5 are able to solve the 6 remaining VCs.

# 4 Accuracy of Sum-Log-Sum-Exp

The goal of this section is to analyse the accuracy of a floating-point implementation of the function SLSE that we introduced in Section 1. To obtain a bound on its accuracy, we proceed in a more or less systematic way by stating lemmas of propagation of errors, in the same vein as the Lemma 5 we stated for the natural logarithm.

## 4.1 Error Propagation Lemmas

**Lemma 7 (Propagation of errors through addition)** Assuming $|\hat{x} - x| \le A\tilde{x}$ with $|x| \le \tilde{x}$ and $|\hat{y} - y| \le B\tilde{y}$ with $|y| \le \tilde{y}$, for any non-negative parameters A, B, we have

$$|(\hat{x} \oplus \hat{y}) - (x + y)| \le (A + B + \varepsilon)(\tilde{x} + \tilde{y})$$

*Proof* We first prove the following inequality:

$$|(\hat{x} \oplus \hat{y}) - (x + y)| \le (\varepsilon + B)\tilde{x} + (\varepsilon + A)\tilde{y} \quad (7)$$

We distinguish three cases:
• if $\tilde{x} \le \varepsilon\tilde{y}$ then by Property (3) we get $|(\hat{x} \oplus \hat{y}) - (\hat{x} + \hat{y})| \le |\hat{x}| \le A\tilde{x} + |x| \le A\tilde{x} + \tilde{x} \le A\varepsilon\tilde{y} + \varepsilon\tilde{y} \le (\varepsilon + B)\tilde{x} + (\varepsilon + A)\tilde{y}$ (because $\varepsilon \le 1$) which proves (7).

• if $\tilde{y} \le \varepsilon\tilde{x}$, by using Property (4) and the same reasoning as above.

• otherwise, we have $\varepsilon\tilde{x} \le \tilde{y}$ and $\varepsilon\tilde{y} \le \tilde{x}$. We have, by Property (2), $|(\hat{x} \oplus \hat{y}) - (\hat{x} + \hat{y})| \le \varepsilon|\hat{x} + \hat{y}| \le \varepsilon|(\hat{x} - x + x) + (\hat{y} - y + y)| \le \varepsilon(|\hat{x} - x| + |x| + |\hat{y} - y| + |y|) \le \varepsilon(A\tilde{x} + \tilde{x} + B\tilde{y} + \tilde{y}) \le A\tilde{y} + B\tilde{x} + \varepsilon(\tilde{x} + \tilde{y}) \le (\varepsilon + B)\tilde{x} + (\varepsilon + A)\tilde{y}$

We now use Property (7) to prove our lemma as follows. For the upper bound we have $\hat{x} \oplus \hat{y} \le (\hat{x} + \hat{y}) + (\varepsilon + B)\tilde{x} + (\varepsilon + A)\tilde{y} \le x + A\tilde{x} + y + B\tilde{y} + (\varepsilon + B)\tilde{x} + (\varepsilon + A)\tilde{y} = (x+y)+(A+B+\varepsilon)(\tilde{x}+\tilde{y})$ and similarly for the lower bound $\hat{x} \oplus \hat{y} \ge (\hat{x}+\hat{y})-(\varepsilon+B)\tilde{x}-(\varepsilon+A)\tilde{y} \ge x - A\tilde{x} + y - B\tilde{y} - (\varepsilon + B)\tilde{x} - (\varepsilon + A)\tilde{y} = (x + y) - (A + B + \varepsilon)(\tilde{x} + \tilde{y})$. $\square$

We devise a similar propagation lemma for multiplication, but first need to state and prove an auxiliary lemma.

**Lemma 8** *Assuming $|\hat{x} - x| \le A\tilde{x}$ with $|x| \le \tilde{x}$ and $0 \le A$, for any real $z$ we have $|\hat{x}z - xz| \le A\tilde{x}|z|$.*

*Proof* We have $|\hat{x}z - xz| = |(\hat{x} - x)z| = |(\hat{x} - x)| \times |z| \le A\tilde{x}|z|$. $\square$

**Lemma 9** (Propagation lemma through multiplication) *Assuming $|\hat{x} - x| \le A\tilde{x}$ with $|x| \le \tilde{x}$ and $|\hat{y} - y| \le B|x|$ with $|y| \le \tilde{y}$, for any non-negative parameters $A, B$, we have*

$$|(\hat{x} \otimes \hat{y}) - xy| \le (\varepsilon + (A + B + AB)(1 + \varepsilon))\tilde{x}\tilde{y} + \eta$$

*Proof* We have

$$\hat{x}\hat{y} \le x\hat{y} + A\tilde{x}|\hat{y}| \quad \text{(by Lemma 8 with } z = \hat{y}\text{)}$$

$$\le xy + B|x|\tilde{y} + A\tilde{x}|\hat{y}| \quad \text{(by Lemma 8 with } z = x\text{)}$$

$$\le xy + B\tilde{x}\tilde{y} + A\tilde{x}|\hat{y}|$$

We can use Lemma 8 one last time with $z = \tilde{x}$ to obtain $\tilde{x}\hat{y} \le \tilde{x}y + B\tilde{x}\tilde{y}$ hence $\tilde{x}|\hat{y}| \le \tilde{x}|y| + B\tilde{x}\tilde{y} \le (1 + B)\tilde{x}\tilde{y}$. Combining it with inequality above, we get

$$\hat{x}\hat{y} \le xy + (A + B + AB)\tilde{x}\tilde{y} \quad (8)$$

Finally, we have

$$\hat{x} \otimes \hat{y} - xy = (\hat{x} \otimes \hat{y} - \hat{x}\hat{y}) + (\hat{x}\hat{y} - xy) \quad \text{(by Property (1))}$$

$$\le (\varepsilon|\hat{x}\hat{y}| + \eta) + (\hat{x}\hat{y} - xy) \quad \text{(by (8))}$$

and since $|\hat{x}| \le (1 + A)\tilde{x}$ and $|\hat{x}| \le (1 + B)\tilde{y}$ we get

$$\hat{x} \otimes \hat{y} - xy \le \varepsilon(1 + A)(1 + B)\tilde{x}\tilde{y} + \eta + (A + B + AB)\tilde{x}\tilde{y} \\ = (\varepsilon + (1 + \varepsilon)(A + B + AB))\tilde{x}\tilde{y} + \eta$$

We proceed similarly to prove the lower bound for $\hat{x} \otimes \hat{y} - xy$. $\square$
We have now error propagation lemmas for addition and multiplication, we continue by stating propagation lemmas for the exponential and the logarithm in base 2.

**Lemma 10 (Propagation lemma through exponential)** Assuming $|\hat{x} - x| \le A\tilde{x} + C$ with $|x| \le \tilde{x}$ and $|\hat{x}| \le M_{\text{exp}}$, for any non-negative parameters A, C, we have

$$|\widehat{\exp}(\hat{x}) - \exp(x)| \le \exp(x)(E_{\text{exp}} + (\exp(A\tilde{x} + C) - 1)(1 + E_{\text{exp}}))$$

*Proof* We first show the formula

$$|\exp(\hat{x}) - \exp(x)| \le \exp(x)(\exp(A\tilde{x} + C) - 1) \quad (9)$$

by distinguishing two cases for the upper bound and the lower bound. For the upper bound, by assumption we have

$$\hat{x} \le x + A\tilde{x} + C$$

hence, since exp is monotonic,

$$\exp(\hat{x}) \le \exp(x + A\tilde{x} + C) = \exp(x) \exp(A\tilde{x} + C)$$

which entails Formula 9. For the lower bound, by assumption we have $\hat{x} \ge x - A\tilde{x} - C$ therefore $\exp(\hat{x}) \ge \exp(x - A\tilde{x} - C) = \exp(x) \exp(-A\tilde{x} - C)$ thus

$$\exp(\hat{x}) - \exp(x) \ge \exp(x)(\exp(-A\tilde{x} - C) - 1) \quad (10)$$

Furthermore, it is a known mathematical fact<fnref n="3" /> that for any y, $\exp(y) + \exp(-y) \ge 2$, therefore $\exp(A\tilde{x} + C) + \exp(-A\tilde{x} - C) \ge 2$ hence $\exp(-A\tilde{x} - C) - 1 \ge -\exp(A\tilde{x} + C) + 1$ so $\exp(x)(\exp(-A\tilde{x} - C) - 1) \ge \exp(x)(-\exp(A\tilde{x} + C) + 1) = -\exp(x)(\exp(A\tilde{x} + C) - 1)$. We have thus proved Formula 9.

To obtain the desired formula we again reason with two cases for the upper bound and the lower bound. For the upper bound we have $\widehat{\exp}(\hat{x}) \le \exp(\hat{x})(1 + E_{\text{exp}})$ so from Formula 9

$$\widehat{\exp}(\hat{x}) - \exp(x) \le \exp(x)(\exp(A\tilde{x} + C)(1 + E_{\text{exp}}) - 1)$$

and the factor of $\exp(x)$ is

$$\exp(A\tilde{x} + C)(1 + E_{\text{exp}}) - 1 = (\exp(A\tilde{x} + C) - 1)(1 + E_{\text{exp}}) + (1 + E_{\text{exp}}) - 1 \\ = E_{\text{exp}} + (\exp(A\tilde{x} + C) - 1)(1 + E_{\text{exp}})$$

Similarly, for the lower bound, we have $\widehat{\exp}(\hat{x}) \ge \exp(\hat{x})(1 - E_{\text{exp}})$ so from Formula 9

$$\widehat{\exp}(\hat{x}) - \exp(x) \ge (\exp(x) - \exp(x)(\exp(A\tilde{x} + C) - 1))(1 - E_{\text{exp}}) - \exp(x) \\ = -\exp(x)(1 + ((\exp(A\tilde{x} + C) - 1) - 1)(1 - E_{\text{exp}}))$$

and the factor of $-\exp(x)$ is

$$1 + ((\exp(A\tilde{x} + C) - 1) - 1)(1 - E_{\text{exp}}) = 1 + (\exp(A\tilde{x} + C) - 1)(1 - E_{\text{exp}}) - (1 - E_{\text{exp}}) \\ = E_{\text{exp}} + (\exp(A\tilde{x} + C) - 1)(1 - E_{\text{exp}}) \\ \le E_{\text{exp}} + (\exp(A\tilde{x} + C) - 1)(1 + E_{\text{exp}})$$

□

<footnote n="3">³a short proof is $\exp(y) - 2 + \exp(-y) = \exp(y)(1 - 2\exp(-y) + (\exp(-y))^2) = \exp(y)(1 - \exp(-y))^2 \ge 0$.</footnote>
As we did in Section 3.1 for functions exp and log, we assume given an implementation of the base 2 logarithm, with the following property: for all $x$ such that $0 < x \le M_{\log_2}$, $|\log_2(x) - \log_2(x)| \le E_{\log_2}|\log_2(x)|$ where the parameter $M_{\log_2}$ is positive and the parameter $E_{\log_2}$ is positive and smaller than or equal to 1.

**Lemma 11 (Propagation lemma for $\log_2$)** For any positive real numbers $x$ and $\hat{x}$ such that $|\hat{x} - x| \le Ax$ with $0 \le A < 1$, and $|\hat{x}| \le M_{\log_2}$, we have:

$$|\widehat{\log_2(\hat{x})} - \log_2(x)| \le E_{\log_2}|\log_2(x)| - \log_2(1 - A)(1 + E_{\log_2})$$

*Proof* We have $\log_2(x) = \frac{\log(x)}{\log(2)}$, therefore by applying Lemma 5 we get $|\log_2(\hat{x}) - \log_2(x)| \le -\log_2(1 - A)$. We have $|\widehat{\log_2(\hat{x})} - \log_2(\hat{x})| \le E_{\log_2}|\log_2(\hat{x})|$, therefore, we can apply the triangular inequality to get

$$\begin{aligned}
|\widehat{\log_2(\hat{x})} - \log_2(x)| &\le |\widehat{\log_2(\hat{x})} - \log_2(\hat{x})| + |\log_2(\hat{x}) - \log_2(x)| \\
&\le E_{\log_2}|\log_2(\hat{x})| + |\log_2(\hat{x}) - \log_2(x)| \\
&= E_{\log_2}|\log_2(\hat{x}) - \log_2(x) + \log_2(x)| + |\log_2(\hat{x}) - \log_2(x)| \\
&\le E_{\log_2}|\log_2(\hat{x})| + (1 + E_{\log_2})|\log_2(\hat{x}) - \log_2(x)| \\
&\le E_{\log_2}|\log_2(\hat{x})| - (1 + E_{\log_2})\log_2(1 - A)
\end{aligned}$$

□

## 4.2 Analysis of Accuracy of SLSE function

For readability, let us denote $g(x, y, \rho) = \exp\left(-\frac{(x+\rho-y)^2}{2}\right)$ so that the definition of the SLSE function can be rewritten as

$$\text{SLSE}(a, \rho) = \sum_{0 \le i < n} \log_2 \left( \sum_{0 \le j < n} g(a_i, a_j, \rho) \right)$$

Similarly, the floating-point approximation of $g$ is $\hat{g}(x, y, \rho) = \widehat{\exp}(\ominus((((x \oplus \rho) \ominus y) \otimes ((x \oplus \rho) \ominus y)) \otimes 2))$ and the floating-point approximation of the SLSE function is $\widehat{\text{SLSE}}(a, \rho) = \bigoplus_{0 \le i < n} \widehat{\log_2}\left(\bigoplus_{0 \le j < n} \hat{g}(a_i, a_j, \rho)\right)$.

To achieve our main result on the accuracy of SLSE function, we need a few more restrictions on the input parameters. We assume that the components of the vector $a$ are upper bounded with some constant $M_a$. We pose the following constraint which is issued from the need of calling the exponential function on arguments which must be smaller than $M_{\exp}$:

$$(2M_a + |\rho|)^2(\frac{1}{2} + A_\epsilon) + \frac{3}{2}\eta \le M_{\exp} \quad (11)$$

where $A_\epsilon = \frac{\epsilon}{2}(1 + 4(1 + \epsilon)^2)$.
We are now ready to state and prove, as a first step, a result on the accuracy of the g auxiliary function. Notice that trivially $0 < g(x, y, \rho) \le 1$ for any $x, y, \rho$.

**Lemma 12 (Accuracy of the g function)** *Assuming x and y satisfying $|x| \le M_a$ and $|y| \le M_a$, we have*

$$|\hat{g}(x, y, \rho) - g(x, y, \rho)| \le g(x, y, \rho) \left( E_{\text{exp}} + \left( \exp \left( A_{\varepsilon} \left( |x| + |\rho| + |y| \right)^2 + \frac{3}{2} \eta \right) - 1 \right) (1 + E_{\text{exp}}) \right)$$

*Proof* Let $f(x, y, \rho) = x + \rho - y$ and $\hat{f}(x, y, \rho) = (x \oplus \rho) \ominus y$. From Property (2) we have $|(x \oplus \rho) - (x + \rho)| \le \varepsilon|x + \rho|$ hence from Lemma 7, with $\tilde{x} = |x + \rho|$, $\tilde{y} = |y|$, $A = \varepsilon$, $B = C = D = 0$, we have $|((x \oplus \rho) \ominus y) - ((x + \rho) - y)| \le 2\varepsilon(|x + \rho| + |y|)$ then, since $|(x + \rho) - y| \le |x| + |\rho| + |y|$, using Lemma 9, with $\tilde{x} = \tilde{y} = |x| + |\rho| + |y|$, $A = B = 2\varepsilon$ we get

$$|(\hat{f}(x, y, \rho) \otimes \hat{f}(x, y, \rho)) - f(x, y, \rho)^2| \le (\varepsilon + (2\varepsilon + 2\varepsilon + 4\varepsilon^2)(1 + \varepsilon))(|x| + |\rho| + |y|)^2 + \eta \\ \le 2A_{\varepsilon}(|x| + |\rho| + |y|)^2 + \eta$$

Dividing by 2 is an exact operation except for de-normalised numbers for which a constant error of $\eta$ can occur. Taking the opposite is an exact operation, therefore

$$| \ominus ((\hat{f}(x, y, \rho) \otimes \hat{f}(x, y, \rho)) \ominus 2) - (-f(x, y, \rho)^2/2) | \le A_{\varepsilon} (|x| + |\rho| + |y|)^2 + \frac{\eta}{2} + \eta$$

We need now to call the exponential function on the floating-point argument so we need to show that its absolute value is bounded by $M_{\text{exp}}$. From the previous formula we have

$$\begin{aligned} | \ominus ((\hat{f}(x, y, \rho) \otimes \hat{f}(x, y, \rho)) \ominus 2) | &\le |f(x, y, \rho)^2/2| + A_{\varepsilon} (|x| + |\rho| + |y|)^2 + \frac{3}{2}\eta \\ &\le (|x| + |\rho| + |y|)^2 / 2 + A_{\varepsilon} (|x| + |\rho| + |y|)^2 + \frac{3}{2}\eta \\ &= (|x| + |\rho| + |y|)^2 \left( \frac{1}{2} + A_{\varepsilon} \right) + \frac{3}{2}\eta \\ &\le (2M_a + |\rho|)^2 \left( \frac{1}{2} + A_{\varepsilon} \right) + \frac{3}{2}\eta \\ &\le M_{\text{exp}} \end{aligned} \quad (\text{by Hypothesis (11)})$$

thus justifying the choice of Hypothesis (11). We can now apply Lemma 10 for error propagation by the exponential with $\tilde{x} = (|x| + |\rho| + |y|)^2$, $A = A_{\varepsilon}$ and $C = \frac{3}{2}\eta$ to get the expected result. □

To ease the reading of the remaining, let us introduce the constant $G = E_{\text{exp}} + (\exp(A_{\varepsilon}(2M_a + |\rho|)^2 + \frac{3}{2}\eta) - 1)(1 + E_{\text{exp}})$. We now state a simple corollary of Lemma 12 and two auxiliary lemmas about G.

**Corollary 13** *Assuming x and y satisfying $|x| \le M_a$ and $|y| \le M_a$, we have $|\hat{g}(x, y, \rho) - g(x, y, \rho)| \le G \times g(x, y, \rho)$.

**Lemma 14 (bound on G)** *Assume that, for the chosen floating-point format and implementation of exponential, $2A_{\varepsilon}M_{\text{exp}} + \frac{3}{2}\eta \le \frac{1}{64}$, and assume also that Hypothesis (11) holds, then $G \le \frac{9}{16}$.*
**Proof** From Hypothesis (11) we have $(2M_a + |\rho|)^2 \le \frac{M_{\text{exp}} - \frac{3}{2}\eta}{\frac{1}{2} + A_{\varepsilon}} \le \frac{M_{\text{exp}}}{\frac{1}{2}} = 2M_{\text{exp}}$ hence $A_{\varepsilon}(2M_a + |\rho|)^2 + \frac{3}{2}\eta \le 2A_{\varepsilon}M_{\text{exp}} + \frac{3}{2}\eta$ hence $G \le E_{\text{exp}} + (\exp(2A_{\varepsilon}M_{\text{exp}} + \frac{3}{2}\eta) - 1)(1 + E_{\text{exp}})$. Moreover, we assumed $E_{\text{exp}} \le \frac{1}{2}$ so $G \le \frac{1}{2} + 2(\exp(2A_{\varepsilon}M_{\text{exp}} + \frac{3}{2}\eta) - 1)$. Since we assumed $2A_{\varepsilon}M_{\text{exp}} + \frac{3}{2}\eta \le \frac{1}{64}$ we have$^4$ $\exp(2A_{\varepsilon}M_{\text{exp}} + \frac{3}{2}\eta) - 1 \le 4A_{\varepsilon}M_{\text{exp}} + 3\eta \le \frac{1}{32}$ hence $G \le \frac{1}{2} + 2\frac{1}{32} = \frac{9}{16}$. $\square$

Note that the extra hypothesis above is satisfied for the double format where $A_{\varepsilon} = \frac{\varepsilon}{2}(1 + 4(1 + \varepsilon)^2) \le 2^{-54} \times 6$ and we assumed $M_{\text{exp}}$ no larger than 708 so $2A_{\varepsilon}M_{\text{exp}} + \frac{3}{2}\eta \le 2^{-53} \times 4248 + 2^{-1073} \le 2^{-53} \times 2^{12} = 2^{-41}$. It is also satisfied for the single format where $A_{\varepsilon} = \frac{\varepsilon}{2}(1 + 4(1 + \varepsilon)^2) \le 2^{-25} \times 6$ and, assuming $M_{\text{exp}}$ no larger than 88, $2A_{\varepsilon}M_{\text{exp}} + \frac{3}{2}\eta \le 2^{-24} \times 528 + 2^{-149} \le 2^{-24} \times 2^{10} = 2^{-14}$.

**Lemma 15** (constant bounds on function $g$) For any $x$ and $y$ such that $|x| \le M_a$ and $|y| \le M_a$, we have $0 < \hat{g}(x, y, \rho) \le 1 + G$.

**Proof** By Corollary 13 we have $|\hat{g}(x, y, \rho) - g(x, y, \rho)| \le Gg(x, y, \rho)$ hence $\hat{g}(x, y, \rho) \ge (1 - G)g(x, y, \rho) > 0$ and $\hat{g}(x, y, \rho) \le (1 + G)g(x, y, \rho) \le 1 + G$. $\square$

Our main theorem for the accuracy of the floating point implementation of SLSE is then as follows. To ease the reading of the formula, we introduce the term

$$|\text{SLSE}|(a, \rho) = \sum_{0 \le i < n} \left| \log_2 \left( \sum_{0 \le j < n} g(a_i, a_j, \rho) \right) \right|$$

which is actually a formula similar to SLSE$(a, \rho)$ but with taking the absolute values for the elements of the outer sum.

**Theorem 16** (Accuracy of SLSE) For any $n \ge 1$ such that $(n - 1)\varepsilon \le \frac{1}{4}$ and any vector $a$ such that for all $a_i$, $|a_i| \le M_a$, assuming Hypothesis (11), and finally assuming

$$n(1 + (G + (n - 1)\varepsilon(1 + G))) \le M_{\log_2} \quad (12)$$

we have

$$\begin{aligned} \left| \widehat{\text{SLSE}}(a, \rho) - \text{SLSE}(a, \rho) \right| &\le (E_{\log_2} + (n - 1)\varepsilon(1 + E_{\log_2})) |\text{SLSE}|(a, \rho) \\ &\quad - \log_2 (1 - (G + (n - 1)\varepsilon(1 + G))) (1 + E_{\log_2}) n (1 + (n - 1)\varepsilon) \end{aligned}$$

**Proof** For any $i$ and $j$, $0 \le i < n$ and $0 \le j < n$, we invoke Corollary 13 to get $|\hat{g}(a_i, a_j, \rho) - g(a_i, a_j, \rho)| \le g(a_i, a_j, \rho) \times G$ so that we can invoke Lemma 4 with $A = G$ and $B = 0$, to get

$$\left| \bigoplus_{0 \le j < n} \hat{g}(a_i, a_j, \rho) - \sum_{0 \le j < n} g(a_i, a_j, \rho) \right| \le (G + (n - 1)\varepsilon(1 + G)) \sum_{0 \le j < n} g(a_i, a_j, \rho) \quad (13)$$

We would now like to invoke Lemma 11 with $\hat{x} = \bigoplus_{0 \le j < n} \hat{g}(a_j, a_i, \rho)$, $x = \sum_{0 \le j < n} g(a_j, a_i, \rho)$ and $A = G + (n - 1)\varepsilon(1 + G)$. We need to prove that $\hat{x}$ is positive and smaller than $M_{\log_2}$,

<footnote>$^4$because for any real $x$, if $0 \le x \le 1$ then $\exp x \le 1 + 2x$.</footnote>
and to prove the hypothesis $A < 1$ which here is equivalent to $G + (n-1)\varepsilon(1+G) < 1$. From Lemma 15 we know that $\hat{x}$ is positive, and together with Property (13) we have $\hat{x} \le (1 + (G + (n-1)\varepsilon(1+G))) \sum_{0 \le j < n} g(a_i, a_j, \rho) \le n(1 + (G + (n-1)\varepsilon(1+G))) \le M_{\log_2}$ from Hypothesis (12). Also, from Lemma 14 and since we assumed $(n-1)\varepsilon \le \frac{1}{4}$, $G + (n-1)\varepsilon(1+G) \le \frac{9}{16} + \frac{1}{4} \times \frac{25}{16} = \frac{61}{64} < 1$. That's it, we can apply Lemma 11 and get

$$ \left| \widehat{\log_2} \left( \bigoplus_{0 \le j < n} \hat{g}(a_i, a_j, \rho) \right) - \log_2 \left( \sum_{0 \le j < n} g(a_i, a_j, \rho) \right) \right| \le E_{\log_2} \left| \log_2 \left( \sum_{0 \le j < n} g(a_i, a_j, \rho) \right) \right| - \log_2 (1 - (G + (n-1)\varepsilon(1+G))) (1 + E_{\log_2}) $$

Finally, we apply Lemma 4 with $\hat{a}_i = \widehat{\log_2} \left( \bigoplus_{0 \le j < n} \hat{g}(a_i, a_j, \rho) \right)$, $a_i = \log_2 \left( \sum_{0 \le i < n} g(a_i, a_j, \rho) \right)$, $A = E_{\log_2}$ and $B = -\log_2 (1 - (G + (n-1)\varepsilon(1+G))) (1 + E_{\log_2})$ which gives us

$$ \left| \bigoplus_{0 \le i < n} \widehat{\log_2} \left( \bigoplus_{0 \le j < n} \hat{g}(a_i, a_j, \rho) \right) - \sum_{0 \le i < n} \log_2 \left( \sum_{0 \le j < n} g(a_i, a_j, \rho) \right) \right| \le (E_{\log_2} + (n-1)\varepsilon(1+E_{\log_2})) \sum_{0 \le i < n} \left| \log_2 \left( \sum_{0 \le j < n} g(a_i, a_j, \rho) \right) \right| - \log_2 (1 - (G + (n-1)\varepsilon(1+G))) (1 + E_{\log_2}) n (1 + (n-1)\varepsilon) $$

□

## 4.3 Variations of the Bound on Accuracy of SLSE

As we did in Section 3.3 for the LSE function, let us analyse how the bound given by Theorem 16 varies depending of the various parameters involved. That error bound comes in two parts. First, there is a (more or less) relative part $E_{\log_2} + (n-1)\varepsilon(1+E_{\log_2})$. This is not properly speaking a relative error because it is a coefficient of $|\text{SLSE}|(a, \rho)$ and not $|\text{SLSE}(a, \rho)|$. Indeed, since the argument passed to $\log_2$ can be smaller or larger than 1, there is no way to avoid taking the absolute value of the $\log_2$. Unlike for the LSE function, this relative error depends on the size $n$ of the input vector, whereas for LSE it was only depending on the relative error $E_{\log}$ of the log. That relative error coefficient is approximately $E_{\log_2} + n\varepsilon$, and grows linearly in $E_{\log_2}$ and $n$, but none of the two terms necessarily dominates the other. Second, the error bound comes with an absolute part, which is $-\log_2(1 - (G + (n-1)\varepsilon(1+G))) (1 + E_{\log_2}) n (1 + (n-1)\varepsilon)$ and can be approximated to $2(G + (n-1)\varepsilon(1+G)) \times 2n \times 2 \le 8n(G + 2n\varepsilon)$ so the bound essentially depends linearly on the size of the vector. The factor of linearity is $8G$ so for the absolute error to remain small, the constant $G$ should be kept small, meaning that $E_{\exp}$ should be kept small. See below for examples of values that can be chosen.

To give concrete values for the bound, let us assume, as in Section 3.3, a reasonable bound on $n$ in practice, that is 1024. We assume the input numbers are smaller than $M_a = 3$ in absolute value, and we assume say that $|\rho| \le 1$, so that $|a_i + \rho - a_j| \le 7$ and then the argument of the exponential is between $-\frac{49}{2}$ and 0, thus fitting with the assumption $M_{\exp} \le 25$ taken in Section 3.3. Indeed, with such values the hypothesis 11 requires $M_{\exp}$ to satisfy $(2M_a + |\rho|)^2(\frac{1}{2} + A_\varepsilon) \le M_{\exp} - \frac{3}{2}\eta$
that is $M_{\text{exp}} \ge \frac{3}{2}\eta + 49(\frac{1}{2} + A_{\varepsilon})$. Let us assume first that we have some implementations of exponential and logarithm with correct rounding, that is $E_{\text{exp}} = 2^{-53}$ and $E_{\text{log}} = 2^{-53}$. With these values we have $A_{\varepsilon} = \frac{\varepsilon}{2}(1 + 4(1 + \varepsilon)^2) \le 3 \times 2^{-53}$, the constant $G$ satisfy $G = E_{\text{exp}} + (\exp(A_{\varepsilon}(2M_a + |\rho|)^2 + \frac{3}{2}\eta) - 1)(1 + E_{\text{exp}}) \le 2^{-53} + (\exp(49 \times 3 \times 2^{-53} + \frac{3}{2}2^{-1074}) - 1)(1 + 2^{-53}) \le 2^{-53}(1 + 2 \times 49 \times 3) \le 2^{-44}$. To ensure hypothesis 12 of Theorem 16 holds, it suffices to have $M_{\text{log}_2}$ larger than $n(1 + G + (n-1)\varepsilon(1+G)) \le 2^{10}(1 + 2^{-44}) + 2^{10} \times 2^{-53}(1 + 2^{-44}) \le 1025$. Indeed, there is nothing surprising here, as the sum taken as argument of $\log_2$ is a sum of at most 1024 numbers between 0 and 1. Thus, assuming that the implementation of $\log_2$ is correctly rounded on the domain given by the bound $M_{\text{log}_2}$ above, the relative error on SLSE is bounded by $E_{\text{log}_2} + (n-1)\varepsilon(1+E_{\text{log}_2}) \le 2^{10} \times 2\varepsilon = 2^{-43} = 1.14 \times 10^{-13}$ and the absolute error is bounded by $8n(G+2n\varepsilon) \le 2^{13}(2^{-44} + 2^{-42}) = 2.33 \times 10^{-9}$. Let us assume less precise implementations of exponential and logarithm with $E_{\text{exp}} = 2^{-40}$ and $E_{\text{log}_2} = 2^{-36}$. With these values the constant $G$ satisfy $G = E_{\text{exp}} + (\exp(A_{\varepsilon}(2M_a + |\rho|)^2 + \frac{3}{2}\eta) - 1)(1 + E_{\text{exp}}) \le 2^{-40} + (\exp(49 \times 3 \times 2^{-53} + \frac{3}{2}2^{-1074}) - 1)(1 + 2^{-40}) \le 2^{-40}(1 + 2 \times 49 \times 3) \le 2^{-31}$. To ensure hypothesis 12 of Theorem 16 holds, it suffices to have $M_{\text{log}_2}$ larger than $n(1+G+(n-1)\varepsilon(1+G)) \le 2^{10}(1+2^{-31})+2^{10} \times 2^{-53}(1+2^{-31}) \le 1025$. It should again not be surprising. The relative error on SLSE is bounded by $E_{\text{log}_2}+(n-1)\varepsilon(1+E_{\text{log}_2}) \le 2^{10} \times 2\varepsilon = 1.14 \times 10^{-13}$ so essentially the same as for correctly rounded exponential and logarithm, and the absolute error is bounded by $8n(G+2n\varepsilon) \le 2^{13}(2^{-31} + 2^{-42}) = 3.82 \times 10^{-6}$ which is this time very significantly larger than with the correctly rounded implementation of exponential.

A general conclusion is that the dominant part of the accuracy of SLSE is coming from the accuracy of the implementation of the exponential. It is also more or less linear in the size of the input vector.

## 4.4 Formalization in WhyML

As is the pen-and-paper proof, the WhyML formalization of the accuracy of SLSE is significantly longer and more complex than those of LSE. We present here only an excerpt of that formalization, focusing on the important parts.

Figure 10 presents the formalization up to Lemma 12. The proof is 61 lines long. It generates 50 VCs. Among them, 6 are proved by Alt-Ergo, 3 by CVC4, 39 by cvc5 and 2 by Z3.

Finally, Figure 11 presents the formalization of Theorem 16. The text of the WhyML proof of g\_accuracy is 61 lines long. The proof amounts to discharge 50 VCs, 6 are proved by Alt-Ergo, 3 by CVC4, 39 by cvc5 and 2 by Z3. The text of the WhyML proof of slse\_accuracy is pages long. The proof amounts to discharge 52 VCs, 19 are proved by Alt-Ergo, 31 by cvc5 and 2 by Z3. Notice additionally that we need an auxiliary mathematical lemma to establish that $\exp(x) \le 1 + 2x$ when $0 \le x \le 1$, that lemma is proved using Rocq and its interval tactic [52]. As before, see the details on the dedicated page of the Toccata gallery [51].
let function sqr (x:real) : real = x × x
constant a.eps : real = eps × (1.0 + 4.0 × sqr (1.0 + eps)) / 2.0

predicate g.hypothesis (rho:udouble) (a_max:real) =
sqr (2.0 × a_max + abs (to_real rho)) × (0.5 + a.eps) + 1.5 × eta
≤ exp_max_value

let ghost g accuracy (x y rho : udouble) (a_max : real)
requires { abs (to_real x) ≤ a_max }
requires { abs (to_real y) ≤ a_max }
requires { g.hypothesis rho a_max }
ensures {
  let t1 = abs (to_real x) + abs (to_real y) + abs (to_real rho) in
  let t2 = a.eps × sqr t1 + 1.5 × eta in
  let error = exp_error + (exp t2 - 1.0) × (1. + exp_error) in
  abs (to_real (u_g x y rho) - g x y rho) ≤ error × abs (g x y rho)
}

<caption>Fig. 10 Accuracy result on function *g*, formalizing Lemma 12.</caption>

function g_max_error (rho:udouble) (a_max:real) : real =
  exp_error
+ (exp (a.eps × sqr (2.0 × a_max + abs(to_real rho)) + 1.5 × eta) − 1.0)
  × (1.0 + exp_error)

predicate slse_hypothesis (rho:udouble) (a_max:real) (n:int) =
  let gg = g_max_error rho a_max in
  from_int n × (1.0 + (gg + from_int (n-1) × eps × (1.0 + gg)))
    ≤ log2_max_value

let ghost slseAccuracy (a:int → udouble) (n:int) (rho:udouble) (a_max:real)
  requires { 1 <= n }
  requires { from_int (n-1) × eps ≤ 0x1p-4 }
  requires {forall i. 0 <= i < n → abs (to_real (a i)) ≤ a_max }
  requires { g_hypothesis rho a_max }
  requires { slse_hypothesis rho a_max n }
  ensures {
    let gg = g_max_error rho a_max in
    abs (to_real (u_slse a n rho) - slse a n rho) ≤
      (log_error + from_int (n-1) × eps × (1.0 + log2_error)) × abs_slse a n rho
    - log2 (1.0 - (gg + from_int (n-1) × eps × (1.0 + gg)))
      × (1.0 + log2_error) × from_int n × (1.0 + from_int (n-1) × eps)
  }

<caption>Fig. 11 Accuracy result on SLSE.</caption>
# 5 Conclusions

We presented bounds on the accuracy of implementations of the LSE and the SLSE functions. The resulting expressions for the bounds are parametric in the precision of the underlying implementations of exp and log, and also parameterized by bounds on the size of the argument vectors, and bounds on the values of the vectors components. The given bounds are proved on paper and then in WhyML, using Why3 constructs such as lemma functions to provide proofs that follow more or less the paper proofs. We also proved some C implementation by reusing the results proved in WhyML. The proofs are made simpler by using a theory of unbounded floating-point numbers, allowing us to separate the reasoning on accuracy from the reasoning on absence of overflow.

## 5.1 Related Work

As far as we know, the only contribution focusing on roundoff errors of LSE-based algorithms is the work of Blanchard et al. [10]. They bound the roundoff errors of the LSE function and its gradient, namely the softmax function. In this work, the authors assume that the exponential and logarithm functions are implemented with correct rounding, i.e., $E_{\text{exp}} = E_{\text{log}} = \varepsilon$. In contrast, we provide a proof which is parameterized with arbitrary error bounds for the called functions. It means that we can rely on any implementations of these functions without invalidating the bounds. Going back to Theorem 6, if we take $E_{\text{exp}} = \varepsilon$, we get a bound in which the relative error term is $E_{\text{log}} = \varepsilon$ and the constant term is $-\log(1-(\varepsilon+(n-1)\varepsilon)(1+\varepsilon)) (1+\varepsilon)$, that is, about $-\log(1-n\varepsilon)+O(u^2)$. Blanchard et al.'s relative error term is identical. Their constant term is about $(n+1)\varepsilon+O(u^2)$, which is actually very slightly tighter than ours, but of comparable order of magnitude. Yet, a strength of our accuracy result is that is parametric in the accuracies of exp and log, instead of assuming ideally precise implementations. Moreover, another main strength compared to this work is the fact that we made formal proofs.

Our work can also be compared to other contributions targeting the end-to-end formal proof of numerical software written in C. For instance, Appel and Bertot [3] have combined the Verifiable Software Toolchain (VST) [2, 4], Flocq [19, 20] and Gappa [31] to formally-verify an example of square root implementation using the Newton method. VST ensures the correctness of the C code, while Flocq and Gappa are used as backend tools to check numerical accuracy facts. Kellison et al. [47] propose a Rocq formal proofs library, called LAProof, for roundoff error analysis of basic linear algebra operations, e.g. inner product or matrix-matrix multiplication. As an application example, the authors prove a C program computing a sparse matrix-vector multiplication using the VST [2] approach and the LAProof library.

Boldo et al. [13, 14] formalized a numerical integration scheme for a wave partial differential equation in Rocq. They not only formally proved a bound on the mathematical errors [14], but also a bound on the roundoff errors [13]. These results have been used as a basis for the formal verification of a wave equation resolution C program [15], based on the Jessie [57] plug-in of Frama-C. The most complex proof obligations related to numerical errors are discharged to Rocq.
Becker et al. [8] developed a CakeML extension for optimizing floating-point arithmetic in Standard ML. Their approach relies on an end-to-end soundness proof linking a real-number specification of the initial code with the code obtained after optimization. The approach is entirely automated (code and proof generation) and targets the optimization of floating-point kernels, which are essentially blocks of floating-point computations free of control flow instructions. The roundoff errors are obtained and proved by using the prover FloVer (interval-based prover in HOL4).

## 5.2 Future Work

The bounds exhibited by Blanchard et al. [10] show that naive implementations of these functions are relatively well-behaved regarding numerical accuracy, but prone to spurious overflow. The authors then study an alternative implementation to bypass this issue. The principle is to find the maximal value $\alpha = \max(a_i)$ among the components of the input vector and to rewrite the LSE expression as $\text{LSE}(a) = \alpha + \log\left(\sum_{1 \le i \le n} \exp(a_i - \alpha)\right)$. Since for all $i$, $a_i - \alpha \le 0$, the exponential takes a reasonable value whatever is the magnitude of the components of $a$, therefore, the risk of overflow is limited. They also prove that the accuracy of this alternative evaluation is not only as good as with the standard evaluation, but even slightly better. In practice, most applications using the LSE function rely on the shifted version. The proof of the error bounds associated to this alternative evaluation which is presented by Blanchard et al. uses rather sophisticated mathematical arguments, e.g. Taylor series expansions of the $\log(1+x)$ quantity. Providing a formal proof of this result could be an interesting perspective.

In the presented work, some parts of the proof are performed with a high level of automation. Making formal verification processes automatic and push-button has a strong impact on their industrial applicability. While the studied example is rather intricate and would be difficult to fully automate, there are plenty of applicable source code in which simple combinations of mathematical functions calls appear. For instance, we could try to apply our methodology on benchmarks from the FPBench [27] or COPRIN projects [55]. Most examples from these benchmarks are loop-free which strongly eases the verification process. We could try to handle these examples in a fully automatic way, providing bounds depending on error bounds for the mathematical functions implementations appearing in the source code.

For now, we assume error bounds on the implementations of the mathematical functions $\log$ and $\exp$. Our formally proved bounds apply only when implementations satisfying the assumptions are provided. It is known in the literature that correctly-rounded implementations of these functions can be achieved, e.g. in the recent CORE-MATH library [60]. To complement our work, we envision the formal verification of the errors induced by such implementations. In the late 20th century, Harrison [41, 42] formally verified implementations of specific floating-point exponential and trigonometric functions implementations in HOL, but the proofs were *ad hoc*, low-level and far from automatic. More recently, the Gappa tool [31] has been partly used to bound roundoff errors of floating-point implementations of functions from the CR-LIBM library [28, 33, 32]. However, proofs are not fully automatic and are
devoted to specific implementations. In addition, these works only focus on roundoff
errors, without taking the mathematical approximation errors into account. Geneau
de Lamarlière et al. [39] provide a methodology and tooling to ease the formal proofs
of low-level floating-point components with a minimal user effort. Their approach
relies on a framework for modeling and reasoning on floating-point expressions with
some facilities, without neglecting potential exceptional behaviors. For that purpose,
they offer tools in the Rocq proof assistant to automate the proof of the absence of
exceptional behaviors, so that the user can reason on real numbers representation of
floating-point expressions. This work is typically applied on mathematical functions
implementations, e.g. exp and log. Combining our approach with their methodology
would be valuable to complete the toolchain. A longer-term goal could be the develop-
ment of a formally verified implementation synthesis tool, in the spirit of the Metalibm
tools [49, 25]. Metalibm already enables the generation of Gappa scripts certifying the
synthesized implementations, but this feature is limited to some pieces of code.

Our work is not limited to the LSE and SLSE functions. We already applied our
methodology to others, including some involving other mathematical functions like
sine and cosine. Achieving formal proofs becomes quickly quite a burden, and we seek
to improve our methodology but introducing more automation. We started a work
aimed at automating the application of the “error propagation lemmas” as they are
introduced in Section 4.1 above [21]. We believe that automating the application of
such lemmas would naturally be reused for proving any code proceeding by composing
numerical functions and operations. Another concern is related to the methodology
to deal directly with C code. Our current methodology is far from being usable by
non-expert users. This is illustrated for example by the numerous assertions that we
had to add in our C code for LSE, on lines 31–39 of Figure 8. There are constructs
available in WhyML that would be nice to have at the C level: we think in particular,
on one hand, about arbitrary lambda-expressions, and on the other hand the ability
to call ghost functions. In fact, ghost functions are in principle present in ACSL [9],
but they are limited to ghost C programs, whereas we would need to have ghost logic
functions that would accept logic types are parameters: these include real numbers,
unbounded floats, functions (lambda-expressions), etc.

Acknowledgments

This research was supported partly by the bilateral contract ProofInUse-MERCE and
partly by the Décysif project funded by the Île-de-France region and by the French
government in the context of “Plan France 2030”.

References

[1] Behzad Akbarpour and Lawrence C. Paulson. Metitarski: An automatic theorem
prover for real-valued special functions. Journal of Automated Reasoning,
44(3):175–205, 2010. Tool page at http://www.cl.cam.ac.uk/~lp15/papers/
Arith/. doi:10.1007/s10817-009-9149-2.

[2] Andrew Appel. Verified software toolchain. In European Symposium on Programming, volume 6602 of Lecture Notes in Computer Science, pages 1–17.
Springer, 2011. Tool page at https://vst.cs.princeton.edu/. doi:10.5555/1987211.1987212.

[3] Andrew Appel and Yves Bertot. C-language floating-point proofs layered with VST and Flocq. *Journal of Formalized Reasoning*, 13(1):1–16, 2020. URL: https://inria.hal.science/hal-03130704/.

[4] Andrew Appel and Ariel Kellison. VCFloat2: Floating-point error analysis in Coq. In *13th ACM SIGPLAN International Conference on Certified Programs and Proofs*, pages 14–29, 2024. doi:10.1145/3636501.3636953.

[5] Ali Ayad and Claude Marché. Multi-prover verification of floating-point programs. In Jürgen Giesl and Reiner Hähnle, editors, *Fifth International Joint Conference on Automated Reasoning*, volume 6173 of *Lecture Notes in Artificial Intelligence*, pages 127–141, Edinburgh, Scotland, July 2010. Springer. URL: http://hal.inria.fr/inria-00534333.

[6] Haniel Barbosa, Clark W. Barrett, Martin Brain, Gereon Kremer, Hanna Lachnitt, Makai Mann, Abdalrhman Mohamed, Mudathir Mohamed, Aina Niemetz, Andres Nötzli, Alex Ozdemir, Mathias Preiner, Andrew Reynolds, Ying Sheng, Cesare Tinelli, and Yoni Zohar. cvc5: A versatile and industrial-strength SMT solver. In Dana Fisman and Grigore Rosu, editors, *Tools and Algorithms for the Construction and Analysis of Systems*, volume 13243 of *Lecture Notes in Computer Science*, pages 415–442. Springer, 2022. doi:10.1007/978-3-030-99524-9_24.

[7] Clark Barrett, Christopher L. Conway, Morgan Deters, Liana Hadarean, Dejan Jovanović, Tim King, Andrew Reynolds, and Cesare Tinelli. CVC4. In *Computer Aided Verification*, volume 6806 of *Lecture Notes in Computer Science*, pages 171–177. Springer, 2011. http://cvc4.cs.stanford.edu/web/. doi:10.1007/978-3-642-22110-1_14.

[8] Heiko Becker, Robert Rabe, Eva Darulova, Magnus O. Myreen, Zachary Tatlock, Ramana Kumar, Yong Kiam Tan, and Anthony C. J. Fox. Verified compilation and optimization of floating-point programs in cakeml. In Karim Ali and Jan Vitek, editors, *36th European Conference on Object-Oriented Programming*, volume 222 of *LIPIcs*, pages 1:1–1:28, Berlin, Germany, 2022. Schloss Dagstuhl - Leibniz-Zentrum für Informatik. doi:10.4230/LIPICS.ECOOP.2022.1.

[9] Allan Blanchard, Claude Marché, and Virgile Prevosto. *Guide to Software Verification with Frama-C — Core Components, Usages, and Applications*, chapter Formally Expressing what a Program Should Do: the ACSL Language, pages 3–80. Springer-Verlag, 2024. URL: https://inria.hal.science/hal-04265707, doi:10.1007/978-3-031-55608-1_1.

[10] Pierre Blanchard, Desmond J Higham, and Nicholas J Higham. Accurately computing the log-sum-exp and softmax functions. *IMA Journal of Numerical Analysis*, 41(4):2311–2330, 2021. doi:10.1093/imanum/draa038.

[11] François Bobot, Jean-Christophe Filliâtre, Claude Marché, and Andrei Paskevich. Let's verify this with Why3. *International Journal on Software Tools for Technology Transfer (STTT)*, 17(6):709–727, 2015. See also http://toccata.gitlabpages.inria.fr/toccata/gallery/fm2012comp.en.html. URL: http://hal.inria.fr/hal-00967132/en, doi:10.1007/s10009-014-0314-5.
[12] Sylvie Boldo. Floats & Ropes: a case study for formal numerical program verification. In *36th International Colloquium on Automata, Languages and Programming*, volume 5556 of *Lecture Notes in Computer Science*, pages 91–102, Rhodos, Greece, July 2009. Springer. doi:10.1007/978-3-642-02930-1_8.

[13] Sylvie Boldo. Floats and ropes: A case study for formal numerical program verification. In *36th International Colloquium on Automata, Languages and Programming*, volume 5556 of *Lecture Notes in Computer Science*, pages 91–102. Springer, 2009. doi:10.1007/978-3-642-02930-1_8.

[14] Sylvie Boldo, François Clément, Jean-Christophe Filliâtre, Micaela Mayero, Guillaume Melquiond, and Pierre Weis. Formal proof of a wave equation resolution scheme: the method error. In *Interactive Theorem Proving*, volume 6172 of *Lecture Notes in Computer Science*, pages 147–162. Springer, 2010. URL: http://hal.inria.fr/inria-00450789/en, doi:10.1007/978-3-642-14052-5_12/.

[15] Sylvie Boldo, François Clément, Jean-Christophe Filliâtre, Micaela Mayero, Guillaume Melquiond, and Pierre Weis. Wave equation numerical resolution: a comprehensive mechanized proof of a C program. *Journal of Automated Reasoning*, 50(4):423–456, April 2013. URL: http://hal.inria.fr/hal-00649240/en/, doi:10.1007/s10817-012-9255-4.

[16] Sylvie Boldo and Jean-Christophe Filliâtre. Formal verification of floating-point programs. In *18th IEEE International Symposium on Computer Arithmetic*, pages 187–194, 2007. URL: https://usr.lmf.cnrs.fr/~jcf/publis/caduceus-floats.pdf, doi:10.1109/ARITH.2007.20.

[17] Sylvie Boldo, Claude-Pierre Jeannerod, Guillaume Melquiond, and Jean-Michel Muller. Floating-point arithmetic. *Acta Numerica*, 32:203–290, 2023. URL: https://hal.science/hal-04095151, doi:10.1017/S0962492922000101.

[18] Sylvie Boldo and Claude Marché. Formal verification of numerical programs: from C annotated programs to mechanical proofs. *Mathematics in Computer Science*, 5:377–393, 2011. URL: http://hal.inria.fr/hal-00777605, doi:10.1007/s11786-011-0099-9.

[19] Sylvie Boldo and Guillaume Melquiond. *Computer Arithmetic and Formal Proofs: Verifying Floating-point Algorithms with the Coq System*. ISTE Press - Elsevier, December 2017. URL: https://hal.inria.fr/hal-01632617.

[20] Sylvie Boldo and Guillaume Melquiond. Some formal tools for computer arithmetic: Flocq and Gappa. In Mioara Joldes and Fabrizio Lamberti, editors, *28th IEEE International Symposium on Computer Arithmetic*, 2021. URL: https://hal.inria.fr/hal-03233227.

[21] Paul Bonnot, Benoît Boyer, Florian Faissole, and Claude Marché. Generating and certifying accuracy properties of floating-point programs. Technical Report RR-9564, Inria, 2024. URL: https://inria.hal.science/hal-04820735.

[22] Paul Bonnot, Benoît Boyer, Florian Faissole, Claude Marché, and Raphaël Rieu-Helft. Formally verified bounds on rounding errors in concrete implementations of logarithm-sum-exponential functions. Research Report 9531, Inria, 2023. URL: https://inria.hal.science/hal-04343157.

[23] Paul Bonnot, Benoît Boyer, Florian Faissole, Claude Marché, and Raphaël
Rieu-Helft. Formally verified rounding errors of the logarithm-sum-exponential function. In *Formal Methods in Computer-Aided Design*. IEEE, 2024. URL: https://inria.hal.science/hal-04674600.

[24] Sven Brüggemann and Corrado Possieri. On the use of difference of log-sum-exp neural networks to solve data-driven model predictive control tracking problems. *IEEE Control Systems Letters*, 5(4):1267-1272, 2020. doi:10.1109/LCSYS. 2020.3032083.

[25] Nicolas Brunie, Christoph Lauter, and Guillaume Revy. Precision adaptation for fast and accurate polynomial evaluation generation. In *30th International Conference on Application-specific Systems, Architectures and Processors*, volume 2160-052X, pages 41–41, New York, NY, USA, 2019. IEEE. doi:10.1109/ASAP. 2019.00-32.

[26] Sylvain Conchon, Mohamed Iguernlala, Kailiang Ji, Guillaume Melquiond, and Clément Fumex. A three-tier strategy for reasoning about floating-point numbers in SMT. In *Computer Aided Verification*, volume 10427 of *Lecture Notes in Computer Science*, pages 419–435, 2017. URL: https://hal.inria.fr/hal-01522770, doi:10.1007/978-3-319-63390-9_22.

[27] Nasrine Damouche, Matthieu Martel, Pavel Panchekha, Chen Qiu, Alexander Sanchez-Stern, and Zachary Tatlock. Toward a standard benchmark format and suite for floating-point analysis. In *Numerical Software Verification*, pages 63–77. Springer, 2017. doi:10.1007/978-3-319-54292-8_6.

[28] Catherine Daramy, David Defour, Florent de Dinechin, and Jean-Michel Muller. Cr-libm: a correctly rounded elementary function library. In *Advanced Signal Processing Algorithms, Architectures, and Implementations XIII*, volume 5205, pages 458–464. SPIE, 2003.

[29] Catherine Daramy-Loirat, David Defour, Florent de Dinechin, Matthieu Gallet, Nicolas Gast, Christoph Lauter, and Jean-Michel Muller. CR-LIBM: a library of correctly rounded elementary functions in double-precision. Research report, LIP, 2006. URL: https://ens-lyon.hal.science/ensl-01529804.

[30] Marc Daumas and Guillaume Melquiond. Certification of bounds on expressions involving rounded operators. *Transactions on Mathematical Software*, 37(1):1–20, 2010. URL: http://www.lri.fr/~melquion/doc/09-toms.pdf, doi:10.1145/ 1644001.1644003.

[31] Marc Daumas and Guillaume Melquiond. Certification of bounds on expressions involving rounded operators. *ACM Transactions on Mathematical Software (TOMS)*, 37(1):1–20, 2010.

[32] Florent De Dinechin, Christoph Lauter, and Guillaume Melquiond. Certifying the floating-point implementation of an elementary function using gappa. *IEEE Transactions on Computers*, 60(2):242–253, 2010.

[33] Florent De Dinechin, Christoph Quirin Lauter, and Guillaume Melquiond. Assisted verification of elementary functions using gappa. In *Proceedings of the 2006 ACM symposium on Applied computing*, pages 1318–1322, 2006.

[34] Leonardo de Moura and Nikolaj Bjørner. Z3, an efficient SMT solver. In *TACAS*, volume 4963 of *Lecture Notes in Computer Science*, pages 337–340. Springer, 2008. doi:10.1007/978-3-540-78800-3_24.
[35] Xavier Denis, Jacques-Henri Jourdan, and Claude Marché. Creusot: a foundry for the deductive verification of Rust programs. In *International Conference on Formal Engineering Methods - ICFEM*, Lecture Notes in Computer Science, Madrid, Spain, 2022. Springer. URL: https://hal.inria.fr/hal-03737878.

[36] Clément Fumex, Claude Marché, and Yannick Moy. Automating the verification of floating-point programs. In Andrei Paskevich and Thomas Wies, editors, *Verified Software: Theories, Tools, and Experiments. Revised Selected Papers Presented at the 9th International Conference VSTTE*, number 10712 in Lecture Notes in Computer Science, Heidelberg, Germany, December 2017. Springer. URL: https://hal.inria.fr/hal-01534533/.

[37] Bolin Gao and Lacra Pavel. On the properties of the softmax function with application in game theory and reinforcement learning. *arXiv preprint arXiv:1704.00805*, 2017.

[38] Sicun Gao, Jeremy Avigad, and Edmund M. Clarke. $\delta$-complete decision procedures for satisfiability over the reals. In Bernhard Gramlich, Dale Miller, and Uli Sattler, editors, *Automated Reasoning*, pages 286–300. Springer, 2012. doi:978-3-642-31365-3_23.

[39] Paul Geneau de Lamarlière, Guillaume Melquiond, and Florian Faissole. Slimmer formal proofs for mathematical libraries. In Theo Drane and Anastasia Volkova, editors, *Proceedings of the 30th IEEE International Symposium on Computer Arithmetic*, Portland, OR, USA, September 2023.

[40] Cédric Gernigon, Silviu-Ioan Filip, Olivier Sentieys, Clément Coggiola, and Mickaël Bruno. Low-precision floating-point for efficient on-board deep neural network processing, 2023. arXiv:2311.11172.

[41] John Harrison. Floating point verification in hol light: the exponential function. In *International Conference on Algebraic Methodology and Software Technology*, pages 246–260. Springer, 1997.

[42] John Harrison. Formal verification of floating point trigonometric functions. In *International conference on formal methods in computer-aided design*, pages 254–270. Springer, 2000.

[43] Nicholas J Higham. *Accuracy and stability of numerical algorithms*. SIAM, 2002. doi:10.1137/1.9780898718027.

[44] Taocheng Hu and Jinhui Yu. LogSumExp for unlabeled data processing. In *15th International Conference on Software Engineering Research, Management and Applications*, pages 63–69. IEEE, 2017. doi:10.1109/SERA.2017.7965708.

[45] IEEE standard for floating-point arithmetic, 2008. doi:10.1109/IEEESTD.2008.4610935.

[46] Claude-Pierre Jeannerod and Siegfried M. Rump. On relative errors of floating-point operations: optimal bounds and applications. *Mathematics of Computation*, 87:803–819, 2018. URL: https://hal.inria.fr/hal-00934443, doi:10.1090/mcom/3234.

[47] Ariel E. Kellison, Andrew W. Appel, Mohit Tekriwal, and David S. Bindel. LAProof: A library of formal proofs of accuracy and correctness for linear algebra programs. In *30th IEEE International Symposium on Computer Arithmetic*, 2023. doi:10.1109/ARITH58626.2023.00021.
[48] Florent Kirchner, Nikolai Kosmatov, Virgile Prevosto, Julien Signoles, and Boris Yakobowski. Frama-c: A software analysis perspective. *Formal Aspects of Computing*, 27(3):573–609, May 2015. doi:10.1007/s00165-014-0326-7.

[49] Olga Kupriianova and Christoph Lauter. Metalibm: A mathematical functions code generator. In *International Congress on Mathematical Software*, volume 8592 of *Lecture Notes in Computer Science*, pages 713–717, Berlin, Heidelberg, 2014. Springer. doi:10.1007/978-3-662-44199-2_106.

[50] Radek Mackowiak, Lynton Ardizzone, Ullrich Kothe, and Carsten Rother. Generative classifiers as a basis for trustworthy image classification. In *Computer Vision and Pattern Recognition*, pages 2971–2981, 2021. doi:10.1109/CVPR46437. 2021.00299.

[51] Claude Marché. Toccata gallery of verified programs, example “LogSumExp-based computations”. https://toccata.gitlabpages.inria.fr/toccata/gallery/log_sum_exp_computations.en.html, 2025.

[52] Érik Martin-Dorel and Guillaume Melquiond. Proving tight bounds on univariate expressions with elementary functions in Coq. *Journal of Automated Reasoning*, 2016. URL: https://hal.inria.fr/hal-01086460, doi:10.1007/s10817-015-9350-4.

[53] John W. McCormick and Peter C. Chapin. *Building High Integrity Applications with SPARK*. Cambridge University Press, 2015. doi:10.1017/CB09781139629294.

[54] Guillaume Melquiond. *Formal Verification for Numerical Computations, and the Other Way Around*. Habilitation à diriger des recherches, Université Paris Sud, April 2019. URL: https://tel.archives-ouvertes.fr/tel-02194683.

[55] Jean-Pierre Merlet. *Parallel Robots*, chapter Structural synthesis and architectures, pages 19–94. Springer, 2006. doi:10.1007/1-4020-4133-0_2.

[56] Taiki Miyagawa and Akinori F Ebihara. The power of log-sum-exp: Sequential density ratio matrix estimation for speed-accuracy optimization. In *International Conference on Machine Learning*, volume 139 of *Proceedings of Machine Learning Research*, pages 7792–7804, 2021. URL: https://proceedings.mlr.press/v139/miyagawa21a.html.

[57] Yannick Moy and Claude Marché. Modular inference of subprogram contracts for safety checking. *Journal of Symbolic Computation*, 45:1184–1211, 2010. URL: http://hal.inria.fr/inria-00534331/en/, doi:10.1016/j.jsc.2010.06.004.

[58] Jean-Michel Muller, Nicolas Brunie, Florent de Dinechin, Claude-Pierre Jéan-nerod, Mioara Joldes, Vincent Lefèvre, Guillaume Melquiond, Nathalie Revol, and Serge Torres. *Handbook of Floating-point Arithmetic (2nd edition)*. Birkhäuser Basel, July 2018. URL: https://hal.inria.fr/hal-01766584, doi:10.1007/978-3-319-76526-6.

[59] Pierre Roux. Formal Proofs of Rounding Error Bounds - With Application to an Automatic Positive Definiteness Check. *Journal of Automated Reasoning*, 57(2):135–156, 2016. doi:10.1007/s10817-015-9339-z.

[60] Alexei Sibidanov, Paul Zimmermann, and Stéphane Glondu. The CORE-MATH project. In *29th IEEE Symposium on Computer Arithmetic*, pages 26–34, 2022. URL: https://inria.hal.science/hal-03721525, doi:10.1109/ARITH54963.2022.00014.
[61] US Government Accountability Office. Defense patriot missile: Software problem led to system failure at dhahran, saudi arabia. *US Government Accountability Office Reports*, *rapport no. GAO/IMTEC-92-26*, 1992.