Gaëtan Leurent
email: gaetan.leurent@inria.fr

Lei Wang
email: wang.lei@ntu.edu.sg

The Sum Can Be Weaker Than Each Part

Keywords: Hash functions, combiners, XOR combiner, preimage attack

In this paper we study the security of summing the outputs of two independent hash functions, in an effort to increase the security of the resulting design, or to hedge against the failure of one of the hash functions. The exclusive-or (XOR) combiner H1(M )⊕H2(M ) is one of the two most classical combiners, together with the concatenation combiner H1(M ) H2(M ). While the security of the concatenation of two hash functions is well understood since Joux's seminal work on multicollisions, the security of the sum of two hash functions has been much less studied. The XOR combiner is well known as a good PRF and MAC combiner, and is used in practice in TLS versions 1.0 and 1.1. In a hash function setting, Hoch and Shamir have shown that if the compression functions are modeled as random oracles, or even weak random oracles (i.e. they can easily be inverted -in particular H1 and H2 offer no security), H1 ⊕ H2 is indifferentiable from a random oracle up to the birthday bound. In this work, we focus on the preimage resistance of the sum of two narrowpipe n-bit hash functions, following the Merkle-Damgård or HAIFA structure (the internal state size and the output size are both n bits). We show a rather surprising result: the sum of two such hash functions, e.g. SHA-512 ⊕ Whirlpool, can never provide n-bit security for preimage resistance. More precisely, we present a generic preimage attack with a complexity of Õ(2 5n/6 ). While it is already known that the XOR combiner is not preserving for preimage resistance (i.e. there might be some instantiations where the hash functions are secure but the sum is not), our result is much stronger: for any narrow-pipe functions, the sum is not preimage resistant. Besides, we also provide concrete preimage attacks on the XOR combiner (and the concatenation combiner) when one or both of the compression functions are weak; this complements Hoch and Shamir's proof by showing its tightness for preimage resistance. Of independent interests, one of our main technical contributions is a novel structure to control simultaneously the behavior of independent hash computations which share the same input message. We hope that breaking the pairwise relationship between their internal states will have applications in related settings.

Introduction

Hash functions are a very important class of primitive in modern cryptography, used in almost every secure system. A hash function H : {0, 1} * → {0, 1} n 1 takes an arbitrary length input and produces an n-bit output or digest. Hash functions are used in many setting with various security requirements; the general expectation is that a hash function should behave like a random function from {0, 1} * to {0, 1} n . More concretely, the main security notions expected from a hash function are:

Collision resistance. It should be hard to find two messages M = M with H(M ) = H(M ). Second-preimage resistance. Given a message M , it should be hard to find M = M with H(M ) = H(M ). Preimage resistance. Given a target hash value H, it should be hard to find M with H(M ) = H.

Since generic collision attacks require 2 n/2 work, and generic preimage attacks require 2 n work, a secure hash function should have the same level of resistance.

In order to build more secure hash functions, or to protect oneself against future cryptanalysis advances, such as the devastating attacks of Wang et al. against the SHA family [START_REF] Wang | How to Break MD5 and Other Hash Functions[END_REF][START_REF] Wang | Finding Collisions in the Full SHA-1[END_REF], a practical countermeasure might be to combine two different hash functions. The goal is that the combined hash function can only be broken when both components are weak. In particular, this reasoning was used by the designers of SSL [START_REF] Freier | The Secure Sockets Layer (SSL) Protocol Version 3.0[END_REF] and TLS [START_REF] Dierks | The TLS Protocol Version 1.0. RFC 2246[END_REF], who combined MD5 and SHA-1 in various ways. More precisely, the Key Derivation Function of TLS v1.0/v1.1 uses a sum of HMAC-MD5 and HMAC-SHA-1. 3 The designers explain [START_REF] Dierks | The TLS Protocol Version 1.0. RFC 2246[END_REF]: "In order to make the PRF as secure as possible, it uses two hash algorithms in a way which should guarantee its security if either algorithm remains secure."

There are two classical hash function combiners: the concatenation combiner H 1 (M ) H 2 (M ) and the XOR combiner H 1 (M ) ⊕ H 2 (M ). In a seminal work [START_REF] Joux | Multicollisions in Iterated Hash Functions. Application to Cascaded Constructions[END_REF], Joux showed that the concatenation combiner with narrow-pipe hash functions offers much less security than could be expected: it has 2n bits of output, but essentially offers the same security as an n-bit hash function. In this work, we carry a similar analysis for the XOR combiner. Previous work has shown that it is indifferentiable from a random oracle up to the birthday bound [START_REF] Hoch | On the Strength of the Concatenated Hash Combiner When All the Hash Functions Are Weak[END_REF], even if the initial functions are weak; in particular, it has optimal collision resistance of n/2 bits. However, we show that the preimage security is much less than one might expect, with a generic preimage attack with complexity Õ(2 5n/6 ).

Since the goal of a combiner is to keep some security even if one of the functions is found to be weak, it is natural that the two hash functions H 1 and H 2 are independent in practice. Throughout this paper the two hash functions used in a combiner are always assumed to be independent without specifying it explicitly 4 . Iterated hash function. In this paper we consider iterated hash functions, following the Merkle-Damgård construction [START_REF] Merkle | One Way Hash Functions and DES[END_REF][START_REF] Damgård | A Design Principle for Hash Functions[END_REF] or the more general HAIFA construction [START_REF] Biham | A Framework for Iterative Hash Functions -HAIFA[END_REF], as shown in Figure 1. We focus on narrow-pipe designs, i.e. we assume that the internal state size is the same as the output size n. The message M is first split into blocks m 0 , . . . m , and the hash function iterates a series of compression functions h i over an internal state x, with the initial value denoted as IV . Finally, the hash value is computed with a finalization function g:

x 0 = IV x i+1 = h i (x i , m i ) H(M ) = g(x +1 , |M |)
In the following, we assume that the compression function is the same at every step (∀i, h i = h) in order to simplify the notations, but it is straightforward to apply the attack with the corresponding function at each step.

|M |

IV g

Fig. 1. Iterated hash function.

Related works

Combiners have been studied in several different settings. For generic attacks, the compression functions are modeled as random functions, in order to devise attacks that don't use any weakness of the compression functions. On the other hand, some work assumes that the compression function is a weak random oracle (that can easily be inverted), and prove that some constructions are still secure in this model. Finally, more theoretical work focus on the notion of robustness, i.e. the existence of a reduction from attacks on the hash functions to attacks on the combiner.

Analysis of the concatenation combiner. The concatenation combiner H 1 (M ) H 2 (M ) is probably the most studied one. In 2004, Joux [START_REF] Joux | Multicollisions in Iterated Hash Functions. Application to Cascaded Constructions[END_REF] described surprising attacks on the concatenation of two narrow-pipe hash functions using multicollisions: while the output size is 2n bits, the concatenation can at most provide n/2-bit security for collision resistance and n-bit security for preimage resistance 5 . In particular, the concatenation is not security-amplifying. On the other hand, the concatenation combiner is robust for preimages and collisions, which gives a matching lower bound for generic attacks.

Later, Hoch and Shamir [START_REF] Hoch | On the Strength of the Concatenated Hash Combiner When All the Hash Functions Are Weak[END_REF] evaluated the security of the concatenation combiner with two weak hash functions. More precisely, the two hash functions are narrow-pipe Merkle-Damgård, and the compression functions are modeled as weak random oracles (as defined by Liskov [START_REF] Liskov | Constructing an Ideal Hash Function from Weak Ideal Compression Functions[END_REF]), i.e., the adversary is given additional interfaces to receive (random) preimages of the compression functions. They have proven that in this model, the concatenation combiner is still indifferentiable from a random oracle with n/2-bit security, implying (at least) the same security bound for collision resistance and preimage resistance. The bound is matched by Joux's attack for collisions, but there is a gap with Joux's attack for preimages, with complexity 2 n , which might be interesting to investigate further.

Mendel et al. analyzed some dedicated instantiations of the concatenation combiner [START_REF] Mendel | MD5 Is Weaker Than Weak: Attacks on Concatenated Combiners[END_REF], in particular using the hash function MD5. We omit the details and refer interested readers to [START_REF] Mendel | MD5 Is Weaker Than Weak: Attacks on Concatenated Combiners[END_REF].

Analysis of the XOR combiner. The XOR combiner has received less analysis. The work of Hoch and Shamir [START_REF] Hoch | On the Strength of the Concatenated Hash Combiner When All the Hash Functions Are Weak[END_REF] actually proves the security of the XOR combiner as an intermediate result: it is also indifferentiable from a random oracle up to 2 n/2 queries in the weak random oracle model. In particular, this proves that there are no generic attacks with complexity smaller than 2 n/2 . For collision resistance, the bound is tight, since it is matched with the generic birthday attack bound. On the other hand, for preimage resistance, there exists a gap between the n/2-bit proven bound and the n-bit expected ideal security bound.

To the best of our knowledge, no preimage attacks have been shown against the XOR combiner. Therefore, the preimage security of the XOR combiner against generic attacks is still an open problem, and will be the main topic of our work. We will also consider the preimage security of the XOR combiner with weak hash functions, and study the tightness of Hoch and Shamir's bound.

Robust combiners. In the last years, the general problem of combining two (or more) hash functions H 1 and H 2 has been extensively studied from a theoretical point of view. These works focus on the notion of a robust combiner: a robust combiner is secure with respect to property α as long as one of the underlying hash functions is secure for α. It can be shown that the concatenation combiner is a robust combiner for collision-resistant hash functions and for MACs, while the XOR combiner is robust for PRFs and for MACs [START_REF] Lehmann | On the Security of Hash Function Combiners[END_REF]. More advanced combiners have been constructed in order to be robust for multiple properties simultaneously [START_REF] Fischlin | Multi-property Preserving Combiners for Hash Functions[END_REF][START_REF] Fischlin | Robust Multi-property Combiners for Hash Functions Revisited[END_REF][START_REF] Fischlin | Robust Multi-Property Combiners for Hash Functions[END_REF]. The notion was mostly studied via the black-box reduction model. A series of results have showed that robust combiners for collision resistance and preimage resistance cannot have an output length significantly shorter than the sum of the output length of the underlying hash functions [START_REF] Boneh | On the Impossibility of Efficiently Combining Collision Resistant Hash Functions[END_REF][START_REF] Pietrzak | Non-trivial Black-Box Combiners for Collision-Resistant Hash-Functions Don't Exist[END_REF][START_REF] Rjasko | On Existence of Robust Combiners for Cryptographic Hash Functions[END_REF][START_REF] Mittelbach | Hash Combiners for Second Pre-image Resistance, Target Collision Resistance and Pre-image Resistance Have Long Output[END_REF]. Since the XOR combiner is length preserving, this shows that it is not robust for collision resistance and preimage resistance.

Actually, the impossibility results are in part due to the stringent requirement from the black-box reduction model. In order to overcome this limitation, Mittel-bach introduced the idealized random oracle model [START_REF] Mittelbach | Cryptophia's Short Combiner for Collision-Resistant Hash Functions[END_REF]. He gives a construction of a short output combiner with optimal collision and preimage security in this model 6 (assuming that one of the functions is ideal): Cryptophia's short combiner uses the sum of two hash functions with some pre-processing of the messages (to allow non-independent functions).

More generally, we point out that a combiner being non-robust does not necessarily mean there is an attack. The non-robustness results only show that the security of the combiner cannot be proved with a reduction from the security of the hash functions. In particular, the XOR combiner is not robust for collisionresistance, or even collision-resistance preserving. However, Hoch and Shamir's work proves that there are no generic collision attacks on this construction, either with ideal compression function, or even with weak compression functions. This arguably makes the XOR a useful combiner for collision resistance. Regarding preimage security, the non-robustness result does not imply that the XOR of two concrete hash functions is weak, and the simplicity and short output of this construction still make it quite attractive.

Our results

In this work, we study the preimage security of the XOR combiner, the main remaining open problem for classical combiners. We show that, surprisingly, the sum of two narrow-pipe n-bit hash functions can never achieve n-bit security for preimage resistance. More precisely, we find a generic preimage attack with a complexity of Õ(2 5n/6 ). It does not exploit any structural weakness of the compression functions and hence is applicable even if the compression functions are two ideal random oracles. Thus, even if the two hash functions are n-bit secure for preimage resistance, the XOR combiner is at most 5n/6-bit secure for primage resistance. In other words, the sum can be weaker than each part.

The attack is based on a novel technique to break the pairwise relationship between the internal states of the two hash functions. More precisely, the two hash functions H 1 and H 2 share the same input message, and hence the internal states of their iterative compression function computations are related. We control the computation chains of H 1 and H 2 simultaneously by constructing a new message structure M, and two sets of internal states A for H 1 and B for H 2 such that: for any value A from A and any value B from B, we can derive a message M A,B from M such that H 1 (M A,B ) produces A and H 2 (M A,B ) produces B. Hence we can select states from A and B independently. After that, we use a birthday match to find a message block m, a value A from A and a value B from B such that h 1 (A, m) ⊕ h 2 (B, m) is equal to the target hash digest, where h 1 and h 2 are the compression functions of H 1 and H 2 respectively. Finally we derive the message M A,B from M, and output M A,B m as a preimage of the target hash digest.

Our preimage attack is also applicable to Cryptophia's short combiner [START_REF] Mittelbach | Cryptophia's Short Combiner for Collision-Resistant Hash Functions[END_REF][START_REF] Mennink | Breaking and Fixing Cryptophia's Short Combiner[END_REF]. This construction has been proven to provide optimal collision and preimage resistance, assuming that at least one of the initial functions is a monolithic random oracle, but our attack does not violate the security proof, because we use the fact that both functions have an iterative structures with an n-bit internal state. Still, this shows that with many practical hash functions, the combiner will be weaker that the initial functions. Our results also show that the XOR combiner and Cryptophia's combiner are not robust in the semi-black-box model introduced by Mittelbach [29]7 , even with independant hash functions H 1 and H 2 .

Our analysis on the XOR combiner is also interesting for dedicated hash function design. The hash function family RIPEMD [START_REF] Dobbertin | RIPEMD-160: A Strengthened Version of RIPEMD[END_REF] is based on a compression function with two parallel lanes, added together at the end of each compression function. Interestingly, RIPEMD-160 has been quite resilient to cryptanalysis [START_REF] Mendel | On the collision resistance of RIPEMD-160[END_REF][START_REF] Sasaki | Distinguishers beyond Three Rounds of the RIPEMD-128/-160 Compression Functions[END_REF][START_REF] Mendel | Differential attacks on reduced RIPEMD-160[END_REF][START_REF] Mendel | Improved Cryptanalysis of Reduced RIPEMD-160[END_REF], and are still considered secure. Several more recent designs use parallel lanes in a similar way (combining them at the end of each compression function call), such as HAS-V [START_REF] Park | HAS-V: A New Hash Function with Variable Output Length[END_REF], FORK [START_REF] Hong | A New Dedicated 256-Bit Hash Function: FORK-256[END_REF] and LANE [START_REF] Indesteege | The lane hash function[END_REF]. It might be tempting to use parallel lanes during the full iteration, and to combine them only at the end. Indeed, the designers of SHA-V [START_REF] Her | A Design of Cryptographic Hash Function Group with Variable Output-Length Based on SHA-1[END_REF] used this approach: the 160-bit version of SHA-V has two parallel lanes, combined at the end with a modular sum. This is equivalent to summing two different hash functions, and hence our attack can be applied to SHA-V.

Another contribution of this paper is to present concrete preimage attacks on the XOR combiner with one or both weak hash functions (defined in [START_REF] Liskov | Constructing an Ideal Hash Function from Weak Ideal Compression Functions[END_REF]). The complexity of our attacks is Õ(2 n/2 ). Furthermore, the attack can be extended to the concatenation combiner with two weak hash functions under the same complexity. It can be seen that these attacks match the bound of Hoch and Shamir's security proof [START_REF] Hoch | On the Strength of the Concatenated Hash Combiner When All the Hash Functions Are Weak[END_REF], and hence fulfill the gaps pointed out in Section 1.1. It implies the tightness of Hoch and Shamir's proof on the classical combiners with weak hash functions for preimage resistance.

Finally, we would like to highlight the technical interests of this paper. We devise a novel structure named interchange structure to simultaneously control two (or more) hash lanes with the same input message, and succeed in further relaxing the pairwise relation between the internal states of lanes. It is indeed a step of technical advance compared with previous extensive studies on this topic, and hence will hopefully have applications or lead to new technical development in related settings. We refer to the open discussions in Section 7 for more details.

Notations and roadmap in the rest of paper

We use the following notations: Finally we conclude the paper and discuss future directions in Section 7.

H 1 ,

Overview of the attack

We first give an overview of the techniques and the structures used in the attack, while more detailed descriptions will be given in the following sections. The main idea is to consider several chains of internal states reached by processing a common message M from different starting points (note that the message M is not fixed in advance, but will be determined when building the structure). More precisely, the message M is denoted as the primary message, and divided in several chunks: M = M 0 M 1 . . . (as discussed later, a chunk will consist of several message blocks). We denote chains of internal states for H 1 as a j , and the individual states of the chain as a i j , with h * 1 (a i j , M i ) = a i+1 j . Similarly, we denote chains for

H 2 as b k , with h * 2 (b i k , M i ) = b i+1 k .
When considering both hash functions, message block M i leads from the pair of states (a i j , b i k ) to (a i+1 j , b i+1 k ), which is denoted:

(a i j , b i k ) Mi (a i+1 j , b i+1 k ).
Switch structure. Next we build special structures called switches in order to jump between chains in a controlled way. A switch allows to jump from a specific pair of chains (a j0 , b k0 ) to a different pair of chains (a j0 , b k1 ) using a secondary message chunk M i , in addition to the normal transitions using chunk M i of the primary message M :

(a i j , b i k ) Mi (a i+1 j , b i+1 k ) : normal transition for each chain (a i j0 , b i k0 ) Mi (a i+1 j0 , b i+1 k1 ) : jump from chains (a j0 , b k0 ) to (a j0 , b k1 )
In order to simplify the notations, we often omit the chunk index, in order to show only the chains that are affected by the switch. The main message chunk M i and the secondary message chunk M i are determined when building the switch, and the main message defines the next state of all the chains. We note that the secondary message chunk M i should only be used when the state is (a i j0 , b i k0 ). A simple example is depicted in Figure 2. Alternatively, a switch can be designed to jump from (a j0 , b k0 ) to (a j1 , b k0 ). We defer the details of the construction to Section 3; it can be built with a complexity of Õ(2 n/2 ).

Interchange structure. By combining several simple switches, we can build an interchange structure with starting points IV 1 and IV 2 and ending points A j , j = 0 . . . 2 t -1 and B k , k = 0 . . . 2 t -1 , so that we can select a message ending in any state (A j , B k ). Figure 3 shows one possible way to build such a structure, and Figure 4 shows how to select a given message in the structure. An interchange structure with 2 t chains for each function requires about 2 2t switches. Since we can build a switch for a cost of Õ(2 n/2 ), the total structure is built with Õ(2 2t+n/2 ) operations.

Preimage search. Finally, we can use an interchange structure with ending points A j , j = 0 . . . 2 t -1 and B k , k = 0 . . . 2 t -1 , to build a preimage attack as follows. Let H denote the target value. We select a random message block m, and we compute two lists by evaluating the compression functions after the interchange structure:

A j = h 1 (A j , m), j = 0 . . . 2 t -1 and B k = H ⊕ h 2 (B k , m), k = 0 . . . 2 t -1 .
We expect a match between the lists with probability 2 2t-n . After about 2 n-2t random choices of m, we get a match (j * , k * ): Therefore, we can construct a preimage of H by concatenating the message leading to (A j * , B k * ) in the interchange structure, and the block m (we ignore the finalization function in this section).

h 1 (A j * , m) = H ⊕ h 2 (B k * , m) i.e. h 1 (A j * , m) ⊕ h 2 (B k * , m) = H.
The complexity of the preimage search is about 2 n-t evaluations of the compression function, using an interchange structure with 2 t end-points.

Complexity analysis. Building the interchange structures requires about 2 2t+n/2 evaluations of the compression function, while the preimage search requires about 2 n-t . The optimal complexity is reached when both steps take the same time, i.e. t = n/6. This gives a complexity of Õ(2 5n/6 ).

The switch structure

We now explain how to build the switch structure at the core of our attack. This construction is strongly based on the multicollision technique of Joux [START_REF] Joux | Multicollisions in Iterated Hash Functions. Application to Cascaded Constructions[END_REF].

Given states a i j0 , b i k0 and b i k1 , we want to build message chunks M i and M i in order to have the following transitions:

(a i j0 , b i k0 ) Mi (a i+1 j0 , b i+1 k0 ) (a i j0 , b i k1 ) Mi (a i+1 j0 , b i+1 k1 ) (a i j0 , b i k0 ) Mi (a i+1 j0 , b i+1 k1 ).
The main message chunk M i is used to define the next state of all the remaining chains, while the secondary message chunk M i will be used to jump from chains (a j0 , b k0 ) to (a j0 , b k1 ). We note that M i will only be used when the state is (a i j0 , b i k0 ). In particular, M i and M i must satisfy:

a i+1 j0 = h * 1 (a i j0 , M i ) = h * 1 (a i j0 , M i ) b i+1 k1 = h * 2 (b i k1 , M i ) = h * 2 (b i k0 , M i ) b i+1 k0 = h * 2 (b i k0 , M i ) = b i+1 k1
We first build a multicollision for h * 1 , starting from state a i j0 , i.e. a large set M of 2 n/2 messages that all reach the same state a i+1 j0 (∀M ∈ M, h * 1 (a i j0 , M ) = a i+1 j0 ). As shown by Joux, this can be done efficiently by sequentially building n/2 collisions.

Next, we evaluate h * 2 (b i k0 , M ) and h * 2 (b i k1 , M ) for all the messages M in the set M. With high probability there is match between the sets of values8 . We denote the colliding messages as M i and M i , so that we have h * 2 (b i k0 , M i ) = h * 2 (b i k1 , M i ). Finally we compute the missing chains using the message m i :

a i+1 j = h * 1 (a i j , m i ), b i+1 k = h * 2 (b i k , m i ).
With high probability all the chains reach distinct values; if this is not the case, we restart the construction with a new multicollision. The full algorithm is shown as Algorithm 1, and illustrated by Figure 5; it requires about n/2 • 2 n/2 evaluations of the compression functions.

The interchange structure

Let us now describe the combination of switch structures into an interchange structure. The goal of this structure is to select the final value of the H 1 computation and the H 2 computation independently. More precisely, the structure defines two sets of final values A j and B k , and a set of messages M jk such that:

(IV 1 , IV 2 ) M jk (A j , B k ).
In order to build this structure, we initialize the first chains with a 0 0 = IV 1 , b 0 0 = IV 2 , and set the other starting points randomly. Then, we use switches to

H1 H2 a i j 0 a i+1 j 0 M b i k 0 b i k 1 M M b i+1 k 1 M M b i+1 k 0 M Fig. 5
. Building a switch structure. First, M and M are selected from M to generate a collision (defining the new b k 1 ), then b k 0 is evaluated using M .

jump for an already reachable pair (a j0 , b k0 ) to a different pair (a j0 , b k1 ) (or to (a j1 , b k0 ), respectively). By using 2 2t -1 switches, we can make all pairs reachable. There are many way to combine the switches; a simple one can be described as follow:

1. first, build switches from (a 0 , b 0 ) to each of the (a 0 , b k )'s; 2. then for each k, build a series of switches from (a 0 , b k ) to all the (a j , b k )'s.

In order to reach the chains (a j , b k ), one would activate the k-th switch in the first part to jump from (a 0 , b 0 ) to (a 0 , b k ), and then the j-th switch in the k-th series of the second part to jump from (a 0 , b k ) to (a j , b k ). This structure is shown in Figure 3 and a pseudo-code description is given by Algorithm 2, where the Interchange functions builds the structure, and the SelectMessage function extracts the message reaching (a j , b k ). The structure can be somewhat optimized using the fact that the extra chains have no prespecified initial values. We show how to take advantage of this in Appendix A, using multicollision structures in addition to the switch structures. However, this doesn't change significantly the complexity: we need (2 t -1)(2 t -1) switches instead of 2 2t -1. In total, we need about n/2 • 2 2t+n/2 evaluations of the compression functions to build a 2 t -interchange structure.

We believe that a 2 t -interchange structure based on switches will need at least Θ(2 2t ) switches, because every switch can only increase the number of reachable pairs (a j , b k ) by one. As shown in Appendix A some switches can be saved in the beginning but it seems that new ideas would be needed to reduce the total complexity below Θ(2 2t+n/2 ).

Preimage Attack

Finally, we describe the full preimage attack. We first build an interchange structure with 2 t chains for each of H 1 and H 2 . We denote the ending points as A j , j = 0 . . . 2 t -1 and B k , k = 0 . . . 2 t -1 , and we know how to select a Algorithm 1 Building a single switch message M jk to reach any state (A j , B k ). When adding a message block m to one of the messages M jk in the interchange structure, the output of the combiner can be written as:

function Switch(h1, h2, a, b, b ) x ← a M ← ∅ for 0 ≤ i < n/2 do (m, m ) ← Collision(h1, x) M ← (M m) ∪ (M m ) x ← h1(x, m) end for H ← {} for M ∈ M do y ← h * 2 (b, M ) H[y] ← M end for for M ∈ M do y ← h * 2 (b , M ) if H[y]
H 1 (M jk m) ⊕ H 2 (M jk m) = g 1 (h 1 (A j , m), + 1) ⊕ g 2 (h 2 (B k , m), + 1),
where g 1 and g 2 are the finalization functions of H 1 and H 2 , respectively, and is the length of the messages in the structure.

In order to reach a target value H, we select a random block m, and we evaluate

A j = g 1 (h 1 (A j , m), + 1), j = 0 . . . 2 t -1 and B k = H ⊕ g 2 (h 2 (B k , m), +
2 n-2t times on average, therefore the total cost of the preimage search is about 2 n-t evaluations of h 1 and h 2 .

As explained in the previous section, building a 2 t -interchange structure requires about n/2 • 2 2t+n/2 operations. Using t = n/6 we balance the two steps of the attack, and reach the optimal complexity of about n/2 • 2 5n/6 operations for the preimage attack.

Message length and memory complexity

The attack uses messages of length n/2 • 2 2t , and the memory complexity of the attack10 is also n/2 • 2 2t . The optimal choice t = n/6 gives messages of length n/2 • 2 n/3 . The memory requirement is probably not an issue 11 for an attacker that can spend time 2 5n/6 , but the message length can be a problem with some hash functions that don't accept long inputs. For instance SHA-256 is only defined for message with less than 2 64 bits (i.e. 2 55 blocks).

In this case, one can apply the attack with a smaller value of t: this reduces the length of the messages, at the cost of more time spent in the preimage search step. For instance, we can mount a preimage attack against SHA-256 ⊕ BLAKE-256 with complexity 2 232 using t = 24, while the optimal attack with n = 256 would cost only 2 220.3 . Similarly, our attack applied to SHA-512 ⊕ Whirlpool has a complexity of 2 461 , rather than 2 434.7 .

Applications and Extensions

The attack works identically if the hash functions use the HAIFA mode rather than the plain Merkle-Damgård iteration. Also it can easily be extended to H 1 (M ) H 2 (M ) where denotes an easy to invert group operation (for instance, a modular addition rather than the exclusive or). The attack can also be extended to hash functions H 1 and/or H 2 using an internal check-sum, such as the GOST family of hash functions, using pairs of blocks with a constant sum.

Application to the sum of wide-pipe hash functions

The attack can also be used when the internal state size is larger than the output size n. The complexity of building a 2 t -interchange structure is related to as /2 • 2 2t+ /2 . On the other hand, the complexity of the meet-in-the-middle preimage search is related to n as 2 n-t . The optimal complexity is /2 • 2 2n/3+ /6 by matching the two complexities with t = n/3 -/6. Therefore our attack can be applied as long as + 6 log( ) ≤ 2n holds. For instance, we can compute preimages of SHA-224 ⊕ BLAKE-224 with complexity roughly 2 199 .

Application to Cryptophia's short combiner

Our attack can also be applied to Cryptophia's short combiner, as proposed by Mittelbach [START_REF] Mittelbach | Cryptophia's Short Combiner for Collision-Resistant Hash Functions[END_REF], and to the revised version of Mennink and Preneel [START_REF] Mennink | Breaking and Fixing Cryptophia's Short Combiner[END_REF]. This combiner computes the sum of two hash functions with some pre-processing of the message, to allow non-independent functions:

C(M ) = H 1 m1 1 . . . m1 ⊕ H 2 m2 1 . . . m2 m1 j = H 1 (0 l 1 m j ⊕ k 1 ) ⊕ H 2 (0 l 2 m j ⊕ k 2 ) m2 j = H 1 (1 l 1 m j ⊕ k 1 ) ⊕ H 2 (1 l 2 m j ⊕ k 2 )
where k 1 , k 2 , l 1 , l 2 is a randomly chosen key. The security proof in the ideal model shows that C is optimally preimage resistant if at least one of the hash functions is ideal. However, if both H 1 and H 2 are narrow-pipe, we can apply our preimage attack with complexity Õ(2 5n/6 ). This does not violate the security proof because we need both functions to be narrow-pipe, hence not n-bit ideal 12 . From a practical point of view, though, it shows that in many cases (e.g. using SHA-512 and Whirlpool) the combiner is weaker than the initial functions.

Improvements using weaknesses of the hash functions

If H 1 or H 2 has known cryptographic weaknesses, more efficient attacks are possible. More precisely, if the compression function of one of the hash functions can be inverted 13 in time 2 t , then we can find a preimage of H 1 ⊕ H 2 with complexity only Õ(2 (n+t)/2 ).

The attack is presented using the case, where (at least) one compression function is modeled as a weak compression function defined in [START_REF] Liskov | Constructing an Ideal Hash Function from Weak Ideal Compression Functions[END_REF], as an example. Without loss of the generality, we assume the compression function of H 2 , h 2 (x, y) = z, is such a weak compression function, which is a random oracle with two additional interfaces as below.

• Backward interface. On a query (?, y, z), it returns either a value x uniformly chosen from all the values satisfying h 2 (x, y) = z, or ⊥ if no such x exists.

• Bridging interface. On a query (x, ?, z), it return either a value y uniformly chosen from all the values satisfying h 2 (x, y) = z, or ⊥ if no such y exists.

Note that the inversion of compression function h 2 takes unit time and hence the attack against H 1 ⊕ H 2 takes time Õ(2 n/2 ). The procedure is detailed as follows, which is also illustrated in Figure 6.

Let the target hash digest be denoted as H. We firstly build an n-block long multicollision on H 1 following Joux's approach [START_REF] Joux | Multicollisions in Iterated Hash Functions. Application to Cascaded Constructions[END_REF]. Let the final output be denoted as Z. It contains a set of up to 2 n messages that link IV 1 to Z on H 1 . We split every n-block long multicolliding message into halves, and collect the first half n/2-block long messages as a set M 1 and the second half as another set M 2 . Hence for any message M 1 ∈ M 1 and any M 2 ∈ M 2 , the concatenated message M 1 M 2 links IV 1 to Z on H 1 . Secondly, for the messages M 2 ∈ M 2 , we use the additional backward interface of h 2 to carry out backward computations from H ⊕ Z, and get values X such that h * 2 (X, M 2 ) = H ⊕ Z. We store (X, M 2 ) in a table T X . On average T X contains 2 n/2 elements. 14 Finally, for each message M 1 ∈ M 1 , we iteratively compute h 2 forward from IV 2 , and get an internal state Y . We match Y to the elements in T X . A match implies that concatenated M 1 M 2 links IV 2 to H ⊕ Z, and in turn is a preimage of

H 1 (M 1 M 2 ) ⊕ H 2 (M 1 M 2 ) = Z ⊕ H ⊕ Z = H.
The success probability of finding such a match is not negligible since there are 2 n/2 X's and Y 's. The complexity is around n • 2 n/2 that is Õ(2 n/2 ) by ignoring the polynomial factors. Moreover, the above preimage attack can be extended to the concatenation combiner H 1 H 2 if both H 1 and H 2 are weak by a minor modification. This shows that the proof of Hoch and Shamir [START_REF] Hoch | On the Strength of the Concatenated Hash Combiner When All the Hash Functions Are Weak[END_REF] is tight for preimage resistance. Here we mainly highlight the modifications. Let the target hash digest be H 1 H 2 , where H 1 is from H 1 and H 2 from H 2 . After we build the multicollision on H 1 and let the output internal state be denoted as Z, we link Z to H 1 by using the bridging interface of h 1 to receive a message m such that h 1 (Z, m) = H 1 . This gives us a set of messages linking IV 1 to H 1 on H 1 . Also note that for the backward computations on H 2 , the starting value should be H 2 . The rest of the attack procedure remains the same. Hence it is easy to get that the complexity is also Õ(2 n/2 ).

Extension to the sum of three or more hash functions

The attack can be extended to the sum of three or more hash functions. In order to attack the sum of k functions, two different strategies are possible: either we use a simpler structure that only gives two degrees of freedom, and fixes k -2 functions to a constant value, or we build an interchange structure to control all the k functions independently.

Controlling only two functions. The easiest way to extend the attack is to use a single chain in the k -2 extra hash functions. The procedure to build a switch is modified in order to use multicollisions for k -1 functions instead a simple multicollisions for one function; this costs O(n k-1 • 2 n/2 ) using Joux's method [START_REF] Joux | Multicollisions in Iterated Hash Functions. Application to Cascaded Constructions[END_REF].

As in the basic attack, we need O(t 2 ) switches to generate a 2 t -interchange for two functions, and the preimage search costs O(2 n-t ); the optimal complexity is therefore O(n k-1 • 2 5n/6 ) with t = n/6. Controlling all the functions. Alternatively, we can build a more complex interchange structure in order to control all the functions independently. When attacking three functions, we will use the switch structure to jump from chains (a j0 , b k0 , c l0 ) to (a j0 , b k0 , c l1 ) (or (a j0 , b k1 , c l0 ) or (a j1 , b k0 , c l0 ), respectively). We need 2 3t -1 switches in the interchange structure to reach all the 2 3t triplets of chains (a switch makes only one new triplet reachable). Each switch is built using a 2 n/2 -multicollision on two functions, which can be built for a cost of O(n 2 • 2 n/2 ) following Joux's technique [START_REF] Joux | Multicollisions in Iterated Hash Functions. Application to Cascaded Constructions[END_REF]. Therefore we can build a 2 t -interchange for a cost of O(n 2 • 2 3t+n/2 ). More generally, for the sum of k hash functions, we can build an interchange structure for k functions for a cost of O(n k-1 • 2 kt+n/2 ).

In the preimage search phase, we generate k lists of size 2 t , and we need to detect efficiently whether we can combine then to generate a zero sum. This problem can be solved using a algorithm similar to Wagner's generalized birthday algorithm [START_REF] Wagner | A Generalized Birthday Problem[END_REF]. If k = 2 κ , we find a solution with probability O(2 n-(κ+1)•t ) for a cost of O(k • 2 t ). Therefore the preimage search costs O(k • 2 n-κt ). With k = 4 (i.e. κ = 2), this yields a complexity of O(n 3 • 2 5n/6 ). However, this approach is less efficient than the previous one for k = 3 and for k > 4.

To summarize, attacking the sum of k hash functions (k ≥ 2) costs O(n k-1 • 2 5n/6 ). Controlling chains independently in more than two hash function might be useful for further work, but it doesn't improve the preimage attack on the sum of k hash functions.

Conclusion and open discussions

In this work, we gave the first generic attack on the XOR combiner. Our result is rather surprising: the sum of two ideal narrow-pipe hash functions only has about 5n/6 bits of security against preimage attacks. In particular, the sum is easier to break than the initial functions. Since most practical hash functions are narrow-pipe (e.g. SHA-1, SHA-256, SHA-512, Whirlpool, RIPEMD, GOST, BLAKE, Skein...), the XOR combiner will usually provide a weaker security than the component hash functions.

Moreover, we would like to discuss a few directions for future work.

On controlling multiple hash lanes. Since 2004, several generic attacks have been found against narrow-pipe hash functions, such as the multicollision attack [START_REF] Joux | Multicollisions in Iterated Hash Functions. Application to Cascaded Constructions[END_REF], the long-message second preimage attack [START_REF] Kelsey | Second Preimages on n-Bit Hash Functions for Much Less than 2 n Work[END_REF] and the herding attack [START_REF] Kelsey | Herding Hash Functions and the Nostradamus Attack[END_REF]. There has been extensive work to extend these attacks to more complex constructions with several computation chains, such as the concatenation H 1 (M ) H 2 (M ) and the cascade H 2 (H 1 (M ), M ).

As in our present work, the essential difficulty in those attacks comes from the fact that several lanes of computation share the same input message, and hence their outputs are related. If an adversary considers a naïve set of messages, the set of outputs gives random pairs of n-bit values {(A i , B i ) : i ∈ I}: selecting a value for the first entry of the pair gives a single candidate for the second entry, and the adversary is essentially working with a 2n-bit state. Previous works [START_REF] Nandi | Multicollision Attacks on Some Generalized Sequential Hash Functions[END_REF][START_REF] Hoch | Breaking the ICE -Finding Multicollisions in Iterated Concatenated and Expanded (ICE) Hash Functions[END_REF][START_REF] Andreeva | Herding, Second Preimage and Trojan Message Attacks beyond Merkle-Damgård[END_REF] have developed various message structures (mostly based on multicollision and diamond structures), in order to relax this relation. They mainly result in a set of messages M such that the corresponding outputs are in a more structured set {(A, B i ) : i ∈ I}: the first entry is a constant value A, but several options B i can be selected for the second entry. For any value (A, B i ), it is then possible to select a message in M so that the first lane reaches A, while the second lane reaches B i . This allows to modify the value of the second lane without affecting the first lane.

Our result is quite stronger: with the interchange structure we have a set of message such that the corresponding outputs are a set {(A i , B j ) : i ∈ I, j ∈ J }, where both lanes have several options that can be selected independently. We hope that our technique will have applications or lead to new technical development in related settings, e.g., the open problem of generic second preimage attacks (with long messages) on the concatenation hash or on the Zipper hash [START_REF] Liskov | Constructing an Ideal Hash Function from Weak Ideal Compression Functions[END_REF].

On extending to practical hash function. Several practical hash functions such as RIPEMD [START_REF] Dobbertin | RIPEMD-160: A Strengthened Version of RIPEMD[END_REF] and HAS-V [START_REF] Park | HAS-V: A New Hash Function with Variable Output Length[END_REF] are based on a compression function with more than one independent lanes, which interacts with each other at the end of each compression function call. It is very interesting to investigate if our attack can be further modified to attack these hash functions in future. Again the obstacle comes from the relations between internal states of lanes. Particularly, the internal states of the lanes interact with each other, which makes the relation even tighter and in turn harder to attack.

Fig. 2 .

 2 Fig.2. A single switch: jump from (a0, b0) to (a0, b1) by using M (dashed lines) instead of M (solid lines).

Fig. 3 .Fig. 4 .

 34 Fig. 3. Overview of an interchange structure

Fig. 6 .

 6 Fig. 6. Preimage attack on the XOR combiner with a weak H2

  H 2 : hash functions IV 1 , IV 2 : initial values for H 1 and H 2 , respectively h 1 , h 2 : compression functions of H 1 and H 2 , respectively , b k : chains for H 1 and H 2 , respectively a j denotes a generic chain, while a j0 denotes a particular chain A j , B k : end points of the chains n : hash function output size Roadmap. Section 2 provides an overview of our generic preimage attack on the XOR combiner. Sections 3, 4, and 5 elaborate the attack procedure step by step in details. Section 6 presents the applications and extensions of the attack.

	h * 1 , h * 2 : compression functions iterated over several blocks (in particular, H i (M ) = h * i (IV i , M ))
	m	: message block
	M	: message chunk (n/2 blocks)
	M	: long message (several chunks)
	a j

We note that this MD5/SHA-1 combiner has been replaced by primitives based on single hash function (e.g., SHA-256) since TLS v1.2[START_REF] Dierks | The Transport Layer Security (TLS) Protocol Version 1.2. RFC 5246[END_REF].

If the two hash functions can be related, it is trivial that the XOR combiner is not security-preserving. For instance, let H2(M ) := H1(M ) ⊕ const, where const is a constant. If H1 is ideally secure, then H2 is also ideally secure. However, the XOR combiner H1(M ) ⊕ H2(M ) = const for any message M .

The attacks actually require only one of the functions to be narrow-pipe.

A mistake in the initial proof and construction was later fixed by Mennink and Preneel[START_REF] Mennink | Breaking and Fixing Cryptophia's Short Combiner[END_REF].

Loosely speaking, a combiner is robust with respect to property α if it is (at least) as secure as the stronger underlying hash function for α.

If this is not the case, we build a new multicollision.

It takes O(t • 2 t ) operations by sorting the lists, but only 2 • 2 t using a hash table.

We only need to store the messages M and M

For instance, the attack is on the verge of practicality with n = 64; the time complexity is 2 58.3 and the memory complexity is 226.3 .

A large multi-collisions can be built with a cost of roughly 2 n/2 in a narrow-pipe function, but costs almost 2 n for an ideal hash function.

finding an input chaining value from the output chaining value and the message

The backward interface may output ⊥ for some message block. To compensate it, for the other message blocks, we make multiple queries, since they may have more than one preimages. On average, the backward interface should produce one preimage for each query.

Acknowledgements

The authors would like to thank the anonymous referees for their helpful comments. Lei Wang is supported by the Singapore National Research Foundation Fellowship 2012 (NRF-NRFF2012-06).

Algorithm 2 Building and using a T -interchange structure function Interchange(h1, h2, IV1, IV2) a0 ← IV1, b0 ← IV2 for 1 ≤ k < T do a k ← $, b k ← $ end for for 1 ≤ j < T do (M, M ) ← Switch(h1, h2, a0, b0, bj)

If there is a match (j * , k * ) between the two lists, we have:

For a random choice of m, we expect that a match exists with probability 2 2t-n , and testing it requires about 2 t operations 9 . We will have to repeat this procedure

A Optimized Interchange Structure

We now describe an optimized attack using only (2 t -1)(2 t -1) switches rather than 2 2t -1. The attack also requires multicollision structures, as introduced by Joux [START_REF] Joux | Multicollisions in Iterated Hash Functions. Application to Cascaded Constructions[END_REF].

We replace the first 2 t -1 switches with a 2 t -multicollision in H 1 , and we use those messages to initialize all the b k chains in H 2 . We can also optimize the first series of switches in H 2 in the same way: we build a 2 t -multicollision in H 2 starting from b 0 , and we use those messages to initialize the a j chains in H 1 . This is illustrated by Figure 7, and the detailed attack is given as Algorithm 3.


Julien Alexandre

Sandretto

Alexandre Chapoutot
email: alexandre.chapoutot@ensta-paristech.fr

Validated Simulation of Differential Algebraic Equations with Runge-Kutta Methods

Keywords: Differential algebraic equations, Validated simulation. AMS subject classifications: 65L80, 65G20, 65G40

Differential Algebraic Equations (DAEs) are a general and implicit form of differential equations. This mathematical object is often used to represent physical systems such as dynamics of solid or chemical interactions. These equations are different from Ordinary Differential Equations (ODEs) in sense that some of the dependent variables occur without their derivatives. These variables are called "algebraic variables", which means free of derivatives and not with respect to abstract algebra. Validated simulation of ODEs has recently known different developments such as guaranteed Runge-Kutta integration schemes, explicit and implicit ones. Not so far from an ODE, solving a DAE consists of searching a consistent initial value and computing a trajectory. Nevertheless, DAEs are in generally much more difficult to solve than ODEs. In this paper, we focus on the semi-explicit form of index one, called Hessenberg index-1 form. We propose a validated way to simulate this kind of differential equations. Finally, our method is applied to different examples in order to show its efficiency.

Chapter 1

Motivations

Our recent results on validated simulation of Ordinary Differential Equations (ODEs) with implicit Runge-Kutta schemes [START_REF] Alexandre | Validated Solution of Initial Value Problem for Ordinary Differential Equations based on Explicit and Implicit Runge-Kutta Schemes[END_REF] lead us to go up in complexity of the kind of differential equations we want to solve. Indeed, we are able to simulate ODEs with interval parameters which is one of the requirements for our solver of Differential Algebraic Equation (DAEs). Indeed, a DAE can be seen as a parametrized ODE. We currently focus on DAE in Hessenberg index-1 form, that is ẏ = f (t, x, y) 0 = g(t, x, y)

(

with f : R n × R m × R → R n and g : R n × R m × R → R m . In Equation (1.1), y ∈ R n is the vector of state variables and x ∈ R m is the vector of algebraic variables (without an expression for its derivative) and ẏ stands for the time derivative of y. This kind of DAE is common and used by a majority of simulation tools as Simulink and Modelica-like software. In this paper, we will present a validated method to solve initial value problem given in the form of a DAE.

The article is organized as followed. In Section 2, a short description of our ODE solver and by the same way a state of the art of validated simulation of ODEs will be done. In Section 3, we will perform a review of literature about few existing validated approaches for DAEs. In Section 4, we present in details our method. In Section 5, we apply our solver on three different examples: a basic problem, a problem with exact solution known and the classic pendulum problem. In Section 6, we present how to take into account additional constraints in IVP for DAE before concluding in Section 7.

Notations

ẏ denotes the time derivative of y, i.e., dy dt . a denotes a real values while a represents a vector of real values. [a] represents an interval values and [a] represents a vector of interval values. The midpoint of an interval [x] is denoted by m([x]). The variables y is assigned to the state variables of the system, when x is reserved for the algebraic variables and t obviously kept for the time. The functions f and g are used to represent a differential algebraic equation system, with f the function given the time derivative of state variable and g the constraint on the algebraic variable. We also denote y(t, yj) the exact solution of a differential equation at time t with a known value yj at time tj, in our case it implies that t > tj. Sets will be represented by calligraphic letter such as X or Y.

Chapter 2

Validated simulation of Ordinary differential equations

A simulation of an ordinary differential equation consists on a discretization of time and an iterative and approximative computation of the state of the system, by the help of an integration scheme. In details, an integration scheme, starting from an initial value yn at time tn and a finite time horizon h, the stepsize, produces an approximation yn+1 at time tn+1, with tn+1 -tn = h, of the solution y(tn+1; yn). In this section, we briefly recall how a validated simulation of an ODE is computed. We refer to [START_REF] Nedialkov | Validated solutions of initial value problems for ordinary differential equations[END_REF] for a more detailled presentation. The notations follow [START_REF] Nedialkov | Validated solutions of initial value problems for ordinary differential equations[END_REF].

Ordinary differential equations

An initial value problem (IVP) defined through an ODE is defined by: ẏ = f (t, y) with y(0) = y0, y0 ∈ R n and t ∈ [0, t end ] .

(2.1)

In the classical approach [START_REF] Lohner | Enclosing the solutions of ordinary initial and boundary value problems[END_REF][START_REF] Nedialkov | Validated solutions of initial value problems for ordinary differential equations[END_REF] to define validated method for IVP, each step of an integration scheme consists in two steps: a priori enclosure and solution tightening. Starting from a valid enclosure [yj] at time tj, the two following steps are applied

Step 1. Compute an a priori enclosure [ỹj] of the solution using Banach's theorem and the Picard-Lindelöf operator. This enclosure has the three major properties:

• y(t, [yj]
) is guaranteed to exist for all t ∈ [tj, tj+1], i.e., along the current step, and for all yj ∈ [yj].

• y(t, [yj]) ⊆ [ỹj] for all t ∈ [tj, tj+1].
• the step-size hj = tj+1 -tj is as larger as possible in terms of accuracy and existence proof for the IVP solution.

Step 2. Compute a tighter enclosure of [yj+1] such that y(tj+1, [yj]) ⊆ [yj+1]. The main issue in this phase is how to counteract the well known wrapping effect 5 [START_REF] Moore | Interval Analysis[END_REF][START_REF] Lohner | Enclosing the solutions of ordinary initial and boundary value problems[END_REF][START_REF] Nedialkov | Validated solutions of initial value problems for ordinary differential equations[END_REF]. This phenomenon appears when we try to enclose a set with an interval vector (geometrically a box). The arising overestimation creates a false dynamic for the next step, and, with accumulation, can lead to intervals with an unacceptably large width.

The different enclosures computed during one step are shown on Figure 2.1.

t y(t)

[ỹ j ] [y j ] [y j+1 ] t j t j+1 h j

Figure 2.1: Enclosures appeared during one step

Some algorithms useful to perform these two steps are described in the following.

A priori solution enclosure

In order to compute the a priori enclosure, we use our interval version of Picard-Lindelöf operator. This operator is based on the following theorem.

Theorem 2.2.1 (Banach fixed-point theorem) Let (K, d) a complete metric space and let g :

K → K a contraction that is for all x, y in K there exists c ∈]0, 1[ such that d (g(x), g(y)) ≤ c • d(x, y) ,
then g has a unique fixed-point in K.

In context of IVP, we consider the space of continuously differentiable functions C 0 ([tj, tj+1], R n ) and the Picard-Lindelöf operator

P f (y) = t → yj + t tn f (s, y(s))ds . (2.2)
Note that this operator is associated to the integral form of Equation (2.1). So the solution of this operator is also the solution of Equation (2.1).

The Picard-Lindelöf operator is used to check the contraction of the solution on an integration step in order to prove the existence and the uniqueness of the solution of Equation (2.1) as stated by the Banach's fixed-point theorem. Furthermore, this operator is used to compute an enclosure of the solution of IVP over a time interval [tj, tj+1].

Rectangular method for the a priori enclosure

Using interval analysis and with a first order integration scheme we can define a simple interval Picard-Lindelöf operator such that Remark that the operator defined in Equation (2.3) can also define a contractor (in a sens of interval analysis [START_REF] Chabert | Contractor programming[END_REF]) on [r] after the contraction proved for

P f ([r]) = [yj] + [0, h] • f ([r]), (2.3
P f such that [r] ← [r] ∩ [yj] + [0, h].f ([r]) . (2.4)
Hence, we can reduce the width of the a priori enclosure in order to increase the accuracy of the integration. The operator defined in Equation (2.3) and its associated contractor defined in Equation (2.4) can be defined over a more accurate integration scheme (at the condition that it is a guaranteed scheme like the interval rectangle rule). For example, the evaluation of t t j f (s)ds can be easily improved with a Taylor or a Runge-Kutta scheme (see [START_REF] Alexandre | Validated Solution of Initial Value Problem for Ordinary Differential Equations based on Explicit and Implicit Runge-Kutta Schemes[END_REF]).

A priori enclosure with Taylor series

Interval version of Taylor series for ODE integration gives

[yj+1] ⊂ N k=0 f [k] ([yj])h k + f [N +1] ([ỹj])h N +1 , (2.5)
with

f [0] = [yj], f [1] = f ([yj]),. . . , f [k] = 1 k ( ∂f [k-1] ∂y f )([yj]
). By replacing h with interval [0, h], this scheme becomes an efficient Taylor Picard-Lindelöf operator, with a parametric order N such that

yj+1([tj,tj+1]; [r]) = yj + N k=0 f [k] ([yj])[0, h k ] + f [N +1] ([r])[0, h N +1 ] . (2.6) In consequence, if [r] ⊇ yj+1 ([tj, tj+1], [r]
), then Equation (2.6) defined a contraction map and Theorem 2.2.1 can be applied.

In our tool, we use it at order 3 by default, it seems to be a good compromise between efficiency and computation quickness.

Note that the scheme defined in Equation (2.5) is usually evaluated in a centered form for a more accurate result

[yj+1] ⊂ N k=0 f [k] (ŷj)h k + f [N +1] ([ỹj])h N +1 + N k=0 J(f [k] , [yj])h i )([yj] -ŷj , (2.7) with ŷj ∈ [yj] J(f [k] , [yj]
) is the Jacobian of f [k] evaluated at [yj]. This scheme can also be combined with a QR-factorization to increase stability and counteract wrapping effect [START_REF] Nedialkov | Validated solutions of initial value problems for ordinary differential equations[END_REF]. These two "tricks", with a strong computational cost, can be avoided by using the affine arithmetic.

Picard-Lindelöf operator, as defined in Equation (2.6), gives an a priori enclosure [r], using Theorem 2.2.1. Picard-Lindelöf operator is proven to be contracting on [r], we can then use this operator to contract the box [r] till a fixpoint is reached

In our tool, the default contractor uses a Taylor expansion as follow

[r] ∩ yj + N k=0 f [k] ([yj])[0, h k ] + f [N +1] ([r])[0, h N +1 ]
It is very important to contract as much as possible this box [r] because the Taylor remainder is function of [r] and the step-size is function of the Taylor remainder.

Tighter enclosure and truncation error

Suppose that Step 1 has been done for the current integration step and that we dispose of the enclosure [ỹj] such that

y(t, tj, [yj]) ⊆ [ỹj] ∀t ∈ [tj, tj+1] .
In particular, we have y(tj+1, tj, One way to do that consists in computing an approximate solution yj+1 ≈ y(tj+1, tj, [yj]) with an integration scheme Φ(tj+1, tj, [yj]), and then the associated local truncation error LT EΦ(t, tj, [yj]). Indeed, a guaranteed integration scheme has the property that there exists a time ξ ∈ [tj, tj+1] such that

y(tj+1, tj, [yj]) ⊆ Φ(tj+1, tj, [yj]) + LT EΦ(ξ, tj, [yj]) ⊆ [ỹj] . So [yj+1] = Φ(tj+1, tj, [yj]) + LT EΦ(ξ, tj, [yj]) is an acceptable tight enclosure.
The guaranteed solution of IVP using interval arithmetic is mainly based on two kinds of methods to compute Φ:

1. Interval Taylor series methods [START_REF] Moore | Interval Analysis[END_REF][START_REF] Lohner | Enclosing the solutions of ordinary initial and boundary value problems[END_REF][START_REF] Berz | Verified integration of odes and flows using differential algebraic methods on high-order taylor models[END_REF][START_REF] Nedialkov | Validated solutions of initial value problems for ordinary differential equations[END_REF][START_REF] Kapela | A lohner-type algorithm for control systems and ordinary differential inclusions[END_REF][START_REF] Wilczak | Cr-lohner algorithm[END_REF][START_REF] Chen | Taylor model flowpipe construction for non-linear hybrid systems[END_REF][START_REF] Lin | Validated solutions of initial value problems for parametric odes[END_REF],

2. Interval Runge-Kutta methods [START_REF] Gajda | A survey of interval runge-kutta and multistep methods for solving the initial value problem[END_REF][START_REF] Bouissou | GRKLib: a Guaranteed Runge Kutta Library[END_REF][START_REF] Bouissou | Enclosing temporal evolution of dynamical systems using numerical methods[END_REF].

The former is the oldest method used in this context. Indeed, R. Moore [START_REF] Moore | Interval Analysis[END_REF] already applied this method in the sixties and until now it is the most used method to solve Equation (2.1). The latter is more recent, see in particular [START_REF] Bouissou | GRKLib: a Guaranteed Runge Kutta Library[END_REF][START_REF] Bouissou | Enclosing temporal evolution of dynamical systems using numerical methods[END_REF], but Runge-Kutta methods have many interesting properties, as strong stability, that we would like to exploit in the context of validated solution of DAEs. The challenge lies in the computation of LT EΦ.

Runge-Kutta methods

To obtain yn+1, a Runge-Kutta method computes s evaluations of f at predetermined time instants. The number s is known as the number of stages of a Runge-Kutta method. More precisely, a Runge-Kutta method is defined by

yn+1 = yn + h s i=1 biki , (2.8)
with ki defined by

ki = f tn + cih, yn + h s j=1 aijkj .
(2.9)

The coefficient ci, aij and bi, for i, j = 1, 2

. bs

In function of the form of the matrix A, made of the coefficients aij, a Runge-Kutta method can be

• explicit, e.g., the classical Runge-Kutta method of order 4 given in Figure 2.2(a).

In other words, the computation of an intermediate ki only depends on the previous steps kj for j < i. These methods are often the faster for a given order, but with less properties than the followers;

• diagonally implicit, e.g., a diagonally implicit method of order 4 given in Figure 2.2(b). In this case, the computation of an intermediate step ki involves the value ki and so non-linear systems in ki must be solved. Nevertheless, the structure of this kind of method leads to a faster solving;

• fully implicit, e.g., the Runge-Kutta method with a Lobatto quadrature formula of order 4 given in Figure 2.2(c). In this last case, the computation of intermediate steps involves the solution of a non-linear system of equations in all the values ki for i = 1, 2, • • • , s. In this class of implicit methods, some of them have strong properties such as A-stability and stiffly accurate capability.

Local truncature error

Recently, we defined an unified approach to express LT E for explicit and implicit Runge-Kutta methods. More precisely, for a Runge-Kutta of order p we have LTE(t, y(ξ)) := y(tn; yn-1) -yn =

h p+1 (p + 1)! r(τ )=p+1 α(τ ) [1 -γ(τ )ψ(τ )] F (τ )(y(ξ)) ξ ∈ [tn, tn+1] (2.10) with 0 0 0 0 0 1 2 1 2 0 0 0 1 2 0 1 2 0 0 1 0 0 1 0 1 6 1 3 1 3 1 6 (a) RK4
• τ is a rooted tree • F (τ ) is the elementary differential associated to τ • r(τ ) is the order of τ (number of nodes) • γ(τ ) is the density • α(τ ) is the number of equivalent trees • ψ(τ )
Note that y(ξ) is a particular, and a priori unknown, solution of Equation (2.1) at a time instant ξ. This solution can be over-approximated using Picard-Lindelöf operator presented in Section 2.2.

Wrapping effect

Runge-Kutta methods also suffer from the wrapping effect, such as discussed before. The problem of reducing this latter has been studied in many different ways [START_REF] Neumaier | The wrapping effect, ellipsoid arithmetic, stability and confidence regions[END_REF][START_REF] Nedialko (ned | Computing Rigorous Bounds on the Solution of an Initial Value Problem for an Ordinary Differential Equation[END_REF]. One of the most known and effective is the QR-factorization [START_REF] Lohner | Enclosing the solutions of ordinary initial and boundary value problems[END_REF]. This method improves the stability of the Taylor series in the Vnode-LP tool [START_REF] Nedialkov | Validated solutions of initial value problems for ordinary differential equations[END_REF]. An other way is to modify the geometry of the enclosing set (parallelepipeds [START_REF] Eijgenraam | The solution of initial value problems using interval arithmetic[END_REF][START_REF] Moore | Interval Analysis[END_REF], ellipsoids [START_REF] Neumaier | The wrapping effect, ellipsoid arithmetic, stability and confidence regions[END_REF], convex polygons [START_REF] Rihm | Interval methods for initial value problems in odes[END_REF] and zonotopes [START_REF] Stewart | A heuristic to reduce the wrapping effect in the numerical solution ofx= f (t, x)[END_REF][START_REF] Bouissou | HySon: Set-based simulation of hybrid systems[END_REF]).

In our work, an efficient affine arithmetic allows us to counteract the wrapping effect [START_REF] De Figueiredo | Self-Validated Numerical Methods and Applications[END_REF][START_REF] Alexandre | Validated Solution of Initial Value Problem for Ordinary Differential Equations based on Explicit and Implicit Runge-Kutta Schemes[END_REF].

Chapter 3

Differential Algebraic Equations

We can distinguish at least two families of DAEs, the fully implicit ones and the semi-explicit ones, the fully explicit are similar to ODEs.

Fully implicit DAEs

The first class of differential algebraic equations is the fully implicit ones. It is the most general representation of a differential system such as:

f (t, y, ẏ, ...) = 0, t0 ≤ t ≤ t end (3.1)
The order of this kind of DAE is defined with respect to the highest order of present derivatives. For example, if f is function of ẏ and no higher derivatives then f is a 1 st order DAE. If ÿ occurs then f is a 2 nd order DAE, and so on. Nevertheless, all DAEs can be rewritten in DAE of 1 st order by increasing the size of the system.

Semi-explicit DAEs or in Hessenberg form

The second, and the most used DAE form in science and engineering, is the semiexplicit DAEs, or also called DAEs in Hessenberg form. In this formalism, the index is a differentiation index [START_REF] Eleda Brenan | Numerical solution of initial-value problems in differential-algebraic equations[END_REF], that is to say the distance to the related ODE. For example, the index 1 Hessenberg form is described by:

ẏ = f (t, x, y) 0 = g(t, x, y) (3.2)
while the index 2 Hessenberg form is described by:

ẏ = f (t, x, y) 0 = g(t, x) (3.3)
In the way that some of dependent variables occur without their derivatives, these differential systems are different from an ODE with an additive constraint (which are independent). 11

Review of literature

If there exists many tools for solving DAEs in a numerical way (DAETS, Sundials, Mathematica) and some particular implementation embedded in simulation softwares (Simulink, Dymola), only few attempts have been done in a guaranteed approach. We can notice an extension of Valencia-IVP [START_REF] Rauh | A novel interval arithmetic approach for solving differential-algebraic equations with ValEncIA-IVP[END_REF], a tool for ODE validated simulation. The chosen approach starts with an approximation of the solution obtained by a numerical method (DASSL, DAETS) and try a posteriori to enclose the solution in a guaranteed way (with a Krawczyk iteration). The validity of this approach is not proven, in particular the fact that the existence test of the algebraic variable is separated from the state variable should not lead to a correct solution.

Another method [START_REF] Althoff | Reachability analysis of nonlinear differentialalgebraic systems[END_REF] consists to compute the reachable set of DAEs by an error linearization and the help of zonotopes. This approach seems close to [START_REF] Rauh | A novel interval arithmetic approach for solving differential-algebraic equations with ValEncIA-IVP[END_REF], and it is not clear on some points, such as existence and uniqueness proofs. Moreover, the results are not sufficiently explained to judge the quality of the method.

An alternative paper dealing with validated solution of DAEs available in the literature is based on Taylor models [START_REF] Hoefkens | Computing validated solutions of implicit differential equations[END_REF]. This method starts by computing a highorder polynomial approximating the solution of the ODE, and after that attempts to inflate it till validating an inclusion test. The approach presented is then different than our method. Unfortunately, none of these two latter approaches seems to be continued.

Finally, and more recently, one important work has been done in [START_REF] Josephk | Interval bounds on the solutions of semiexplicit index-one daes. part 1: analysis[END_REF][START_REF] Josephk | Interval bounds on the solutions of semiexplicit index-one daes. part 2: computation[END_REF]. In these papers, authors present the analysis and the computation of bounds by using i) an existence and uniqueness test for the solutions of DAEs, based on Hansen-Sengupta, and ii) sufficient conditions for two functions to enclose lower and upper bounds of the solution. Another method is presented unifying these two steps. Regrettably, the algorithms are not sufficiently clearly described to allow a real analysis and comparison. The results presented are also difficult to estimate. Nevertheless, we are able to make few remarks. Firstly, the presented approach is not completely validated (as honestly said in the paper). Secondly, if it is well known that Hansen-Sengupta test is more efficient than Krawczyk one [START_REF] Goldsztejn | Comparison of the hansen-sengupta and the frommer-lang-schnurr existence tests[END_REF], in [START_REF] Josephk | Interval bounds on the solutions of semiexplicit index-one daes. part 1: analysis[END_REF] it is used in non-parametric version while we advocate to use Krawczyk in parametric version. Indeed, Hansen-Sengupta requires the solving of a linear interval system, which can be time consuming function of problem size and done many times at each integration step (till obtain conditions of Theorem 2.2.1). Moreover, the operator is applied to a non linear function with state variable as interval parameter, which can be large due to initial state, the use of a parametric version is then necessary to obtain a sufficiently sharp solution (see Section 4.1).

It is important to remark that the problem of existence and uniqueness of the solution is a common issue considered in the literature, whether in the validated approaches -presented above -or in the numerical field, for which consistent initialization is one of the main issue [START_REF] Vieira | An overview of initialization approaches for differential algebraic equations[END_REF].

Chapter 4 Our Method for Validated Simulation of DAE

In this work, we present a method to solve the initial value problems written in index-1 Hessenberg form:

ẏ = f (t, x, y) 0 = g(t, x, y) with y(0) ∈ [y0] and x(0) ∈ [x0] . (4.1)
In Section 2.1, we recalled the classical two steps method presented before by Lohner [START_REF] Lohner | Enclosing the solutions of ordinary initial and boundary value problems[END_REF], and used by the community of ODE validated integration. We used the same approach for our DAE integration method, as [START_REF] Josephk | Interval bounds on the solutions of semiexplicit index-one daes. part 2: computation[END_REF]. The two steps are:

• Compute an a priori enclosure of all the solution of the DAE on an integration step [t, t + h] with a novel Picard-like operator;

• Refine this enclosure at t + h with integration Runge-Kutta scheme and contractors.

In the next subsections, we will present these two steps in details.

New Picard-like operator

For an ordinary differential equation with interval parameter (similar to differential inclusion) described by

ẏ = f (t,

Analysis of existence and uniqueness

If we assume that ∂g ∂x is locally reversal, we are theoretically able to find the unique x = ψ(y) (with the help of the implicit function theorem), and then we can write ẏ = f (t, ψ(y), y). Finally, we could apply Picard-Lindelöf to f in order to prove existence and uniqueness of the solution. Of course it is not feasible nor realistic in general, because ψ is unknown and uncomputable. Anyway, this approach leads us to the following solution.

We propose a theorem based on Frobenius theorem with a set membership view. Frobenius theorem is too long to be completely recalled here, but it is available in [START_REF] Dieudonné | Foundations of modern analysis[END_REF]. Proof: If ∀y ∈ Y, ∃! x ∈ X : g(x, y) = 0, then the Implicit Function Theorem prove that there is a unique function x = ψ(y) for all (x, y) ∈ A × B. Then the differential Equation (4.4) becomes ẏ = f (t, ψ(y), y) with ψ(y) maps Y into X . Finally, the statement ∀x ∈ X , ∃! y ∈ Y : ẏ = f (t, y, x) implies, with the help of Frobenius theorem, that there exists a unique mapping solution of Equation (4.4) in Y.

Adapted to DAEs

The operator

Finally, with the last Proposition 4.1.2, we define a novel operator Picard-Krawczyk PK:

If P([ỹ], [x]) K([ỹ], [x]) ⊂ Int [ỹ] [x] then ∃! solution of DAE in [ỹ] (4.6)
with

• P a Picard-Lindelöf operator for the differential inclusion ẏ ∈ f ([x], y)

• K a parametric preconditioned Krawczyk operator for the constraint g(x, y) = 0, ∀y ∈ [ỹ] In details, the operators are given by: Picard-Lindelöf operator with Taylor (N = 3):

P f ([yj], [xj], [r], [x], h) = [yj] + N k=0 f [k] ([xj], [yj])[0, h k ] + f [N +1] ([x], [r])[0, h N +1 ] . (4.7) If P f ([yj], [xj], [r], [x], h) ⊂ Int([r]) then f is integrable and [yj+1] ⊂ [r].
Parametric preconditioned Krawczyk operator:

Kg([ỹ], [r]) = m([r]) -Cg(m([r]), m([ỹ]))- (C ∂g ∂x ([r], [ỹ]) -I)([r] -m([r]))- C ∂g ∂y (m([r]), [ỹ])([ỹ] -m([ỹ])) (4.8) If Kg([ỹ], [r]) ⊂ Int([r]
) then for each y ∈ [ỹ] there exists one and only one x ∈ [r] : g(x, y) = 0.

Parametric preconditioned Krawczyk operator in hybrid form:

An hybrid form [START_REF] Goldsztejn | Sensitivity analysis using a fixed point interval iteration[END_REF] is also available:

[s] = Cg(m([r]), m([ỹ])) + C ∂g ∂y (m([r]), [ỹ])([ỹ] -m([ỹ])) [s] = [s] ∩ (Cg(m([r]), [ỹ])) Kg([ỹ], [r]) = m([r]) -[s] -(C ∂g ∂x ([r], [ỹ]) -I)([r] -m([r])) If Kg([ỹ], [r]) ⊂ Int([r]
) then for each y ∈ [ỹ] there exists one and only one x ∈ [r] : g(x, y) = 0.

Algorithm implementing the Picard-Krawczyk operator

We propose an algorithm implementing this Picard-Krawczyk operator to compute the enclosures of the algebraic and state variables simultaneously in Algorithm 1. The inputs are the initial estimation for enclosures 0 [x], 0 [ỹ], a tolerancy on LT E tol and a guessed step size h. The outputs are wheither a success flag with the validated enclosures 1 [x], 1 [ỹ] and the already computed lte, either a fail flag. We iterate two times the size of function f in order to prapagate enoughly the inflation.

Algorithm 1 Compute the a priori enclosures Require:

0 [x], 0 [ỹ], tol, h, iter = 0 1 [x] = K( 0 [ỹ], 0 [x]) // 4.1.3 1 [ỹ] = P( 0 [ỹ], 0 [x]) // 4.1.3 while ( 1 [x] ⊂ 0 [x]) and ( 1 [ỹ] ⊂ 0 [ỹ]) and (iter < 2 * size(f )) do iter = iter + 1 0 [ỹ] = 1 [ỹ] ± 1% 0 [x] = 1 [x] ± 1% 1 [x] = K( 0 [ỹ], 0 [x]) // 4.1.3 1 [ỹ] = P( 0 [ỹ], 0 [x]) // 4.1.3 end while if ( 1 [x] ⊂ 0 [x]) and ( 1 [ỹ] ⊂ 0 [ỹ]) then // conditions of 4.1.2 obtained lte = LTE( 1 [x], 1 [ỹ]) // 2.4 if lte < tol then // Acceptable lte return SUCCESS( 1 [x], 1 [ỹ], lte) end if end if return FAILED
// No enclosures found with the inputs

Contractors

After the guaranteed enclosures on [t, t + h] obtained, we can contract these enclosures around t + h, because obviously x(t + h) ⊂ [x] and y(t + h) ⊂ [ỹ]. Firstly, it is more efficient to start by the contraction of the state variable y, which can be strongly refined by the well chosen integration scheme. The differential inclusion to integrate is given by: ẏ ∈ f (t, [x], y). This differential equation is often stiff (this assumption will be verified in section 5.1) and has interval coefficient. As demonstrated in [START_REF] Alexandre | Validated Solution of Initial Value Problem for Ordinary Differential Equations based on Explicit and Implicit Runge-Kutta Schemes[END_REF], an implicit Runge-Kutta scheme is efficient for this kind of interval parametrized differential equation. In the family of Implicit Runge-Kutta (IRK), Radau IIA is one of the more powerful method at order 3 for stiff problem. Its Butcher tableau is given in Fig. 4.1. This IRK method attains an order three with two stages, it is fully implicit as shown in its tableau (Fig. 4.1) and A-stable. The corresponding local truncature error (lte) can be computed by using the Butcher trees as shown in [START_REF] Alexandre | Validated Solution of Initial Value Problem for Ordinary Differential Equations based on Explicit and Implicit Runge-Kutta Schemes[END_REF].

Secondly, with a tightened state variable, we are able to refine the algebraic variable x under the constraint g(t, x, y) = 0, ∀y ∈ [ỹ]. For this task, we combine the Krawczyk

Complete algorithm

The algorithm for a validated simulation of a DAE in index-1 Hessenberg form is given in a simplified way, without the step size controller (available in [START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF]), in Algorithm 3. The inputs are initial states, a time to reach, a given minimal step size and a tolerancy. The algorithm builds a list of data computed at each integration step, as described in Section 2.1.

Example with exact solution

The second example is a little more complex. It is described by the differential equation:

y =    -y2y1 -(1 + y2)x0 y2y0 -(1 + y2)x1 1
and the algebraic relation: This problem is interesting because it has a known exact solution:

(y0 -x1)/5 -cos(y 2
       y0 = sin(t) + 5 * cos((t 2 )/2.0) y1 = cos(t) + 5 * sin((t 2 )/2.0) x0 = -cos(t) x1 = sin(t)
We perform a simulation till T f inal = 2 seconds with a desired tolerancy tol = 1e -22. The initial value consistency is checked with Krawczyk which leads to x(0) ∈ [-1, -1] -1.e -18 , 1.e -18 .

The results at final time is given in Tab. 5.2. These results allow us to verify the enclosure of exact solution by the solution provided by our tool. In fact, we perform a verification of inclusion on the fly, that is to say, at each step of simulation. Moreover, at t = 2 seconds, the diameter of our solution is lower than 0.00056 when one of the best numerical method (without guarantee) [START_REF] Akinfenwa | Solving semi-explicit index-1 dae systems using l-stable extended block backward differentiation formula with continuous coefficients[END_REF], Extended Block Backward Differentiation Formula (EBBDF) at order 4, provides a result at 0.0002 from the exact value, which is comparable to our results (see Tab. 5.2). A ball of mass m, suspended to a massless rod of length l, under gravity g is defined by coordinates (p, q) such that:

Our method

       p = u q = v mu = -pλ mv = -qλ -g
with u and v the velocity of pendulum, and λ the force in the rod. And the algebraic relation: m(u 2 + v 2 ) -gq -l 2 λ = 0 with (p, q, u, v)0 = (1, 0, 0, 0) and λ0 ∈ [-0.1, 0.1]

We perform a simulation till T f inal = 1 second with a desired tolerancy tol = 1e -10. The initial value consistency is checked with Krawczyk which leads to λ(0) = 0. The computation takes less than 20 seconds on a dual core computer. With the files containing state and algebraic values w.r.t. time, we are able to plot three figures, given in Fig. 5.3.

Discussion

We see through these three examples that our method is correct (Ex. 5.2), efficient in comparison to the existing non guaranteed methods (Ex. 5.2), able to solve a classical problem (Ex. 5.3), while providing statistics on the computation progress (Ex. 5.1). Nevertheless, we can also conclude that the diameter of solution growth quickly and that time spent for the computation is sometimes significant. In the following section, we present an add-on in order to improve our results.

Additive contractors

Our tool is based on contractor programming [START_REF] Chabert | Contractor programming[END_REF], it means that we are able to add any contractors to the Alg. 2 to take into account some other constraints than g(t, x, y) = 0.

In the case of physical systems, constraints can appear from the context such as energy conservation, mechanical constraints, or whatever. This approach is substantially similar to [START_REF] Braun | Simulation of constrained mechanical systems part i: An equation of motion[END_REF]. In this latter, authors use the energy conservation as constraint in order to improve the precision of a numerical simulation. More interesting, some constraints come directly from the Pantelides algorithm [START_REF] Pantelides | The consistent initialization of differential-algebraic systems[END_REF] which is used to reduce the order of DAEs.

For example, on the case of Pendulum (Ex. 5.3), Pantelides algorithm is used to obtained m(u 2 + v 2 ) -gq -l 2 λ = 0 by differentiation of circle equation p 2 + q 2 -l 2 = 0:

   p 2 + q 2 -l 2 = 0 pu + qv = 0 m(u 2 + v 2 ) -gq 2 -l 2 p = 0
These three additives constraints can be used to generate a forward/backward contractor (with a propagation principle between constraints) for both state and algebraic variables. This contractor improve the result of our algorithm by changing Alg. 2 to Alg. 4.

Algorithm 4 Contract the a priori enclosures around solution

Require: [x], [ỹ], lte [x j+1 ] = [x], [y j+1 ] = [ỹ] while ([x j+1 ] or [y j+1 ] is improved) do [y j+1 ]∩ = (RADAU 3 f ([x j+1 ], [y j+1 ]) + lte) [x j+1 ]∩ = (F wdBwd g ([x j+1 ], [y j+1 ]) ∩ K g ([x j+1 ], [y j+1 ])) ([x

Conclusions and future work

To conclude, we present in this paper our method for the validated simulation of differential algebraic equations. The main issue being the existence and unicity of solution, we develop a novel operator which compute simultaneously the a priori enclosures of algebraic and state variables. For that, our operator mixes a classical Picard-Lindelöf (based on Taylor series) and a parametric Krawczyk operator. Then these enclosures is refined with a contractor programming approach by combining an integration scheme (Runge-Kutta RadauIIA) and an algebraic contractor (made with Krawczyk and forward/backward). Our complete algorithm is applied to three examples that we thought sufficient to prove the concept and the efficiency of the method. Moreover, a third contribution is also proposed in order to improve the result obtained. Indeed, it is often possible to obtain additive constraints on a physical system and these constraints can be easily used as contractors for both state and algebraic variables. Finally and to conclude on the contributions, our approach can naturally solve the initial consistency, which is one of the main issue in DAE community.

Several potential improvements to our method can be considered. The first one and the most important is the integration method used. Currently we have tried only Radau IIA order 3 method. A serious improvement can be probably obtained with an higher order Runge-Kutta method such as Radau IIA order 5 (but its coefficients are not exact) or Gauss order 6 (which has good properties but its local truncature error is computation time consuming). In addition, many improvement can be done on the global algorithm such as a better stepsize control and a better estimation for algebraic variable in the first trial of Picard-Krawczyk in order to reduce the number of rejected steps that is very important for the computation time.

  [yj]) ⊆ [ỹj]. The goal of Step 2 is thus to compute the tighter enclosure [yj+1] such that y(tj+1, tj, [yj]) ⊆ [yj+1] ⊆ [ỹj] .

Figure 2 . 2 :

 22 Figure 2.2: Different kinds of Runge-Kutta methods

13

 13

Theorem 4 . 1 . 1 ( 5 )

 4115 Part of Frobenius theorem) Let X and Y be Banach spaces, and A ⊂ X , B ⊂ Y a pair of open sets. Let F : A × B → L(X , Y) be a continuously differentiable function of the Cartesian product (which inherits a differentiable structure from its inclusion into X × Y) into the space L(X , Y) of continuous linear transformations of X into Y. A differentiable mapping u : A → B is a solution of the differential equationẏ = F (x, y) (4.3) if u(x) = F (x, u(x)) for all x ∈ A. Equation (4.3) is completely integrable if for each (x0, y0) ∈ A × B, there is a neighborhood U of x0 such that Equation (4.3) has a unique solution u(x) defined on U such that u(x0) = y0.Theorem 4.1.2 (Banach space version of Implicit Function Theorem) Let X , Y, Z be Banach spaces. Let the mapping f : X ×Y → Z be continuously Fréchet differentiable. If (x0, y0) ∈ X ×Y, f (x0, y0) = 0, and y → Df (x0, y0)(0, y) is a Banach space isomorphism from Y onto Z, then there exist neighborhoods U of x0 and V of y0 and a Fréchet differentiable function g : U → V such that f (x, g(x)) = 0 and f (x, y) = 0 if and only if y = g(x), for all (x, y) ∈ U × V. Proposition 4.1.1 (Set membership view of Frobenius) Let A and B two Banach spaces and g, f two continuously differentiable function of Cartesian product A×B into the set L(A, B) of continuous linear transformations of A into B, defined by: ẏ = f (t, x, y) Let X ⊂ A, and Y ⊂ B be a pairs of set. If ∀y ∈ Y, ∃! x ∈ X satisfaying Equation (4.5) and ∀x ∈ X , ∃! y ∈ Y satisfaying Equation (4.4), simultaneously, then Equation (4.4) is completely integrable and the solution exists and is unique in Y.

A

  direct translation of Proposition 4.1.1 to our problem gives: Proposition 4.1.2 Let [ỹ] be a Picard solution of ẏ ∈ f (t, [x], y) and [x] defined such that for each y ∈ [ỹ], ∃! x ∈ [x] : g(x, y) = 0, then ∃! ψ on the neighborhood of [x], and the solution of DAE exists and is unique in [ỹ]. In addition, the algebraic variable is enclosed by [x]. Proof: Direct application of Proposition 4.1.1.

Algorithm 2

 2 Contract the a priori enclosures around solutionRequire: [x], [ỹ], lte [x j+1 ] = [x], [y j+1 ] = [ỹ] while ([x j+1 ] or [y j+1 ] is improved) do // fix point [y j+1 ]∩ = (RADAU 3 f ([x j+1 ], [y j+1 ]) + lte) // 2.4 and 4.1 [x j+1 ]∩ = (F wdBwd g ([x j+1 ], [y j+1 ]) ∩ K g ([x j+1 ], [y j+1 ])) // [4] and 4.1.3 end while

  e -14 , 1.e-14

Finally, the lastFigure 5 . 2 :

 52 Figure 5.2: Problem statement of Ex. 5.3

Figure 5 . 3 :

 53 Figure 5.3: Plot of state variable p w.r.t. q (a), algebraic variable w.r.t. time (b) and algebraic variable w.r.t. state variable p (c) (for prob. 5.3)

  j+1 ]; [y j+1 ])∩ = F wdBwd ctc (([x j+1 ]; [y j+1 ])) end whileOur tool is applied with Alg. 2 and Alg. 4 to the Pendulum problem for a simulation till T f inal = 1.6 seconds and for a desired tolerancy tol = 10 -18 . The results are shown on Fig.6.1 and Fig.6.2. It is obvious that the addition of contractors has improved the simulation by reducing the size (and thus the pessimism) of the enclosures, approximately 50%, and even the time computation.

Figure 6 . 1 :

 61 Figure 6.1: Plot state variable p w.r.t. q (a), algebraic variable w.r.t. time (b) and algebraic variable w.r.t. state variable p (c) (for prob. 5.3) -with Alg. 2: computation time 28 minutes

Figure 6 . 2 :

 62 Figure 6.2: Plot state variable p w.r.t. q (a), algebraic variable w.r.t. time (b) and algebraic variable w.r.t. state variable p (c) (for prob. 5.3) -with Alg. 4: computation time 27 minutes

  Partially funded by the Academic and Research Chair:"Complex Systems Engineering"-Ecole polytechnique ∼ THALES ∼ FX ∼ DGA ∼ DASSAULT AVIATION ∼ DCNS Research ∼ ENSTA ParisTech ∼ Telecom ParisTech ∼ Fondation ParisTech ∼ FDO ENSTA

  w.r.t. a single constraint such that no solution of the constraint is lost in the box. Using a tree representation of the constraint for accelerating the contraction, this contractor isolates every occurrence xi in the expression and performs a natural evaluation of the corresponding function to contract[xi]. This combination can be easily done by the contractor programming view of our tool.Finally, these two obtained contractors are embedded in a fixpoint presented in Algorithm 2. The inputs of this algotithm are the enclosures and the lte computed with Algorithm 1. The outputs are the guaranteed enclosures of the state and algebraic solutions at the end of current time step.

	1/3 5/12 -1/12
	1	3/4	1/4
		3/4	1/4
	Figure 4.1: Butcher Tableau of Radau IIA order 3 (an Implicit Runge-Kutta
	method)
	operator from Section 4.1 and a forward/backward contractor coming from constraint
	programming. The Forward/Backward contractor (also called HC4-Revise [4]) can
	contract a box

Validated Simulation of

Algorithm 3 Simulation of DAE Require: y(0), x(0), T f inal , h min , tol // Problem statement t = 0 while (t < T f inal ) do // Till desired time

Examples

We apply our algorithm to three examples. The first one is quite basic but allows us to show the issue appearing in DAE simulation and the results obtained by our tool. This example is sufficient to point out the problem of stiffness generated by DAE problems. The second example is interesting because it has a known exact solution and allows us to verify the exactness of our method and judge the results with respect to the best existing numerical method. Finally, the third example is the classical pendulum problem, which permits to highlight the efficiency of our method.

Basic example

We start our experimentations with a basic example in one dimension for the state variable and one dimension for the algebraic variable:

We perform a simulation till T f inal = 4 seconds with a desired tolerancy tol = 1e -16. The initial value consistency is checked with Krawczyk which leads to x(0) = -1. The computation takes between 16 and 30 seconds depending of computer. Our tool provides three files containing:

• The values of state variable enclosure w.r. With the files containing state and algebraic values w.r.t. time, we are able to plot three figures, given in Fig. 5.1. On the figure (a) and (b), it is apparent that even if state variable and algebraic variable evolves exponentially but quite slowly, algebraic variable evolves in a stiff way w.r.t. state variable (c). In general, DAEs leads to a stiff problem, which is the motivation for using RADAU IIA method. This method is known for its efficiency in front of stiff problems.