Sylvie Boldo
email: sylvie.boldo@inria.fr

Jacques-Henri Jourdan
email: jacques-henri.jourdan@inria.fr

Xavier Leroy
email: xavier.leroy@inria.fr

Guillaume Melquiond
email: guillaume.melquiond@inria.fr

A Formally-Verified C Compiler Supporting Floating-Point Arithmetic

Keywords: floating-point arithmetic, verified compilation, formal proof, floating-point semantic preservation

Floating-point arithmetic is known to be tricky: roundings, formats, exceptional values. The IEEE-754 standard was a push towards straightening the field and made formal reasoning about floating-point computations easier and flourishing. Unfortunately, this is not sufficient to guarantee the final result of a program, as several other actors are involved: programming language, compiler, architecture. The CompCert formally-verified compiler provides a solution to this problem: this compiler comes with a mathematical specification of the semantics of its source language (a large subset of ISO C90) and target platforms (ARM, PowerPC, x86-SSE2), and with a proof that compilation preserves semantics. In this paper, we report on our recent success in formally specifying and proving correct CompCert's compilation of floating-point arithmetic. Since CompCert is verified using the Coq proof assistant, this effort required a suitable Coq formalization of the IEEE-754 standard; we extended the Flocq library for this purpose. As a result, we obtain the first formally verified compiler that provably preserves the semantics of floating-point programs.

I. INTRODUCTION

Use and study of floating-point (FP) arithmetic have intensified since the 70s [START_REF] Sterbenz | Floating point computation[END_REF], [START_REF] Dekker | A floating point technique for extending the available precision[END_REF]. At that time, computations were not standardized and, due to the differences between architectures, the use of the same source program in different contexts gave different results. Since the IEEE-754 standard of 1985 and its revision in 2008 [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF], things should have changed as reproducibility was a keyword. Each basic operation is guaranteed to be computed as if the computation was done with infinite precision and then rounded. The goal was that the same program could be run on various platforms and give the same result. It allowed the development of many algorithms coming with mathematical proofs based on the fact that operations were correctly rounded. Since the 2000s, this was even pushed to formal proofs of algorithms or hardware components: in PVS [START_REF] Carreño | Specification of the IEEE-854 floatingpoint standard in HOL and PVS[END_REF], in ACL2 [START_REF] Russinoff | A mechanically checked proof of IEEE compliance of the floating point multiplication, division and square root algorithms of the AMD-K7 processor[END_REF], in HOL-light [START_REF] Harrison | Formal verification of floating point trigonometric functions[END_REF] and in Coq [START_REF] Boldo | Preuves formelles en arithmétiques à virgule flottante[END_REF], [START_REF] Boldo | Flocq: A unified library for proving floating-point algorithms in Coq[END_REF]. The basic axiom for algorithms and the basic goal for hardware components was still that all the operations are correctly rounded.

To complicate matters further, the processor architecture is not the only party responsible for the computed results. Stand also accused the programming language and the compiler used. We will focus on the compiler, as it can deviate from what the programmer wants or what was proved from the written code.

To illustrate what the compiler can change, here is a small example in C: int main () { double y, z; y = 0x1p-53 + 0x1p-78; // y = 2 -53 + 2 -78 z = 1. + y -1. -y; printf("%a\n", z); return 1; } Experts may have recognized a Fast-Two-Sum [START_REF] Dekker | A floating point technique for extending the available precision[END_REF] that computes the round-off error of a FP addition by ((a+b)-a)-b for |a| ≥ |b|. This very simple program compiled with GCC 4.6.3 gives three different answers on an x86 32-bit architecture depending on the chosen level of optimization.

Optimization level Program result -O0

-0x1p-78 -O1, -O2, -O3 0x1.ffffffp-54 -Ofast 0x0p+0

How can we explain the three results? For the first two results, the answer lies in the x86 architecture: it may compute with double precision (64 bits, 53 bits of precision) or with extended precision (80 bits, 64 bits of precision). For each operation, the compiler may choose to round the infinitelyprecise result either to extended precision, or to double precision, or first to extended and then to double precision. The latter is called a double rounding. In all cases, y is computed exactly: y = 2 -53 + 2 -78 .

With the -O0 optimization, all the computations are performed with extended precision and rounded in double precision only once at the end. With -O1 and higher, the result (1 + y) -1 is pre-computed by the compiler and the program only computes the last subtraction and prints the value. With -Ofast, there is no computation at all in the program but only the output of the constant 0. This optimization level turns on -funsafe-math-optimizations which allows the reordering of FP operations. It is explicitly stated in GCC documentation that this option "can result in incorrect output for programs which depend on an exact implementation of IEEE or ISO rules/specifications for math functions".

Another possible discrepancy comes from the use of the fused-multiply-add operator (FMA). For example, consider a× b + c × d and let us denote the rounding operator by •. When a FMA is available, the compiler may choose either

•(a × b + •(c × d)), or •(•(a × b) + c × d), or •(•(a × b) + •(c × d))
which may give different results. A wide set of examples of strange FP behaviors can be found in [START_REF] Monniaux | The pitfalls of verifying floating-point computations[END_REF], [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF].

As surprising as it may seem, all the discrepancies described so far are allowed by the ISO C standard [START_REF] Iso | International standard ISO/IEC 9899:2011, Programming languages -C[END_REF], which leaves much freedom to the compiler in the way it implements FP computations. Sometimes, optimizing compilers take additional liberties with the source programs, generating executable code that exhibits behaviors not allowed by the specification of the source language. This is called miscompilation. Consider the following example, adapted from GCC's infamous "bug #323":

void test(double x, double y) { const double y2 = x + 1.0; if (y != y2) printf("error\n"); } int main() {
const double x = .012; const double y = x + 1.0; test(x, y); return 0; } For an x86 32-bit target at optimization level -O1, all versions of GCC prior to 4.5 miscompile this code as follows: the expression x + 1.0 in function test is computed in extended precision, as allowed by C, but the compiler omits to round it back to double precision when assigning to y2, as prescribed by the C standard. Consequently, y and y2 compare different, while they must be equal according to the C standard. Miscompilation happens more often than one may think: Yang et al [START_REF] Yang | Finding and understanding bugs in C compilers[END_REF] tested many production-quality C compilers using differential random testing, and found hundreds of cases where the compiler either crashes at compile-time or-much worsesilently generates an incorrect executable from a correct source program.

As the compiler gives so few guarantees on how it implements FP arithmetic, it therefore seems impossible to guarantee the result of a program. In fact, most analysis of FP programs assume correct compilation and a strict application of the IEEE-754 standard where no extended registers nor FMA are used. This assumption is correct for embedded software such as those used in avionics. For the automatic analysis of C programs, a successful approach is based on abstract interpretation, and tools include Astre [START_REF] Cousot | The ASTR ÉE analyzer[END_REF], [START_REF] Monniaux | Analyse statique : de la théorie à la pratique[END_REF] and Fluctuat [START_REF] Delmas | Towards an industrial use of FLUCTUAT on safety-critical avionics software[END_REF]. Another method to specify and prove behavioral properties of FP programs is deductive verification system: specification languages have to take into account FP arithmetic. This has been done for Java in JML [START_REF] Leavens | Not a number of floating point problems[END_REF],

for C in ACSL [START_REF] Boldo | Formal Verification of Floating-Point Programs[END_REF], [START_REF] Ayad | Multi-prover verification of floating-point programs[END_REF]. However, all these works only follow strictly the IEEE-754 standard, with neither FMA, nor extended registers, nor considering optimization aspects. Recently, several possibilities have been offered to take these aspects into account. One approach is to cover all the ways a compiler may have compiled each FP operation and to compute an error bound that stands correct whatever the compiler choices [START_REF] Boldo | Proofs of numerical programs when the compiler optimizes[END_REF]. Another approach is to analyze the assembly code to get all the precision information [START_REF] Nguyen | Hardware-dependent proofs of numerical programs[END_REF].

Our approach is different: rather than trying to account for all the changes a compiler may have silently introduced in a FP program, we have focused on getting a correct and predictable compiler that supports FP arithmetic. Concerning compilers and how to make them more trustworthy, Milner and Weyrauch [START_REF] Milner | Proving compiler correctness in a mechanized logic[END_REF] were the first to mechanically prove the correctness of a compiler, although for a very simple language of expressions. Moore [START_REF] Moore | A mechanically verified language implementation[END_REF] extended this approach to an implementation of the Piton programming language. Li et al [START_REF] Li | Structure of a proof-producing compiler for a subset of higher order logic[END_REF] showed that one can compile programs with proof, directly from the logic of the HOL4 theorem prover. A year later, Myreen [START_REF] Myreen | Formal verification of machine-code programs[END_REF] made contributions both to approaches for verification of programs and methods for automatically constructing correct code.

To build our compiler, we started from CompCert [START_REF] Leroy | Formal verification of a realistic compiler[END_REF], a formally-verified compiler described in Section III and extended it with FP arithmetic. As CompCert is developed using the Coq proof assistant, we had to build on a Coq library formalizing FP arithmetic: we relied on the Flocq library [START_REF] Boldo | Flocq: A unified library for proving floating-point algorithms in Coq[END_REF] and extended it to serve the needs of a verified compiler. With all these components, we were able to get a correct, predictable compiler that conforms strictly to the IEEE-754 standard.

In this article, we present in Section II the semantics of FP arithmetic in programs, depending in particular on the programming language. In Section III, we describe the CompCert certified compiler. We will explain in Section IV the required additions to Flocq to represent all IEEE-754 FP numbers. In Section V, we detail what modifications to CompCert were needed to handle FP arithmetic.

II. SEMANTICS OF FLOATING-POINT ARITHMETIC

Starting from an algorithm using FP arithmetic, there is a long road until one gets some machine code running on a processor. First, there is the question of what the original algorithm is supposed to compute. Hopefully, the programmer has used the same semantic as the IEEE-754 standard for the operations, the goal being to get portable code and reproducible results. Then the programmer chooses a high-level programming language, since assembly languages would defeat the point of portability. Unfortunately, high-level language semantics are often rather vague with respect to FP operations, so as to account for as many execution environments as possible, even non-IEEE-754-compliant ones. So the programmer has to make some assumptions on how compilers will interpret the program. Unfortunately, compilers might have made different assumptions while still being compliant with the language standard, or they might depart from the standard for the sake of execution speed (possibly controlled by a compilation flag). Finally, the operating system and various libraries play a role too, as they might modify the default behavior of FP units or emulate features not supported in hardware, e.g. subnormal numbers.

A. Java

Let us have an overview of some of the possible semantics through the lens of three major programming languages. Java, being a relatively recent language, started with the most specified description of FP arithmetic. It proposed two data types that match the binary32 and binary64 formats of IEEE-754. Moreover, arithmetic operators are mapped to the corresponding operators from IEEE-754, but rounding modes other than default are not supported, and neither are the override of exceptional behaviors. The latter is hardly ever supported by languages so we will not focus on it in the remaining of this paper.

Unfortunately, a non-negligible part of the architectures the Java language was targeting had only access to x87-like FP units, which allow to set the precision of computation but not the allowed range of exponents. Thus, they behave as if they were working with exotic FP formats that have the usual IEEE-754 precision but an extended exponent range. On such architectures, complying with the Java semantics was therefore highly inefficient. As a consequence, the language later evolved and the FP semantics were relaxed to account for a potential extended exponent range:

Within an expression that is not FP-strict, some leeway is granted for an implementation to use an extended exponent range to represent intermediate results. (15.4 FP-strict expressions, Java SE 7) The Java language specification, however, introduced a strictfp keyword for reinstating the early IEEE-754compliant behavior.

B. C

The C language comes from a time where FP units were more exotic, so the wording of the standard leaves much more liberty to the compiler. Intermediate results can not only be computed with an extended range, they can also have an extended precision.

The values of operations with floating operands [. . . ] are evaluated to a format whose range and precision may be greater than required by the type. (5.2.4.2.2 Characteristics of floating types, C11) In fact, most compilers interpret the standard in an even more relaxed way: values of local variables that are not spilled to memory might preserve their extended range and precision.

Note that this optimization opportunity also applies to the use of a FMA operator for computing the expression a × b + c, as the intermediate product is then performed with a much greater precision.

While Annex F of the C standard allows a compiler to advertise compliance with IEEE-754 FP arithmetic if it supports a specified set of features, none of these features reduces the leeway compilers have in choosing intermediate formats. Moreover, features of Annex F are optional anyway.

C. Fortran

The Fortran language gives even more leeway to compilers, allowing them to rewrite expressions as long as they do not change the value that would be obtained if the computations were to be infinitely-precise.

Two expressions of a numeric type are mathematically equivalent if, for all possible values of their primaries, their mathematical values are equal. (7.1.5.2.4 Evaluation of numeric intrinsic operations, Fortran 2008) The standard, however, forbids such transformations when they would violate the "integrity of parentheses". For instance, (a+b)-a-b can be rewritten as 0, but ((a+b)-a)-b cannot, since it would break the integrity of the outer parentheses.

This allowance for assuming FP operations to be associative and distributive has unfortunately leaked to compilers for other languages, which do not even have the provision about preserving parentheses. For instance, the seemingly innocuous -Ofast option of GCC will enable this optimization for the sake of speed, at the expense of the conformance with the C standard.

D. Stricter Semantics

Fortunately, thanks to the IEEE-754 standard and to hardware makers willing to design strictly-compliant FP units [START_REF] Nickolls | The GPU computing era[END_REF], the situation is improving. It is now possible to specify programming languages without having to keep the FP semantic vague and obscure so that vastly incompatible architectures can be supported. Moreover, even if the original description of a language was purposely unhelpful, compilers can now document precisely how they interpret FP arithmetic for several architectures at once. In fact, in this work, we are going further: not only are we documenting what the expected semantic of our compiler is, but we are formally proving that the compiler follows it for all the architectures it supports.

III. FORMALLY-VERIFIED COMPILATION

As mentioned in Introduction, ordinary compilers sometimes miscompile source programs: starting with a correct source, they can produce executable machine code that crashes or computes the wrong results. Formally-verified compilers such as CompCert C come with a mathematical proof of semantic preservation that rules out all possibilities of miscompilation. Intuitively, the semantic preservation theorem says that the executable code produced by the compiler always executes as prescribed by the semantics of the source program.

Before proving a semantic preservation theorem, we must make its statement mathematically precise. This entails (1) specifying precisely the program transformations (compiler passes) performed by the compiler, and (2) giving mathematical semantics to the source and target languages of the compiler (in the case of CompCert, the CompCert C subset of ISO C90 and ARM/PowerPC/x86 assembly languages, respectively). The semantics used in CompCert associate observable behaviors to every program. Observable behaviors include normal termination, divergence (the program runs forever), and abnormal termination on an undefined behavior (such as an out-of-bounds array access). They also include traces of all input/output operations performed by the program: calls to I/O library functions (such as printf) and accesses to volatile memory locations.

Equipped with these formal semantics, we can state precisely the desired semantic preservation results. Here is one such result that is proved in CompCert:

Theorem 1 (Semantic preservation) Let S be a source C program. Assume that S is free of undefined behaviors. Further assume that the CompCert compiler, invoked on S, does not report a compile-time error, but instead produces executable code E. Then, any observable behavior B of E is one of the possible observable behaviors of S.

The statement of the theorem leaves two important degrees of freedom to the compiler. First, a C program can have several legal behaviors, owing to underspecification in expression evaluation order, and the compiler is allowed to pick any one of them. Second, undefined C behaviors need not be preserved during compilation, as the compiler can optimize them away. This is not the only possible statement of semantic preservation: indeed, CompCert proves additional, stronger statements that imply the theorem above. The bottom line, however, is that the correctness of a compiler can be characterized in a mathematically-precise, yet intuitively understandable way, as soon as the semantics of the source and target languages are specified.

Concerning arithmetic operations in C and in assembly languages, their semantics are specified in terms of two Coq libraries, Int and Float, which provide Coq types for integer and FP values, and Coq functions for the basic arithmetic and logical operations, for conversions between these types, and for comparisons. The CompCert semantics map C language constructs to these basic operations, making fully precise a number of points that the C standards (ISO C 90, 99, and 2011) leave to the discretion of the implementation. For example, the C standards do not specify the precision and range of the float and double types; CompCert C maps them to IEEE-754 binary32 and binary64 numbers, respectively. Likewise, C does not specify the precision of intermediate results during expression evaluation, requiring only that each FP operation is evaluated with a precision greater or equal to that of each operand; CompCert specifies that all intermediate results are computed in double precision. Finally, the C standards allow the compiler to "contract" several FP operations, such as a multiplication and an addition, in a single operation, such as FMA; the CompCert semantics disallow this contraction. 1These choices of implementation are somewhat arbitrary, but they provide programmers with a completely specified, easy-to-understand model of FP arithmetic, which is guaranteed to be implemented faithfully by the compiler. For example, as a consequence of this choice of C semantics and of the semantic preservation theorem, the x86 code generator of CompCert is guaranteed not to generate x87 FP instructions (which operate in extended precision and cannot implement IEEE-754 double precision exactly), generating SSE2 "scalar double" operations instead.

In early versions of CompCert (up to and including 1.11), the formalization of FP arithmetic is, however, less complete and less satisfactory than that of integer arithmetic. The Int library defines machine integers and their operations in a fully constructive manner, as Coq mathematical integers (type Z) modulo 2 32 . In turn, Coq's mathematical integers are defined from first principles, essentially as lists of bits plus a sign. As a consequence of these constructive definitions, all the algebraic identities over machine integers used to justify optimizations and code generation idioms are proved correct in Coq, such as the equivalence between left-shift by n ≥ 0 bits and multiplication by 2 n .

In contrast, in early versions of CompCert, the Float library was not constructed, but only axiomatized: the type of FP numbers is an abstract type, the arithmetic operations are just declared as functions but not realized, and the algebraic identities exploited during code generation are not proved to be true, but only asserted as axioms. (Section V-B shows examples of these identities.) Consequently, conformance to IEEE-754 could not be guaranteed, and the validity of the axioms could not be machine-checked. Moreover, this introduced a regrettable dependency on the host platform (the platform that runs the CompCert compiler), as we now explain.

The Int and Float Coq libraries are used not only to give semantics to the CompCert languages, modeling runtime computations, but also to specify the CompCert passes that perform numerical computations at compile-time. For instance, the constant propagation pass transforms the expression 2.0 * 3.0 into the constant 6.0 obtained by evaluating Float.mul(2.0,3.0) at compile-time. All the verified passes of the CompCert compiler are specified in executable style, as Coq recursive functions, from which an executable compiler is automatically generated by Coq's extraction mechanism, which produces equivalent OCaml code that is then compiled to an executable. For a fully-constructive library such as Int, this process produces an implementation of machine integers that is provably correct and entirely independent from the host platform, and can therefore safely be used during compilation. 2In contrast, for an axiomatized library such as the early versions of Float, there is no other choice than to map FP operations of the library onto those of the host, namely the FP operations provided by OCaml. However, OCaml's FP arithmetic is not guaranteed to implement IEEE-754 double precision: on the x86 architecture running in 32-bit mode, OCaml compiles FP operations to x87 machine instructions, resulting in excess precision and double-rounding issues. Likewise, conversion of decimal FP literals to binary32 or binary64 during lexing and parsing was achieved by calling into the corresponding OCaml library functions, which then call into the strtod and strtof C library functions, which are known to produce incorrectly-rounded results in several C standard libraries.

The discussion above points to a strong need for a fullyconstructive Coq formalization of IEEE-754 arithmetic, providing implementations of FP arithmetic and conversions that are proved correct against the IEEE-754 standard, and can be invoked during compilation to perform constant propagation and other optimizations without being dependent on the host platform. We now describe how we extended the Flocq library to reach these goals.

IV. A BIT-LEVEL COQ FORMALIZATION OF IEEE-754

BINARY FLOATING-POINT ARITHMETIC Flocq (Floats for Coq) is a formalization for the Coq system [START_REF] Boldo | Flocq: A unified library for proving floating-point algorithms in Coq[END_REF]. It provides a comprehensive library of theorems on a multi-radix multi-precision arithmetic. In particular, it encompasses radix-2 and 10 arithmetics, all the standard rounding modes, and it supports fixed-and floating-point arithmetics. The latter comes in two flavors depending on whether underflow is gradual or abrupt. The core of Flocq does not comply with IEEE-754 though, as it only sees FP numbers as subsets of real numbers, that is, it neither distinguishes the sign of zero nor handles special values. We therefore had to extend it to fully support IEEE-754 binary arithmetic. Moreover, this extension had to come with some effective computability so that it could be used in CompCert.

A. Formats and Numbers

Binary FP data with numeric values can be seen as rational numbers m • 2 e , that is, pairs of integers (m, e). This is the generic representation that Flocq manipulates. Support for exceptional values is built upon this representation by using a dependent sum. The above Coq code says that a value of type binary_float can be obtained in four different ways (depending on whether one wants a zero, an infinity, a NaN, or a finite number), and that, for instance, to build a finite number, one has to provide a boolean s, a positive integer m, an integer e, and a proof of the property bounded m e = true.

This property ensures that both m and e are integers that fit into the represented format. This format is described by two variables (precision and exponent range) that are implicit in the above definition. By setting these variables later, one gets specific instances of binary_float, for instance the traditional formats binary32 and binary64. The bounded predicate also checks that m is normalized whenever e is not the minimal exponent. This constraint does not come from the IEEE-754 standard: any element of a FP cohort could be used, but it helps in some proofs to know that this element is unique.

In addition to finite numbers (both normal and subnormal), the binary_float type also supports signed zeros and signed infinities. Notice that there is a single datum NaN though and no way to distinguish between the numerous bit-level encodings the IEEE-754 standard mandates. We chose to abstract NaNs because the IEEE-754 standard underspecifies what happens to their sign and their payload. By ignoring them, we get a data type that encompasses all the compliant architectures. Note that this introduces an inconsistency with the actual architecture: if a program were to create a NaN by an operation between two constants (say 1.0 / 0.0) and test the non-exponent bits of its encoding, CompCert would presumably optimize the program incorrectly.

The function B2R converts a binary_float value to a real number. For finite values, it returns (-1) s ×m×2 e . Otherwise it returns zero. The sign of a value can be obtained by applying the Bsign function.

B. Executable Operations

Once the types are defined, the next step is to implement FP operators and prove their usual properties. An operator takes one or more binary_float inputs and a rounding mode, which tells which FP value to choose when the infinitelyprecise result cannot be represented. The code of these operators always has the same structure. First, they perform a pattern matching on the inputs and handle all the special cases. Only finite numbers are left.

There are two different approaches for defining arithmetic operations. The first one is to have a round function that takes a rounding mode and a real number and return the closest FP number (according to the rounding mode m). For instance, the sum of two finite FP numbers can be defined as a ⊕ b = round(m, B2R(a) + B2R(b)), assuming it does not overflow. The upside is that this operation trivially matches the IEEE-754 standard, since that is the way the standard defines arithmetic operations. The downside is that it depends on an abstract addition and an abstract rounding function, and thus it does not carry any computable content. As such, it cannot be used in a compiler that needs to perform FP operations to propagate constant values. This approach is used in the Pff [START_REF] Boldo | Preuves formelles en arithmétiques à virgule flottante[END_REF] library and in the Flocq core library [START_REF] Boldo | Flocq: A unified library for proving floating-point algorithms in Coq[END_REF].

The second approach is to define arithmetic operators that actually perform computations on integers to construct a FP result. This time, the code of these operators can be used by a compiler for emulating FP operations, which is what we want. The downside is that, not only are these functions complicated to write, but there is no longer any guarantee that they are compliant with the IEEE-754 standard. So one also has to formally prove such theorems. This approach is used in the FP formalization for ACL2 [START_REF] Russinoff | A mechanically checked proof of IEEE compliance of the floating point multiplication, division and square root algorithms of the AMD-K7 processor[END_REF].

As done in HOL Light [START_REF] Harrison | A machine-checked theory of floating point arithmetic[END_REF], [START_REF] Harrison | Formal verification of floating point trigonometric functions[END_REF], we have mixed both approaches for our purpose: the second one offers effective computability, while stating and proving that the first one is equivalent provides concise specifications for our operations. Currently supported operations are opposite, addition, subtraction, multiplication, division, and square root. Since other operations like FMA, remainder, or square root are standard library functions, they are not needed in our compiler formalization, as there are no specific inlining optimizations for them. As an example of our approach, here is the correctness theorem for the FP multiplication Bmult.

Theorem 2 (Bmult correct) Given x and y two binary_float numbers, a rounding mode m, and denoting round(m, B2R(x) × B2R(y)) by z, we have

   B2R(Bmult(m, x, y)) = z if |z| < 2 E , Bmult(m, x, y) = overflow(m, Bsign(x) × Bsign(y)) otherwise.
Note that Flocq's round function returns a real number that would be representable by a FP number if the format had no upper bound on the exponents. In particular, if the product overflows, then z is a number larger than the largest representable FP number (1-2 -p )•2 E . In that case, the overflow function is used to select the proper result depending on the rounding mode (either an infinity or the largest representable number) according to the IEEE-754 standard.

Notice that the theorem works even for exceptional inputs since B2R maps them to zero. This makes it a bit simpler to apply. This simplicity is even more sensible for the square root, since it is proved to never overflow. Moreover, Coq's square root is a total function that returns zero for negative inputs, while the FP operator Bsqrt returns NaN. So the input x does not even have to be nonnegative for the theorem to hold.

Theorem 3 (Bsqrt correct) Given x a binary_float number and m a rounding mode, we have

B2R(Bsqrt(m, x)) = round m, B2R(x) .
These correctness theorems specify fully only the case when both inputs of the operators are finite numbers. Indeed this is the difficult case. When one or both inputs are exceptional, no theorems are needed since one can simply execute the operators to recover their values.

C. Bit-Level Representation

Finally, the last part needed to build a compiler is the ability to go from and to the representation of FP numbers as integer words. We provide two functions for this purpose and a few theorems about them. Again, it is important that these functions are effectively computable.

The binary_float_of_bits function takes an integer, splits it into the three parts of a FP datum, looks whether the biased exponent is minimal (meaning the number is zero or subnormal), maximal (meaning infinity or NaN), or in between (meaning a normal number with an implicit bit), and constructs the resulting FP number of type binary_float. The bits_of_binary_float performs the converse operation. Note that it always returns the same NaN (all bits set to 1, except for the sign bit).

Both functions have been proved to be inverse of each other (except for NaNs) for bounded integers. This property also guarantees that we did not get these conversion functions too wrong. Indeed, it ensures that all the bits of the memory representation are accounted for and that there is no overlap between the three fields of the binary representation.

V. A VERIFIED COMPILER FOR FLOATING-POINT COMPUTATIONS

We integrated the Coq formalization of IEEE-754 arithmetic described in Section IV into the CompCert compiler, version 1.12, effectively replacing the axiomatization of FP arithmetic used in earlier versions (see Section III) by a provably-correct, executable implementation.

As a first benefit, we obtain more precise semantic specifications for the source and target languages of CompCert. The semantics for the source CompCert C language now guarantee that FP arithmetic is performed as prescribed by IEEE-754, a guarantee that programmers can rely on. Symmetrically, the semantics for the target assembly languages (ARM, PowerPC, x86) now require that the hardware implements IEEE-754 correctly. Two of CompCert's target architectures have several FP instruction sets, with different characteristics. Our semantics only model the instructions actually generated by CompCert: for ARM, the scalar VFD instruction set, omitting vector instructions; for x86, the scalar SSE2 instruction set, leaving aside vector instructions and x87 extended-precision instructions.

As another benefit of building on a Coq formalization of IEEE-754 arithmetic, we can now prove, as Coq theorems, the axioms about the float abstract type previously used by CompCert. As we explain in the following, these theorems prove the correctness of CompCert's compile-time handling of FP arithmetic: first, FP computations performed at compiletime by the compiler (such as FP literal parsing or constant propagation); second, the code generation strategies used to implement C's FP operations in terms of the instructions provided by the target architectures.

A. Verifying Compile-Time Computations

The CompCert compiler performs FP computations at different stages of compilation: (1) parsing of FP literals, (2) the constant propagation optimization, and (3) conversion of FP numbers to their bit-level representation when generating the final executable code. For conducting these operations, we need an implementation of FP arithmetic that is proved correct in Coq, executable via extraction from Coq to OCaml, and reasonably efficient. As shown in Section IV, our extension to the Flocq library provides such an implementation. In particular, the bits_of_binary_float function described in Section IV-C directly answers usage (3) above. We now discuss the use of Flocq for purposes (1) and [START_REF] Dekker | A floating point technique for extending the available precision[END_REF].

Constant propagation is a basic but important optimization in compilers. It consists in evaluating, at compile-time, arithmetic and logical operations whose arguments can be statically determined. For instance, the Fast-Two-Sum example of the introduction is reduced to the printing of a single constant; no FP operations are performed by the executable code. For another example, consider the following C code fragment:

inline double f(double x) { if (x < 1.0) return 1.0; else return 1.0 / x; } double g(void) { return f(3.0); }

Combining

constant propagation with function inlining, the body of function g is optimized into return 0x1.5555555555555p-2. Not only the division 1.0 / x but also the conditional statement x < 1.0 have been evaluated at compile-time. These evaluations are performed by the executable operations provided by the Flocq library, making them independent from the FP arithmetic of the host platform running the compiler, and guaranteeing that the constant propagation optimization preserves the semantics of the source program. 3The evaluation of FP literals is delicate: literals are often written in decimal, requiring nontrivial conversion to IEEE-754 binary format; moreover, correct rounding must be guaranteed [START_REF] Clinger | How to read floating-point numbers accurately[END_REF]. For example, until recently, the strtod and strtof functions of the GNU C standard library incorrectly rounded the result in some corner cases. 4 To avoid these pitfalls, we use a simple but correct Flocq-based algorithm for evaluating these literals.

In C, a FP literal consists of an integral part, a fractional part, an exponent part, and a precision suffix (which indicates at which precision the literal should be evaluated). Each of these parts can be omitted, in which case 0 is used as default value for the first three parts. (This operation is done in an early stage of parsing in our compiler.) The integral and fractional parts may be written in either decimal or hexadecimal notation; the use of hexadecimal (in both parts) is indicated if the integral part begins with the prefix "0x". The exponent is given as a power of 2 if hexadecimal is used or as a power of 10 if decimal is used. To summarize, a literal number always has the form I.F × b E with b = 2 or b = 10.

The first part of our algorithm consists in shifting the point to the right, while modifying the exponent in order to transfer the fractional part F into the integral part I. Then, it parses both the exponent and the new integral part as arbitraryprecision integers. The last part consists in actually evaluating the FP number, using Flocq with the precision specified by the precision suffix. When E ≥ 0, we compute I × b E using exact integer arithmetic, then round the result to the nearest representable FP number. When E < 0, we first compute b -E using exact integer arithmetic, then perform the FP division •(I/b -E ), using the proved division of Flocq. Notice that, since Flocq formalizes a multi-precision arithmetic, numbers I and b -E do not have to fit into the target format; the division can cope with arbitrarily large numbers.

It is clear that the result is evaluated as in the reals before being rounded at the very last step. We believe this implementation is one of the simplest one could give, and we would use it as a specification to a more complicated algorithm if better performance is needed.

B. Verifying Code Generation for Floating-Point Operations

Most FP operations of the C language map directly to hardware-implemented instructions of the target platforms. However, some operations, such as certain comparisons and conversions between integers and FP numbers, are not directly supported by some target platforms, forcing the compiler to implement these operations by sometimes convoluted combinations of other instructions. The correctness of these code generation strategies depends on the validity of algebraic identities over FP operations, identities that we were able to verify in Coq using the theorems provided by Flocq.

A first example is FP comparisons on the PowerPC architecture. The PowerPC provides an fcmp instruction that produces 4 bits of output: "less than", "equal", "greater", and "uncomparable", and conditional branch instructions that test any one of these bits. To compile a large inequality test such as "less than or equal", CompCert produces code that performs the logical "or" of the "less than" and "equal" bits, then conditionally branches on the resulting bit. Semantically, this is justified by the identity (x ≤ y) ≡ (x < y) ∨ (x = y), which holds for any two FP numbers x and y. Note that, even if two NaNs are equal from the mathematical point of view of Coq, the comparison operators defined by the compiler still know that NaNs shall be unordered [START_REF]IEEE Standard for Floating-Point Arithmetic[END_REF].

Another example is conversions between integers and FP numbers, which come in 4 variants, depending on the direction of the conversion (from an integer or to an integer) and on the type of the integer: either signed, in the range [-2 31 , 2 31 ), or unsigned, in the range [0, 2 32 ). Of the 3 CompCert target platforms, only ARM provides hardware implementations of all 4 conversions. The x86-SSE2 instruction set only provides conversions to and from signed integers, requiring the unsigned integer conversions to be implemented by case analysis: if n is an unsigned 32-bit integer variable, the C conversion (double)n is compiled like the C conditional expression n < 0x80000000 ? (double)((int) n) : (double)((int)(n -0x80000000)) + 0x1.p31

Likewise, if d is a binary64 variable, the C conversion (unsigned int)d is compiled like

d < 0x1.p31 ? (int) d : (int)(d -0x1.p31) + 0x8000000
We proved the correctness of this code generation strategy, using the fact that all 32-bit signed and unsigned integers are exactly representable as binary64, and that conversions from FP numbers to integers are undefined if the argument falls outside the range of the destination integer type. The PowerPC 32-bit architecture is even more problematic, as the only conversion it implements in hardware is binary64 to signed integer. Conversion to an unsigned integer is implemented as shown above. Conversions from integers to FP numbers are synthesized via bit-level manipulations over the binary64 numbers, as suggested by IBM [START_REF] Ibm | The PowerPC Compiler Writer's Guide[END_REF]. If n is an unsigned 32-bit integer, the conversion (double)n is compiled as fmake(0x43000000, n) -fmake(0x43000000, 0)

where fmake(hi,lo) is a compiler built-in function that returns a double whose 64-bit binary representation is the concatenation of the 32 bits of integer hi followed by the 32 bits of integer lo. If, instead, n is a signed integer, (double)n is compiled as fmake(0x43000000, n ˆ0x80000000) -fmake(0x43000000, 0x80000000)

The correctness of this implementation technique is far from obvious. Taking the unsigned case as example, we first note that fmake(0x43000000, n) is equal to 2 52 + n, and that fmake(0x43000000, 0) is 2 52 . We then prove that the binary64 subtraction between these two numbers is exact and produces the binary64 number equal to n. Mechanizing this proof in Coq brings much confidence in this implementation. The discussion above focuses on conversions between 32bit integers and binary64 FP numbers, as these are the only conversions between integers and FP numbers used by Comp-Cert's implementation and formal semantics. Conversions to and from binary32 FP numbers go through an intermediate binary64 result, with no risk of double rounding errors. Support for 64-bit integers is being considered for future versions of CompCert. Synthesizing conversions between 64-bit integers and binary64 FP numbers from the 32-bit conversion instructions provided by the target processors raises interesting issues that deserve further formal verification.

VI. CONCLUSIONS

In this article, we have presented a formally-verified compiler that supports FP computations. Producing such a compiler required us to define the FP semantics for the C language and for the target architectures, and to prove that the compiler preserves the semantics between a C program and the produced executable code. Flocq has been extended with a formalization of the IEEE-754 standard; this formalization is used by CompCert to define the semantics, parse literal FP constants, and perform constant propagation at compiletime. This development has been integrated into version 1.12 of CompCert available at http://compcert.inria.fr/. This work required to add about 3000 lines of new Coq proofs to both CompCert and Flocq.

This approach gives a correct and predictable compiler that conforms to the IEEE-754 standard. This means that, among the several possibilities allowed by the ISO C standard, we have chosen a single way to compile and we have formally proved its correctness. This compilation choice can be discussed: for example all intermediate results are computed in double precision, therefore with (usually) less accuracy than with extended registers. The first reason is that this is sorely needed to be able to prove algorithms or programs. The second reason is that we favored reproducibility over possible higher accuracy. The actual interpretation of FP operations can be seen in the Float module of CompCert; one does not have to wade through all the optimization passes to understand what happens to them, since their semantics is provably preserved. Another advantage is that having strict semantics paves the way to simpler, more precise, and even verified, static analyzers.

For the sake of completeness, one should note that Comp-Cert's formal semantics does not support certain features of the IEEE-754 standard. First, CompCert does not know about directed rounding modes and assumes that all the FP operations are performed with the default rounding mode. As a consequence, on architectures that have dynamic rounding modes, changing the mode prevents CompCert's semantics from being preserved. For instance, constant propagation might give a different result from actual execution. CompCert could be extended to support a dynamic mode, e.g. by representing it as a pseudo global variable. Constant propagation would then only happen if either the rounding mode is statically known, or if the result would be the same whatever the mode.

Another peculiarity of CompCert is that all the intermediate computations are performed in double precision, as allowed by the C standard. It is still possible to achieve binary32 computations by following each operation by a store to a binary32 variable. Double rounding occurs but is known to produce the correctly-rounded result [START_REF] Figueroa | When is double rounding innocuous?[END_REF].

The integration of Flocq and CompCert opens the way to adding more optimizations specific to FP arithmetic, and to prove them correct. FP identities such as x -0.0 ≡ x or x × 2 ≡ x + x or x/2 n ≡ x × 2 -n can be exploited to generate shorter or cheaper instruction sequences. For the semantic preservation theorem to hold, however, only identities that hold for all representable FP numbers can be used, and there are few of them. Other tempting simplifications are incorrect for some values of their arguments: for example, x + 0.0 ≡ x does not hold if x = -0.0, and Brisebarre et al's technique to replace a FP division by a constant with a multiplication and an FMA [START_REF] Brisebarre | Accelerating correctly rounded floating-point division when the divisor is known in advance[END_REF] is not always correct for subnormal arguments. The only way to exploit these simplifications while preserving semantics would be to apply them conditionally, based on the results of a static analysis (such as FP interval analysis) that can exclude the problematic cases.

The problem is even more acute for aggressive loop optimizations such as vectorization, which often entail reassociating FP operations, and cannot be guaranteed to preserve semantics except in very special cases. We conclude that the compiler is probably the wrong place to perform aggressive program transformations over FP operations, because it lacks much of the information necessary for this endeavor. Automatic code generation tools, however, are in a more favorable position to preserve or improve precision by reassociation and other aggressive transformations [START_REF] Ioualalen | A new abstract domain for the representation of mathematically equivalent expressions[END_REF].

Inductive

  binary_float := | B754_zero : bool -> binary_float | B754_infinity : bool -> binary_float | B754_nan : binary_float | B754_finite : forall (s : bool) (m : positive) (e : Z), bounded m e = true -> binary_float.

This work was supported by the VERASCO project (ANR-11-INSE-003) of the French National Agency for Research (ANR).

On target platforms that support them, CompCert makes FMA instructions available as compiler built-in functions, but they must be explicitly used by the programmer.

This is similar in spirit to GCC's use of exact, GMP-based integer arithmetic during compilation, to avoid dependencies on the integer types of its host platform.

The CompCert C semantics gives programmers no way to change the FP rounding mode during program execution, therefore guaranteeing that all FP arithmetic rounds to nearest even. Programs that need other rounding modes fall outside the perimeter of CompCert's semantic preservation results. They can, however, be supported via a compiler option, -ffloat-const-prop 0, which turns FP constant propagation off.

Bug 3479 -Incorrect rounding in strtod(), http://sourceware.org/bugzilla/ show bug.cgi?id=3479


T Goudon

P Lafitte

The lovebirds problem: why solve Hamilton-Jacobi-Bellman equations matters in love affairs

The lovebirds problem consists in finding the compromise between well-being and efforts that are necessary to sustain a sentimental relationship. According to a modeling introduced by J.-M. Rey, the problem can be described as finding the initial data for a certain dynamical system, guaranteeing that the associated trajectory belongs to the stable manifold. We further discuss this interpretation by means of the Dynamical Programming Principle and the Hamilton-Jacobi-Bellman framework. It allows us to propose an algorithm that determines numerically the solution of the lovebirds problem.

Introduction

In a recent paper [START_REF] Rey | A mathematical model of sentimental dynamics accounting for marital dissolution[END_REF], J.-M. Rey introduces a convincing mathematical description of the "sentimental dynamics". The couple is seen as an entity, and the model is based on a contest between the quest of a common well-being, embodied into a quantity t → x(t), depending only on the time variable, and the cost t → c(t) of the efforts necessary to maintain a satisfactory well-being. Indeed, the well-being has a natural tendency to fade, and efforts can counter-balance the erosion and sustain the well-being. The model is completed by a utility structure which accounts for the balance between the valuation of the well-being and sacrifices induced by the efforts. Too small values of the well-being as well as too high values of the efforts are not tolerable, which leads to the dissolution of the relationship. Despite its simplicity, the model proposed in [START_REF] Rey | A mathematical model of sentimental dynamics accounting for marital dissolution[END_REF] provides interesting information on the dynamics of sentimental partnerships. In particular it brings an explanation to the so-called "failure paradox": while partners consistently plan a firm and eternal sentimental relationship, a large proportion of unions are eventually going to break up. As we shall recall below, according to J.-M. Rey's modeling, the sentimental dynamics can be mathematically understood in terms of dynamical systems and analysis of the phase portrait t → (x(t), c(t)). It turns out that the underlying dynamical system has a single equilibrium point, denoted hereafter (x, c), and the corresponding linearized system admits eigenvalues of opposite signs. The instability of the equilibrium state explains the difficulty in maintaining a durable relationship. Overall, the partners have to determine their efforts policy so that the trajectory t → (x(t), c(t)) remains on the stable manifold of the dynamical system. In this paper we are thus interested in the "lovebirds problem": given an initial well-being state x 0 , which is hopefully supposed to be quite high, how can we find the corresponding effort rate c 0 which guarantees an unbreakable relationship, with a loving trajectory staying on the stable manifold? We shall see that the problem can be interpreted by means of the Dynamic Programming Principle that leads to the Hamilton-Jacobi-Bellman formalism. We describe a numerical procedure based on this interpretation which allows to find a fair approximation of the searched initial state.

The paper is organized as follows. At first we review the basis of J.-M. Rey's modeling. In particular we make the dynamical system which governs the sentimental dynamics appear. In Section 3 we detail several questions of mathematical analysis on the model to make the role of the stable manifold clear. Finally, we discuss the problem in the framework of Dynamic Programming, which allows us to design a numerical method that gives an approximate solution to the lovebirds problem.

The dynamics of sentimental relationships

Let us remind the basis of J.-M. Rey's model. Readers interested in further mathematical modeling of sentimental dynamics, based on different ideas, can consult [START_REF] Gottman | The mathematics of marriage: dynamic nonlinear models[END_REF] or [START_REF] Murray | Mathematical Biology: I. An Introduction[END_REF]Chapter 5]. Also an interesting model of parental care relying on game theory and with common technical features with the problem under consideration here can be found in [START_REF] Ewald | Parental care as a differential game: A dynamic extension of the houston-davies game[END_REF]. The evolution of the well-being is governed by the ODE

d dt x = -rx + c, x(0) = x 0 . ( 1
)
Given x 0 and t → c(t), we denote by X x 0 ,c the corresponding solution of (1):

X x 0 ,c : t -→ X x 0 ,c (t) = x 0 e -rt + t 0 e -r(t-s) c(s) ds.
The coefficient r > 0 contains the natural erosion of the sentimental feeling. The value of the erosion parameter r depends on the considered relationship, as well as the threshold x min below which the relationship becomes unsatisfactory and cannot survive. Next, we introduce the utility structure. The quest for happiness then results in a competition between feeling and efforts, which can be expressed as the optimization of the following cost function, the so-called "total satisfaction"

J : (x, c) -→ ∞ 0 e -ρt (U(x(t)) -D(c(t))) dt,
the functional framework being clarified later on. This quantity involves the so-called impatience parameter ρ > 0, which enters in the definition of the weight in the history of the relationship evaluated with J(x, c). We refer to [START_REF] Rey | A mathematical model of sentimental dynamics accounting for marital dissolution[END_REF] for further details on the modeling issues. We seek to maximize J under the constraint [START_REF] Barles | Solutions de viscosité des équations de Hamilton-Jacobi[END_REF]. Therefore we wish to define the mapping

J : x 0 -→ max c {J (x 0 , c)} where J : (x 0 , c) -→ ∞ 0 e -ρt (U(X x 0 ,c (t)) -D(c(t)) dt.
As a matter of fact, we deduce from i)-v) that D is non-negative, convex and satisfies 0 ≤ D(c

) + α|c -c | 2 ≤ D(c) ≤ D(c ) + A|c -c | 2
while U is concave and satisfies, for any x ≥ 0

0 ≤ U(x) ≤ U(0) + U (0)x.
This remark allows to make the functional framework precise: given x 0 > 0, J (x 0 , •) makes sense on the set

L 2 ρ,+ = c ≥ 0, ∞ 0 e -ρt |c(t)| 2 dt < ∞ .
Indeed, it is clear that, for any x 0 > 0 and c ∈ L 2 ρ,+ , we have X x 0 ,c ∈ L 1 ρ,+ , and

X x 0 ,c ∈ L 2 κ,+ for any κ > ρ. Therefore J (x 0 , •) is α-concave on L 2 ρ,+
, which is a weakly closed convex subset of a Banach space. Computing the derivative with respect to c (denoted here with the prime symbol) we get

J (x 0 , c)(h) = ∞ 0 e -ρt U (X x 0 ,c (t))e -rt t 0 e rs h(s) ds -D (c(t))h(t) dt.
Integrating by parts yields

J (x 0 , c)(h) = ∞ 0 h(t) e rt ∞ t U (X x 0 ,c (s))e -(r+ρ)s ds -e -ρt D (c(t)) dt.
We conclude with the following statement. Proposition 2.1 The functional J (x 0 , •) admits a unique maximizer t → c ♥ (t) in L 2 ρ,+ , characterized by the property:

For any h ∈ L 2 ρ,+ , we have J (x 0 , c ♥ )(c ♥ -h) ≥ 0.
Assuming that c ♥ remains positive, we can choose the trial function

h = c ♥ ± ζ, with ζ ≥ 0 in C ∞ c ((0, ∞))
and we end up with the following relation, which holds for any t ≥ 0, c ♥ (t) = (D ) -1 e (ρ+r)t ∞ t e -(r+ρ)s U x 0 e -rs + s 0 e -r(s-σ) c ♥ (σ)dσ ds .

(

) 2
In turns out that the optimal solution can be investigated by means of dynamical systems, as a consequence of the following facts (which are already presented in [START_REF] Rey | A mathematical model of sentimental dynamics accounting for marital dissolution[END_REF]). c) The pair

(x ♥ = X x 0 ,c ♥ , c ♥ ) is a solution of the differential system d dt x = -rx + c, D (c) d dt c = -U (x) + (r + ρ)D (c). (3)
d) Furthermore, we have lim t→∞ (x ♥ (t), c ♥ (t)) = (x, c), where (x, c) is the (unique) equilibrium of system (3); in other words (x ♥ , c ♥ ) belongs to the stable manifold of (3).

For the forthcoming discussion, it can be helpful to have in mind the typical phase portrait of the dynamical system (3), as displayed in Figure 1. Before detailing the analysis, it is worth explaining the question by means of the Lagrangian framework. We denote by L the following Lagrangian functional

L : (x, c; λ) -→ J(x, c) - ∞ 0 λ(t)e -ρt d dt x + rx -c (t) dt.
For the functional framework, L(x, c; λ) is well defined

with c ∈ L 2 ρ , x ∈ L 2 κ with d dt x ∈ L 2 κ , where ρ < κ < 2ρ so that L 2 κ ⊂ L 1 ρ , and λ ∈ L 2 2ρ-κ . Theorem 2.3 If (x ♥ , c ♥ , λ ♥ ) is a saddle point of L then we have J(x ♥ , c ♥ ) = max J(x, c), d dt x = -rx + c, x(0) = x 0 = max c J (x 0 , c) = J(x 0 ). ( 4
)
Figure 1: Typical phase portrait of the sentimental dynamics.

Furthermore, as long as c ♥ (t) > 0, (x ♥ , c ♥ , λ ♥ ) is a solution of the following differential system

d dt x ♥ (t) = -rx ♥ (t) + c ♥ (t), (5)
d dt c ♥ (t) = (r + ρ)D (c ♥ (t)) -U (x ♥ (t)) D (c ♥ (t)) , ( 6
) d dt λ ♥ (t) = (r + ρ)λ ♥ (t) -U (x ♥ (t)). (7)
Proof. For any trial function t → µ(t) we observe that

L(x ♥ , c ♥ ; λ ♥ + µ) -L(x ♥ , c ♥ ; λ ♥ ) = - ∞ 0 e -ρt µ(t) d dt x ♥ (t) + rx ♥ (t) -c ♥ (t) dt ≥ 0
holds, which allows to justify [START_REF] Ewald | Parental care as a differential game: A dynamic extension of the houston-davies game[END_REF]. Then, we deduce [START_REF] Evans | Partial Differential Equations[END_REF]. By the same token, we study the behavior of L with respect to perturbations of x ♥ and c ♥ . Owing to the convexity assumptions on U and D we show that for any t → ξ(t) and t → κ(t) such that c ♥ (t) + κ(t) ≥ 0, the quantity

∞ 0 e -ρt U (x ♥ (t))ξ(t) -D (c ♥ (t))κ(t) -λ ♥ (t) d dt ξ(t) + rξ(t) -κ(t) dt = ∞ 0 e -ρt ξ(t) U (x ♥ (t)) -(ρ + r)λ ♥ (t) + d dt λ ♥ (t) -κ(t) D (c ♥ (t)) -λ ♥ (t) dt
is non positive. It implies the relations ( 6) and ( 7).

Theorem 2.2 establishes the connection between the optimization problem and the determination of the stable manifold of the associated differential system: given x 0 ≥ 0, finding the optimal strategy c ♥ reduces to finding c 0 ≥ 0 such that (x 0 , c 0 ) lies in the stable manifold of (3). However, this observation, which has a key role in the discussion in [START_REF] Rey | A mathematical model of sentimental dynamics accounting for marital dissolution[END_REF], left open how the optimal state c 0 could be found. Formula (2) defines c 0 and the trajectory c ♥ implicitly, but the formula is of reduced practical use: how to derive a numerical method based on (2) to find the optimal strategy is far from clear. We shall see that c 0 can be obtained as a function of x 0 by means of the resolution of a HJB equation. This viewpoint which involves now the resolution of a PDE might look complicated. However, it turns out that efficient numerical procedures can be used to solve this PDE, thus determining the stable manifold of (3).

Proof of Theorem 2.2

We start by discussing the existence-uniqueness of the equilibrium point of (3). For x ≥ 0 given, since D and U are monotone, their asymptotic behavior being prescribed in Assumptions ii) and iv) above, there exists a unique γ(x)

∈ [c , +∞) such that (r +ρ)D (γ(x)) = U (x). We note that γ (x) = U (x) (r+ρ)D (γ(x)) < 0 and lim x→∞ γ(x) = c . Hence, in the phase plane (x, c), the curve x → γ(x) intersects the line c = rx at a unique point (x, c) ∈ (0, ∞) × (c , ∞).
It defines the unique equilibrium point of (3). As observed in [START_REF] Rey | A mathematical model of sentimental dynamics accounting for marital dissolution[END_REF], there is an effort gap: the rate of effort c at the equilibrium is larger than the rate of least effort c . Note also that the equilibrium values x and c = rx are non increasing functions of ρ, because

r + ρ = U (x) D (rx)
where the right-hand side is a non increasing function of x, which means that the longer the memory lasts, the better the feeling.

The goal is to show the existence of a curve c 0 = ϕ(x 0 ) such that the points (x 0 , ϕ(x 0 )) are zeroes of the mapping 3) associated to the initial data (x 0 , c 0 ). The mapping under consideration corresponds to the evaluation of (2) at time t = 0. To start with, we observe that (x, c) is such a zero:

Φ(x 0 , c 0 ) = D (c 0 ) - ∞ 0 e -(r+ρ)s U (x(x 0 , c 0 ; s)) ds with t → (x, c)(x 0 , c 0 ; t) solution of (
Φ(x, c) = D (c) - U (x) r + ρ = 0.
We are going to consider Φ over certain perturbations of this specific solution.

To this end, we need to introduce a couple of definitions. The linearization of (3)

yields d dt y = Ay with A =   -r 1 - U (x) D (c) r + ρ   .
The characteristic polynomial reads

P A (λ) = λ 2 -ρλ + U (x) D (c) -r(r + ρ) = (λ -λ + )(λ -λ -).
The discriminant ∆ = (2r + ρ) 2 -4 U (x) D (c) is positive. Hence the system admits two eigenvalues of opposite signs

λ ± = 1 2 (ρ ± √ ∆). ( 8
)
Some associated eigenvectors are

V ± = 1 N ± 1 r + λ ± , N ± = 1 + (r + λ ± ) 2 .
For further purposes we remark that

r + λ -= r + ρ -λ + = 2r + ρ 2 1 -1 - 4U (x) (2r + ρ) 2 D (c) < 0, λ + -λ -= √ ∆ > 0, (9)
holds. Let P = (V + V -) be the passage matrix from the canonical basis (e 1 , e 2 ) to the eigenbasis (V + , V -). For the inverse matrix we find

P -1 = N -N + λ --λ +     r + λ - N - - 1 N - - r + λ + N + 1 N +     .
We introduce the projection Π -on the negative eigenspace Span{V -} following the direction V + , namely in the canonical basis

Π -= P 0 0 0 1 P -1 .
We set

(x, c) = x c + Π -x c and Ψ(x, c) = D ( (x, c) • e 2 ) - ∞ 0 e -(r+ρ)s U (x( (x, c); s)) ds = Φ( (x, c)).
Of course, we have Ψ(0, 0) = 0. We rephrase the problem of finding zeroes (x 0 , c 0 ) of Φ as searching for perturbations ỹ = (x, c) that make Ψ vanish. To this end, we are going to use the Implicit Function Theorem.

We compute the derivative

∇ ỹΨ(0) T = D (c)(0 1)Π -- ∞ 0 e -(ρ+r)s U (x) X(s) ds
where s → X(s) is given by the first row of the matrix Y solution of the linear system

d dt Y = AY, Y (0) = Π -.
In other words, since A = P diag(λ + , λ -) P -1 , we have

Y (t) = P e λ + t 0 0 e λ -t 0 0 0 1 P -1 = P 0 0 0 e λ -t P -1 .
It follows that

∇ ỹΨ(0) T = D (c)(0 1) - ∞ 0 e -(ρ+r)s e λ -s U (x)(1 0) ds P 0 0 0 1 P -1 = D (c)(0 1) - U (x) ρ + r -λ - (1 0) ds P 0 0 0 1 P -1 = ∂ xΨ(0) ∂ cΨ(0) .
We compute

P 0 0 0 1 P -1 =     r + λ + λ + -λ - N + N -(λ + -λ -) - (r + λ + )(r + λ -) λ + -λ - - r + λ - λ + -λ -     .
Then, coming back to (9), we obtain

∂ cΨ(0) = -(r + λ -) λ + -λ - D (c) + -U (x) r + ρ -λ - N + N -(λ + -λ -) > 0.
We can thus apply the Implicit Function Theorem and we are led to the following statement.

Lemma 3.1 There exists a neighborhood V of (0, 0) and an application x → ψ(x)

such that (x, ψ(x)) ∈ V and Ψ(x, ψ(x)) = 0 = Φ( (x, ψ(x))). The set {(x, ψ(x))} characterizes the zeroes of Ψ in V.
It remains to establish clearly the connection with the optimization problem and claims in Theorem 2.2. This is the purpose of the following sequence of lemmas.

Lemma 3.2 Assume that c ♥ satisfies (2). Then, (x ♥ = X x 0 ,c ♥ , c ♥ ) fulfills (3), c ♥ (t) > 0 for any t ≥ 0 and t → c ♥ (t) is uniformly bounded. Proof. By definition x ♥ satisfies d dt x ♥ = -rx ♥ +c ♥ . Next, we differentiate (2) written as D (c ♥ (t)) = e (r+ρ)t ∞ t U (x ♥ (s))e (r+ρ)s ds.
Finally, (2) implies that D (c ♥ (t)) > 0 = D (c ). Since D is increasing, we conclude that c ♥ (t) > c ≥ 0. We also observe that

0 ≤ e (r+ρ)t ∞ t U (x ♥ (s))e (r+ρ)s ds ≤ U (0) r + ρ , since U is non increasing. Thus, (2) leads to c ♥ (t) ≤ D U (0) r+ρ . Lemma 3.3 Assume that (x, c) is a solution of (3) associated to an initial data (x 0 , c 0 ) satisfying Φ(x 0 , c 0 ) = 0. Then, t → c(t) fulfills (2).
Proof. We have x = X x 0 ,c and

d dt e -(ρ+r)t D (c(t)) = e -(ρ+r)t U (X x 0 ,c (t)).
Integrating this relation and using Φ(x 0 , c 0 ) = 0 yield

D (c(t)) = e (ρ+r)t D (c 0 ) - t 0 e -(ρ+r)s U (X x 0 ,c (s)) ds = e (ρ+r)t ∞ 0 e -(ρ+r)s U (X x 0 ,c (s)) ds - t 0 e -(ρ+r)s U (X x 0 ,c (s)) ds = e (ρ+r)t ∞ t e -(ρ+r)s U (X x 0 ,c (s)) ds which is (2).
Lemma 3. [START_REF] Evans | Partial Differential Equations[END_REF] The set {(x, c) ∈ (0, ∞) × (0, ∞), Φ(x, c) = 0} is invariant under the flow associated to [START_REF] Michael | Viscosity solutions of Hamilton-Jacobi equations[END_REF].

Proof. This is a consequence of the fact that (3) is autonomous. Let us denote (x(x 0 , c 0 ; t), c(x 0 , c 0 ; t)) the evaluation at time t ≥ 0 of the solution of (3) issued from (x 0 , c 0 ) at time 0. Therefore we have

x(x 0 , c 0 ; t+s), c(x 0 , c 0 ; t+s) = x(x(x 0 , c 0 ; s), c(x 0 , c 0 ; s); t), c(x(x 0 , c 0 ; s), c(x 0 , c 0 ; s); t) .

By Lemma 3.3, when Φ(x 0 , c 0 ) = 0, we can write

c(x 0 , c 0 ; s) = D -1 e (ρ+r)s ∞ s e -(ρ+r)σ U (x(x 0 , c 0 ; σ)) dσ = D -1 ∞ 0 e -(ρ+r)τ U (x(x 0 , c 0 ; s + τ )) dτ = D -1 ∞ 0 e -(ρ+r)τ U (x(x(x 0 , c 0 ; s), c(x 0 , c 0 ; s); τ )) dτ .
We deduce that Φ x(x 0 , c 0 ; s), c(x 0 , c 0 ; s) = 0.

Lemma 3.5 Assume that (x, c) is a solution of (3) associated to an initial data (x 0 , c 0 ) satisfying Φ(x 0 , c 0 ) = 0. Then, for any t ≥ 0, c(t) lies in the interval defined by c and the nullcline {(r + ρ)D (c) = U (x)}.

Proof. The result follows from the analysis of the phase portrait of (3). We justify the result for c 0 , the general case being a consequence of Lemma 3.4. Assume that x 0 > x.

If c 0 > c, then lim t→∞ c(t) = ∞; if c 0 is below the nullcline, then t → c(t) vanishes in finite time. In both case it contradicts Lemma 3.2. Similar reasoning applies when 0 < x 0 < x. In particular, it means that if either the initial well-being or the initial effort is too low a lasting relationship is hopeless.

By Lemma 3.1, we have constructed a solution of (3), associated to an initial data which satisfies Φ(x 0 , c 0 ) = 0. Then, by Lemma 3.3, c(t) satisfies (2) for any t ≥ 0 and we have c(t) > 0. Therefore the obtained c realizes the optimum of the functional J (x 0 , •). By Lemma 3.5 we deduce that t → c(t) and t → x(t) are monotone (for instance t → c(t) is increasing and t → x(t) is decreasing when x 0 > x). Accordingly x(t) and c(t) admit limits, denoted x ∞ and c ∞ respectively, as time goes to ∞. Now letting t → ∞ in (2) we are led to

c(t) -→ t→∞ c ∞ = (D ) -1 U (x ∞ ) r + ρ .
The cluster point is nothing but the expected equilibrium point (x, c). It shows that the obtained trajectory lies in the stable manifold of (3).

Relation with the HJB equation; numerical approach of the lovebirds problem

We shall reinterpret the problem by means of optimal control theory, which in turn will lead to define the stable manifold through the resolution of the Hamilton-Jacobi-Bellman equation. Let us denote

f : (x, c) -→ -rx + c, g : (x, c) -→ U(x) -D(c).
Given the function c, the evolution equation for

X x 0 ,c reads d dt X x 0 ,c = f (X x 0 ,c , c) X x 0 ,c (0) = x 0 .
We want to find

J : x 0 -→ sup c J (x 0 , c) = sup c ∞ 0 e -ρt g X x 0 ,c (t), c(t) dt.
Let us apply the Dynamic Programming Principle [1, Section 3.1.2, Th. 3.1]: for any t > 0, we have

J(x 0 ) = sup c t 0 e -ρs g(X x 0 ,c (s), c(s)) ds + e -ρt J(X x 0 ,c (t)) . ( 10
)
Accordingly, we deduce that J solves the following stationary HJB equation [1, Section 3.1.4, Th.

3.3] ρJ(x 0 ) = H(x 0 , J (x 0 )) (11)
with

H : (x 0 , p) → sup κ≥0 g(x 0 , κ) + pf (x 0 , κ) = g x 0 , (D ) -1 (p) + pf x 0 , (D ) -1 (p)
since p = D (κ) at the optimum. For the sake of completeness, the proof of (10) and ( 11) is sketched in the Appendix. We can therefore rewrite (11) as follows

ρJ(x 0 ) = g x 0 , (D ) -1 (J (x 0 ) + J (x 0 )f x 0 , (D ) -1 (J (x 0 )) .
Equation ( 11) is completed by J(x) = c. Unfortunately solving (11) is at least as complicated as solving the integral relation [START_REF] Michael | Some properties of viscosity solutions of Hamilton-Jacobi equations[END_REF]. Nevertheless, we bear in mind that what we are really interested in is to determine the initial control c 0 (x 0 ) such that the trajectory issued from (x 0 , c 0 (x 0 )) lies on the stable manifold. Knowing (x 0 , c 0 (x 0 )), the whole trajectory can be obtained by coming back to the differential system (3), which can be solved by standard numerical routines for ODEs. To this end we go back to the Dynamic Programming Principle. We shall use (10) to find an approximation of the initial control. The method we propose works as follows.

i) We need a discrete set of initial well-beings

P wb = {x 1 0 , ...., x J 0 }, with step size 0 < ∆x = x j+1 0 -x j 0 1.
ii) We consider a discrete set of controls P control = {c 1 0 , ..., c K 0 }, with K J and step size 0

< ∆c = c k+1 0 -c k 0 1.
Of course, we use as far as possible the a priori knowledge of the dynamical system (e. g. provided by statements like Lemma 3.2, Lemma 3.5 and the phase portrait analysis) so that the initial controls associated to the elements of P wb likely belong to [c 1 0 , c K 0 ]. iii) We pick 0 < τ 1, which plays the role of time step. The controls are assumed to remain constant on the time interval [0, τ ]. iv) Given c k 0 and x j 0 , we approximate the well being at time τ by using the simple forward Euler approximation of (1); it defines

X j,k (τ ) = x j 0 + τ (-rx j 0 + c k 0 ).
v) Starting from J 0 = 0, we construct a sequence of functions (J m ) m≥1 , defined on P wb as follows

J m+1 (x j 0 ) = max k∈{1,...,K} 1 2 τ g((x j 0 , c k 0 ) + e -ρτ g(X j,k (τ ), c k 0 ) + e -ρτ J m (X j,k (τ )) .
(12) Namely, we use the trapezoidal rule for replacing the integral in the right hand side of (10). The definition of the last term in (12) is slightly misleading because J m is not defined on a continuous set, but on the discrete set P wb only. Hence, J m (X j,k (τ )) needs to be defined through a suitable interpolation procedure using the available values of J m . This is the reason why we cannot consider a single point in P wb . We denote by c k m+1 j 0 the obtained maximizers, k m+1 j ∈ {1, . . . , K} being the index of the discrete control that solves (12). vi) We wish that J m (x j 0 ) resembles J(x j 0 ) as m → ∞, while the c k m j 0 's approach the controls c 0 (x j 0 ) that we seek. In practice we stop the algorithm when the relative error between two consecutive iterates remains under a (small) given threshold. vii) This procedure allows to associate to each x j 0 ∈ P wb a control c k j 0 , depending on ∆x, ∆c, which is intended to approximate c 0 (x j 0 ) with (x j 0 , c 0 (x j 0 )) on the stable manifold. Having at hand the pairs (x j 0 , c k j 0 ) we solve the dynamical system (3) with a standard high-order scheme, for instance the Runge-Kutta algorithm. Actually we can complete this approach with a dichotomy procedure from the guess ck j 0 in order to determine a fair evaluation of the control c 0 (x j 0 ). We point out that the computational effort essentially relies on the iteration procedure in Steps i)-vi), while in comparison the resolution of the differential system is almost cost-free. Then, having obtained a fair estimate of c 0 in Steps i)-vi), it is worth refining the result by computing a few trajectories.

We illustrate the behavior of the method with the following simple example:

U : x → 5 log(1 + x) and D : c → (c -0.2) 2 2 , ρ = 0.1, r = 2.
The equilibrium point is x 0.772, c 1.544 and the phase portrait can be found in Figure 1. We work with 256 points equidistributed in the range 0 ≤ x ≤ 3; it defines P wb . For the set of controls, we work with 4 × 128 points equidistributed in the range 1 ≤ c ≤ 5; it defines P control while the time step is τ = 0.01. The dynamical system (3) is approached with the 4th order Runge-Kutta scheme with the constant time step ∆t = 10 -2 . In the figures below, the trajectories are all represented up to the same final time T = 2.

In Figure 2, the snapshot represents the trajectories computed for the initial state x 0 = 3 for several values of initial control. The algorithm described above returns c 0 2.750 and we compute trajectories for a couple of controls around this value. Based on the numerical experiment represented on the picture, we assert that we have found c 0 up to 10 -3 . In Figure 3, the snapshot corresponds to the same data with the initial state x 0 = 2.53. We obtain c 0 2.81. Based on the numerical experiment represented on the picture, we assert that we have found c 0 up to 10 -3 . In Figure 4, we change the value of the impatience factor: ρ = 1 here. The positive eigenvalue λ + defined in ( 8) is a non decreasing function of ρ and thus, the larger ρ, the more unstable the system. The problem thus should be more demanding for numerics. However, the method is still able to produce a decent approximation of the control. The snapshot corresponds to the initial state x 0 = 2.53 and we find c 0 2.57. Based on the numerical experiment represented on the picture, we assert that we have found c 0 up to 10 -2 .

Remark 4.1

The discussion uses the fact that the dynamical system is autonomous. This is a questionable aspect of the model which does not take into account external events which can impact, positively or negatively, either the well-being or the effort policy. As mentioned in [START_REF] Rey | A mathematical model of sentimental dynamics accounting for marital dissolution[END_REF] it could be relevant also to insert in the modeling alerts and corrections strategies intended to restore from time to time a deteriorated effort policy at a safe level.

Remark 4.2

The model can be slightly generalized. For instance it can be relevant to deal with a non linear well-being equation, generalizing (1) to d dt x = φ(x) + c, with φ a suitable non increasing function. On the mathematical viewpoint it could be interesting also to think the other way around: given a 2 × 2 differential system with stable and unstable directions at equilibrium, can we design a numerical method for finding points on the stable manifold? The difficulty to handle such a general case relies on the definition of the cost function, analog of J.

A Proof of (10) and (11)

We follow [START_REF] Barles | Solutions de viscosité des équations de Hamilton-Jacobi[END_REF] and we also refer to [START_REF] Evans | Partial Differential Equations[END_REF]Chap. 10] for further details and different arguments. Let us denote K(x 0 , t) = sup c t 0 e -ρs g(X x 0 ,c (s), c(s)) ds + e -ρt J(X x 0 ,c (t)) the rhs of (10). A priori this quantity depends on 0 ≤ t < ∞. Let c be a fixed function and denote X = X x 0 ,c (t). By definition of J(X x 0 ,c (t)) = J( X) as a supremum, for any > 0, there exists t → c (t) such that

J( X) -≤ ∞ 0 e -ρs g(X X,c (s), c (s)) ds = J ( X, c ) ≤ J( X).
We set c(s) = c(s)1 0≤s≤t + c (s -t)1 s>t . We have

J(x 0 ) ≥ ∞ 0 e -ρs g(X x 0 ,c (s), c(s)) ds ≥ t 0 e -ρs g(X x 0 ,c (s), c(s)) ds + ∞ t e -ρs g(X x 0 ,c (s), c (s -t)) ds ≥ t 0 e -ρs g(X x 0 ,c (s), c(s)) ds + e -ρt ∞ 0 e -ρs g(X x 0 ,c (t + s), c (s)) ds.
because for 0 ≤ s ≤ t, c(s) = c(s) and thus X x 0 ,c , X x 0 ,c satisfy the same ODE and are issued form the same initial state x 0 . Furthermore, we get X = X x 0 ,c (t) = X x 0 ,c (t), which is the data at s = 0 of s → X x 0 ,c (t + s) and s → X X,c (s), both being solution of the ODE d ds y(s) = f (y(s), c (s)) for s ≥ 0. We deduce that X x 0 ,c (t + s) = X X,c (s) for any s ≥ 0. We thus arrive at

J(x 0 ) ≥ t 0 e -ρs g(X x 0 ,c (s), c(s)) ds + e -ρt ∞ 0 e -ρs g(X X,c (s), c (s)) ds ≥ t 0 e -ρs g(X x 0 ,c (s), c(s)) ds + e -ρt J ( X, c ) ≥ t 0 e -ρs g(X x 0 ,c (s), c(s)) ds + e -ρt (J( X) -).
It holds for any > 0 and any function c so that we infer J(x 0 ) ≥ K(x 0 , t).

Similarly, we can find t → c (t) such that

J(x 0 ) -≤ ∞ 0 e -ρs g(X x 0 ,c (s), c (s)) ds = J (x 0 , c ) ≤ J(x 0 ).
Proceeding the same way, we are led to J(x 0 ) -≤ t 0 e -ρs g(X x 0 ,c (s), c (s)) ds + e -ρt ∞ 0 e -ρs g(X x 0 ,c (t + s), c (s)) ds.

The last integral is dominated by J(X x 0 ,c (t)) and thus we obtain

J(x 0 ) -≤ t 0
e -ρs g(X x 0 ,c (s), c (s)) ds + e -ρt J(X x 0 ,c (t)) ≤ K(x 0 , t).

Letting go to 0 ends the proof of (10).

We turn to the proof of (11). A crucial ingredient is the following simple remark: for any ϕ ∈ C 1 (R), we have

d dt e -ρt ϕ(X x 0 ,c (t)) = e -ρt -ρϕ(X x 0 ,c (t)) + ϕ (X x 0 ,c (t))f (X x 0 ,c (t), c(t)) (13)
owing to the chain rule. Let us explain formally the formula (11): assume that the supremum in (10) is attained by a certain function t → c (t) and that J is smooth. Then, for any t > 0 and κ ≥ 0, (10) implies

1 t t 0 e ρs g(X x 0 ,c (s), c (s)) ds + e -ρt J(X x 0 ,c (t)) -J(x 0 ) t = 0 ≥ 1 t t 0 e ρs g(X x 0 ,κ (s), κ) ds + e -ρt J(X x 0 ,κ (t)) -J(x 0 ) t .
Letting t → 0 yields g(x 0 , c (0)) + J (x 0 )f (x 0 , c (0)) -ρJ(x 0 ) = 0 ≥ g(x 0 , κ) + J (x 0 )f (x 0 , κ) -ρJ(x 0 ). by using (13). It proves that ρJ(x 0 ) = H(x 0 , J (x 0 )) holds.

The technical difficulty is related to a lack of regularity of the function J, which make these manipulations questionable. To circumvent the difficulty, one needs to introduce the notion of viscosity solutions of (11). Definition A. [START_REF] Barles | Solutions de viscosité des équations de Hamilton-Jacobi[END_REF] We say that J is a viscosity solution of (11) if for any trial function ϕ ∈ C 1 such that x 0 is a local minimum (resp. maximum) of J -ϕ, we have H(x 0 , φ (x 0 )) -ρφ(x 0 ) ≤ 0 (resp. ≥ 0). Clearly if J is a viscosity solution which is known to be smooth at x 0 , it suffices to set ϕ = J to obtain (11). Conversely, as far J is when x 0 is an extremum of J -ϕ we have J (x 0 ) = ϕ (x 0 ) so that a classical solution of (11) is a viscosity solution. The name "viscosity solution" comes from the possible approximation of such solutions by the solutions u of the elliptic equations ρu (x 0 ) -H(x 0 , u (x 0 )) = u (x 0 ) as the regularizing parameter > 0 goes to 0, and the analogies with the theory of conservation laws in gas dynamics. We refer on these aspects to the seminal works [START_REF] Michael | Some properties of viscosity solutions of Hamilton-Jacobi equations[END_REF][START_REF] Michael | Viscosity solutions of Hamilton-Jacobi equations[END_REF][START_REF] Lions | Generalized solutions of Hamilton-Jacobi equations[END_REF].

Let us discuss (11) in the framework of viscosity solutions. Let us first assume that x 0 is a local minimum of x → J(x) -ϕ(x). Adding a constant to ϕ if necessary, we thus have: J(x 0 ) = ϕ(x 0 ), for some δ > 0, J(x) ≥ ϕ(x) holds for any |x -x 0 | ≤ δ.

Let κ > 0; using (10) yields ϕ(x 0 ) = J(x 0 ) ≥ t 0 e -ρs g(X x 0 ,κ (s), κ) ds + e -ρt J(X x 0 ,κ (t)).

By continuity, as far as t > 0 is small enough, |X x 0 ,κ (t) -x 0 | ≤ δ so that J(X x 0 ,κ (t)) ≥ ϕ(X x 0 ,κ (t)) and therefore ϕ(x 0 ) -e -ρt ϕ(X x 0 ,κ (t)) t ≥ 1 t t 0 e -ρs g(X x 0 ,κ (s), κ) ds.

By using (13), letting t → 0 leads to +ρϕ(x 0 ) -ϕ (x 0 )f (x 0 , κ) ≥ g(x 0 , κ).

This relation holds for any κ ≥ 0 and we conclude that ρϕ(x 0 ) -H(x 0 , ϕ (x 0 )) ≥ 0 holds. Next, we assume that x 0 is a local maximum of x → J(x) -ϕ(x):

J(x 0 ) = ϕ(x 0 ), for some δ > 0, J(x) ≤ ϕ(x) holds for any |x -x 0 | ≤ δ.

By using (10) as above, we obtain ϕ(x 0 ) = J(x 0 ) ≤ sup e -ρs -ρϕ(X x 0 ,c (s)) + H X x 0 ,c (s), ϕ (X x 0 ,c (s)) ds .

Finally, we divide by t > 0 and let t go to 0; by using a continuity argument we arrive at 0 ≤ -ρϕ(x 0 ) + H(x 0 , ϕ (x 0 )).

  A satisfaction function x → U(x) returns a valuation of the feeling state x. While the shape of the function obeys general rules, that are listed below, the details depend again on the individual situation: different couples might give different valuations to the same common feeling. Producing efforts at a certain level reduces the satisfaction: the function c → D(c) evaluates the dissatisfaction due to the effort rate c. The two functions U (utility) and D (disutility) are defined on [0, ∞), and they are required to satisfy the following properties: i) U and D are (at least) C 2 non negative functions, ii) U (x) > 0 for any x ≥ 0 and lim x→∞ U (x) = 0, iii) U (x) ≤ 0 for any x ≥ 0, iv) 0 < α ≤ D (c) ≤ A < ∞ for any c ≥ 0 and lim c→∞ D (c) = ∞, v) there exists c ≥ 0 such that D(c ) ≥ 0, D (c ) = 0.

Theorem 2 . 2

 22 Let x 0 ≥ 0. Let c ♥ be the maximizer of J (x 0 , •). The following assertions hold true: a) The function t → c ♥ (t) does not vanish and remains positive a. e.. b) It is thus characterized by (2).

Figure 2 :

 2 Figure 2: Trajectories issued from x 0 = 3 and several values of c 0 around the estimated value c 0 = 2.75.

Figure 3 :

 3 Figure 3: Trajectories issued from x 0 = 2.53 and several values of c 0 around the estimated value c 0 = 2.81.

Figure 4 :

 4 Figure 4: Trajectories issued from x 0 = 2.53 and several values of c 0 around the estimated value c 0 = 2.57.

c t 0 e≤ sup c t 0 esup c t 0 e

 000 -ρs g(X x 0 ,c (s), c(s)) ds + e -ρt ϕ(X x 0 ,c (t)) But (13) allows us to writee -ρt ϕ(X x 0 ,c (t)) -ϕ(x 0 ) = t 0 e -ρt -ρϕ(X x 0 ,c (s)) + ϕ (X x 0 ,c (s))f (X x 0 ,c (s), c(s)) ds which leads to 0 -ρs g(X x 0 ,c (s), c(s)) -ρϕ(X x 0 ,c (s)) + ϕ (X x 0 ,c (s))f (X x 0 ,c (s), c(s)) ds ≤ -ρs sup κ≥0 g(X x 0 ,c (s), κ) -ρϕ(X x 0 ,c (s)) + ϕ (X x 0 ,c (s))f (X x 0 ,c (s), κ) ds ≤ sup c t 0